Index: /branches/rel10_ifa/psModules/src/objects/pmModel.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmModel.c (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmModel.c (revision 6537) @@ -0,0 +1,216 @@ +/** @file pmObjects.c + * + * This file will ... + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmObjects.h" +#include "pmModelGroup.h" + +static void modelFree(pmModel *tmp) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + psFree(tmp->params); + psFree(tmp->dparams); + psTrace(__func__, 4, "---- %s() end ----\n", __func__); +} + +/****************************************************************************** +pmModelAlloc(): Allocate the pmModel structure, along with its parameters, +and initialize the type member. Initialize the params to 0.0. +XXX EAM: simplifying code with pmModelParameterCount +*****************************************************************************/ +pmModel *pmModelAlloc(pmModelType type) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel)); + + tmp->type = type; + tmp->chisq = 0.0; + tmp->nIter = 0; + tmp->radius = 0; + tmp->status = PM_MODEL_UNTRIED; + + psS32 Nparams = pmModelParameterCount(type); + if (Nparams == 0) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return(NULL); + } + + tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32); + + for (psS32 i = 0; i < tmp->params->n; i++) { + tmp->params->data.F32[i] = 0.0; + tmp->dparams->data.F32[i] = 0.0; + } + + psMemSetDeallocator(tmp, (psFreeFunc) modelFree); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + +// copy model to a new structure +pmModel *pmModelCopy (pmModel *model) +{ + + pmModel *new = pmModelAlloc (model->type); + + new->chisq = model->chisq; + new->nDOF = model->nDOF; + new->nIter = model->nIter; + new->status = model->status; + new->radius = model->radius; + + for (int i = 0; i < new->params->n; i++) { + new->params->data.F32[i] = model->params->data.F32[i]; + new->dparams->data.F32[i] = model->dparams->data.F32[i]; + } + + return (new); +} + +/****************************************************************************** + pmModelEval(source, level, row): evaluates the model function at the specified coords. + + NOTE: The coords are in subImage source->pixel coords, not image coords. + + XXX: Use static vectors for x (NO: needs to be thread safe) +*****************************************************************************/ +psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(image, false); + PS_ASSERT_PTR_NON_NULL(model, false); + PS_ASSERT_PTR_NON_NULL(model->params, false); + + // Allocate the x coordinate structure and convert row/col to image space. + // + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + x->data.F32[0] = (psF32) (col + image->col0); + x->data.F32[1] = (psF32) (row + image->row0); + psF32 tmpF; + pmModelFunc modelFunc; + + modelFunc = pmModelFunc_GetFunction (model->type); + tmpF = modelFunc (NULL, model->params, x); + psFree(x); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpF); +} + +bool AddOrSubModel(psImage *image, + psImage *mask, + pmModel *model, + bool center, + bool sky, + bool add + ) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + PS_ASSERT_PTR_NON_NULL(model, false); + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); + + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + psVector *params = model->params; + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + psS32 imageCol; + psS32 imageRow; + psF32 skyValue = params->data.F32[0]; + psF32 pixelValue; + + for (psS32 i = 0; i < image->numRows; i++) { + for (psS32 j = 0; j < image->numCols; j++) { + if ((mask != NULL) && mask->data.U8[i][j]) + continue; + + // Convert i/j to image coord space: + // 'center' option changes meaning of i,j: + if (center) { + imageCol = j - 0.5*image->numCols + model->params->data.F32[2]; + imageRow = i - 0.5*image->numRows + model->params->data.F32[3]; + } else { + imageCol = j + image->col0; + imageRow = i + image->row0; + } + + x->data.F32[0] = (float) imageCol; + x->data.F32[1] = (float) imageRow; + + // set the appropriate pixel value for this coordinate + if (sky) { + pixelValue = modelFunc (NULL, params, x); + } else { + pixelValue = modelFunc (NULL, params, x) - skyValue; + } + + + // add or subtract the value + if (add + ) { + image->data.F32[i][j] += pixelValue; + } + else { + image->data.F32[i][j] -= pixelValue; + } + } + } + psFree(x); + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); +} + +/****************************************************************************** + *****************************************************************************/ +bool pmModelAdd(psImage *image, + psImage *mask, + pmModel *model, + bool center, + bool sky) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psBool rc = AddOrSubModel(image, mask, model, center, sky, true); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +/****************************************************************************** + *****************************************************************************/ +bool pmModelSub(psImage *image, + psImage *mask, + pmModel *model, + bool center, + bool sky) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psBool rc = AddOrSubModel(image, mask, model, center, sky, false); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +// check for success of model fit +bool pmModelFitStatus (pmModel *model) +{ + + bool status; + + pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type); + status = statusFunc (model); + + return (status); +} + Index: /branches/rel10_ifa/psModules/src/objects/pmModel.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmModel.h (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmModel.h (revision 6537) @@ -0,0 +1,116 @@ +/** @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.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_MODEL_H +# define PM_MODEL_H + +// type of model carried by the pmModel structure +typedef int pmModelType; + +typedef enum { + PM_MODEL_UNTRIED, ///< model fit not yet attempted + PM_MODEL_SUCCESS, ///< model fit succeeded + PM_MODEL_NONCONVERGE, ///< model fit did not converge + PM_MODEL_OFFIMAGE, ///< model fit drove out of range + PM_MODEL_BADARGS ///< model fit called with invalid args +} pmModelStatus; + +/** pmModel data structure + * + * Every source may have two types of models: a PSF model and a EXT (extended-source) + * 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 EXT model represents the best fit of + * the given model to the object, with all shape 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 + pmModelStatus status; ///< fit status + float radius; ///< fit radius actually used +} +pmModel; + +/** pmModelAlloc() + * + */ +pmModel *pmModelAlloc(pmModelType type); + +// copy model to a new structure +pmModel *pmModelCopy (pmModel *model); + +psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) + +/** pmModelAdd() + * + * 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 pmModelAdd( + 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 + bool sky ///< A boolean flag that determines if the sky is subtracted +); + + +/** pmModelSub() + * + * 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 pmModelSub( + 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 + bool sky ///< A boolean flag that determines if the sky is subtracted +); + +/** 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 ///< Model to be used +); + +# endif /* PM_MODEL_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmMoments.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmMoments.c (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmMoments.c (revision 6537) @@ -0,0 +1,218 @@ +/** @file pmObjects.c + * + * This file will ... + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmObjects.h" +#include "pmModelGroup.h" + +/****************************************************************************** +pmMomentsAlloc(): Allocate the pmMoments structure and initialize the members +to zero. +*****************************************************************************/ +pmMoments *pmMomentsAlloc() +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments)); + tmp->x = 0.0; + tmp->y = 0.0; + tmp->Sx = 0.0; + tmp->Sy = 0.0; + tmp->Sxy = 0.0; + tmp->Sum = 0.0; + tmp->Peak = 0.0; + tmp->Sky = 0.0; + tmp->nPixels = 0; + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + +/****************************************************************************** +pmSourceMoments(source, radius): this function takes a subImage defined in the +pmSource data structure, along with the peak location, and determines the +various moments associated with that peak. + +Requires the following to have been created: + pmSource + pmSource->peak + pmSource->pixels + pmSource->weight + pmSource->mask + +XXX: The peak calculations are done in image coords, not subImage coords. + +XXX EAM : this version clips input pixels on S/N +XXX EAM : this version returns false for several reasons +*****************************************************************************/ +# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) + +bool pmSourceMoments(pmSource *source, + psF32 radius) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); + + // + // XXX: Verify the setting for sky if source->moments == NULL. + // + psF32 sky = 0.0; + if (source->moments == NULL) { + source->moments = pmMomentsAlloc(); + } else { + sky = source->moments->Sky; + } + + // + // Sum = SUM (z - sky) + // X1 = SUM (x - xc)*(z - sky) + // X2 = SUM (x - xc)^2 * (z - sky) + // XY = SUM (x - xc)*(y - yc)*(z - sky) + // + psF32 peakPixel = -PS_MAX_F32; + psS32 numPixels = 0; + psF32 Sum = 0.0; + psF32 Var = 0.0; + psF32 X1 = 0.0; + psF32 Y1 = 0.0; + psF32 X2 = 0.0; + psF32 Y2 = 0.0; + psF32 XY = 0.0; + psF32 x = 0; + psF32 y = 0; + psF32 R2 = PS_SQR(radius); + + psF32 xPeak = source->peak->x; + psF32 yPeak = source->peak->y; + psF32 xOff = source->pixels->col0 - source->peak->x; + psF32 yOff = source->pixels->row0 - source->peak->y; + + // XXX why do I get different results for these two methods of finding Sx? + // XXX Sx, Sy would be better measured if we clip pixels close to sky + // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? + // We loop through all pixels in this subimage (source->pixels), and for each + // pixel that is not masked, AND within the radius of the peak pixel, we + // proceed with the moments calculation. need to do two loops for a + // numerically stable result. first loop: get the sums. + // XXX EAM : mask == 0 is valid + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 *vPix = source->pixels->data.F32[row]; + psF32 *vWgt = source->weight->data.F32[row]; + psU8 *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row]; + + for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { + if ((vMsk != NULL) && *vMsk) { + vMsk++; + continue; + } + + psF32 xDiff = col + xOff; + psF32 yDiff = row + yOff; + + // radius is just a function of (xDiff, yDiff) + if (!VALID_RADIUS(xDiff, yDiff, R2)) { + if (vMsk != NULL) + vMsk++; + continue; + } + + psF32 pDiff = *vPix - sky; + psF32 wDiff = *vWgt; + + // XXX EAM : check for valid S/N in pixel + // XXX EAM : should this limit be user-defined? + if (PS_SQR(pDiff) < wDiff) { + if (vMsk != NULL) + vMsk++; + continue; + } + + Var += wDiff; + Sum += pDiff; + + psF32 xWght = xDiff * pDiff; + psF32 yWght = yDiff * pDiff; + + X1 += xWght; + Y1 += yWght; + + XY += xDiff * yWght; + X2 += xDiff * xWght; + Y2 += yDiff * yWght; + + peakPixel = PS_MAX (*vPix, peakPixel); + numPixels++; + if (vMsk != NULL) + vMsk++; + } + } + + // if we have less than (1/4) of the possible pixels, force a retry + // XXX EAM - the limit is a bit arbitrary. make it user defined? + if ((numPixels < 0.75*R2) || (Sum <= 0)) { + psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return (false); + } + + psTrace (".psModules.pmSourceMoments", 5, + "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", + sky, Sum, X1, Y1, X2, Y2, XY, numPixels); + + // + // first moment X = X1/Sum + xc + // second moment X = sqrt (X2/Sum - (X1/Sum)^2) + // Sxy = XY / Sum + // + x = X1/Sum; + y = Y1/Sum; + if ((fabs(x) > radius) || (fabs(y) > radius)) { + psTrace (".psModules.pmSourceMoments", 3, + "large centroid swing; invalid peak %d, %d\n", + source->peak->x, source->peak->y); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return (false); + } + + source->moments->x = x + xPeak; + source->moments->y = y + yPeak; + + // XXX EAM : Sxy needs to have x*y subtracted + source->moments->Sxy = XY/Sum - x*y; + source->moments->Sum = Sum; + source->moments->SN = Sum / sqrt(Var); + source->moments->Peak = peakPixel; + source->moments->nPixels = numPixels; + + // XXX EAM : these values can be negative, so we need to limit the range + source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); + source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); + + psTrace (".psModules.pmSourceMoments", 4, + "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", + sky, Sum, source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, source->moments->Sxy); + + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); +} + Index: /branches/rel10_ifa/psModules/src/objects/pmMoments.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmMoments.h (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmMoments.h (revision 6537) @@ -0,0 +1,60 @@ +/** @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.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +/** 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 dSky; ///< local Sky variance + float SN; ///< approx signal-to-noise + int nPixels; ///< Number of pixels used. +} +pmMoments; + +/** pmMomentsAlloc() + * + */ +pmMoments *pmMomentsAlloc(); + +/** 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 +); Index: /branches/rel10_ifa/psModules/src/objects/pmPeaks.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPeaks.c (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmPeaks.c (revision 6537) @@ -0,0 +1,544 @@ +/** @file pmObjects.c + * + * This file will ... + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmObjects.h" +#include "pmModelGroup.h" + +/****************************************************************************** +myListAddPeak(): A private function which allocates a psArray, if the list +argument is NULL, otherwise it adds the peak to that list. +XXX EAM : changed the output to psArray +XXX EAM : Switched row, col args +XXX EAM : NOTE: this was changed in the call, so the new code is consistent +*****************************************************************************/ +static psArray *myListAddPeak(psArray *list, + psS32 row, + psS32 col, + psF32 counts, + pmPeakType type) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); + + if (list == NULL) { + list = psArrayAlloc(100); + list->n = 0; + } + psArrayAdd(list, 100, tmpPeak); + psFree (tmpPeak); + // XXX EAM : is this free appropriate? (does psArrayAdd increment memory counter?) + + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return(list); +} + +/****************************************************************************** +getRowVectorFromImage(): a private function which simply returns a +psVector containing the specified row of data from the psImage. + +XXX: Is there a better way to do this? +XXX EAM: does this really need to alloc a new vector??? +*****************************************************************************/ +static psVector *getRowVectorFromImage(psImage *image, + psU32 row) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); + + psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); + for (psU32 col = 0; col < image->numCols ; col++) { + tmpVector->data.F32[col] = image->data.F32[row][col]; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return(tmpVector); +} + +/****************************************************************************** +isItInThisRegion(): a private function which simply returns a +boolean denoting if specified coordinate is in the region. +XXX: Macro this. +*****************************************************************************/ +static bool isItInThisRegion(const psRegion valid, + psS32 x, + psS32 y) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + if ((x >= valid.x0) && + (x <= valid.x1) && + (y >= valid.y0) && + (y <= valid.y1)) { + psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); + return(true); + } + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return(false); +} + +/****************************************************************************** +pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members. +*****************************************************************************/ +pmPeak *pmPeakAlloc(psS32 x, + psS32 y, + psF32 counts, + pmPeakType type) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak)); + tmp->x = x; + tmp->y = y; + tmp->counts = counts; + tmp->type = type; + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + +// XXX are these private or public? +// psSort comparison function for peaks +int pmComparePeakAscend (const void **a, const void **b) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 diff; + + diff = A->counts - B->counts; + if (diff < FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); + return (-1); + } else if (diff > FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); + return (+1); + } + psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); + return (0); +} + +// psSort comparison function for peaks +int pmComparePeakDescend (const void **a, const void **b) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 diff; + + diff = A->counts - B->counts; + if (diff < FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); + return (+1); + } else if (diff > FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); + return (-1); + } + psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); + return (0); +} + +/****************************************************************************** +pmFindVectorPeaks(vector, threshold): Find all local peaks in the given vector +above the given threshold. Returns a vector of type PS_TYPE_U32 containing +the location (x value) of all peaks. + +XXX: What types should be supported? Only F32 is implemented. + +XXX: We currently step through the input vector twice; once to determine the +size of the output vector, then to set the values of the output vector. +Depending upon actual use, this may need to be optimized. +*****************************************************************************/ +psVector *pmFindVectorPeaks(const psVector *vector, + psF32 threshold) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_VECTOR_NON_NULL(vector, NULL); + PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL); + PS_ASSERT_VECTOR_TYPE(vector, PS_TYPE_F32, NULL); + int count = 0; + int n = vector->n; + + // + // Special case: the input vector has a single element. + // + if (n == 1) { + psVector *tmpVector = NULL; + ; + if (vector->data.F32[0] > threshold) { + tmpVector = psVectorAlloc(1, PS_TYPE_U32); + tmpVector->data.U32[0] = 0; + } else { + tmpVector = psVectorAlloc(0, PS_TYPE_U32); + } + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpVector); + } + + // + // Determine if first pixel is a peak + // + if ((vector->data.F32[0] > vector->data.F32[1]) && + (vector->data.F32[0] > threshold)) { + count++; + } + + // + // Determine if interior pixels are peaks + // + for (psU32 i = 1; i < n-1 ; i++) { + if ((vector->data.F32[i] > vector->data.F32[i-1]) && + (vector->data.F32[i] > vector->data.F32[i+1]) && + (vector->data.F32[i] > threshold)) { + count++; + } + } + + // + // Determine if last pixel is a peak + // + if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && + (vector->data.F32[n-1] > threshold)) { + count++; + } + + // + // We know how many peaks exist, so we now allocate a psVector to store + // those peaks. + // + psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32); + count = 0; + + // + // Determine if first pixel is a peak + // + if ((vector->data.F32[0] > vector->data.F32[1]) && + (vector->data.F32[0] > threshold)) { + tmpVector->data.U32[count++] = 0; + } + + // + // Determine if interior pixels are peaks + // + for (psU32 i = 1; i < (n-1) ; i++) { + if ((vector->data.F32[i] > vector->data.F32[i-1]) && + (vector->data.F32[i] > vector->data.F32[i+1]) && + (vector->data.F32[i] > threshold)) { + tmpVector->data.U32[count++] = i; + } + } + + // + // Determine if last pixel is a peak + // + if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && + (vector->data.F32[n-1] > threshold)) { + tmpVector->data.U32[count++] = n-1; + } + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpVector); +} + + +/****************************************************************************** +pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage +above the given threshold. Returns a psArray containing location (x/y value) +of all peaks. + +XXX: I'm not convinced the peak type definition in the SDRS is mutually +exclusive. Some peaks can have multiple types. Edges for sure. Also, a +digonal line with the same value at each point will have a peak for every +point on that line. + +XXX: This does not work if image has either a single row, or a single column. + +XXX: In the output psArray elements, should we use the image row/column offsets? + Currently, we do not. +XXX EAM : this function needs to return peaks in *parent* coords + +XXX: Merge with CVS 1.20. This had the proper code for images with a single +row or column. + +*****************************************************************************/ +psArray *pmFindImagePeaks(const psImage *image, + psF32 threshold) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); + if ((image->numRows == 1) || (image->numCols == 1)) { + psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns."); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + } + psVector *tmpRow = NULL; + psU32 col = 0; + psU32 row = 0; + psArray *list = NULL; + + psU32 col0 = image->col0; + psU32 row0 = image->row0; + + // + // Find peaks in row 0 only. + // + row = 0; + tmpRow = getRowVectorFromImage((psImage *) image, row); + psVector *row1 = pmFindVectorPeaks(tmpRow, threshold); + // pmFindVectorPeaks returns coords in the vector, not corrected for col0 + + for (psU32 i = 0 ; i < row1->n ; i++ ) { + col = row1->data.U32[i]; + // + // Determine if pixel (0,0) is a peak. + // + if (col == 0) { + if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + } else if (col < (image->numCols - 1)) { + if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && + (image->data.F32[row][col] > image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + + } else if (col == (image->numCols - 1)) { + if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) { + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + + } else { + psError(PS_ERR_UNKNOWN, true, "peak specified valid column range."); + } + } + psFree (tmpRow); + psFree (row1); + + // + // Exit if this image has a single row. + // + if (image->numRows == 1) { + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(list); + } + + // + // Find peaks in interior rows only. + // + for (row = 1 ; row < (image->numRows - 1) ; row++) { + tmpRow = getRowVectorFromImage((psImage *) image, row); + row1 = pmFindVectorPeaks(tmpRow, threshold); + + // Step through all local peaks in this row. + for (psU32 i = 0 ; i < row1->n ; i++ ) { + pmPeakType myType = PM_PEAK_UNDEF; + col = row1->data.U32[i]; + + if (col == 0) { + // If col==0, then we can not read col-1 pixels + if ((image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + myType = PM_PEAK_EDGE; + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); + } + } else if (col < (image->numCols - 1)) { + // This is an interior pixel + if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + if (image->data.F32[row][col] > threshold) { + if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] > image->data.F32[row][col+1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { + myType = PM_PEAK_LONE; + } + + if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) || + (image->data.F32[row][col] == image->data.F32[row-1][col]) || + (image->data.F32[row][col] == image->data.F32[row-1][col+1]) || + (image->data.F32[row][col] == image->data.F32[row][col-1]) || + (image->data.F32[row][col] == image->data.F32[row][col+1]) || + (image->data.F32[row][col] == image->data.F32[row+1][col-1]) || + (image->data.F32[row][col] == image->data.F32[row+1][col]) || + (image->data.F32[row][col] == image->data.F32[row+1][col+1])) { + myType = PM_PEAK_FLAT; + } + + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); + } + } + } else if (col == (image->numCols - 1)) { + // If col==numCols - 1, then we can not read col+1 pixels + if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col])) { + myType = PM_PEAK_EDGE; + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); + } + } else { + psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); + } + + } + psFree (tmpRow); + psFree (row1); + } + + // + // Find peaks in the last row only. + // + row = image->numRows - 1; + tmpRow = getRowVectorFromImage((psImage *) image, row); + row1 = pmFindVectorPeaks(tmpRow, threshold); + for (psU32 i = 0 ; i < row1->n ; i++ ) { + col = row1->data.U32[i]; + if (col == 0) { + if ( (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col+1])) { + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + } else if (col < (image->numCols - 1)) { + if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1])) { + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + + } else if (col == (image->numCols - 1)) { + if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] > image->data.F32[row][col-1])) { + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + } else { + psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); + } + } + psFree (tmpRow); + psFree (row1); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(list); +} + + +/****************************************************************************** +psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have +a peak value above the given maximum, or fall outside the valid region. + +XXX: Should the sky value be used when comparing the maximum? + +XXX: warning message if valid is NULL? + +XXX: changed API to create a NEW output psArray (should change name as well) + +XXX: Do we free the psList elements of those culled peaks? + +XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset? +*****************************************************************************/ +psList *pmCullPeaks(psList *peaks, + psF32 maxValue, + const psRegion valid) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(peaks, NULL); + + psListElem *tmpListElem = (psListElem *) peaks->head; + psS32 indexNum = 0; + + // printf("pmCullPeaks(): list size is %d\n", peaks->size); + while (tmpListElem != NULL) { + pmPeak *tmpPeak = (pmPeak *) tmpListElem->data; + if ((tmpPeak->counts > maxValue) || + (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))) { + psListRemoveData(peaks, (psPtr) tmpPeak); + } + + indexNum++; + tmpListElem = tmpListElem->next; + } + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(peaks); +} + +// XXX EAM: I changed this to return a new, subset array +// rather than alter the existing one +// XXX: Fix the *valid pointer. +psArray *pmPeaksSubset( + psArray *peaks, + psF32 maxValue, + const psRegion valid) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(peaks, NULL); + + psArray *output = psArrayAlloc (200); + output->n = 0; + + psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n); + + for (int i = 0; i < peaks->n; i++) { + pmPeak *tmpPeak = (pmPeak *) peaks->data[i]; + if (tmpPeak->counts > maxValue) + continue; + if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)) + continue; + psArrayAdd (output, 200, tmpPeak); + } + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(output); +} Index: /branches/rel10_ifa/psModules/src/objects/pmPeaks.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPeaks.h (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmPeaks.h (revision 6537) @@ -0,0 +1,140 @@ +/** @file pmPeaks.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.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_PEAKS_H +# define PM_PEAKS_H + +/** 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 type; ///< Description of peak. +} +pmPeak; + + +/** 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 type ///< The type of peak pixel +); + + +/** 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. +); + +# endif /* PM_PEAKS_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSource.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSource.c (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSource.c (revision 6537) @@ -0,0 +1,410 @@ +/** @file pmObjects.c + * + * This file will ... + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmObjects.h" +#include "pmModelGroup.h" + +static void sourceFree(pmSource *tmp) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + psFree(tmp->peak); + psFree(tmp->pixels); + psFree(tmp->weight); + psFree(tmp->mask); + psFree(tmp->moments); + psFree(tmp->modelPSF); + psFree(tmp->modelEXT); + psFree(tmp->blends); + psTrace(__func__, 4, "---- %s() end ----\n", __func__); +} + +/****************************************************************************** +pmSourceAlloc(): Allocate the pmSource structure and initialize its members +to NULL. +*****************************************************************************/ +pmSource *pmSourceAlloc() +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource)); + tmp->peak = NULL; + tmp->pixels = NULL; + tmp->weight = NULL; + tmp->mask = NULL; + tmp->moments = NULL; + tmp->blends = NULL; + tmp->modelPSF = NULL; + tmp->modelEXT = NULL; + tmp->type = PM_SOURCE_UNKNOWN; + tmp->mode = PM_SOURCE_DEFAULT; + psMemSetDeallocator(tmp, (psFreeFunc) sourceFree); + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + +// x,y are defined in the parent image coords of readout->image +bool pmSourceDefinePixels(pmSource *mySource, + const pmReadout *readout, + psF32 x, + psF32 y, + psF32 Radius) +{ + psRegion srcRegion; + + // Grab a subimage of the original image of size (2 * outerRadius). + srcRegion = psRegionForSquare (x, y, Radius); + srcRegion = psRegionForImage (readout->image, srcRegion); + + mySource->pixels = psImageSubset(readout->image, srcRegion); + mySource->weight = psImageSubset(readout->weight, srcRegion); + mySource->mask = psImageSubset(readout->mask, srcRegion); + mySource->region = srcRegion; + + return true; +} + +bool pmSourceRedefinePixels(pmSource *mySource, + const pmReadout *readout, + psF32 x, + psF32 y, + psF32 Radius) +{ + bool extend; + psRegion newRegion; + + if (Radius == 0) + return false; + + // check to see if new region is completely contained within old region + newRegion = psRegionForSquare (x, y, Radius); + newRegion = psRegionForImage (readout->image, newRegion); + + extend = false; + extend |= (int)(newRegion.x0) < (int)(mySource->region.x0); + extend |= (int)(newRegion.x1) > (int)(mySource->region.x1); + extend |= (int)(newRegion.y0) < (int)(mySource->region.y0); + extend |= (int)(newRegion.y1) > (int)(mySource->region.y1); + + extend |= (mySource->pixels == NULL); + extend |= (mySource->weight == NULL); + extend |= (mySource->mask == NULL); + + // extend = true; + if (extend) { + // Grab a new subimage + psFree (mySource->pixels); + psFree (mySource->weight); + psFree (mySource->mask); + + mySource->pixels = psImageSubset(readout->image, newRegion); + mySource->weight = psImageSubset(readout->weight, newRegion); + mySource->mask = psImageSubset(readout->mask, newRegion); + mySource->region = newRegion; + } + + return extend; +} + +/****************************************************************************** + pmSourcePSFClump(source, metadata): Find the likely PSF clump in the + sigma-x, sigma-y plane. return 0,0 clump in case of error. +*****************************************************************************/ + +// XXX EAM include a S/N cutoff in selecting the sources? +pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + # define NPIX 10 + # define SCALE 0.1 + + psArray *peaks = NULL; + pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; + pmPSFClump psfClump = emptyClump; + + PS_ASSERT_PTR_NON_NULL(sources, emptyClump); + PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); + + // find the sigmaX, sigmaY clump + { + psStats *stats = NULL; + psImage *splane = NULL; + int binX, binY; + bool status; + + psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX"); + if (!status) + SX_MAX = 10.0; + psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX"); + if (!status) + SY_MAX = 10.0; + + // construct a sigma-plane image + // psImageAlloc does zero the data + splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32); + for (int i = 0; i < splane->numRows; i++) + { + memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); + } + + // place the sources in the sigma-plane image (ignore 0,0 values?) + for (psS32 i = 0 ; i < sources->n ; i++) + { + pmSource *tmpSrc = (pmSource *) sources->data[i]; + if (tmpSrc == NULL) { + continue; + } + if (tmpSrc->moments == NULL) { + continue; + } + + // Sx,Sy are limited at 0. a peak at 0,0 is artificial + if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { + continue; + } + + // for the moment, force splane dimensions to be 10x10 image pix + binX = tmpSrc->moments->Sx/SCALE; + if (binX < 0) + continue; + if (binX >= splane->numCols) + continue; + + binY = tmpSrc->moments->Sy/SCALE; + if (binY < 0) + continue; + if (binY >= splane->numRows) + continue; + + splane->data.F32[binY][binX] += 1.0; + } + + // find the peak in this image + stats = psStatsAlloc (PS_STAT_MAX); + stats = psImageStats (stats, splane, NULL, 0); + peaks = pmFindImagePeaks (splane, stats[0].max / 2); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); + + psFree (splane); + psFree (stats); + + } + // XXX EAM : possible errors: + // 1) no peak in splane + // 2) no significant peak in splane + + // measure statistics on Sx, Sy if Sx, Sy within range of clump + { + pmPeak *clump; + psF32 minSx, maxSx; + psF32 minSy, maxSy; + psVector *tmpSx = NULL; + psVector *tmpSy = NULL; + psStats *stats = NULL; + + // XXX EAM : this lets us takes the single highest peak + psArraySort (peaks, pmComparePeakDescend); + clump = peaks->data[0]; + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); + + // define section window for clump + minSx = clump->x * SCALE - 0.2; + maxSx = clump->x * SCALE + 0.2; + minSy = clump->y * SCALE - 0.2; + maxSy = clump->y * SCALE + 0.2; + + tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); + tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); + tmpSx->n = 0; + tmpSy->n = 0; + + // XXX clip sources based on flux? + // create vectors with Sx, Sy values in window + for (psS32 i = 0 ; i < sources->n ; i++) + { + pmSource *tmpSrc = (pmSource *) sources->data[i]; + + if (tmpSrc->moments->Sx < minSx) + continue; + if (tmpSrc->moments->Sx > maxSx) + continue; + if (tmpSrc->moments->Sy < minSy) + continue; + if (tmpSrc->moments->Sy > maxSy) + continue; + tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; + tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; + tmpSx->n++; + tmpSy->n++; + if (tmpSx->n == tmpSx->nalloc) { + psVectorRealloc (tmpSx, tmpSx->nalloc + 100); + psVectorRealloc (tmpSy, tmpSy->nalloc + 100); + } + } + + // measures stats of Sx, Sy + stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + + stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); + psfClump.X = stats->clippedMean; + psfClump.dX = stats->clippedStdev; + + stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); + psfClump.Y = stats->clippedMean; + psfClump.dY = stats->clippedStdev; + + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); + // these values should be pushed on the metadata somewhere + + psFree (stats); + psFree (peaks); + psFree (tmpSx); + psFree (tmpSy); + } + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return (psfClump); +} + +/****************************************************************************** + pmSourceRoughClass(source, metadata): make a guess at the source + classification. + + XXX: push the clump info into the metadata? + + XXX: How can this function ever return FALSE? + + EAM: I moved S/N calculation to pmSourceMoments, using weight image +*****************************************************************************/ + +bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + psBool rc = true; + + int Nsat = 0; + int Next = 0; + int Nstar = 0; + int Npsf = 0; + int Ncr = 0; + int Nsatstar = 0; + // psRegion allArray = psRegionSet (0, 0, 0, 0); + psRegion inner; + + // report stats on S/N values for star-like objects + psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); + starsn->n = 0; + + // check return status value (do these exist?) + bool status; + psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); + + // XXX allow clump size to be scaled relative to sigmas? + // make rough IDs based on clumpX,Y,DX,DY + for (psS32 i = 0 ; i < sources->n ; i++) { + + pmSource *tmpSrc = (pmSource *) sources->data[i]; + + tmpSrc->peak->type = 0; + + psF32 sigX = tmpSrc->moments->Sx; + psF32 sigY = tmpSrc->moments->Sy; + + // XXX EAM : can we use the value of SATURATE if mask is NULL? + // inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2); + inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2); + int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED); + + // saturated star (size consistent with PSF or larger) + // Nsigma should be user-configured parameter + bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); + big = true; + if ((Nsatpix > 1) && big) { + tmpSrc->type = PM_SOURCE_STAR; + tmpSrc->mode = PM_SOURCE_SATSTAR; + Nsatstar ++; + continue; + } + + // saturated object (not a star, eg bleed trails, hot pixels) + if (Nsatpix > 1) { + tmpSrc->type = PM_SOURCE_SATURATED; + tmpSrc->mode = PM_SOURCE_DEFAULT; + Nsat ++; + continue; + } + + // likely defect (too small to be stellar) (push out to 3 sigma) + // low S/N objects which are small are probably stellar + // only set candidate defects if + if ((sigX < 0.05) || (sigY < 0.05)) { + tmpSrc->type = PM_SOURCE_DEFECT; + tmpSrc->mode = PM_SOURCE_DEFAULT; + Ncr ++; + continue; + } + + // likely unsaturated extended source (too large to be stellar) + if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { + tmpSrc->type = PM_SOURCE_EXTENDED; + tmpSrc->mode = PM_SOURCE_DEFAULT; + Next ++; + continue; + } + + // the rest are probable stellar objects + starsn->data.F32[starsn->n] = tmpSrc->moments->SN; + starsn->n ++; + Nstar ++; + + // PSF star (within 1.5 sigma of clump center, S/N > limit) + psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); + if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) { + tmpSrc->type = PM_SOURCE_STAR; + tmpSrc->mode = PM_SOURCE_PSFSTAR; + Npsf ++; + continue; + } + + // random type of star + tmpSrc->type = PM_SOURCE_STAR; + tmpSrc->mode = PM_SOURCE_DEFAULT; + } + + { + psStats *stats = NULL; + stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); + stats = psVectorStats (stats, starsn, NULL, NULL, 0); + psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); + psFree (stats); + psFree (starsn); + } + + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Next: %3d\n", Next); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); + + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + Index: /branches/rel10_ifa/psModules/src/objects/pmSource.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSource.h (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSource.h (revision 6537) @@ -0,0 +1,190 @@ +/** @file pmSource.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_H +# define PM_SOURCE_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 { + PM_SOURCE_MASK_CLEAR = 0x00, + PM_SOURCE_MASK_INVALID = 0x01, + PM_SOURCE_MASK_SATURATED = 0x02, + PM_SOURCE_MASK_MARKED = 0x08, +} psphotMaskValues; + +/** 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_TYPE_UNKNOWN, ///< a cosmic-ray + PM_SOURCE_TYPE_DEFECT, ///< a cosmic-ray + PM_SOURCE_TYPE_SATURATED, ///< random saturated pixels + PM_SOURCE_TYPE_STAR, ///< a good-quality star + PM_SOURCE_TYPE_EXTENDED, ///< an extended object (eg, galaxy) +} pmSourceType; + +typedef enum { + PM_SOURCE_MODE_DEFAULT = 0x0000, ///< + PM_SOURCE_MODE_PSFMODEL = 0x0001, ///< + PM_SOURCE_MODE_EXTMODEL = 0x0002, ///< + PM_SOURCE_MODE_SUBTRACTED = 0x0004, ///< + PM_SOURCE_MODE_FITTED = 0x0008, ///< + PM_SOURCE_MODE_FAIL = 0x0010, ///< + PM_SOURCE_MODE_POOR = 0x0020, ///< + PM_SOURCE_MODE_PAIR = 0x0040, ///< + PM_SOURCE_MODE_PSFSTAR = 0x0080, ///< + PM_SOURCE_MODE_SATSTAR = 0x0100, ///< + PM_SOURCE_MODE_BLEND = 0x0200, ///< + PM_SOURCE_MODE_LINEAR = 0x0400, ///< + PM_SOURCE_MODE_TEMPSUB = 0x0800, ///< XXX get me a better name! +} pmSourceMode; + +/** 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 *modelEXT; ///< EXT (floating) Model fit (parameters and type). + pmSourceType type; ///< Best identification of object. + pmSourceMode mode; ///< Best identification of object. + psArray *blends; + float fitMag; + float apMag; + psRegion region; // area on image covered by selected pixels +} +pmSource; + +/** 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; + + +/** pmSourceAlloc() + * + */ +pmSource *pmSourceAlloc(); + + +/** 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: 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. +); + +bool pmSourceDefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius); +bool pmSourceRedefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius); + +/** 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 +); + + +# endif /* PM_SOURCE_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSourceContour.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceContour.c (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceContour.c (revision 6537) @@ -0,0 +1,422 @@ +# include "psphot.h" + +/****************************************************************************** +findValue(source, level, row, col, dir): a private function which determines +the column coordinate of the model function which has the value "level". If +dir equals 0, then you loop leftwards from the peak pixel, otherwise, +rightwards. + +XXX: reverse order of row,col args? + +XXX: Input row/col are in image coords. + +XXX: The result is returned in image coords. +*****************************************************************************/ +# define LEFT false +# define RIGHT true + +// return the first coordinate at or below the threshold in the requested direction +static int findContourNeg( + psImage *image, + float threshold, + int x, + int y, + bool right) +{ + + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + + int incr; + int subCol; + int lastColumn; + if (right) { + incr = 1; + lastColumn = image->numCols - 1; + } else { + incr = -1; + lastColumn = 0; + } + + subCol = x; + + while (subCol != lastColumn) { + float value = image->data.F32[y][subCol]; + if (value <= threshold) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (subCol); + } + subCol += incr; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (lastColumn); +} + +// return the last coordinate at or below the threshold in the requested direction +static int findContourPos( + psImage *image, + float threshold, + int x, + int y, + bool right, + int xEnd) +{ + + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + + int incr; + int subCol; + int lastColumn; + if (right) { + incr = 1; + lastColumn = PS_MIN (image->numCols - 1, xEnd); + } else { + incr = -1; + lastColumn = PS_MAX (0, xEnd); + } + + subCol = x; + while (subCol != lastColumn) { + float value = image->data.F32[y][subCol]; + if (value >= threshold) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + if (subCol == x) { + return (subCol); + } + return (subCol); + } + subCol += incr; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (lastColumn); +} + +/****************************************************************************** +findValue(source, level, row, col, dir): a private function which determines +the column coordinate of the model function which has the value "level". If +dir equals 0, then you loop leftwards from the peak pixel, otherwise, +rightwards. + +XXX: reverse order of row,col args? + +XXX: Input row/col are in image coords. + +XXX: The result is returned in image coords. +*****************************************************************************/ +static psF32 findValue(pmSource *source, + psF32 level, + psU32 row, + psU32 col, + psU32 dir) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + // + // Convert coords to subImage space. + // + psU32 subRow = row - source->pixels->row0; + psU32 subCol = col - source->pixels->col0; + + // Ensure that the starting column is allowable. + if (!((0 <= subCol) && (subCol < source->pixels->numCols))) { + psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); + return(NAN); + } + if (!((0 <= subRow) && (subRow < source->pixels->numRows))) { + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); + psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range"); + return(NAN); + } + + // XXX EAM : i changed this to match pmModelEval above, but see + // XXX EAM the note below in pmSourceContour + psF32 oldValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow); + if (oldValue == level) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return(((psF32) (subCol + source->pixels->col0))); + } + + // + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + // + psS32 incr; + psS32 lastColumn; + if (dir == 0) { + incr = -1; + lastColumn = -1; + } else { + incr = 1; + lastColumn = source->pixels->numCols; + } + subCol+=incr; + + while (subCol != lastColumn) { + psF32 newValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow); + if (oldValue == level) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return((psF32) (subCol + source->pixels->col0)); + } + + if ((newValue <= level) && (level <= oldValue)) { + // This is simple linear interpolation. + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); + } + + if ((oldValue <= level) && (level <= newValue)) { + // This is simple linear interpolation. + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); + } + + subCol+=incr; + } + + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); + return(NAN); +} + +/****************************************************************************** +new implementation of source contour function +*****************************************************************************/ +psArray *pmSourceContour (psImage *image, int xc, int yc, float threshold) +{ + + int xR, yR, x0, x1, x0s, x1s; + + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(image, false); + + int x = xc - image->col0; + int y = yc - image->row0; + + // the requested point must be within the contour + if (image->data.F32[y][x] < threshold) + return NULL; + + // Ensure that the starting column is allowable. + if (x < 0) + return NULL; + if (y < 0) + return NULL; + if (x >= image->numCols) + return NULL; + if (y >= image->numRows) + return NULL; + + // Allocate data for x/y pairs. + psVector *xVec = psVectorAlloc(100, PS_TYPE_F32); + psVector *yVec = psVectorAlloc(100, PS_TYPE_F32); + + // First row: find the left and right end-points + int Npt = 0; + + x0 = findContourNeg (image, threshold, x, y, LEFT); + x1 = findContourNeg (image, threshold, x, y, RIGHT); + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + y; + yVec->data.F32[Npt + 1] = image->row0 + y; + Npt += 2; + + x0s = x0; + x1s = x1; + + // look for contour outline above row + xR = x0s; + yR = y + 1; + while (yR < image->numRows) { + if (image->data.F32[yR][xR] < threshold) { + x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); + if (x0 == x1) { + // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + goto pt1; + } + x1 = findContourNeg (image, threshold, x0, yR, RIGHT); + x0--; + } else { + x0 = findContourNeg (image, threshold, xR, yR, LEFT); + x1 = findContourNeg (image, threshold, xR, yR, RIGHT); + } + // fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1); + + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + yR; + yVec->data.F32[Npt + 1] = image->row0 + yR; + Npt += 2; + + if (Npt >= xVec->nalloc - 1) { + psVectorRealloc (xVec, xVec->nalloc + 100); + psVectorRealloc (yVec, yVec->nalloc + 100); + } + yR ++; + xR = x0; + } + +pt1: + // look for contour outline below row + xR = x0s; + x1 = x1s; + yR = y - 1; + while (yR >= 0) { + if (image->data.F32[yR][xR] < threshold) { + x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); + if (x0 == x1) { + // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + goto pt2; + } + x1 = findContourNeg (image, threshold, x0, yR, RIGHT); + x0--; + } else { + x0 = findContourNeg (image, threshold, xR, yR, LEFT); + x1 = findContourNeg (image, threshold, xR, yR, RIGHT); + } + // fprintf (stderr, "neg: %d (%d - %d)\n", yR, x0, x1); + + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + yR; + yVec->data.F32[Npt + 1] = image->row0 + yR; + Npt += 2; + + if (Npt >= xVec->nalloc - 1) { + psVectorRealloc (xVec, xVec->nalloc + 100); + psVectorRealloc (yVec, yVec->nalloc + 100); + } + yR --; + } +pt2: + xVec->n = Npt; + yVec->n = Npt; + + // fprintf (stderr, "done\n"); + psArray *tmpArray = psArrayAlloc(2); + tmpArray->data[0] = (psPtr *) xVec; + tmpArray->data[1] = (psPtr *) yVec; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpArray); +} + +/****************************************************************************** + pmSourceContour(src, img, level, mode): For an input subImage, and model, this + routine returns a psArray of coordinates that evaluate to the specified level. + + XXX: Probably should remove the "image" argument. + XXX: What type should the output coordinate vectors consist of? col,row? + XXX: Why a pmArray output? + XXX: doex x,y correspond with col,row or row/col? + XXX: What is mode? + XXX: The top, bottom of the contour is not correctly determined. + XXX EAM : this function is using the model for the contour, but it should + be using only the image counts +*****************************************************************************/ +psArray *pmSourceContour_Crude(pmSource *source, + const psImage *image, + psF32 level) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(image, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->modelEXT, false); + // XXX EAM : what is the purpose of modelPSF/modelEXT? + + // + // Allocate data for x/y pairs. + // + psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); + psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); + + // + // Start at the row with peak pixel, then decrement. + // + psS32 col = source->peak->x; + for (psS32 row = source->peak->y; row>= 0 ; row--) { + // XXX: yVec contain no real information. Do we really need it? + yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); + yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); + + // Starting at peak pixel, search leftwards for the column intercept. + psF32 leftIntercept = findValue(source, level, row, col, 0); + if (isnan(leftIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; + + // Starting at peak pixel, search rightwards for the column intercept. + + psF32 rightIntercept = findValue(source, level, row, col, 1); + if (isnan(rightIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); + xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; + + // Set starting column for next row + col = (psS32) ((leftIntercept + rightIntercept) / 2.0); + } + // + // Start at the row (+1) with peak pixel, then increment. + // + col = source->peak->x; + for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) { + // XXX: yVec contain no real information. Do we really need it? + yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); + yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); + + // Starting at peak pixel, search leftwards for the column intercept. + psF32 leftIntercept = findValue(source, level, row, col, 0); + if (isnan(leftIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; + + // Starting at peak pixel, search rightwards for the column intercept. + psF32 rightIntercept = findValue(source, level, row, col, 1); + if (isnan(rightIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; + + // Set starting column for next row + col = (psS32) ((leftIntercept + rightIntercept) / 2.0); + } + + // + // Allocate an array for result, store coord vectors there. + // + psArray *tmpArray = psArrayAlloc(2); + tmpArray->data[0] = (psPtr *) yVec; + tmpArray->data[1] = (psPtr *) xVec; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpArray); +} Index: /branches/rel10_ifa/psModules/src/objects/pmSourceContour.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceContour.h (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceContour.h (revision 6537) @@ -0,0 +1,37 @@ +/** @file pmSourceContour.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_CONTOUR_H +# define PM_SOURCE_CONTOUR_H + +psArray *pmSourceContour (psImage *image, int xc, int yc, float threshold) +{ + + /** 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_Crude( + pmSource *source, ///< The input pmSource + const psImage *image, ///< The input image (float) (this arg should be removed) + float level ///< The level of the contour + ); + + # endif /* PM_SOURCE_PHOTOMETRY_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.c (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.c (revision 6537) @@ -0,0 +1,228 @@ +/** @file pmObjects.c + * + * This file will ... + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmObjects.h" +#include "pmModelGroup.h" + +// save a static values so they may be set externally +static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; +static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; + +bool pmSourceFitModelInit (float nIter, float tol) +{ + + PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter; + PM_SOURCE_FIT_MODEL_TOLERANCE = tol; + return true; +} + +bool pmSourceFitModel (pmSource *source, + pmModel *model, + const bool PSF) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->weight, false); + + // XXX EAM : is it necessary for the mask & weight to exist? the + // tests below could be conditions (!NULL) + + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + psVector *params = model->params; + psVector *dparams = model->dparams; + psVector *paramMask = NULL; + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + + int nParams = PSF ? 4 : params->n; + + // maximum number of valid pixels + psS32 nPix = source->pixels->numRows * source->pixels->numCols; + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(nPix); + psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); + + nPix = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + // skip masked points + if (source->mask->data.U8[i][j]) { + continue; + } + // skip zero-weight points + if (source->weight->data.F32[i][j] == 0) { + continue; + } + + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + // Convert i/j to image space: + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[nPix] = (psPtr *) coord; + y->data.F32[nPix] = source->pixels->data.F32[i][j]; + // psMinimizeLMChi2 takes wt = 1/dY^2 + yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + nPix++; + } + } + x->n = nPix; + y->n = nPix; + yErr->n = nPix; + if (nPix < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + model->status = PM_MODEL_BADARGS; + psFree (x); + psFree (y); + psFree (yErr); + return(false); + } + + // XXX EAM : the new minimization API supplies the constraints as a struct + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + psMinConstrain *constrain = psMinConstrainAlloc(); + + // PSF model only fits first 4 parameters, EXT model fits all + if (PSF) { + paramMask = psVectorAlloc (params->n, PS_TYPE_U8); + for (int i = 0; i < 4; i++) { + paramMask->data.U8[i] = 0; + } + for (int i = 4; i < paramMask->n; i++) { + paramMask->data.U8[i] = 1; + } + } + constrain->paramMask = paramMask; + + // Set the parameter range checks + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); + modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax); + + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); + + psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); + + fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if ((paramMask != NULL) && paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // save the resulting chisq, nDOF, nIter + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // get the Gauss-Newton distance for fixed model parameters + if (paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if (!paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + psFree (delta); + } + + // set the model success or failure status + if (!fitStatus) { + model->status = PM_MODEL_NONCONVERGE; + } else { + model->status = PM_MODEL_SUCCESS; + } + + // models can go insane: reject these + onPic &= (params->data.F32[2] >= source->pixels->col0); + onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (params->data.F32[3] >= source->pixels->row0); + onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + if (!onPic) { + model->status = PM_MODEL_OFFIMAGE; + } + + source->mode |= PM_SOURCE_FITTED; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(constrain->paramMask); + psFree(constrain->paramMin); + psFree(constrain->paramMax); + psFree(constrain->paramDelta); + psFree(constrain); + + rc = (onPic && fitStatus); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +pmModel *pmSourceSelectModel (pmSource *source) +{ + switch (source->type) { + case PM_SOURCE_STAR: + return source->modelPSF; + + case PM_SOURCE_EXTENDED: + return source->modelEXT; + + default: + return NULL; + } + return NULL; +} + +/****************************************************************************** + pmSourceModelGuess(source, model): This function allocates a new + pmModel structure based on the given modelType specified in the argument list. + The corresponding pmModelGuess function is returned, and used to + supply the values of the params array in the pmModel structure. + + XXX: Many parameters are based on the src->moments structure, which is in + image, not subImage coords. Therefore, the calls to the model evaluation + functions will be in image, not subImage coords. Remember this. +*****************************************************************************/ +pmModel *pmSourceModelGuess(pmSource *source, + pmModelType modelType) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + + pmModel *model = pmModelAlloc(modelType); + + pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); + modelGuessFunc(model, source); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(model); +} + Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.h (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.h (revision 6537) @@ -0,0 +1,52 @@ +/** @file pmSourceFitModel.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_FIT_MODEL_H +# define PM_SOURCE_FIT_MODEL_H + +bool pmSourceFitModelInit( + float nIter, ///< max number of allowed iterations + float tol ///< convergence criterion +); + +/** 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 EXT? +); + + +pmModel *pmSourceSelectModel (pmSource *source); + +/** 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. +); + +# endif /* PM_SOURCE_FIT_MODEL_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.c (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.c (revision 6537) @@ -0,0 +1,297 @@ +/** @file pmSourceFitSet.h + * + * @author EAM, IfA + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2, +// nPar = nSrc*(nOnePar - 1) + 1 + +# include +# include +# include +# include "pmSource.h" + +static pmModelFunc mFunc; +static int nPar; +static psVector *onePar; +static psVector *oneDeriv; + +bool pmModelFitSetInit (pmModelType type) +{ + + mFunc = pmModelFunc_GetFunction (type); + nPar = pmModelParameterCount (type); + + onePar = psVectorAlloc (nPar, PS_DATA_F32); + oneDeriv = psVectorAlloc (nPar, PS_DATA_F32); + + return true; +} + +void pmModelFitSetClear (void) +{ + + psFree (onePar); + psFree (oneDeriv); + return; +} + +psF32 pmModelFitSet(psVector *deriv, + const psVector *params, + const psVector *x) +{ + + psF32 value; + psF32 model; + + psF32 *PAR = onePar->data.F32; + psF32 *dPAR = oneDeriv->data.F32; + + psF32 *pars = params->data.F32; + psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32; + + int nSrc = (params->n - 1) / (nPar - 1); + + PAR[0] = model = pars[0]; + for (int i = 0; i < nSrc; i++) { + int nOff = i*nPar - i; + for (int n = 1; n < nPar; n++) { + PAR[n] = pars[nOff + n]; + } + if (deriv == NULL) { + value = mFunc (NULL, onePar, x); + } else { + value = mFunc (oneDeriv, onePar, x); + for (int n = 1; n < nPar; n++) { + dpars[nOff + n] = dPAR[n]; + } + } + model += value; + } + if (deriv != NULL) { + dpars[0] = dPAR[0]*2.0; + } + return (model); +} + +/* +i: 0 1 2 +n: 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 +i*6 + n: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 +*/ + +# define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 +# define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 + +bool pmSourceFitSet (pmSource *source, + psArray *modelSet, + const bool PSF) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->weight, false); + + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + // base values on first model + pmModel *model = modelSet->data[0]; + + // set the static variables + pmSourceFitSetInit (model->type); + + int nSrc = modelSet->n; + int nPar = model->params->n - 1; // number of object parameters (excluding sky) + + // define parameter vectors for source set + psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); + + // define limits for single-source model + psVector *oneDelta; + psVector *oneParMin; + psVector *oneParMax; + modelLimits (&oneDelta, &oneParMin, &oneParMax); + + psMinConstrain *constrain = psMinConstrainAlloc(); + constrain->paramDelta = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + constrain->paramMin = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + constrain->paramMax = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + constrain->paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL; + + // all but the sky are allowed to vary independently (subject to PSF) + // set the parameters from the multiple models + // also, set limits based on single-source limits + params->data.F32[0] = model->params->data.F32[0]; + constrain->paramDelta->data.F32[0] = oneDelta->data.F32[0]; + constrain->paramMin->data.F32[0] = oneParMin->data.F32[0]; + constrain->paramMax->data.F32[0] = oneParMax->data.F32[0]; + constrain->paramMask->data.U8[0] = 0; + for (int i = 0; i < nSrc; i++) { + model = modelSet->data[i]; + for (int n = 1; n < nPar + 1; n++) { + params->data.F32[i*nPar + n] = model->params->data.F32[n]; + dparams->data.F32[i*nPar + n] = model->dparams->data.F32[n]; + constrain->paramDelta->data.F32[i*nPar + n] = oneDelta->data.F32[n]; + constrain->paramMin->data.F32[i*nPar + n] = oneParMin->data.F32[n]; + constrain->paramMax->data.F32[i*nPar + n] = oneParMax->data.F32[n]; + if (PSF) { + constrain->paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1; + } + } + } + psFree (oneDelta); + psFree (oneParMin); + psFree (oneParMax); + + if (psTraceGetLevel(__func__) >= 5) { + for (int i = 0; i < params->n; i++) { + fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constrain->paramMask->data.U8[i]); + } + } + + // PSF model only fits first 4 parameters, EXT model fits all + int nParams = PSF ? nSrc*3 + 1 : nSrc*nPar + 1; + + // maximum number of valid pixels + psS32 nPix = source->pixels->numRows * source->pixels->numCols; + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(nPix); + psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); + + nPix = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + // skip masked points + if (source->mask->data.U8[i][j]) { + continue; + } + // skip zero-weight points + if (source->weight->data.F32[i][j] == 0) { + continue; + } + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + // Convert i/j to image space: + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[nPix] = (psPtr *) coord; + y->data.F32[nPix] = source->pixels->data.F32[i][j]; + // psMinimizeLMChi2 takes wt = 1/dY^2 + yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + nPix++; + } + } + if (nPix < nParams + 1) { + psTrace (__func__, 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s() end : fail pixels ----\n", __func__); + model->status = PM_MODEL_BADARGS; + psFree (x); + psFree (y); + psFree (yErr); + psFree (params); + psFree (dparams); + psFree(constrain->paramMask); + psFree(constrain->paramMin); + psFree(constrain->paramMax); + psFree(constrain->paramDelta); + psFree(constrain); + return(false); + } + x->n = nPix; + y->n = nPix; + yErr->n = nPix; + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); + + psTrace (__func__, 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmModelFitSet); + + // parameter errors from the covariance matrix + for (int i = 0; i < dparams->n; i++) { + if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // get the Gauss-Newton distance for fixed model parameters + if (constrain->paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, pmModelFitSet); + for (int i = 0; i < dparams->n; i++) { + if (!constrain->paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + psFree (delta); + } + + // assign back the parameters to the models + for (int i = 0; i < nSrc; i++) { + model = modelSet->data[i]; + model->params->data.F32[0] = params->data.F32[0]; + for (int n = 1; n < nPar + 1; n++) { + model->params->data.F32[n] = params->data.F32[i*nPar + n]; + model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n]; + } + // save the resulting chisq, nDOF, nIter + // these are not unique for any one source + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // set the model success or failure status + model->status = fitStatus ? PM_MODEL_SUCCESS : PM_MODEL_NONCONVERGE; + + // models can go insane: reject these + onPic &= (model->params->data.F32[2] >= source->pixels->col0); + onPic &= (model->params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (model->params->data.F32[3] >= source->pixels->row0); + onPic &= (model->params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + if (!onPic) { + model->status = PM_MODEL_OFFIMAGE; + } + } + + source->mode |= PM_SOURCE_FITTED; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(constrain->paramMask); + psFree(constrain->paramMin); + psFree(constrain->paramMax); + psFree(constrain->paramDelta); + psFree(constrain); + psFree(params); + psFree(dparams); + + // free static memory used by pmModelFitSet + pmModelFitSetClear (); + + rc = (onPic && fitStatus); + psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF); + psTrace(__func__, 3, "---- %s end : status %d ----\n", __func__, rc); + return(rc); +} + Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.h (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.h (revision 6537) @@ -0,0 +1,38 @@ +/** @file pmSourceFitset.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_FIT_SET_H +# define PM_SOURCE_FIT_SET_H + +psF32 pmModelFitSet(psVector *deriv, + const psVector *params, + const psVector *x); + +bool pmModelFitSetInit (pmModelType type); + +void pmModelFitSetClear (void); + +/** pmSourceFitSet() + * + * 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 pmSourceFitSet( + pmSource *source, ///< The input pmSource + psArray *modelSet, ///< model to be fitted + const bool PSF ///< Treat model as PSF or EXT? +); + +# endif /* PM_SOURCE_FIT_MODEL_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c (revision 6537) @@ -0,0 +1,217 @@ +/** @file pmSourcePhotometry.C + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +// XXX EAM : the apMag should only be calculated for the brighter sources? +// XXX EAM : SN limit set by user? +// XXX EAM : masked region should be (optionally) elliptical +pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius) +{ + + int status; + bool isPSF; + float x, y; + float rflux; + float radius; + pmModel *model; + + switch (source->type) { + case PM_SOURCE_STAR: + model = source->modelPSF; + if (model == NULL) + return NULL; + radius = (apRadius > 0) ? apRadius : model->radius; + isPSF = true; + break; + + case PM_SOURCE_EXTENDED: + model = source->modelEXT; + if (model == NULL) + return NULL; + radius = model->radius; + isPSF = false; + break; + + default: + return NULL; + } + + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + + // replace source flux + pmModelAdd (source->pixels, source->mask, model, false, false); + + // set aperture mask circle of PSF_FIT_RADIUS + psImageKeepCircle (source->mask, x, y, radius, "OR", PSPHOT_MASK_MARKED); + + // measure object photometry + status = pmSourcePhotometry (&source->fitMag, &source->apMag, model, source->pixels, source->mask); + + // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS + if (isPSF && (psf != NULL)) { + if (psf->growth != NULL) { + source->apMag += pmGrowthCurveCorrect (psf->growth, model->radius); + } + + rflux = pow (10.0, 0.4*source->fitMag); + source->apMag -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux; + source->fitMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0); + } + + // unmask aperture + psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED); + + // subtract object, leave local sky + pmModelSub (source->pixels, source->mask, model, false, false); + + if (!status) + return NULL; + return model; +} + +/* +aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y) +(aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y) +(aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y) + +*/ + +// XXX split into ap, psf, ext mags? +bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) +{ + + float obsSum = 0; + float fitSum = 0; + float sky = model->params->data.F32[0]; + + *fitMag = 99.0; + *obsMag = 99.0; + + pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); + fitSum = modelFluxFunc (model->params); + if (fitSum <= 0) + return false; + if (!isfinite(fitSum)) + return false; + *fitMag = -2.5*log10(fitSum); + + for (int ix = 0; ix < image->numCols; ix++) { + for (int iy = 0; iy < image->numRows; iy++) { + if (mask->data.U8[iy][ix]) + continue; + obsSum += image->data.F32[iy][ix] - sky; + } + } + if (obsSum <= 0) + return false; + *obsMag = -2.5*log10(obsSum); + + return (true); +} + +float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj) +{ + + int Xs, Xe, Ys, Ye; + int xi, xj, yi, yj; + int xIs, xJs, yIs, yJs; + int xIe, yIe; + float flux, wt; + + psImage *Pi = Mi->pixels; + psImage *Pj = Mj->pixels; + + psImage *Wi = Mi->weight; + + psImage *Ti = Mi->mask; + psImage *Tj = Mj->mask; + + Xs = PS_MAX (Pi->col0, Pj->col0); + Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols); + + Ys = PS_MAX (Pi->row0, Pj->row0); + Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows); + + xIs = Xs - Pi->col0; + xJs = Xs - Pj->col0; + yIs = Ys - Pi->row0; + yJs = Ys - Pj->row0; + + xIe = Xe - Pi->col0; + yIe = Ye - Pi->row0; + + // note that this is addressing the same image pixels, + // though only if both are source not model images + flux = 0; + for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { + for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { + if (Ti->data.U8[yi][xi]) + continue; + if (Tj->data.U8[yj][xj]) + continue; + wt = Wi->data.F32[yi][xi]; + if (wt > 0) { + flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt; + } + } + } + return (flux); +} + +float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj) +{ + + int Xs, Xe, Ys, Ye; + int xi, xj, yi, yj; + int xIs, xJs, yIs, yJs; + int xIe, yIe; + float flux, wt; + + psImage *Pi = Mi->pixels; + psImage *Pj = Mj->pixels; + + psImage *Wi = Mi->weight; + + psImage *Ti = Mi->mask; + psImage *Tj = Mj->mask; + + Xs = PS_MAX (Pi->col0, Pj->col0); + Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols); + + Ys = PS_MAX (Pi->row0, Pj->row0); + Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows); + + xIs = Xs - Pi->col0; + xJs = Xs - Pj->col0; + yIs = Ys - Pi->row0; + yJs = Ys - Pj->row0; + + xIe = Xe - Pi->col0; + yIe = Ye - Pi->row0; + + // note that this is addressing the same image pixels, + // though only if both are source not model images + flux = 0; + for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { + for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { + if (Ti->data.U8[yi][xi]) + continue; + if (Tj->data.U8[yj][xj]) + continue; + wt = Wi->data.F32[yi][xi]; + if (wt > 0) { + flux += 1.0 / wt; + } + } + } + return (flux); +} + Index: /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.h (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.h (revision 6537) @@ -0,0 +1,39 @@ +/** @file pmSourcePhotometry.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_PHOTOMETRY_H +# define PM_SOURCE_PHOTOMETRY_H + +/** + * + * 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 +); + +pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); +float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj); +float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj); + +# endif /* PM_SOURCE_PHOTOMETRY_H */ Index: /branches/rel10_ifa/psModules/src/objects/pmSourceSky.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceSky.c (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceSky.c (revision 6537) @@ -0,0 +1,124 @@ + +/****************************************************************************** +pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this +routine creates a new pmSource data structure and sets the following members: + ->pmPeak + ->pmMoments->sky + +The sky value is set from the pixels in the square annulus surrounding the +peak pixel. + +We simply create a subSet image and mask the inner pixels, then call +psImageStats on that subImage+mask. + +XXX: The subImage has width of 1+2*outerRadius. Verify with IfA. + +XXX: Use static data structures for: + subImage + subImageMask + myStats + +XXX: ensure that the inner and out radius fit in the actual image. Should + we generate an error, or warning? Currently an error. + +XXX: Sync with IfA on whether the peak x/y coords are data structure coords, + or they use the image row/column offsets. +XXX EAM : peak->x,y uses parent coordinates + +XXX: Should we simply set pmSource->peak = peak? If so, should we increase +the reference counter? Or, should we copy the data structure? + +XXX: Currently the subimage always has an even number of rows/columns. Is + this correct? Since there is a center pixel, maybe it should have an + odd number of rows/columns. + +XXX: Use psTrace() for the print statements. + +XXX: Don't use separate structs for the subimage and mask. Use the source-> + members. +*****************************************************************************/ + +bool pmSourceLocalSky( + pmSource *source, + psStatsOptions statsOptions, + psF32 Radius) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_IMAGE_NON_NULL(source->pixels, false); + PS_ASSERT_IMAGE_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_INT_POSITIVE(Radius, false); + PS_ASSERT_INT_NONNEGATIVE(Radius, false); + + psImage *image = source->pixels; + psImage *mask = source->mask; + pmPeak *peak = source->peak; + psRegion srcRegion; + + srcRegion = psRegionForSquare(peak->x, peak->y, Radius); + srcRegion = psRegionForImage(mask, srcRegion); + + psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED); + psStats *myStats = psStatsAlloc(statsOptions); + myStats = psImageStats(myStats, image, mask, 0xff); + psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED); + + psF64 tmpF64; + p_psGetStatValue(myStats, &tmpF64); + psFree(myStats); + + if (isnan(tmpF64)) { + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return(false); + } + if (source->moments == NULL) { + source->moments = pmMomentsAlloc(); + } + source->moments->Sky = (psF32) tmpF64; + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return (true); +} + +// A complementary function to pmSourceLocalSky: calculate the local median variance +bool pmSourceLocalSkyVariance( + pmSource *source, + psStatsOptions statsOptions, + psF32 Radius) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_IMAGE_NON_NULL(source->weight, false); + PS_ASSERT_IMAGE_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_INT_POSITIVE(Radius, false); + PS_ASSERT_INT_NONNEGATIVE(Radius, false); + + psImage *image = source->weight; + psImage *mask = source->mask; + pmPeak *peak = source->peak; + psRegion srcRegion; + + srcRegion = psRegionForSquare(peak->x, peak->y, Radius); + srcRegion = psRegionForImage(mask, srcRegion); + + psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED); + psStats *myStats = psStatsAlloc(statsOptions); + myStats = psImageStats(myStats, image, mask, 0xff); + psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED); + + psF64 tmpF64; + p_psGetStatValue(myStats, &tmpF64); + psFree(myStats); + + if (isnan(tmpF64)) { + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return(false); + } + if (source->moments == NULL) { + source->moments = pmMomentsAlloc(); + } + source->moments->dSky = (psF32) tmpF64; + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return (true); +} Index: /branches/rel10_ifa/psModules/src/objects/pmSourceSky.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceSky.h (revision 6537) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceSky.h (revision 6537) @@ -0,0 +1,42 @@ +/** @file pmSourceSky.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-07 06:33:35 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_SKY_H +# define PM_SOURCE_SKY_H + +/** 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 +); + + +// A complementary function to pmSourceLocalSky: calculate the local sky variance +bool pmSourceLocalSkyVariance( + 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 +); + +# endif /* PM_SOURCE_PHOTOMETRY_H */