ImageBeamSet.h

Go to the documentation of this file.

// # Copyright (C) 1996,1997,1998,1999,2000,2001,2003
// # Associated Universities, Inc. Washington DC, USA.
// #
// # This library is free software; you can redistribute it and/or modify it
// # under the terms of the GNU Library General Public License as published by
// # the Free Software Foundation; either version 2 of the License, or (at your
// # option) any later version.
// #
// # This library is distributed in the hope that it will be useful, but WITHOUT
// # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
// # License for more details.
// #
// # You should have received a copy of the GNU Library General Public License
// # along with this library; if not, write to the Free Software Foundation,
// # Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
// #
// # Correspondence concerning AIPS++ should be addressed as follows:
// #        Internet email: casa-feedback@nrao.edu.
// #        Postal address: AIPS++ Project Office
// #                        National Radio Astronomy Observatory
// #                        520 Edgemont Road
// #                        Charlottesville, VA 22903-2475 USA
 
#ifndef IMAGES_IMAGEBEAMSET_H
#define IMAGES_IMAGEBEAMSET_H
 
#include <casacore/casa/aips.h>
#include <casacore/casa/Arrays/Matrix.h>
#include <casacore/scimath/Mathematics/GaussianBeam.h>
// #include <casacore/measures/Measures/Stokes.h>
#include <map>
 
namespace casacore {
 
class SpectralCoordinate;
 
class CoordinateSystem;
 
// <summary>
// Represents a set of restoring beams associated with an image.
// </summary>
 
// <use visibility=export>
 
// <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
// </reviewed>
 
// <prerequisite>
// </prerequisite>
 
// <etymology>
// A Set of Beams associated with an Image.
// </etymology>
 
// <synopsis>
// This class represents a set of restoring beams associated with
// a deconvolved image. Internally, the beams are stored in a Matrix in
// which the first dimension represents the spectral axis and the second
// dimension represents the polarization axis. Methods which take the number
// of channels and stokes as the input parameters will accept 0 for either of
// these, in the cases where the corresponding axis is absent, but internally,
// the associated axis of the storage Matrix will be set to one since a Matrix
// with one dimension of length 0 must be empty. If one (or both) of the axes is
// of length 1, the beam associated with that position is valid for all channels or
// stokes at the position. For example, if one has an image with 10 spectral channels
// and 4 stokes, and one has a beam set of dimensions (1, 4) associated with that image,
// all channels for a given stokes will have the same beam. Similarly, if the beam set
// is of shape (10, 1) in this case, all stokes will have the same beam for a given channel.
// If the axis lengths of the beam set are greater than one, they must be exactly
// the same length of the corresponding axes in the associated image.
// </synopsis>
//
// <example>
 
// </example>
 
// <motivation>
// Restoring beams are used many places in image analysis tasks.
// </motivation>
 
// <todo>
// </todo>
 
class ImageBeamSet {
 public:
  typedef Array<GaussianBeam>::const_iterator BeamIter;
 
  // Construct an empty beam set.
  ImageBeamSet();
 
  // Construct a beam set from an 2-D array of beams representing
  // the frequency and stokes axis.
  // Axis length 1 means it is valid for all channels cq. stokes.
  // If the image has 0 spectral channels or stokes, the corresponding
  // length of the axis in the provided matrix should be 1.
  ImageBeamSet(const Matrix<GaussianBeam>& beams);
 
  // construct an ImageBeamSet representing a single beam which is valid for
  // all channels and stokes
  ImageBeamSet(const GaussianBeam& beam);
 
  // Create an ImageBeamSet of the specified shape with all
  // GaussianBeams initialized to <src>beam</src>.
  ImageBeamSet(uInt nchan, uInt nstokes, const GaussianBeam& beam = GaussianBeam::NULL_BEAM);
 
  // The copy constructor (reference semantics).
  ImageBeamSet(const ImageBeamSet& other);
 
  ~ImageBeamSet();
 
  // Assignment can change the shape (copy semantics).
  ImageBeamSet& operator=(const ImageBeamSet& other);
 
  // Beam sets are equal if the shapes and all corresponding beams are equal.
  Bool operator==(const ImageBeamSet& other) const;
  Bool operator!=(const ImageBeamSet& other) const;
 
  // Beam sets are equivalent if both have no beams or if the
  // expanded sets are equal. Expanded means that an axis can have
  // length 0 or 1 and is (virtually) expanded to the length of the matching
  // axis in the other beam set.
  Bool equivalent(const ImageBeamSet& that) const;
 
  // Get the number of elements in the beam array.
  // <group>
  uInt nelements() const { return _beams.size(); }
  uInt size() const { return _beams.size(); }
  // </group>
 
  Bool hasSingleBeam() const { return _beams.size() == 1; }
 
  // Does this beam set contain multiple beams?
  Bool hasMultiBeam() const { return _beams.size() > 1; }
 
  // Is the beam set empty?
  Bool empty() const { return _beams.empty(); }
 
  // Get the shape of the beam array. The minimum value for
  // a component of the returned IPosition is always 1.
  const IPosition& shape() const { return _beams.shape(); }
 
  // Get the number of channels in the beam array. Note that this will
  // always return a minimum of 1, even if nchan was specified as 0 on construction.
  uInt nchan() const { return _beams.shape()[0]; }
 
  // Get the number of stokes in the beam array. Note that this will always
  // return a minimum of 1, even if nstokes was specified as 0 on construction.
  uInt nstokes() const { return _beams.shape()[1]; }
 
  // Get the single global beam. If there are multiple beams,
  // an exception is thrown.
  const GaussianBeam& getBeam() const;
 
  // Get the beam at the specified location.
  // Note that a single channel or stokes in the beam set is valid for
  // all channels cq. stokes.
  // <group>
  const GaussianBeam& getBeam(Int chan, Int stokes) const;
  const GaussianBeam& operator()(Int chan, Int stokes) const { return getBeam(chan, stokes); }
  // </group>
 
  // Get a beam at the given 2-dim IPosition. It should match exactly,
  // thus a single channel or stokes in the beam set is not valid for all.
  // const GaussianBeam& getBeam(const IPosition& pos) const
  //  { return _beams(pos); }
 
  // Set the beam at the given location.
  // The location must be within the beam set shape.
  // If <src>chan</src> or <src>stokes</src> is negative, then the beam applies
  // to all channels or stokes, respectively. If both are negative, the specified
  // beam becomes the global beam and the beam set is resized to (1, 1).
  void setBeam(Int chan, Int stokes, const GaussianBeam& beam);
 
  // Resize the beam array. <src>nchan</src>=0 or <src>nstokes</src>=0
  // is silently changed to 1.
  void resize(uInt nchan, uInt nstokes);
 
  // Return a subset of the beam array.
  // The slicer is usually the slicer used for a subimage.
  // The slicer can contain multiple stokes or channels, even if the
  // beam set has only one.
  ImageBeamSet subset(const Slicer& imageSlicer, const CoordinateSystem& csys) const;
 
  // Get the beam array.
  const Matrix<GaussianBeam>& getBeams() const { return _beams; }
 
  // Set the beams in this beam set.
  // The shape of the given array must match the beam set.
  // It also matches if an axis in array or beam set has length 1, which
  // means that it expands to the other length.
  void setBeams(const Matrix<GaussianBeam>& beams);
 
  // Set all beams to the same value.
  void set(const GaussianBeam& beam);
 
  // Get the beam in the set which has the smallest area.
  GaussianBeam getMinAreaBeam() const { return _minBeam; }
 
  // Get the beam in the set which has the largest area.
  // Get the beam in the set which has the largest area.
  GaussianBeam getMaxAreaBeam() const { return _maxBeam; }
 
  // Get the beam in the set which has the median area.
  GaussianBeam getMedianAreaBeam() const;
 
  // Get the position of the beam with the minimum area.
  IPosition getMinAreaBeamPosition() const { return _minBeamPos; }
 
  // Get the position of the beam with the maximum area.
  IPosition getMaxAreaBeamPosition() const { return _maxBeamPos; }
 
  // Get the minimal, maximal, and median area beams and positions in the
  // beam set matrix for the given stokes. If the stokes axis has length 1
  // in the beam matrix, it is valid for all stokes and no checking is done
  // that <src>stokes</src> is valid; the requested beam for the entire beam set
  // is simply returned in this case.
  // If the number of stokes in the beam matrix is >1, checking is done that
  // the specified value of <src>stokes</src> is valid and if not, an exception
  // is thrown.
  // <group>
  const GaussianBeam& getMinAreaBeamForPol(IPosition& pos, uInt stokes) const;
 
  const GaussianBeam& getMaxAreaBeamForPol(IPosition& pos, uInt stokes) const;
 
  const GaussianBeam& getMedianAreaBeamForPol(IPosition& pos, uInt stokes) const;
  // </group>
 
  static const String& className();
 
  // Get the beam that has the smallest minor axis. If multiple beams have
  // the smallest minor axis, the beam in this subset with the smallest area
  // will be returned.
  const GaussianBeam getSmallestMinorAxisBeam() const;
 
  // convert ImageBeamSet to and from record
  // <group>
  static ImageBeamSet fromRecord(const Record& rec);
  Record toRecord() const;
  //</group>
 
  // If verbose, log all beams, if not just summarize beam stats.
  void summarize(LogIO& log, Bool verbose, const CoordinateSystem& csys) const;
 
  // Modify the beam set by rotating all beams counterclockwise through the
  // specified angle. If unwrap=True, unwrap the new position angle(s) so that
  // it falls in the range -90 to 90 degrees before setting it.
  void rotate(const Quantity& angle, Bool unwrap = False);
 
  // get all the beam areas in a single quantum matrix.
  const Quantum<Matrix<double>> getAreas() const;
 
  // get all the major axes, minor axes, and pas in quantum matrices.
  // The returned map has keys "major", "minor", and "pa".
  const std::map<String, Quantum<Matrix<double>>> paramMatrices(
      const Unit& majminUnit = Unit("arcsec"), const Unit& paUnit = "deg") const;
 
 private:
  static const String _DEFAULT_AREA_UNIT;
 
  Matrix<GaussianBeam> _beams;
  Matrix<Double> _areas;
  String _areaUnit;
  GaussianBeam _minBeam, _maxBeam;
  IPosition _minBeamPos, _maxBeamPos;
 
  void _calculateAreas();
 
  // common code for replacing a beam in a multi-beam set
  void _replaceBeam(const GaussianBeam& beam, const IPosition& location1,
                    const IPosition& location2, Bool overwriteMaxMin);
 
  // Show the spectral info.
  static void _chanInfoToStream(ostream& os, const SpectralCoordinate* spCoord, const uInt chan,
                                const uInt chanWidth, const uInt freqPrec, const uInt velWidth,
                                const uInt velPrec);
 
  // Show the beam info.
  static void _beamToStream(ostream& os, const GaussianBeam& beam, const Unit& unit);
};
 
// Show the beam set info.
ostream& operator<<(ostream& os, const ImageBeamSet& beamSet);
 
}  // namespace casacore
 
#endif