MomentsBase.h

Go to the documentation of this file.

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

#ifndef IMAGES_MOMENTSBASE_H
#define IMAGES_MOMENTSBASE_H

#include <casa/aips.h>
#include <coordinates/Coordinates/CoordinateSystem.h>
#include <casa/Quanta/Quantum.h>
#include <measures/Measures/MDoppler.h>
#include <casa/Logging/LogIO.h>
#include <casa/Arrays/Vector.h>
#include <casa/iosfwd.h>

namespace casa {

template <class T> class MomentCalcBase;
class IPosition;
class String;
class Unit;

// <summary>
// This class is a base class for generating moments from an image or a spectral data.
// </summary>

// <use visibility=export>

// <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
// </reviewed>

// <prerequisite>
//   <li> <linkto class="ImageMoments">ImageMoments</linkto>   
//   <li> <linkto class="MSMoments">MSMoments</linkto>   
//   <li> <linkto class="LatticeApply">LatticeApply</linkto>   
//   <li> <linkto class="MomentCalcBase">MomentCalcBase</linkto>
// </prerequisite>

// <etymology>
//   This class is an abstract class to compute moments from images or spectral data.
// </etymology>

// <synopsis>
//  The primary goal of MSMoments, ImageMoments, and MSMoments are to help 
//  spectral-line astronomers analyze their multi-dimensional images or 
//  spectral data (in the form of MeasurementSet) by generating moments of 
//  a specified axis.  ImageMoments is a specialized class to generate moments 
//  from images, while MSMoments is designed properly for MeasurementSet input.
//  MSMoments class is an abstract class that is inherited by the above two 
//  concrete classes.
//  MomentsBase provides interface layer to the MomentCalculators so that 
//  functionalities in MomentCalculators can be shared with ImageMoments and 
//  MSMoments.
//  The word "moment" is used loosely here.  It refers to collapsing an axis
//  to one pixel and putting the value of that pixel (for all of the other 
//  non-collapsed axes) to something computed from the data values along
//  the moment axis.  For example, take an RA-DEC-Velocity cube, collapse
//  the velocity axis by computing the mean intensity at each RA-DEC
//  pixel.  This class and its inheritances offer many different moments and a variety of
//  interactive and automatic ways to compute them. 
//

// <motivation>
//  MSMoments is defined so that moments can be generated from both images 
//  and spectral data (in the form of MeasurementSet).  
// </motivation>


template <class T> class MomentsBase {
public:

    // Note that if I don't put MomentCalcBase as a forward declaration
    // and use instead  "friend class MomentCalcBase<T>"  The gnu compiler
    // fails with a typedef problem.  But I can't solve it with say
    // <src>typedef MomentCalcBase<T> gpp_type;</src>  because you can only do a
    // typedef with an actual type, not a <tt>T</tt> !
    friend class MomentCalcBase<T>;

    // The <src>enum MethodTypes</src> is provided for use with the
    // <src>setWinFitMethod</src> function.  It gives the allowed moment
    // methods which are available with this function.  The use of these
    // methods is described further with the description of this function
    // as well as in the general documentation earlier.
    enum MethodTypes {
        // Invokes the spectral windowing method
        WINDOW,
        // Invokes Gaussian fitting
        FIT,
        NMETHODS
    };

    // This <src>enum MomentTypes</src> is provided for use with the
    // <src>setMoments</src> function.  It gives the allowed moment
    // types that you can ask for.

    enum MomentTypes {
        // The average intensity
        AVERAGE,
        // The integrated intensity
        INTEGRATED,
        // The intensity weighted mean coordinate (usually velocity)
        WEIGHTED_MEAN_COORDINATE,
        // The intensity weighted coordinate (usually velocity) dispersion
        WEIGHTED_DISPERSION_COORDINATE,
        // The median intensity
        MEDIAN,
        // The median coordinate (usually velocity). Treat the spectrum as
        // a probability distribution, generate the cumulative distribution,
        // and find the coordinate corresponding to the 50% value.
        MEDIAN_COORDINATE,
        // The standard deviation about the mean of the intensity
        STANDARD_DEVIATION,
        // The rms of the intensity
        RMS,
        // The absolute mean deviation of the intensity
        ABS_MEAN_DEVIATION,
        // The maximum value of the intensity
        MAXIMUM,
        // The coordinate (usually velocity) of the maximum value of the intensity
        MAXIMUM_COORDINATE,
        // The minimum value of the intensity
        MINIMUM,
        // The coordinate (usually velocity) of the minimum value of the intensity
        MINIMUM_COORDINATE,
        // Total number
        NMOMENTS,
        // Default value is the integrated intensity
        DEFAULT = INTEGRATED
    };
    
    // Destructor
    virtual ~MomentsBase();

    // Set the desired moments via an <src>Int</src> array.  Each <src>Int</src>
    // specifies a different moment; the allowed values and their meanings
    // are given by the <src>enum MomentTypes</src>.   A return value
    // of <src>False</src> indicates you asked for an out of range
    // moment.  If you don't call this function, the default state of the
    // class is to request the integrated intensity.
    Bool setMoments (const Vector<Int>& moments);

    // Set the moment axis (0 relative).  A return value of <src>False</src>
    // indicates that the axis was not contained in the image. If you don't
    // call this function, the default state of the class is to set the
    // moment axis to the spectral axis if it can find one.  Otherwise
    // an error will result.
    virtual Bool setMomentAxis (Int) = 0;

    // The method by which you compute the moments is specified by calling
    // (or not calling) the <src>setWinFitMethod</src> and
    // <src>setSmoothMethod</src> functions.  The default state of the class
    // is to compute directly on all (or some according to <src>setInExClude</src>)
    // of the pixels in the spectrum.  Calling these functions modifies the
    // computational state to something more complicated.
    //
    // The <src>setWinMethod</src> function requires an <src>Int</src> array
    // as its argument.  Each <src>Int</src> specifies a different method
    // that you can invoke (either separately or in combination).  The
    // allowed values and their meanings are given by the
    // <src>enum MethodTypes</src>.
    //
    // Both the windowing and fitting methods have interactive modes. The
    // windowing method also has a fitting flavour, so if you set both
    // MomentsBase::WINDOW and MomentsBase::FIT, you would be invoking the
    // windowing method but determining the window by fitting Gaussians
    // automatically (as MomentsBase::INTERACTIVE) was not given.
    //
    // If you don't call this function, then neither the windowing nor fitting
    // methods are invoked.  A return value of <src>False</src> indicates
    // that you asked for an illegal method.
    Bool setWinFitMethod(const Vector<Int>& method);

    // This function invokes smoothing of the input image.  Give <src>Int</src>
    // arrays for the axes (0 relative) to be smoothed and the smoothing kernel
    // types (use the <src>enum KernelTypes</src>) for each axis.  Give a
    // <src>Double</src> array for the widths (full width for BOXCAR and full
    // width at half maximum for GAUSSIAN) in pixels of the smoothing kernels for
    // each axis.  For HANNING smoothing, you always get the quarter-half-quarter
    // kernel (no matter what you might ask for).  A return value of <src>False</src>
    // indicates that you have given an inconsistent or invalid set of smoothing
    // parameters.  If you don't call this function the default state of the
    // class is to do no smoothing.  The kernel types are specified with
    // the VectorKernel::KernelTypes enum
    // <group>
    virtual Bool setSmoothMethod(
            const Vector<Int>&, const Vector<Int>&,
            const Vector<Quantum<Double> >&
    ) = 0;

    Bool setSmoothMethod(
        const Vector<Int>& smoothAxes,
        const Vector<Int>& kernelTypes,
        const Vector<Double>& kernelWidths
    );
    // </group>

    // You may specify a pixel intensity range as either one for which
    // all pixels in that range are included in the moment calculation,
    // or one for which all pixels in that range are excluded from the moment
    // calculations.  One or the other of <src>include</src> and <src>exclude</src>
    // must therefore be a zero length vector if you call this function.
    // A return value of <src>False</src> indicates that you have given both
    // an <src>include</src> and an <src>exclude</src> range.  If you don't call
    // this function, the default state of the class is to include all pixels.
    void setInExCludeRange(
        const Vector<T>& include, const Vector<T>& exclude
    );

    // This function is used to help assess whether a spectrum along the moment
    // axis is all noise or not.  If it is all noise, there is not a lot of point
    // to trying to computing the moment.
    // <src>peakSNR</src> is the signal-to-noise ratio of the peak value in the
    // spectrum below which the spectrum is considered pure noise.
    // <src>stdDeviation</src> is the standard deviation of the noise for the
    // input image.
    //
    // Default values for one or the other parameter are indicated by giving zero.
    //
    // The default state of the class then is to set <src>peakSNR=3</src>
    // and/or to work out the noise level from a Gaussian fit to a histogram
    // (above 25%) of the entire image (it is very hard to get an accurate
    // estimate of the noise a single spectrum).
    void setSnr(const T& peakSNR, const T& stdDeviation);

    // This is the output file name for the smoothed image.   It can be useful
    // to have access this to this image when trying to get the pixel
    // <src>include</src> or <src>exclude</src> range correct for the smooth-clip
    // method.  The default state of the class is to not output the smoothed image.
    Bool setSmoothOutName(const String& smOut);

    // Set Velocity type.  This is used for moments for which the moment axis is
    // a spectral axis for which the coordinate is traditionally presented in
    // km/s   You can select the velocity definition. The default state of the
    // class is to use the radio definition.
    void setVelocityType (MDoppler::Types type);

    // Reset argument error condition.  If you specify invalid arguments to
    // one of the above functions, an internal flag will be set which will
    // prevent the <src>createMoments</src> function, which is defined in its inheritances,
    // from doing anything
    // (should you have chosen to igmore the Boolean return values of the functions).
    // This function allows you to reset that internal state to good.
    void resetError () {goodParameterStatus_p = True; error_p = "";};

    // Recover last error message
    String errorMessage() const {return error_p;};

    // Get CoordinateSystem
    virtual const CoordinateSystem& coordinates() = 0;

    // Helper function to convert a string containing a list of desired methods to
    // the correct <src>Vector<Int></src> required for the <src>setWinFitMethod</src> function.
    // This may be usful if your user interface involves strings rather than integers.
    // A new value is added to the output vector (which is resized appropriately) if any of the
    // substrings "window", "fit" or "interactive" (actually "win", "box" and
    // "inter" will do) is present.
    static Vector<Int> toMethodTypes (const String& methods);

    virtual IPosition getShape() const = 0;

    Bool shouldConvertToVelocity() const { return convertToVelocity_p; }

protected:

    // Constructor takes an image and a <src>LogIO</src> object for logging purposes.
    // You specify whether output images are  automatically overwritten if pre-existing,
    // or whether an intercative choice dialog widget appears (overWriteOutput=F)
    // You may also determine whether a progress meter is displayed or not.
    MomentsBase(
        LogIO &os, Bool overWriteOutput=False,
        Bool showProgress=True
    );

    LogIO os_p;
    Bool showProgress_p;
    Int momentAxisDefault_p {-10};
    T peakSNR_p {T(3)};
    T stdDeviation_p {T(0)};
    T yMin_p {T(0)};
    T yMax_p {T(0)};
    String out_p;
    String smoothOut_p {};
    Bool goodParameterStatus_p {True};
    Bool doWindow_p {False};
    Bool doFit_p {False};
    Bool doSmooth_p {False};
    Bool noInclude_p {True};
    Bool noExclude_p {True};
    Bool fixedYLimits_p {False};

    Int momentAxis_p {momentAxisDefault_p};
    Int worldMomentAxis_p;
    Vector<Int> kernelTypes_p {};
    Vector<Quantity> kernelWidths_p {};
    Vector<Int> moments_p {1, INTEGRATED};
    Vector<T> selectRange_p {};
    Vector<Int> smoothAxes_p {};
    Bool overWriteOutput_p;
    String error_p {};
    Bool convertToVelocity_p {False};
    MDoppler::Types velocityType_p {MDoppler::RADIO};

    // Check that the combination of methods that the user has requested is valid
    // List a handy dandy table if not.
    void _checkMethod();

    // Take the user's data inclusion and exclusion data ranges and
    // generate the range and Booleans to say what sort it is
    void _setIncludeExclude (
        Vector<T>& range, Bool& noInclude,
        Bool& noExclude, const Vector<T>& include,
        const Vector<T>& exclude
    );

    // Set the output image suffixes and units
    Bool _setOutThings(
        String& suffix, Unit& momentUnits,
        const Unit& imageUnits, const String& momentAxisUnits,
        const Int moment, Bool convertToVelocity
    );

    // Make output Coordinates
    CoordinateSystem _makeOutputCoordinates(
        IPosition& outShape, const CoordinateSystem& cSysIn,
        const IPosition& inShape, Int momentAxis,
        Bool removeAxis
    );
};

}

#ifndef AIPS_NO_TEMPLATE_SRC
#include <imageanalysis/ImageAnalysis/MomentsBase.tcc>
#endif

#endif