ImageMoments.h

Classes

ImageMoments -- This class generates moments from an image. (full description)

template <class T> class ImageMoments

Types

enum MomentTypes

AVERAGE

The average intensity

INTEGRATED

The integrated intensity

WEIGHTED_MEAN_COORDINATE

The intensity weighted mean coordinate (usually velocity)

WEIGHTED_DISPERSION_COORDINATE

The intensity weighted coordinate (usually velocity) dispersion

MEDIAN

The median intensity

MEDIAN_COORDINATE

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.

STANDARD_DEVIATION

The standard deviation about the mean of the intensity

RMS

The rms of the intensity

ABS_MEAN_DEVIATION

The absolute mean deviation of the intensity

MAXIMUM

The maximum value of the intensity

MAXIMUM_COORDINATE

The coordinate (usually velocity) of the maximum value of the intensity

MINIMUM

The minimum value of the intensity

MINIMUM_COORDINATE

The coordinate (usually velocity) of the minimum value of the intensity

NMOMENTS

Total number

enum MethodTypes

WINDOW

Invokes the spectral windowing method

FIT

Invokes Gaussian fitting

INTERACTIVE

Invokes interactive methods

NMETHODS

Interface

Public Members

ImageMoments (ImageInterface<T>& image, T &os, Bool overWriteOutput=False, Bool showProgress=True)

ImageMoments(const ImageMoments<T> &other)

ImageMoments(ImageMoments<T> &other)

~ImageMoments()

ImageMoments<T> &operator=(const ImageMoments<T> &other)

Bool setMoments (const Vector<Int>& moments)

Bool setMomentAxis (const Int& momentAxis)

Bool setWinFitMethod(const Vector<Int>& method)

Bool setSmoothMethod(const Vector<Int>& smoothAxes, const Vector<Int>& kernelTypes, const Vector<Quantum<Double> >& kernelWidths)

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

Bool setInExCludeRange(const Vector<T>& include, const Vector<T>& exclude)

Bool setSnr(const T& peakSNR, const T& stdDeviation)

Bool setSmoothOutName(const String& smOut)

Bool setPlotting(PGPlotter& device, const Vector<Int>& nxy, const Bool yInd=False)

void setVelocityType (MDoppler::Types type)

void closePlotting()

void resetError ()

String errorMessage() const

Bool createMoments(PtrBlock<MaskedLattice<T>* >& images, Bool doTemp, const String& outFileName, Bool removeAxes=True)

Bool setNewImage (ImageInterface<T>& image)

CoordinateSystem coordinates() const

static Vector<Int> toMethodTypes (const String& methods)

Private Members

Bool checkMethod()

static void drawHistogram (const Vector<T>& values, const Vector<T>& counts, PGPlotter& plotter)

static void drawLine (const Vector<T>& x, const Vector<T>& y, PGPlotter& plotter)

static void drawVertical (const T& x, const T& yMin, const T& yMax, PGPlotter& plotter)

static Bool getLoc (T& x, T& y, PGPlotter& plotter)

static Float convertT (const T value)

static T convertF (const Float value)

static Bool readCursor (PGPlotter& plotter, Float& x, Float& y, String& ch)

Bool setIncludeExclude (Vector<T>& range, Bool& noInclude, Bool& noExclude, const Vector<T>& include, const Vector<T>& exclude, ostream& os)

Bool setOutThings(String& suffix, Unit& momentUnits, const Unit& imageUnits, const String& momentAxisUnits, const Int moment, Bool convertToVelocity)

Bool smoothImage (PtrHolder<ImageInterface<T> >& pSmoothedImage, T& smoothName)

Bool whatIsTheNoise (T& noise, ImageInterface<T>& image)

CoordinateSystem makeOutputCoordinates (IPosition& outShape, const CoordinateSystem& cSysIn, const IPosition& inShape, Int momentAxis, Bool removeAxis)

Description

Review Status

Date Reviewed:

yyyy/mm/dd

Prerequisite

• ImageInterface
• LatticeApply
• MomentCalcBase

Etymology

Synopsis

This class only accepts images of type Float and Double. This restriction is because of the plotting capabilities which are a bit awkward for other types.

This class makes a distinction between a "moment" and a "method". This boundary is a little blurred, but it claims to refer to the distinction of what you are computing, as to how the pixels that were included in the computation were selected. For example, a "moment" would be the average value of some pixels. A method for selecting those pixels would be a simple range specifying a range for which pixels should be included.

The default state of this class is to do nothing. If you specify an image via the setImage function then invoking the createMoments function will cause it to compute the integrated itensity along the spectral axis of the image (if it can find one). You can change the computational state of the class from this basic form via the remaining set functions. You can call any number of these functions to suit your needs.

Because there are a wide variety of methods available, if you specify an invalid combination, a table showing the available methods is listed. It is reproduced below for convenience.

The basic method is to just compute moments directly from the pixel intensities. This can be modified by applying pixel intensity inclusion or exclusion ranges. You can then also smooth the image and compute a mask based on the inclusion or exclusion ranges applied to the smoothed image. This mask is then applied to the unsmoothed data for moment computation.

The window method does no pixel intensity range selection. Instead a spectral window is found (hopefully surrounding the spectral line feature) and only the pixels in that window are used for computation. This window can be found from the smoothed or unsmoothed data. The moments are always computed from the unsmoothed data. For each spectrum, the window can be found interactively or automatically. There are two interactive methods. Either you just mark the window with the cursor, or you interactively fit a Gaussian to the profile and the +/- 3-sigma window is returned. There are two automatic methods. Either Bosma's converging mean algorithm is used, or an automatically fit Gaussian +/- 3-sigma window is returned.

The fitting method fits Gaussians to spectral features either automatically or interactively. The moments are then computed from the Gaussian fits (not the data themselves).

When an output image is created, it will have N-1 axes, where the input image has N axes. In the output image, the physical axis corresponding to the moment axis will have been removed, but the coordinate information will be retained for future coordinate transformations. For example, if you have a RA-DEC-VELOCITY image and you collapsed axis 2 (the DEC axis) the output images would be RA-VELOCITY with the coordinate information retained for the DEC axis so that the coupled nature of RA/DEC coordinates is preserved.

Output images are created with an all True (good) mask. If, for a given pixel, the moment calculation fails, then the mask is set to F.

When making plots, the order in which the spectra are displayed is determined by the tiling sequence of the image (for optimum speed of access).

                   Allowed Methods
                   ---------------

   Smooth    Window      Fit   in/exclude   Interactive 
   -----------------------------------------------------
     N          N         N        N            N       
     Y/N        N         N        Y            N       
 
     Y/N        Y         N        N            Y       
     Y/N        Y         N        N            N       
     Y/N        Y         Y        N            Y/N     

     N          N         Y        N            Y/N     
   ----------------------------------------------------

Note that the MEDIAN_COORDINATE moment is not very robust. It is very useful for generating quickly a velocity field in a way that is not sensitive to noise. However, it will only give sensible results under certain conditions. It treats the spectrum as a probability distribution, generates the cumulative distribution for the selected pixels (via an include or exclude pixel range, and finds the (interpolated) coordinate coresponding to the 50% value. The generation of the cumulative distribution and the finding of the 50% level really only makes sense if the cumulative distribution is monotonically increasing. This essentially means only selecting pixels which are positive or negative. For this reason, this moment type is *only* supported with the basic method (i.e. no smoothing, no windowing, no fitting) with a pixel selection range that is either all positive, or all negative.

If you ignore return error statuses from the functions that set the state of the class, the internal status of the class is set to bad. This means it will just keep on returning error conditions until you explicitly recover the situation. A message describing the last error condition can be recovered with function errorMessage.

Example

    // Set state function argument values
    
         Vector<Int> moments(2);
         moments(0) = ImageMoments<Float>::AVERAGE;
         moments(1) = ImageMoments<Float>::WEIGHTED_MEAN_COORDINATE;
         Vector<int> methods(2);
         methods(0) = ImageMoments<Float>::WINDOW;
         methods(1) = ImageMoments<Float>::INTERACTIVE;
         Vector<Int> nxy(2);
         nxy(0) = 3;
         nxy(1) = 3;
    
    // Open paged image
        
         PagedImage<Float> inImage(inName);  
    
    // Construct moment helper object
    
         LogOrigin or("myClass", "myFunction(...)", WHERE);
         LogIO os(or);
         ImageMoments<Float> moment(SubImage<Float>(inName), os);
    
    // Specify state via control functions
    
         if (!moment.setMoments(moments)) return 1;
         if (!moment.setWinFitMethod(methods)) return 1;
         if (!moment.setMomentAxis(3)) return 1;
         if (!moment.setPlotting("/xs", nxy)) return 1;
    
    // Create the moments
    
         if (!moment.createMoments()) return 1;
    

Motivation

To Do

• more control over histogram of image noise at start (pixel range and number of bins)
• better algorithm for seeing if spectrum is pure noise
• Make this class extensible so users could add their own method.

Member Description

ImageMoments (ImageInterface<T>& image, T &os, Bool overWriteOutput=False, Bool showProgress=True)

Constructor takes an image and a LogIO 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.

ImageMoments(const ImageMoments<T> &other)

Copy constructor. Uses copy semantics.

ImageMoments(ImageMoments<T> &other)

Copy constructor. Uses copy semantics.

~ImageMoments()

Destructor

ImageMoments<T> &operator=(const ImageMoments<T> &other)

Assignment operator. USes copy semantics.

enum MomentTypes

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

Bool setMoments (const Vector<Int>& moments)

Bool setMomentAxis (const Int& momentAxis)

Set the moment axis (0 relative). A return value of False 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.

enum MethodTypes

The enum MethodTypes is provided for use with the setWinFitMethod 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.

Bool setWinFitMethod(const Vector<Int>& method)

The method by which you compute the moments is specified by calling (or not calling) the setWinFitMethod and setSmoothMethod functions. The default state of the class is to compute directly on all (or some according to setInExClude) of the pixels in the spectrum. Calling these functions modifies the computational state to something more complicated.

The setWinMethod function requires an Int array as its argument. Each Int specifies a different method that you can invoke (either separately or in combination). The allowed values and their meanings are given by the enum MethodTypes.

Both the windowing and fitting methods have interactive modes. The windowing method also has a fitting flavour, so if you set both ImageMoments::WINDOW and ImageMoments::FIT, you would be invoking the windowing method but determining the window by fitting Gaussians automatically (as ImageMoments::INTERACTIVE) was not given.

If you don't call this function, then neither the windowing nor fitting methods are invoked. A return value of False indicates that you asked for an illegal method.

Bool setSmoothMethod(const Vector<Int>& smoothAxes, const Vector<Int>& kernelTypes, const Vector<Quantum<Double> >& kernelWidths)
Bool setSmoothMethod(const Vector<Int>& smoothAxes, const Vector<Int>& kernelTypes, const Vector<Double>& kernelWidths)

This function invokes smoothing of the input image. Give Int arrays for the axes (0 relative) to be smoothed and the smoothing kernel types (use the enum KernelTypes) for each axis. Give a Double 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 False 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

Bool setInExCludeRange(const Vector<T>& include, const Vector<T>& exclude)

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 include and exclude must therefore be a zero length vector if you call this function. A return value of False indicates that you have given both an include and an exclude range. If you don't call this function, the default state of the class is to include all pixels.

Bool setSnr(const T& peakSNR, const T& stdDeviation)

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. This is only needed for the automatic window or fit methods. If you are using an interactive nethod, you assess yourself whether the spectrum contains signal or not.

peakSNR is the signal-to-noise ratio of the peak value in the spectrum below which the spectrum is considered pure noise. stdDeviation 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 peakSNR=3 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). If you have specified a plotting device (see setPlotting) then you get to interact with the fitting procedure if you want to. A return value of False indicates you have set invalid values.

Bool setSmoothOutName(const String& smOut)

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 include or exclude range correct for the smooth-clip method. The default state of the class is to not output the smoothed image.

Bool setPlotting(PGPlotter& device, const Vector<Int>& nxy, const Bool yInd=False)

This sets the name of the PGPLOT plotting device, the number of subplots in x and y per page and whether each spectrum plot is autoscaled individually (yInd=False) or they are plotted with the same range automatically determined from the image. Plotting is not invoked for all states of the class. It is only needed for the interactive methods. If you ask for a method that needs to determine the noise from the image, and you set the plotting device, then this will be done interactively. Similarly, if you invoke the automatic window or fit methods, but set the plotting device, then you will see plots of the spectra and the selected windows and fits, respectively.

The default state of the class is that no plotting characteristics are set. However, if you set device but offer a zero length array for nxy then the latter is set to [1,1]. A return value of False indicates that you gave roo many values in the nxy vector.

void setVelocityType (MDoppler::Types type)

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 closePlotting()

CLose plotter

void resetError ()

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 createMoments function 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.

String errorMessage() const

Bool createMoments(PtrBlock<MaskedLattice<T>* >& images, Bool doTemp, const String& outFileName, Bool removeAxes=True)

The output PtrBlock will hold PagedImages or TempImages (doTemp=True). It is your responsibility to delete the pointers. If doTemp is True, the outFileName is irrelevant.

If you create PagedImages, outFileName is the root name for the output files. Suffixes will be made up internally to append to this root. If you only ask for one moment, this will be the actual name of the output file. If you don't set this variable, the default state of the class is to set the output name root to the name of the input file.

Bool setNewImage (ImageInterface<T>& image)

Set a new image. A return value of False indicates the image had an invalid type (this class only accepts Float or Double images).

CoordinateSystem coordinates() const

Get CoordinateSystem

static Vector<Int> toMethodTypes (const String& methods)

Bool checkMethod()

Check that the combination of methods that the user has requested is valid List a handy dandy table if not.

static void drawHistogram (const Vector<T>& values, const Vector<T>& counts, PGPlotter& plotter)

Plot a histogram

static void drawLine (const Vector<T>& x, const Vector<T>& y, PGPlotter& plotter)

Plot a line

static void drawVertical (const T& x, const T& yMin, const T& yMax, PGPlotter& plotter)

Draw a vertical line of the given length at a given abcissa

static Bool getLoc (T& x, T& y, PGPlotter& plotter)

Read the cursor and return its coordinates

static Float convertT (const T value)

Convert a to a for plotting

static T convertF (const Float value)