MomentCalcBase.h

Go to the documentation of this file.

//# MomentCalcBase.h:
//# Copyright (C) 1997,1999,2000,2001,2002
//# 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: aips2-request@nrao.edu.
//#        Postal address: AIPS++ Project Office
//#                        National Radio Astronomy Observatory
//#                        520 Edgemont Road
//#                        Charlottesville, VA 22903-2475 USA
//#
//# $Id: MomentCalculator.h 20299 2008-04-03 05:56:44Z gervandiepen $

#ifndef IMAGEANALYSIS_MOMENTCALCBASE_H
#define IMAGEANALYSIS_MOMENTCALCBASE_H

namespace casa {

//# Forward Declarations
template <class T> class MomentsBase;
// <summary>
// Abstract base class for moment calculator classes
// </summary>
// <use visibility=export>
// <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
// </reviewed>
//
// <prerequisite>
//   <li> <linkto class="MomentsBase">MomentsBase</linkto>
//   <li> <linkto class="ImageMoments">ImageMoments</linkto>
//   <li> <linkto class="casacore::LatticeApply">casacore::LatticeApply</linkto>
//   <li> <linkto class="casacore::LineCollapser">casacore::LineCollapser</linkto>
// </prerequisite>
//
// <synopsis>
//  This class, its concrete derived classes, and the classes casacore::LineCollapser,
//  ImageMoments, MSMoments, and casacore::LatticeApply are connected as follows.   casacore::LatticeApply offers 
//  functions so that the application programmer does not need to worry about how 
//  to optimally iterate through a casacore::Lattice; it deals with tiling and to a 
//  lesser extent memory.    casacore::LatticeApply functions are used by offering a class 
//  object to them that has a member function with a name and signature 
//  specified by an abstract base class that casacore::LatticeApply uses and the 
//  offered class inherits from.   Specifically, in this case, MomentCalcBase
//  inherits from casacore::LineCollapser and casacore::LatticeApply uses objects and methods of this
//  class (but does not inherit from it).  This defines the functions
//  <src>collapse</src> and <src>multiProcess</src> which operate on a vector
//  extracted from a Lattice.  The former returns one number, the latter a vector
//  of numbers from that profile.  MomentCalcBase is a base class for
//  for moment calculation and the <src>multiProcess</src>
//  functions are used to compute moments  (e.g., mean, sum, sum squared, 
//  intensity weighted velocity etc).
//
//  It is actually the concrete classes derived from MomentCalcBase (call them,
//  as a group, the MomentCalculator classes) that implement the <src>multiProcess</src> 
//  functions.  These derived classes allow different 
//  algorithms to be written with which moments of the vector can be computed. 
//
//  Now, so far, we have a casacore::LatticeApply function which iterates through Lattices,
//  extracts vectors, and offers them up to functions implemented in the derived 
//  MomentCalculator classes to compute the moments.   As well as that, we need some
//  class to actually construct the MomentCalculator classes and to feed them to 
//  LatticeApply.   This is the role of the ImageMoments or MSMoments classes.  
//  They are a high level 
//  class which takes control information from users specifying which moments they 
//  would like to calculate and how. They also provide the ancilliary masking lattice to 
//  the MomentCalculator constructors. The actual computational work is done by the 
//  MomentCalculator classes. So MomentsBase, MomentCalcBase and their derived 
//  MomentCalculator classes are really one unit; none of them are useful without 
//  the others.  The separation of functionality is caused by having the
//  casacore::LatticeApply class that knows all about optimally iterating through Lattices.
//
//  The coupling between these classes is done partly by the "friendship".   MomentsBase and 
//  its inheritances 
//  grant friendship to MomentCalcBase so that the latter has access to the private data and 
//  private functions of the formers.  MomentCalcBase then operates as an interface between 
//  its derived MomentCalculator classes and ImageMoments or MSMoments. It retrieves private data 
//  from these classes, and also activates private functions in these classes, on behalf 
//  of the MomentCalculator classes. The rest of the coupling is done via the constructors 
//  of the derived MomentCalculator classes.  
//
//  Finally, MomentCalcBase also has a number of protected functions that are common to its
//  derived classes (e.g. fitting, accumulating sums etc).  It also has protected
//  data that is common to all the MomentCalculator classes.  This protected data is accessed 
//  directly by name rather than with interface functions as there is too much of it.  Of 
//  course, since MomentCalcBase is an abstract base class, it is up to the MomentCalculator 
//  classes to give the MomentCalcBase protected data objects values.
//
//  For discussion about different moments and algorithms to compute them see the 
//  discussion in <linkto class="MomentsBase">MomentsBase</linkto>, 
//  <linkto class="ImageMoments">ImageMoments</linkto>, 
//  <linkto class="MSMoments">MSMoments</linkto> and also in
//  the derived classes documentation.
// </synopsis>
//
// <example>
//  Since MomentCalcBase is an abstract class, we defer code examples to
//  the derived classes.
// </example>
//
// <motivation>
// We were desirous of writing functions to optimally iterate through Lattices
// so that the application programmer did not have to know anything about tiling
// or memory if possible.   These are the casacore::LatticeApply functions. To incorporate 
// MomentsBase and its inheritances into this scheme required some of it to be shifted into 
// MomentCalcBase and its derived classes.
// </motivation>
//
// <note role=tip>
// Note that there are is assignment operator or copy constructor.
// Do not use the ones the system would generate either.
// </note>
//
// <todo asof="yyyy/mm/dd">
//  <li> Derive more classes !
// </todo>


template <class T> class MomentCalcBase : public casacore::LineCollapser<T,T> {
public:

    using AccumType = typename casacore::NumericTraits<T>::PrecisionType;
    using DataIterator = typename casacore::Vector<T>::const_iterator;
    using MaskIterator = casacore::Vector<casacore::Bool>::const_iterator;

    virtual ~MomentCalcBase();

    // Returns the number of failed fits if doing fitting
    virtual inline casacore::uInt nFailedFits() const { return nFailed_p; }

protected:

    // A number of private data members are kept here in the base class
    // as they are common to the derived classes.  Since this class
    // is abstract, they have to be filled by the derived classes.

    // CoordinateSystem
    casacore::CoordinateSystem cSys_p;

    // This vector is a container for all the possible moments that
    // can be calculated.  They are in the order given by the MomentsBase
    // enum MomentTypes
    casacore::Vector<T> calcMoments_p;
    casacore::Vector<casacore::Bool> calcMomentsMask_p;

    // This vector tells us which elements of the calcMoments_p vector
    // we wish to select
    casacore::Vector<casacore::Int> selectMoments_p;

    // Although the general philosophy of these classes is to compute
    // all the posisble moments and then select the ones we want,
    // some of them are too expensive to calculate unless they are
    // really wanted.  These are the median moments and those that
    // require a second pass.  These control Bools tell us whether
    // we really want to compute the expensive ones.
    casacore::Bool doMedianI_p, doMedianV_p, doAbsDev_p;

    // These vectors are used to transform coordinates between pixel and world
    casacore::Vector<casacore::Double> pixelIn_p, worldOut_p;

    // All computations involving casacore::Coordinate conversions are relatively expensive
    // These Bools signifies whether we need coordinate calculations or not for
    // the full profile, and for some occaisional calculations
    casacore::Bool doCoordProfile_p, doCoordRandom_p;

    // This vector houses the world coordinate values for the profile if it
    // was from a separable axis. This means this vector can be pre computed
    // just once, instead of working out the coordinates for each profile
    // (expensive).  It should only be filled if doCoordCalc_p is true
    casacore::Vector<casacore::Double> sepWorldCoord_p;

    // This vector is used to hold the abscissa values
    casacore::Vector<T> abcissa_p;

    // This string tells us the name of the moment axis (VELO or FREQ etc)
    casacore::String momAxisType_p;

    // This is the number of Gaussian fits that failed.
    casacore::uInt nFailed_p;

    // This scale factor is the increment along the moment axis
    // applied so that units for the Integrated moment are like
    // Jy/beam.km/s (or whatever is needed for the moment axis units)
    // For non-linear velocities (e.g. optical) this is approximate
    // only and is computed at the reference pixel
    casacore::Double integratedScaleFactor_p;

    // Accumulate statistical sums from a vector
    inline void accumSums(
        typename casacore::NumericTraits<T>::PrecisionType& s0,
        typename casacore::NumericTraits<T>::PrecisionType& s0Sq,
        typename casacore::NumericTraits<T>::PrecisionType& s1,
        typename casacore::NumericTraits<T>::PrecisionType& s2,
        casacore::Int& iMin, casacore::Int& iMax, T& dMin, T& dMax,
        casacore::Int i, T datum, casacore::Double coord
    ) const {
        // Accumulate statistical sums from this datum
        //
        // casacore::Input:
        //  i              Index
        //  datum          Pixel value
        //  coord          casacore::Coordinate value on moment axis
        // casacore::Input/output:
        //  iMin,max       index of dMin and dMax
        //  dMin,dMax      minimum and maximum value
        // Output:
        //  s0             sum (I)
        //  s0Sq           sum (I*I)
        //  s1             sum (I*v)
        //  s2             sum (I*v*v)
        typename casacore::NumericTraits<T>::PrecisionType dDatum = datum;
        s0 += dDatum;
        s0Sq += dDatum*dDatum;
        s1 += dDatum*coord;
        s2 += dDatum*coord*coord;
        if (datum < dMin) {
            iMin = i;
            dMin = datum;
        }
        if (datum > dMax) {
            iMax = i;
            dMax = datum;
        }
    }

    // Determine if the spectrum is pure noise
    casacore::uInt allNoise(T& dMean,
        const casacore::Vector<T>& data,
        const casacore::Vector<casacore::Bool>& mask,
        T peakSNR,
        T stdDeviation
    ) const;

    // Check validity of constructor inputs
    void constructorCheck(
        casacore::Vector<T>& calcMoments,
        casacore::Vector<casacore::Bool>& calcMomentsMask,
        const casacore::Vector<casacore::Int>& selectMoments,
        casacore::uInt nLatticeOut
    ) const;

    // Find out from the selectMoments array whether we want
    // to compute the more expensive moments
    void costlyMoments(
        MomentsBase<T>& iMom, casacore::Bool& doMedianI,
        casacore::Bool& doMedianV, casacore::Bool& doAbsDev
    ) const;

    // Return the casacore::Bool saying whether we need to compute coordinates
    // or not for the requested moments
    void doCoordCalc(
        casacore::Bool& doCoordProfile,
        casacore::Bool& doCoordRandom,
        const MomentsBase<T>& iMom
    ) const;

    // Return the casacore::Bool from the ImageMoments or MSMoments object saying whether we
    // are going to fit Gaussians to the profiles or not.
    casacore::Bool doFit(const MomentsBase<T>& iMom) const;

    // Find the next masked or unmasked point in a vector
    casacore::Bool findNextDatum(
        casacore::uInt& iFound, const casacore::uInt& n,
        const casacore::Vector<casacore::Bool>& mask, const casacore::uInt& iStart,
        const casacore::Bool& findGood
    ) const;

    // Fit a Gaussian to x and y arrays given guesses for the gaussian parameters
    casacore::Bool fitGaussian(
        casacore::uInt& nFailed, T& peak, T& pos,
        T& width, T& level, const casacore::Vector<T>& x,
        const casacore::Vector<T>& y, const casacore::Vector<casacore::Bool>& mask,
        T peakGuess, T posGuess, T widthGuess,
        T levelGuess
    ) const;

    // Automatically fit a Gaussian to a spectrum, including finding the
    // starting guesses.
    casacore::Bool getAutoGaussianFit(
        casacore::uInt& nFailed, casacore::Vector<T>& gaussPars,
        const casacore::Vector<T>& x, const casacore::Vector<T>& y,
        const casacore::Vector<casacore::Bool>& mask, T peakSNR,
        T stdDeviation
    ) const;

    // Automatically work out a guess for the Gaussian parameters
    // Returns false if all pixels masked.
    casacore::Bool getAutoGaussianGuess(
        T& peakGuess, T& posGuess,
        T& widthGuess, T& levelGuess,
        const casacore::Vector<T>& x, const casacore::Vector<T>& y,
        const casacore::Vector<casacore::Bool>& mask
    ) const;

    // Compute the world coordinate for the given moment axis pixel
    inline casacore::Double getMomentCoord(
        const MomentsBase<T>& iMom, casacore::Vector<casacore::Double>& pixelIn,
        casacore::Vector<casacore::Double>& worldOut, casacore::Double momentPixel,
        casacore::Bool asVelocity=false
    ) const {
        // Find the value of the world coordinate on the moment axis
        // for the given moment axis pixel value.
        //
        // Input
        //   momentPixel   is the index in the profile extracted from the data
        // casacore::Input/output
        //   pixelIn       Pixels to convert.  Must all be filled in except for
        //                 pixelIn(momentPixel).
        //   worldOut      casacore::Vector to hold result
        //
        // Should really return a casacore::Fallible as I don't check and see
        // if the coordinate transformation fails or not
        //
        // Should really check the result is true, but for speed ...
        pixelIn[iMom.momentAxis_p] = momentPixel;
        cSys_p.toWorld(worldOut, pixelIn);
        if (asVelocity) {
            casacore::Double velocity;
            cSys_p.spectralCoordinate().frequencyToVelocity(
                velocity, worldOut(iMom.worldMomentAxis_p)
            );
            return velocity;
        }
        return worldOut(iMom.worldMomentAxis_p);
    }

    // Examine a mask and determine how many segments of unmasked points
    // it consists of.
    void lineSegments (
        casacore::uInt& nSeg, casacore::Vector<casacore::uInt>& start,
        casacore::Vector<casacore::uInt>& nPts, const casacore::Vector<casacore::Bool>& mask
    ) const;

    // Return the moment axis from the ImageMoments object
    casacore::Int& momentAxis(MomentsBase<T>& iMom) const;

    // Return the name of the moment/profile axis
    casacore::String momentAxisName(
        const casacore::CoordinateSystem&,
        const MomentsBase<T>& iMom
    ) const;

    // Return the peak SNR for determination of all noise spectra from
    // the ImageMoments or MSMoments object
    T& peakSNR(MomentsBase<T>& iMom) const;

    // Return the selected pixel intensity range from the ImageMoments or MSMoments
    // object and the Bools describing whether it is inclusion or exclusion
    void selectRange(
        casacore::Vector<T>& pixelRange,
        casacore::Bool& doInclude,
        casacore::Bool& doExlude,
        MomentsBase<T>& iMom
    ) const;

    // The MomentCalculators compute a vector of all possible moments.
    // This function returns a vector which selects the desired moments from that
    // "all moment" vector.
    casacore::Vector<casacore::Int> selectMoments(MomentsBase<T>& iMom) const;

    // Fill the ouput moments array
    void setCalcMoments(
        const MomentsBase<T>& iMom, casacore::Vector<T>& calcMoments,
        casacore::Vector<casacore::Bool>& calcMomentsMask, casacore::Vector<casacore::Double>& pixelIn,
        casacore::Vector<casacore::Double>& worldOut, casacore::Bool doCoord,
        casacore::Double integratedScaleFactor, T dMedian,
        T vMedian, casacore::Int nPts,
        typename casacore::NumericTraits<T>::PrecisionType s0,
        typename casacore::NumericTraits<T>::PrecisionType s1,
        typename casacore::NumericTraits<T>::PrecisionType s2,
        typename casacore::NumericTraits<T>::PrecisionType s0Sq,
        typename casacore::NumericTraits<T>::PrecisionType sumAbsDev,
        T dMin, T dMax, casacore::Int iMin, casacore::Int iMax
    ) const;

    // Fill a string with the position of the cursor
    void setPosLabel(casacore::String& title, const casacore::IPosition& pos) const;

    // Install casacore::CoordinateSystem and SpectralCoordinate
    // in protected data members
    void setCoordinateSystem(
        casacore::CoordinateSystem& cSys, MomentsBase<T>& iMom
    );

    // Set up separable moment axis coordinate vector and
    // conversion vectors if not separable
    void setUpCoords(
        const MomentsBase<T>& iMom, casacore::Vector<casacore::Double>& pixelIn,
        casacore::Vector<casacore::Double>& worldOut, casacore::Vector<casacore::Double>& sepWorldCoord,
        casacore::LogIO& os, casacore::Double& integratedScaleFactor,
        const casacore::CoordinateSystem& cSys, casacore::Bool doCoordProfile,
        casacore::Bool doCoordRandom
    ) const;

    // Return standard deviation of image from ImageMoments or MSMoments object
    T& stdDeviation(MomentsBase<T>& iMom) const;

protected:
    // Check if #pixels is indeed 1.
    virtual void init (casacore::uInt nOutPixelsPerCollapse);
};

}

#ifndef CASACORE_NO_AUTO_TEMPLATES
#include <imageanalysis/ImageAnalysis/MomentCalcBase.tcc>
#endif
#endif