pmObjects.h

Go to the documentation of this file.

/** @file  pmObjects.h
 *
 * The process of finding, measuring, and classifying astronomical sources on
 * images is one of the critical tasks of the IPP or any astronomical software
 * system. This file will define structures and functions related to the task
 * of source detection and measurement. The elements defined in this section 
 * are generally low-level components which can be connected together to
 * construct a complete object measurement suite.
 *
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.2.8.1 $ $Name: rel9_1 $
 *  @date $Date: 2005/12/12 21:52:22 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#if !defined(PM_OBJECTS_H)
#define PM_OBJECTS_H

#if HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <math.h>
#include "pslib.h"
#include "pmAstrometry.h"
/**
 * In the object analysis process, we will use specific mask values to mark the
 * image pixels. The following structure defines the relevant mask values.
 *
 * XXX: This is probably a bad solution: we will want to set mask values
 * outside of the PSPHOT code.  Perhaps we can set up a registered set of mask
 * values with specific meanings that other functions can add to or define?
 */
typedef enum {
    PSPHOT_MASK_CLEAR     = 0x00,
    PSPHOT_MASK_INVALID   = 0x01,
    PSPHOT_MASK_SATURATED = 0x02,
    PSPHOT_MASK_MARKED    = 0x08,
} psphotMaskValues;


/** pmPeakType
 * 
 *  A peak pixel may have several features which may be determined when the
 *  peak is found or measured. These are specified by the pmPeakType enum.
 *  PM_PEAK_LONE represents a single pixel which is higher than its 8 immediate
 *  neighbors.  The PM_PEAK_EDGE represents a peak pixel which touching the image
 *  edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific
 *  number of neighbors at the same value, within some tolarence:
 * 
 */
typedef enum {
    PM_PEAK_LONE,                       ///< Isolated peak.
    PM_PEAK_EDGE,                       ///< Peak on edge.
    PM_PEAK_FLAT,                       ///< Peak has equal-value neighbors.
    PM_PEAK_UNDEF                       ///< Undefined.
} pmPeakType;


/** pmPeak data structure
 *  
 *  A source has the capacity for several types of measurements. The
 *  simplest measurement of a source is the location and flux of the peak pixel
 *  associated with the source:
 *  
 */
typedef struct
{
    int x;                              ///< X-coordinate of peak pixel.
    int y;                              ///< Y-coordinate of peak pixel.
    float counts;                       ///< Value of peak pixel (above sky?).
    pmPeakType class;                   ///< Description of peak.
}
pmPeak;


/** pmMoments data structure
 *  
 * One of the simplest measurements which can be made quickly for an object
 * are the object moments. We specify a structure to carry the moment information
 * for a specific source:
 *  
 */
typedef struct
{
    float x;    ///< X-coord of centroid.
    float y;    ///< Y-coord of centroid.
    float Sx;    ///< x-second moment.
    float Sy;    ///< y-second moment.
    float Sxy;    ///< xy cross moment.
    float Sum;    ///< Pixel sum above sky (background).
    float Peak;    ///< Peak counts above sky.
    float Sky;    ///< Sky level (background).
    float SN;    ///< approx signal-to-noise
    int nPixels;   ///< Number of pixels used.
}
pmMoments;


/** pmPSFClump data structure
 * 
 * A collection of object moment measurements can be used to determine
 * approximate object classes. The key to this analysis is the location and
 * statistics (in the second-moment plane,
 *  
 */
typedef struct
{
    float X;
    float dX;
    float Y;
    float dY;
}
pmPSFClump;

typedef int pmModelType;
#define PS_MODEL_GAUSS 0
#define PS_MODEL_PGAUSS 1
#define PS_MODEL_QGAUSS 2
#define PS_MODEL_SGAUSS 3


/** pmModel data structure
 * 
 * Every source may have two types of models: a PSF model and a FLT (floating)
 * model. The PSF model represents the best fit of the image PSF to the specific
 * object. In this case, the PSF-dependent parameters are specified for the
 * object by the PSF, not by the fit. The FLT model represents the best fit of
 * the given model to the object, with all parameters floating in the fit.
 *  
 */
typedef struct
{
    pmModelType type;   ///< Model to be used.
    psVector *params;   ///< Paramater values.
    psVector *dparams;   ///< Parameter errors.
    float chisq;   ///< Fit chi-squared.
    int nDOF;    ///< number of degrees of freedom
    int nIter;    ///< number of iterations to reach min
    float radius;   ///< fit radius actually used
}
pmModel;

/** pmSourceType enumeration
 * 
 * A given source may be identified as most-likely to be one of several source
 * types. The pmSource entry pmSourceType defines the current best-guess for this
 * source.
 * 
 * XXX: The values given below are currently illustrative and will require
 * some modification as the source classification code is developed. (TBD)
 * 
 */
typedef enum {
    PM_SOURCE_DEFECT,                   ///< a cosmic-ray
    PM_SOURCE_SATURATED,                ///< random saturated pixels

    PM_SOURCE_SATSTAR,                  ///< a saturated star
    PM_SOURCE_PSFSTAR,                  ///< a PSF star
    PM_SOURCE_GOODSTAR,                 ///< a good-quality star

    PM_SOURCE_POOR_FIT_PSF,             ///< poor quality PSF fit
    PM_SOURCE_FAIL_FIT_PSF,             ///< failed to get a good PSF fit
    PM_SOURCE_FAINTSTAR,                ///< below S/N cutoff

    PM_SOURCE_GALAXY,                   ///< an extended object (galaxy)
    PM_SOURCE_FAINT_GALAXY,             ///< a galaxy below S/N cutoff
    PM_SOURCE_DROP_GALAXY,              ///< ?
    PM_SOURCE_FAIL_FIT_GAL,             ///< failed on the galaxy fit
    PM_SOURCE_POOR_FIT_GAL,             ///< poor quality galaxy fit

    PM_SOURCE_OTHER,                    ///< unidentified
} pmSourceType;

/** pmSource data structure
 *  
 *  This source has the capacity for several types of measurements. The
 *  simplest measurement of a source is the location and flux of the peak pixel
 *  associated with the source:
 *  
 */
typedef struct
{
    pmPeak *peak;   ///< Description of peak pixel.
    psImage *pixels;   ///< Rectangular region including object pixels.
    psImage *weight;   ///< Image variance.
    psImage *mask;   ///< Mask which marks pixels associated with objects.
    pmMoments *moments;   ///< Basic moments measure for the object.
    pmModel *modelPSF;   ///< PSF Model fit (parameters and type)
    pmModel *modelFLT;   ///< FLT (floating) Model fit (parameters and type).
    pmSourceType type;   ///< Best identification of object.
}
pmSource;


/** pmPeakAlloc()
 *
 *  @return pmPeak*    newly allocated pmPeak with all internal pointers set to NULL
 */
pmPeak *pmPeakAlloc(
    int x,    ///< Row-coordinate in image space
    int y,    ///< Col-coordinate in image space
    float counts,   ///< The value of the peak pixel
    pmPeakType class   ///< The type of peak pixel
);


/** pmMomentsAlloc()
 * 
 */
pmMoments *pmMomentsAlloc();


/** pmModelAlloc()
 * 
 */
pmModel *pmModelAlloc(pmModelType type);


/** pmSourceAlloc()
 * 
 */
pmSource  *pmSourceAlloc();


/** pmFindVectorPeaks()
 * 
 * Find all local peaks in the given vector above the given threshold. A peak
 * is defined as any element with a value greater than its two neighbors and with
 * a value above the threshold. Two types of special cases must be addressed.
 * Equal value elements: If an element has the same value as the following
 * element, it is not considered a peak. If an element has the same value as the
 * preceding element (but not the following), then it is considered a peak. Note
 * that this rule (arbitrarily) identifies flat regions by their trailing edge.
 * Edge cases: At start of the vector, the element must be higher than its
 * neighbor. At the end of the vector, the element must be higher or equal to its
 * neighbor. These two rules again places the peak associated with a flat region
 * which touches the image edge at the image edge. The result of this function is
 * a vector containing the coordinates (element number) of the detected peaks
 * (type psU32).
 * 
 */
psVector *pmFindVectorPeaks(
    const psVector *vector,  ///< The input vector (float)
    float threshold   ///< Threshold above which to find a peak
);


/** pmFindImagePeaks()
 * 
 * Find all local peaks in the given image above the given threshold. This
 * function should find all row peaks using pmFindVectorPeaks, then test each row
 * peak and exclude peaks which are not local peaks. A peak is a local peak if it
 * has a higher value than all 8 neighbors. If the peak has the same value as its
 * +y neighbor or +x neighbor, it is NOT a local peak. If any other neighbors
 * have an equal value, the peak is considered a valid peak. Note two points:
 * first, the +x neighbor condition is already enforced by pmFindVectorPeaks.
 * Second, these rules have the effect of making flat-topped regions have single
 * peaks at the (+x,+y) corner. When selecting the peaks, their type must also be
 * set. The result of this function is an array of pmPeak entries.
 * 
 */
psArray *pmFindImagePeaks(
    const psImage *image,  ///< The input image where peaks will be found (float)
    float threshold   ///< Threshold above which to find a peak
);


/** pmCullPeaks()
 * 
 * Eliminate peaks from the psList that have a peak value above the given
 * maximum, or fall outside the valid region.
 * 
 */
psList *pmCullPeaks(
    psList *peaks,   ///< The psList of peaks to be culled
    float maxValue,   ///< Cull peaks above this value
    const psRegion valid                ///< Cull peaks otside this psRegion
);


/** pmPeaksSubset()
 * 
 * Create a new peaks array, removing certain types of peaks from the input
 * array of peaks based on the given criteria. Peaks should be eliminated if they
 * have a peak value above the given maximum value limit or if the fall outside
 * the valid region.  The result of the function is a new array with a reduced
 * number of peaks.
 * 
 */
psArray *pmPeaksSubset(
    psArray *peaks,                     ///< Add comment.
    float maxvalue,                     ///< Add comment.
    const psRegion valid                ///< Add comment.
);


/** pmSourceDefinePixels()
 * 
 * Define psImage subarrays for the source located at coordinates x,y on the
 * image set defined by readout. The pixels defined by this operation consist of
 * a square window (of full width 2Radius+1) centered on the pixel which contains
 * the given coordinate, in the frame of the readout. The window is defined to
 * have limits which are valid within the boundary of the readout image, thus if
 * the radius would fall outside the image pixels, the subimage is truncated to
 * only consist of valid pixels. If readout->mask or readout->weight are not
 * NULL, matching subimages are defined for those images as well. This function
 * fails if no valid pixels can be defined (x or y less than Radius, for
 * example). This function should be used to define a region of interest around a
 * source, including both source and sky pixels.
 * 
 * XXX: must code this.
 * 
 */
// XXX: Uncommenting the pmReadout causes compile errors.
bool pmSourceDefinePixels(
    pmSource *mySource,                 ///< Add comment.
    pmReadout *readout,                 ///< Add comment.
    psF32 x,                            ///< Add comment.
    psF32 y,                            ///< Add comment.
    psF32 Radius                        ///< Add comment.
);


/** pmSourceLocalSky()
 * 
 * Measure the local sky in the vicinity of the given source. The Radius
 * defines the square aperture in which the moments will be measured. This
 * function assumes the source pixels have been defined, and that the value of
 * Radius here is smaller than the value of Radius used to define the pixels. The
 * annular region not contained within the radius defined here is used to measure
 * the local background in the vicinity of the source. The local background
 * measurement uses the specified statistic passed in via the statsOptions entry.
 * This function allocates the pmMoments structure. The resulting sky is used to
 * set the value of the pmMoments.sky element of the provided pmSource structure.
 * 
 */
bool pmSourceLocalSky(
    pmSource *source,   ///< The input image (float)
    psStatsOptions statsOptions, ///< The statistic used in calculating the background sky
    float Radius   ///< The inner radius of the square annulus to exclude
);


/** pmSourceMoments()
 * 
 * Measure source moments for the given source, using the value of
 * source.moments.sky provided as the local background value and the peak
 * coordinates as the initial source location. The resulting moment values are
 * applied to the source.moments entry, and the source is returned. The moments
 * are measured within the given circular radius of the source.peak coordinates.
 * The return value indicates the success (TRUE) of the operation.
 * 
 */
bool pmSourceMoments(
    pmSource *source,   ///< The input pmSource for which moments will be computed
    float radius   ///< Use a circle of pixels around the peak
);


/** pmSourcePSFClump()
 * 
 * We use the source moments to make an initial, approximate source
 * classification, and as part of the information needed to build a PSF model for
 * the image. As long as the PSF shape does not vary excessively across the
 * image, the sources which are represented by a PSF (the start) will have very
 * similar second moments. The function pmSourcePSFClump searches a collection of
 * sources with measured moments for a group with moments which are all very
 * similar. The function returns a pmPSFClump structure, representing the
 * centroid and size of the clump in the sigma_x, sigma_y second-moment plane.
 * 
 * The goal is to identify and characterize the stellar clump within the
 * sigma_x, sigma_y second-moment plane.  To do this, an image is constructed to
 * represent this plane.  The units of sigma_x and sigma_y are in image pixels. A
 * pixel in this analysis image represents 0.1 pixels in the input image. The
 * dimensions of the image need only be 10 pixels. The peak pixel in this image
 * (above a threshold of half of the image maximum) is found. The coordinates of
 * this peak pixel represent the 2D mode of the sigma_x, sigma_y distribution.
 * The sources with sigma_x, sigma_y within 0.2 pixels of this value are then
 *  * used to calculate the median and standard deviation of the sigma_x, sigma_y
 * values. These resulting values are returned via the pmPSFClump structure.
 * 
 * The return value indicates the success (TRUE) of the operation.
 * 
 * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD).
 * XXX: Save the clump parameters on the Metadata (TBD)
 * 
 */
pmPSFClump pmSourcePSFClump(
    psArray *source,   ///< The input pmSource
    psMetadata *metadata  ///< Contains classification parameters
);


/** pmSourceRoughClass()
 * 
 * Based on the specified data values, make a guess at the source
 * classification. The sources are provides as a psArray of pmSource entries.
 * Definable parameters needed to make the classification are provided to the
 * routine with the psMetadata structure. The rules (in SDRS) refer to values which
 * can be extracted from the metadata using the given keywords. Except as noted,
 * the data type for these parameters are psF32.
 * 
 */
bool pmSourceRoughClass(
    psArray *source,   ///< The input pmSource
    psMetadata *metadata,  ///< Contains classification parameters
    pmPSFClump clump   ///< Statistics about the PSF clump
);


/** pmSourceModelGuess()
 * 
 * Convert available data to an initial guess for the given model. This
 * function allocates a pmModel entry for the pmSource structure based on the
 * provided model selection. The method of defining the model parameter guesses
 * are specified for each model below. The guess values are placed in the model
 * parameters. The function returns TRUE on success or FALSE on failure.
 * 
 */
pmModel *pmSourceModelGuess(
    pmSource *source,   ///< The input pmSource
    pmModelType model   ///< The type of model to be created.
);


/** pmContourType
 * 
 * Only one type is defined at present.
 * 
 */
typedef enum {
    PS_CONTOUR_CRUDE,
    PS_CONTOUR_UNKNOWN01,
    PS_CONTOUR_UNKNOWN02
} pmContourType;


/** pmSourceContour()
 * 
 * Find points in a contour for the given source at the given level. If type
 * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak,
 * running along each pixel row until the level is crossed, then interpolating to
 * the level coordinate for that row. This is done for each row, with the
 * starting point determined by the midpoint of the previous row, until the
 * starting point has a value below the contour level. The returned contour
 * consists of two vectors giving the x and y coordinates of the contour levels.
 * This function may be used as part of the model guess inputs.  Other contour
 * types may be specified in the future for more refined contours (TBD)
 * 
 */
psArray *pmSourceContour(
    pmSource *source,   ///< The input pmSource
    const psImage *image,  ///< The input image (float) (this arg should be removed)
    float level,   ///< The level of the contour
    pmContourType mode   ///< Currently this must be PS_CONTOUR_CRUDE
);


/** pmSourceFitModel()
 * 
 * Fit the requested model to the specified source. The starting guess for the
 * model is given by the input source.model parameter values. The pixels of
 * interest are specified by the source.pixelsand source.maskentries. This
 * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
 * on success or FALSE on failure.
 * 
 */
bool pmSourceFitModel(
    pmSource *source,   ///< The input pmSource
    pmModel *model,   ///< model to be fitted
    const bool PSF   ///< Treat model as PSF or FLT?
);


/** pmModelFitStatus()
 * 
 * This function wraps the call to the model-specific function returned by
 * pmModelFitStatusFunc_GetFunction.  The model-specific function examines the
 * model parameters, parameter errors, Chisq, S/N, and other parameters available
 * from model to decide if the particular fit was successful or not.
 * 
 * XXX: Must code this.
 * 
 */
bool pmModelFitStatus(
    pmModel *model                      ///< Add comment.
);


/** pmSourceAddModel()
 * 
 * Add the given source model flux to/from the provided image. The boolean
 * option center selects if the source is re-centered to the image center or if
 * it is placed at its centroid location. The boolean option sky selects if the
 * background sky is applied (TRUE) or not. The pixel range in the target image
 * is at most the pixel range specified by the source.pixels image. The success
 * status is returned.
 * 
 */
bool pmSourceAddModel(
    psImage *image,   ///< The output image (float)
    psImage *mask,   ///< The image pixel mask (valid == 0)
    pmModel *model,   ///< The input pmModel
    bool center    ///< A boolean flag that determines whether pixels are centered
);


/** pmSourceSubModel()
 * 
 * Subtract the given source model flux to/from the provided image. The
 * boolean option center selects if the source is re-centered to the image center
 * or if it is placed at its centroid location. The boolean option sky selects if
 * the background sky is applied (TRUE) or not. The pixel range in the target
 * image is at most the pixel range specified by the source.pixels image. The
 * success status is returned.
 * 
 */
bool pmSourceSubModel(
    psImage *image,   ///< The output image (float)
    psImage *mask,   ///< The image pixel mask (valid == 0)
    pmModel *model,   ///< The input pmModel
    bool center    ///< A boolean flag that determines whether pixels are centered
);


/**
 * 
 * The function returns both the magnitude of the fit, defined as -2.5log(flux),
 * where the flux is integrated under the model, theoretically from a radius of 0
 * to infinity. In practice, we integrate the model beyond 50sigma.  The aperture magnitude is
 * defined as -2.5log(flux) , where the flux is summed for all pixels which are
 * not excluded by the aperture mask. The model flux is calculated by calling the
 * model-specific function provided by pmModelFlux_GetFunction.
 * 
 * XXX: must code this.
 * 
 */
bool pmSourcePhotometry(
    float *fitMag,                      ///< integrated fit magnitude
    float *obsMag,   ///< aperture flux magnitude
    pmModel *model,                     ///< model used for photometry
    psImage *image,                     ///< image pixels to be used
    psImage *mask                       ///< mask of pixels to ignore
);


/**
 * 
 * This function converts the source classification into the closest available
 * approximation to the Dophot classification scheme:
 * 
 * PM_SOURCE_DEFECT: 8
 * PM_SOURCE_SATURATED: 8
 * PM_SOURCE_SATSTAR: 10
 * PM_SOURCE_PSFSTAR: 1
 * PM_SOURCE_GOODSTAR: 1
 * PM_SOURCE_POOR_FIT_PSF: 7
 * PM_SOURCE_FAIL_FIT_PSF: 4
 * PM_SOURCE_FAINTSTAR: 4
 * PM_SOURCE_GALAXY: 2
 * PM_SOURCE_FAINT_GALAXY: 2
 * PM_SOURCE_DROP_GALAXY: 2
 * PM_SOURCE_FAIL_FIT_GAL: 2
 * PM_SOURCE_POOR_FIT_GAL: 2
 * PM_SOURCE_OTHER: ?
 * 
 */
int pmSourceDophotType(
    pmSource *source                    ///< Add comment.
);


/** pmSourceSextractType()
 * 
 * This function converts the source classification into the closest available
 * approximation to the Sextractor classification scheme. the correspondence is
 * not yet defined (TBD) .
 * 
 * XXX: Must code this.
 * 
 */
int pmSourceSextractType(
    pmSource *source                    ///< Add comment.
);

/** pmSourceFitModel_v5()
 * 
 * Fit the requested model to the specified source. The starting guess for the
 * model is given by the input source.model parameter values. The pixels of
 * interest are specified by the source.pixelsand source.maskentries. This
 * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
 * on success or FALSE on failure.
 * 
 */
bool pmSourceFitModel_v5(
    pmSource *source,   ///< The input pmSource
    pmModel *model,   ///< model to be fitted
    const bool PSF   ///< Treat model as PSF or FLT?
);


/** pmSourceFitModel_v7()
 * 
 * Fit the requested model to the specified source. The starting guess for the
 * model is given by the input source.model parameter values. The pixels of
 * interest are specified by the source.pixelsand source.maskentries. This
 * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
 * on success or FALSE on failure.
 * 
 */
bool pmSourceFitModel_v7(
    pmSource *source,   ///< The input pmSource
    pmModel *model,   ///< model to be fitted
    const bool PSF   ///< Treat model as PSF or FLT?
);


/** pmSourcePhotometry()
 * 
 * XXX: Need descriptions
 * 
 */
bool pmSourcePhotometry(
    float *fitMag,
    float *obsMag,
    pmModel *model,
    psImage *image,
    psImage *mask
);

/** pmModelEval()
 *
 *  XXX: Need descriptions
 *
 */
psF32 pmModelEval(
    pmModel *model,
    psImage *image,
    psS32 col,
    psS32 row
);

#endif

---