Index: /branches/eam_rel8_b2/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_rel8_b2/psModules/src/objects/Makefile.am (revision 5606) +++ /branches/eam_rel8_b2/psModules/src/objects/Makefile.am (revision 5606) @@ -0,0 +1,24 @@ +noinst_LTLIBRARIES = libpsmoduleobjects.la + +libpsmoduleobjects_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS) +libpsmoduleobjects_la_LDFLAGS = -release $(PACKAGE_VERSION) +libpsmoduleobjects_la_SOURCES = \ + pmObjects.c \ + pmPSF.c \ + pmPSFtry.c \ + pmModelGroup.c \ + psEllipse.c + +EXTRA_DIST = \ + models/pmModel_GAUSS.c \ + models/pmModel_PGAUSS.c \ + models/pmModel_QGAUSS.c \ + models/pmModel_SGAUSS.c + +psmoduleincludedir = $(includedir) +psmoduleinclude_HEADERS = \ + pmObjects.h \ + pmPSF.h \ + pmPSFtry.h \ + pmModelGroup.h \ + psEllipse.h Index: /branches/eam_rel8_b2/psModules/src/objects/pmObjects.c =================================================================== --- /branches/eam_rel8_b2/psModules/src/objects/pmObjects.c (revision 5606) +++ /branches/eam_rel8_b2/psModules/src/objects/pmObjects.c (revision 5606) @@ -0,0 +1,1934 @@ +/** @file pmObjects.c + * + * This file will ... + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.2.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-26 02:48:02 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmObjects.h" +#include "pmModelGroup.h" +/****************************************************************************** +pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members. +*****************************************************************************/ +pmPeak *pmPeakAlloc(psS32 x, + psS32 y, + psF32 counts, + pmPeakType class) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak)); + tmp->x = x; + tmp->y = y; + tmp->counts = counts; + tmp->class = class; + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + +/****************************************************************************** +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); +} + +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; + 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); +} + +/****************************************************************************** +XXX EAM : we can now free these pixels - memory ref is incremented now +*****************************************************************************/ +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->modelFLT); + 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->modelPSF = NULL; + tmp->modelFLT = NULL; + tmp->type = 0; + psMemSetDeallocator(tmp, (psFreeFunc) sourceFree); + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + +/****************************************************************************** +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); +} + +/****************************************************************************** +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); +} + +/****************************************************************************** +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); +} + +/****************************************************************************** +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: 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; + + // + // Find peaks in row 0 only. + // + row = 0; + tmpRow = getRowVectorFromImage((psImage *) image, row); + psVector *row1 = pmFindVectorPeaks(tmpRow, threshold); + + 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, col, 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, col, 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, col, 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, col, 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, col, 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, col, 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, col, 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, col, 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, col, 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); +} + +// 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); +} + + +/****************************************************************************** +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); +} + +/****************************************************************************** +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: 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); + } + source->moments = pmMomentsAlloc(); + source->moments->Sky = (psF32) tmpF64; + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return (true); +} + +/****************************************************************************** +bool checkRadius2(): private function which simply determines if the (x, y) +point is within the radius of the specified peak. + +XXX: macro this for performance. +XXX: this is rather inefficient - at least compute and compare against radius^2 +*****************************************************************************/ +static bool checkRadius2(psF32 xCenter, + psF32 yCenter, + psF32 radius, + psF32 x, + psF32 y) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + /// XXX EAM should compare with hypot (x,y) for speed + if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { + return(true); + } + + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return(false); +} + +/****************************************************************************** +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 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; + + // 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++) { + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + if ((source->mask != NULL) && (source->mask->data.U8[row][col])) { + continue; + } + + psF32 xDiff = col + source->pixels->col0 - xPeak; + psF32 yDiff = row + source->pixels->row0 - yPeak; + + // XXX EAM : calculate xDiff, yDiff up front; + // radius is just a function of (xDiff, yDiff) + if (!VALID_RADIUS(xDiff, yDiff, R2)) { + continue; + } + + psF32 pDiff = source->pixels->data.F32[row][col] - sky; + + // XXX EAM : check for valid S/N in pixel + // XXX EAM : should this limit be user-defined? + if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) { + continue; + } + + Sum += pDiff; + X1 += xDiff * pDiff; + Y1 += yDiff * pDiff; + XY += xDiff * yDiff * pDiff; + + X2 += PS_SQR(xDiff) * pDiff; + Y2 += PS_SQR(yDiff) * pDiff; + + peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel); + numPixels++; + } + } + // XXX EAM - the limit is a bit arbitrary. make it user defined? + if ((numPixels < 3) || (Sum <= 0)) { + psTrace (".psModules.pmSourceMoments", 5, "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", 5, + "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->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); +} + +// XXX EAM : I used +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); +} + +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); +} + +/****************************************************************************** +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; + + // construct a sigma-plane image + // psImageAlloc does zero the data + splane = psImageAlloc (NPIX/SCALE, NPIX/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? + +XXX EAM : add the saturated mask value to metadata +*****************************************************************************/ + +bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + psBool rc = true; + + int Nsat = 0; + int Ngal = 0; + int Nstar = 0; + int Npsf = 0; + int Ncr = 0; + int Nsatstar = 0; + psRegion allArray = psRegionSet (0, 0, 0, 0); + + psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); + starsn->n = 0; + + // check return status value (do these exist?) + bool status; + psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE"); + psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN"); + psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); + // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE"); + + // 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->class = 0; + + psF32 sigX = tmpSrc->moments->Sx; + psF32 sigY = tmpSrc->moments->Sy; + + // calculate and save signal-to-noise estimates + psF32 S = tmpSrc->moments->Sum; + psF32 A = 4 * M_PI * sigX * sigY; + psF32 B = tmpSrc->moments->Sky; + psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN))); + psF32 SN = (S * sqrt(GAIN) / RT); + tmpSrc->moments->SN = SN; + + // XXX EAM : can we use the value of SATURATE if mask is NULL? + int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, 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)); + if ((Nsatpix > 1) && big) { + tmpSrc->type |= PM_SOURCE_SATSTAR; + Nsatstar ++; + continue; + } + + // saturated object (not a star, eg bleed trails, hot pixels) + if (Nsatpix > 1) { + tmpSrc->type |= PM_SOURCE_SATURATED; + 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; + Ncr ++; + continue; + } + + // likely unsaturated galaxy (too large to be stellar) + if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { + tmpSrc->type |= PM_SOURCE_GALAXY; + Ngal ++; + continue; + } + + // the rest are probable stellar objects + starsn->data.F32[starsn->n] = SN; + starsn->n ++; + Nstar ++; + + // PSF star (within 1.5 sigma of clump center + psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); + if ((SN > PSF_SN_LIM) && (radius < 1.5)) { + tmpSrc->type |= PM_SOURCE_PSFSTAR; + Npsf ++; + continue; + } + + // random type of star + tmpSrc->type |= PM_SOURCE_OTHER; + } + + { + 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, "Ngal: %3d\n", Ngal); + 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); +} + +/** pmSourceDefinePixels() + * + * Define psImage subarrays for the source located at coordinates x,y on the + * image set defined by readout. The pixels defined by this operation consist of + * a square window (of full width 2Radius+1) centered on the pixel which contains + * the given coordinate, in the frame of the readout. The window is defined to + * have limits which are valid within the boundary of the readout image, thus if + * the radius would fall outside the image pixels, the subimage is truncated to + * only consist of valid pixels. If readout->mask or readout->weight are not + * NULL, matching subimages are defined for those images as well. This function + * fails if no valid pixels can be defined (x or y less than Radius, for + * example). This function should be used to define a region of interest around a + * source, including both source and sky pixels. + * + * XXX: must code this. + * + */ +bool pmSourceDefinePixels( + pmSource *mySource, ///< Add comment. + pmReadout *readout, ///< Add comment. + psF32 x, ///< Add comment. + psF32 y, ///< Add comment. + psF32 Radius) ///< Add comment. +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented. Returning FALSE.\n"); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(false); +} + +/****************************************************************************** +pmSourceSetPixelsCircle(source, image, radius) + +XXX: This was replaced by DefinePixels in SDRS. Remove it. +*****************************************************************************/ +bool pmSourceSetPixelsCircle(pmSource *source, + const psImage *image, + psF32 radius) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_FLOAT_COMPARE(0.0, radius, false); + + // + // We define variables for code readability. + // + // XXX: Since the peak->xy coords are in image, not subImage coords, + // these variables should be renamed for clarity (imageCenterRow, etc). + // + psS32 radiusS32 = (psS32) radius; + psS32 SubImageCenterRow = source->peak->y; + psS32 SubImageCenterCol = source->peak->x; + // XXX EAM : for the circle to stay on the image + // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) + psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); + psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1); + psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); + psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); + + // XXX: Must recycle image. + // XXX EAM: this message reflects a programming error we know about. + // i am setting it to a trace message which we can take out + if (source->pixels != NULL) { + psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4, + "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL. Freeing and reallocating.\n"); + psFree(source->pixels); + } + source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, + SubImageStartRow, + SubImageEndCol, + SubImageEndRow)); + + // XXX: Must recycle image. + if (source->mask != NULL) { + psFree(source->mask); + } + source->mask = psImageAlloc(source->pixels->numCols, + source->pixels->numRows, + PS_TYPE_U8); // XXX EAM : type was F32 + + // + // Loop through the subimage mask, initialize mask to 0 or 1. + // XXX EAM: valid pixels should have 0, not 1 + for (psS32 row = 0 ; row < source->mask->numRows; row++) { + for (psS32 col = 0 ; col < source->mask->numCols; col++) { + + if (checkRadius2((psF32) radiusS32, + (psF32) radiusS32, + radius, + (psF32) col, + (psF32) row)) { + source->mask->data.U8[row][col] = 0; + } else { + source->mask->data.U8[row][col] = 1; + } + } + } + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); +} + +/****************************************************************************** +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); +} + +/****************************************************************************** +evalModel(source, level, row): a private function which evaluates the +source->modelPSF function at the specified coords. The coords are subImage, not +image coords. + +NOTE: The coords are in subImage source->pixel coords, not image coords. + +XXX: reverse order of row,col args? + +XXX: rename all coords in this file such that their name defines whether +the coords is in subImage or image space. + +XXX: This should probably be a public pmModules function. + +XXX: Use static vectors for x. + +XXX: Figure out if it's (row, col) or (col, row) for the model functions. + +XXX: For a while, the first psVectorAlloc() was generating a seg fault during +testing. Try to reproduce that and debug. +*****************************************************************************/ + +// XXX EAM : I have made this a public function +// XXX EAM : this now uses a pmModel as the input +// XXX EAM : it was using src->type to find the model, not model->type +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); +} + +/****************************************************************************** +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->modelFLT, 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->modelFLT, 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); +} + +/****************************************************************************** +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 functions is using the model for the contour, but it should + be using only the image counts +*****************************************************************************/ +psArray *pmSourceContour(pmSource *source, + const psImage *image, + psF32 level, + pmContourType mode) +{ + 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->modelFLT, false); + // XXX EAM : what is the purpose of modelPSF/modelFLT? + + // + // 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); +} + +// XXX EAM : these are better starting values, but should be available from metadata? +#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 +#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 +/****************************************************************************** +pmSourceFitModel(source, model): must create the appropiate arguments to the +LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. + +XXX: should there be a mask value? +XXX EAM : fit the specified model (not necessarily the one in source) +*****************************************************************************/ +bool pmSourceFitModel_v5(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); + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + // XXX EAM : is it necessary for the mask & weight to exist? the + // tests below could be conditions (!NULL) + + psVector *params = model->params; + psVector *dparams = model->dparams; + psVector *paramMask = NULL; + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + + int nParams = PSF ? params->n - 4 : params->n; + + // find the number of valid pixels + // XXX EAM : this loop and the loop below could just be one pass + // using the psArrayAdd and psVectorExtend functions + psS32 count = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j] == 0) { + count++; + } + } + } + if (count < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return(false); + } + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(count); + psVector *y = psVectorAlloc(count, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); + psS32 tmpCnt = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j] == 0) { + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + // XXX: Convert i/j to image space: + // XXX EAM: coord order is (x,y) == (col,row) + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[tmpCnt] = (psPtr *) coord; + y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; + yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]); + // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then + // the minimization function calculates sq() + tmpCnt++; + } + } + } + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + // PSF model only fits first 4 parameters, FLT 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; + } + } + + // XXX EAM : covar must be F64? + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); + + psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, 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]); + } + + // XXX EAM: we need to do something (give an error?) if rc is false + // XXX EAM: psMinimizeLMChi2 does not check convergence + + // XXX 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); + + // XXX EAM: save the resulting chisq, nDOF, nIter + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // XXX EAM 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); + } + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(paramMask); + + rc = (onPic && fitStatus); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +// XXX EAM : new version with parameter range limits and weight enhancement +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; + psF32 Ro, ymodel; + + psVector *params = model->params; + psVector *dparams = model->dparams; + psVector *paramMask = NULL; + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + + // XXX EAM : I need to use the sky value to constrain the weight model + int nParams = PSF ? params->n - 4 : params->n; + psF32 So = params->data.F32[0]; + + // find the number of valid pixels + // XXX EAM : this loop and the loop below could just be one pass + // using the psArrayAdd and psVectorExtend functions + psS32 count = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j] == 0) { + count++; + } + } + } + if (count < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return(false); + } + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(count); + psVector *y = psVectorAlloc(count, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); + psS32 tmpCnt = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j] == 0) { + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + // XXX: Convert i/j to image space: + // XXX EAM: coord order is (x,y) == (col,row) + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[tmpCnt] = (psPtr *) coord; + y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; + + // compare observed flux to model flux to adjust weight + ymodel = modelFunc (NULL, model->params, coord); + + // this test modifies the weight based on deviation from the model flux + Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So)); + + // psMinimizeLMChi2 takes wt = 1/dY^2 + if (source->weight->data.F32[i][j] == 0) { + yErr->data.F32[tmpCnt] = 0.0; + } else { + yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro); + } + tmpCnt++; + } + } + } + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + // PSF model only fits first 4 parameters, FLT 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; + } + } + + // XXX EAM : I've added three types of parameter range checks + // XXX EAM : this requires my new psMinimization functions + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); + psVector *beta_lim = NULL; + psVector *params_min = NULL; + psVector *params_max = NULL; + + // XXX EAM : in this implementation, I pass in the limits with the covar matrix. + // in the SDRS, I define a new psMinimization which will take these in + psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); + modelLimits (&beta_lim, ¶ms_min, ¶ms_max); + for (int i = 0; i < params->n; i++) { + covar->data.F64[0][i] = beta_lim->data.F32[i]; + covar->data.F64[1][i] = params_min->data.F32[i]; + covar->data.F64[2][i] = params_max->data.F32[i]; + } + + psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, 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]); + } + + // XXX EAM: we need to do something (give an error?) if rc is false + // XXX EAM: psMinimizeLMChi2 does not check convergence + + // XXX 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); + + // XXX EAM: save the resulting chisq, nDOF, nIter + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // XXX EAM 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(paramMask); + psFree(x); + psFree(y); + psFree(myMin); + + rc = (onPic && fitStatus); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +bool p_pmSourceAddOrSubModel(psImage *image, + psImage *mask, + pmModel *model, + bool center, + psS32 flag) +{ + 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; + + 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; + psF32 pixelValue; + // XXX: Should you be adding the pixels for the entire subImage, + // or a radius of pixels around it? + + // Convert i/j to imace coord space: + // XXX: Make sure you have col/row order correct. + // XXX EAM : 'center' option changes this + // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2] + 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; + pixelValue = modelFunc (NULL, params, x); + // fprintf (stderr, "%f %f %d %d %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue); + + if (flag == 1) { + pixelValue = -pixelValue; + } + + // XXX: Must figure out how to calculate the image coordinates and + // how to use the boolean "center" flag. + + image->data.F32[i][j]+= pixelValue; + } + } + psFree(x); + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); +} + + + +/****************************************************************************** + *****************************************************************************/ +bool pmSourceAddModel(psImage *image, + psImage *mask, + pmModel *model, + bool center) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +/****************************************************************************** + *****************************************************************************/ +bool pmSourceSubModel(psImage *image, + psImage *mask, + pmModel *model, + bool center) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + + Index: /branches/eam_rel8_b2/psModules/src/objects/pmObjects.h =================================================================== --- /branches/eam_rel8_b2/psModules/src/objects/pmObjects.h (revision 5606) +++ /branches/eam_rel8_b2/psModules/src/objects/pmObjects.h (revision 5606) @@ -0,0 +1,625 @@ +/** @file pmObjects.h + * + * The process of finding, measuring, and classifying astronomical sources on + * images is one of the critical tasks of the IPP or any astronomical software + * system. This file will define structures and functions related to the task + * of source detection and measurement. The elements defined in this section + * are generally low-level components which can be connected together to + * construct a complete object measurement suite. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.2.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-26 02:48:42 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#if !defined(PM_OBJECTS_H) +#define PM_OBJECTS_H + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include "pslib.h" +#include "pmAstrometry.h" +/** + * In the object analysis process, we will use specific mask values to mark the + * image pixels. The following structure defines the relevant mask values. + * + * XXX: This is probably a bad solution: we will want to set mask values + * outside of the PSPHOT code. Perhaps we can set up a registered set of mask + * values with specific meanings that other functions can add to or define? + */ +enum { + PSPHOT_MASK_CLEAR = 0x00, + PSPHOT_MASK_INVALID = 0x01, + PSPHOT_MASK_SATURATED = 0x02, + PSPHOT_MASK_MARKED = 0x08, +} psphotMaskValues; + + +/** pmPeakType + * + * A peak pixel may have several features which may be determined when the + * peak is found or measured. These are specified by the pmPeakType enum. + * PM_PEAK_LONE represents a single pixel which is higher than its 8 immediate + * neighbors. The PM_PEAK_EDGE represents a peak pixel which touching the image + * edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific + * number of neighbors at the same value, within some tolarence: + * + */ +typedef enum { + PM_PEAK_LONE, ///< Isolated peak. + PM_PEAK_EDGE, ///< Peak on edge. + PM_PEAK_FLAT, ///< Peak has equal-value neighbors. + PM_PEAK_UNDEF ///< Undefined. +} pmPeakType; + + +/** pmPeak data structure + * + * A source has the capacity for several types of measurements. The + * simplest measurement of a source is the location and flux of the peak pixel + * associated with the source: + * + */ +typedef struct +{ + int x; ///< X-coordinate of peak pixel. + int y; ///< Y-coordinate of peak pixel. + float counts; ///< Value of peak pixel (above sky?). + pmPeakType class; ///< Description of peak. +} +pmPeak; + + +/** pmMoments data structure + * + * One of the simplest measurements which can be made quickly for an object + * are the object moments. We specify a structure to carry the moment information + * for a specific source: + * + */ +typedef struct +{ + float x; ///< X-coord of centroid. + float y; ///< Y-coord of centroid. + float Sx; ///< x-second moment. + float Sy; ///< y-second moment. + float Sxy; ///< xy cross moment. + float Sum; ///< Pixel sum above sky (background). + float Peak; ///< Peak counts above sky. + float Sky; ///< Sky level (background). + float SN; ///< approx signal-to-noise + int nPixels; ///< Number of pixels used. +} +pmMoments; + + +/** pmPSFClump data structure + * + * A collection of object moment measurements can be used to determine + * approximate object classes. The key to this analysis is the location and + * statistics (in the second-moment plane, + * + */ +typedef struct +{ + float X; + float dX; + float Y; + float dY; +} +pmPSFClump; + +typedef int pmModelType; +// XXX EAM : this must not be hard-coded +// #define PS_MODEL_GAUSS 0 +// #define PS_MODEL_PGAUSS 1 +// #define PS_MODEL_QGAUSS 2 +// #define PS_MODEL_SGAUSS 3 + + +/** pmModel data structure + * + * Every source may have two types of models: a PSF model and a FLT (floating) + * model. The PSF model represents the best fit of the image PSF to the specific + * object. In this case, the PSF-dependent parameters are specified for the + * object by the PSF, not by the fit. The FLT model represents the best fit of + * the given model to the object, with all parameters floating in the fit. + * + */ +typedef struct +{ + pmModelType type; ///< Model to be used. + psVector *params; ///< Paramater values. + psVector *dparams; ///< Parameter errors. + float chisq; ///< Fit chi-squared. + int nDOF; ///< number of degrees of freedom + int nIter; ///< number of iterations to reach min + float radius; ///< fit radius actually used +} +pmModel; + +/** pmSourceType enumeration + * + * A given source may be identified as most-likely to be one of several source + * types. The pmSource entry pmSourceType defines the current best-guess for this + * source. + * + * XXX: The values given below are currently illustrative and will require + * some modification as the source classification code is developed. (TBD) + * + */ +typedef enum { + PM_SOURCE_DEFECT, ///< a cosmic-ray + PM_SOURCE_SATURATED, ///< random saturated pixels + + PM_SOURCE_SATSTAR, ///< a saturated star + PM_SOURCE_PSFSTAR, ///< a PSF star + PM_SOURCE_GOODSTAR, ///< a good-quality star + + PM_SOURCE_POOR_FIT_PSF, ///< poor quality PSF fit + PM_SOURCE_FAIL_FIT_PSF, ///< failed to get a good PSF fit + PM_SOURCE_FAINTSTAR, ///< below S/N cutoff + + PM_SOURCE_GALAXY, ///< an extended object (galaxy) + PM_SOURCE_FAINT_GALAXY, ///< a galaxy below S/N cutoff + PM_SOURCE_DROP_GALAXY, ///< ? + PM_SOURCE_FAIL_FIT_GAL, ///< failed on the galaxy fit + PM_SOURCE_POOR_FIT_GAL, ///< poor quality galaxy fit + + PM_SOURCE_OTHER, ///< unidentified +} pmSourceType; + +/** pmSource data structure + * + * This source has the capacity for several types of measurements. The + * simplest measurement of a source is the location and flux of the peak pixel + * associated with the source: + * + */ +typedef struct +{ + pmPeak *peak; ///< Description of peak pixel. + psImage *pixels; ///< Rectangular region including object pixels. + psImage *weight; ///< Image variance. + psImage *mask; ///< Mask which marks pixels associated with objects. + pmMoments *moments; ///< Basic moments measure for the object. + pmModel *modelPSF; ///< PSF Model fit (parameters and type) + pmModel *modelFLT; ///< FLT (floating) Model fit (parameters and type). + pmSourceType type; ///< Best identification of object. + float apMag; + float fitMag; +} +pmSource; + + +/** pmPeakAlloc() + * + * @return pmPeak* newly allocated pmPeak with all internal pointers set to NULL + */ +pmPeak *pmPeakAlloc( + int x, ///< Row-coordinate in image space + int y, ///< Col-coordinate in image space + float counts, ///< The value of the peak pixel + pmPeakType class ///< The type of peak pixel +); + + +/** pmMomentsAlloc() + * + */ +pmMoments *pmMomentsAlloc(); + + +/** pmModelAlloc() + * + */ +pmModel *pmModelAlloc(pmModelType type); + + +/** pmSourceAlloc() + * + */ +pmSource *pmSourceAlloc(); + + +/** pmFindVectorPeaks() + * + * Find all local peaks in the given vector above the given threshold. A peak + * is defined as any element with a value greater than its two neighbors and with + * a value above the threshold. Two types of special cases must be addressed. + * Equal value elements: If an element has the same value as the following + * element, it is not considered a peak. If an element has the same value as the + * preceding element (but not the following), then it is considered a peak. Note + * that this rule (arbitrarily) identifies flat regions by their trailing edge. + * Edge cases: At start of the vector, the element must be higher than its + * neighbor. At the end of the vector, the element must be higher or equal to its + * neighbor. These two rules again places the peak associated with a flat region + * which touches the image edge at the image edge. The result of this function is + * a vector containing the coordinates (element number) of the detected peaks + * (type psU32). + * + */ +psVector *pmFindVectorPeaks( + const psVector *vector, ///< The input vector (float) + float threshold ///< Threshold above which to find a peak +); + + +/** pmFindImagePeaks() + * + * Find all local peaks in the given image above the given threshold. This + * function should find all row peaks using pmFindVectorPeaks, then test each row + * peak and exclude peaks which are not local peaks. A peak is a local peak if it + * has a higher value than all 8 neighbors. If the peak has the same value as its + * +y neighbor or +x neighbor, it is NOT a local peak. If any other neighbors + * have an equal value, the peak is considered a valid peak. Note two points: + * first, the +x neighbor condition is already enforced by pmFindVectorPeaks. + * Second, these rules have the effect of making flat-topped regions have single + * peaks at the (+x,+y) corner. When selecting the peaks, their type must also be + * set. The result of this function is an array of pmPeak entries. + * + */ +psArray *pmFindImagePeaks( + const psImage *image, ///< The input image where peaks will be found (float) + float threshold ///< Threshold above which to find a peak +); + + +/** pmCullPeaks() + * + * Eliminate peaks from the psList that have a peak value above the given + * maximum, or fall outside the valid region. + * + */ +psList *pmCullPeaks( + psList *peaks, ///< The psList of peaks to be culled + float maxValue, ///< Cull peaks above this value + const psRegion valid ///< Cull peaks otside this psRegion +); + + +/** pmPeaksSubset() + * + * Create a new peaks array, removing certain types of peaks from the input + * array of peaks based on the given criteria. Peaks should be eliminated if they + * have a peak value above the given maximum value limit or if the fall outside + * the valid region. The result of the function is a new array with a reduced + * number of peaks. + * + */ +psArray *pmPeaksSubset( + psArray *peaks, ///< Add comment. + float maxvalue, ///< Add comment. + const psRegion valid ///< Add comment. +); + + +/** pmSourceDefinePixels() + * + * Define psImage subarrays for the source located at coordinates x,y on the + * image set defined by readout. The pixels defined by this operation consist of + * a square window (of full width 2Radius+1) centered on the pixel which contains + * the given coordinate, in the frame of the readout. The window is defined to + * have limits which are valid within the boundary of the readout image, thus if + * the radius would fall outside the image pixels, the subimage is truncated to + * only consist of valid pixels. If readout->mask or readout->weight are not + * NULL, matching subimages are defined for those images as well. This function + * fails if no valid pixels can be defined (x or y less than Radius, for + * example). This function should be used to define a region of interest around a + * source, including both source and sky pixels. + * + * XXX: must code this. + * + */ +// XXX: Uncommenting the pmReadout causes compile errors. +bool pmSourceDefinePixels( + pmSource *mySource, ///< Add comment. + pmReadout *readout, ///< Add comment. + psF32 x, ///< Add comment. + psF32 y, ///< Add comment. + psF32 Radius ///< Add comment. +); + + +/** pmSourceLocalSky() + * + * Measure the local sky in the vicinity of the given source. The Radius + * defines the square aperture in which the moments will be measured. This + * function assumes the source pixels have been defined, and that the value of + * Radius here is smaller than the value of Radius used to define the pixels. The + * annular region not contained within the radius defined here is used to measure + * the local background in the vicinity of the source. The local background + * measurement uses the specified statistic passed in via the statsOptions entry. + * This function allocates the pmMoments structure. The resulting sky is used to + * set the value of the pmMoments.sky element of the provided pmSource structure. + * + */ +bool pmSourceLocalSky( + pmSource *source, ///< The input image (float) + psStatsOptions statsOptions, ///< The statistic used in calculating the background sky + float Radius ///< The inner radius of the square annulus to exclude +); + + +/** pmSourceMoments() + * + * Measure source moments for the given source, using the value of + * source.moments.sky provided as the local background value and the peak + * coordinates as the initial source location. The resulting moment values are + * applied to the source.moments entry, and the source is returned. The moments + * are measured within the given circular radius of the source.peak coordinates. + * The return value indicates the success (TRUE) of the operation. + * + */ +bool pmSourceMoments( + pmSource *source, ///< The input pmSource for which moments will be computed + float radius ///< Use a circle of pixels around the peak +); + + +/** pmSourcePSFClump() + * + * We use the source moments to make an initial, approximate source + * classification, and as part of the information needed to build a PSF model for + * the image. As long as the PSF shape does not vary excessively across the + * image, the sources which are represented by a PSF (the start) will have very + * similar second moments. The function pmSourcePSFClump searches a collection of + * sources with measured moments for a group with moments which are all very + * similar. The function returns a pmPSFClump structure, representing the + * centroid and size of the clump in the sigma_x, sigma_y second-moment plane. + * + * The goal is to identify and characterize the stellar clump within the + * sigma_x, sigma_y second-moment plane. To do this, an image is constructed to + * represent this plane. The units of sigma_x and sigma_y are in image pixels. A + * pixel in this analysis image represents 0.1 pixels in the input image. The + * dimensions of the image need only be 10 pixels. The peak pixel in this image + * (above a threshold of half of the image maximum) is found. The coordinates of + * this peak pixel represent the 2D mode of the sigma_x, sigma_y distribution. + * The sources with sigma_x, sigma_y within 0.2 pixels of this value are then + * * used to calculate the median and standard deviation of the sigma_x, sigma_y + * values. These resulting values are returned via the pmPSFClump structure. + * + * The return value indicates the success (TRUE) of the operation. + * + * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD). + * XXX: Save the clump parameters on the Metadata (TBD) + * + */ +pmPSFClump pmSourcePSFClump( + psArray *source, ///< The input pmSource + psMetadata *metadata ///< Contains classification parameters +); + + +/** pmSourceRoughClass() + * + * Based on the specified data values, make a guess at the source + * classification. The sources are provides as a psArray of pmSource entries. + * Definable parameters needed to make the classification are provided to the + * routine with the psMetadata structure. The rules (in SDRS) refer to values which + * can be extracted from the metadata using the given keywords. Except as noted, + * the data type for these parameters are psF32. + * + */ +bool pmSourceRoughClass( + psArray *source, ///< The input pmSource + psMetadata *metadata, ///< Contains classification parameters + pmPSFClump clump ///< Statistics about the PSF clump +); + + +/** pmSourceModelGuess() + * + * Convert available data to an initial guess for the given model. This + * function allocates a pmModel entry for the pmSource structure based on the + * provided model selection. The method of defining the model parameter guesses + * are specified for each model below. The guess values are placed in the model + * parameters. The function returns TRUE on success or FALSE on failure. + * + */ +pmModel *pmSourceModelGuess( + pmSource *source, ///< The input pmSource + pmModelType model ///< The type of model to be created. +); + + +/** pmContourType + * + * Only one type is defined at present. + * + */ +typedef enum { + PS_CONTOUR_CRUDE, + PS_CONTOUR_UNKNOWN01, + PS_CONTOUR_UNKNOWN02 +} pmContourType; + + +/** pmSourceContour() + * + * Find points in a contour for the given source at the given level. If type + * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak, + * running along each pixel row until the level is crossed, then interpolating to + * the level coordinate for that row. This is done for each row, with the + * starting point determined by the midpoint of the previous row, until the + * starting point has a value below the contour level. The returned contour + * consists of two vectors giving the x and y coordinates of the contour levels. + * This function may be used as part of the model guess inputs. Other contour + * types may be specified in the future for more refined contours (TBD) + * + */ +psArray *pmSourceContour( + pmSource *source, ///< The input pmSource + const psImage *image, ///< The input image (float) (this arg should be removed) + float level, ///< The level of the contour + pmContourType mode ///< Currently this must be PS_CONTOUR_CRUDE +); + + +/** pmSourceFitModel() + * + * Fit the requested model to the specified source. The starting guess for the + * model is given by the input source.model parameter values. The pixels of + * interest are specified by the source.pixelsand source.maskentries. This + * function calls psMinimizeLMChi2() on the image data. The function returns TRUE + * on success or FALSE on failure. + * + */ +bool pmSourceFitModel( + pmSource *source, ///< The input pmSource + pmModel *model, ///< model to be fitted + const bool PSF ///< Treat model as PSF or FLT? +); + + +/** pmModelFitStatus() + * + * This function wraps the call to the model-specific function returned by + * pmModelFitStatusFunc_GetFunction. The model-specific function examines the + * model parameters, parameter errors, Chisq, S/N, and other parameters available + * from model to decide if the particular fit was successful or not. + * + * XXX: Must code this. + * + */ +bool pmModelFitStatus( + pmModel *model ///< Add comment. +); + + +/** pmSourceAddModel() + * + * Add the given source model flux to/from the provided image. The boolean + * option center selects if the source is re-centered to the image center or if + * it is placed at its centroid location. The boolean option sky selects if the + * background sky is applied (TRUE) or not. The pixel range in the target image + * is at most the pixel range specified by the source.pixels image. The success + * status is returned. + * + */ +bool pmSourceAddModel( + psImage *image, ///< The output image (float) + psImage *mask, ///< The image pixel mask (valid == 0) + pmModel *model, ///< The input pmModel + bool center ///< A boolean flag that determines whether pixels are centered +); + + +/** pmSourceSubModel() + * + * Subtract the given source model flux to/from the provided image. The + * boolean option center selects if the source is re-centered to the image center + * or if it is placed at its centroid location. The boolean option sky selects if + * the background sky is applied (TRUE) or not. The pixel range in the target + * image is at most the pixel range specified by the source.pixels image. The + * success status is returned. + * + */ +bool pmSourceSubModel( + psImage *image, ///< The output image (float) + psImage *mask, ///< The image pixel mask (valid == 0) + pmModel *model, ///< The input pmModel + bool center ///< A boolean flag that determines whether pixels are centered +); + + +/** + * + * The function returns both the magnitude of the fit, defined as -2.5log(flux), + * where the flux is integrated under the model, theoretically from a radius of 0 + * to infinity. In practice, we integrate the model beyond 50sigma. The aperture magnitude is + * defined as -2.5log(flux) , where the flux is summed for all pixels which are + * not excluded by the aperture mask. The model flux is calculated by calling the + * model-specific function provided by pmModelFlux_GetFunction. + * + * XXX: must code this. + * + */ +bool pmSourcePhotometry( + float *fitMag, ///< integrated fit magnitude + float *obsMag, ///< aperture flux magnitude + pmModel *model, ///< model used for photometry + psImage *image, ///< image pixels to be used + psImage *mask ///< mask of pixels to ignore +); + + +/** + * + * This function converts the source classification into the closest available + * approximation to the Dophot classification scheme: + * + * PM_SOURCE_DEFECT: 8 + * PM_SOURCE_SATURATED: 8 + * PM_SOURCE_SATSTAR: 10 + * PM_SOURCE_PSFSTAR: 1 + * PM_SOURCE_GOODSTAR: 1 + * PM_SOURCE_POOR_FIT_PSF: 7 + * PM_SOURCE_FAIL_FIT_PSF: 4 + * PM_SOURCE_FAINTSTAR: 4 + * PM_SOURCE_GALAXY: 2 + * PM_SOURCE_FAINT_GALAXY: 2 + * PM_SOURCE_DROP_GALAXY: 2 + * PM_SOURCE_FAIL_FIT_GAL: 2 + * PM_SOURCE_POOR_FIT_GAL: 2 + * PM_SOURCE_OTHER: ? + * + */ +int pmSourceDophotType( + pmSource *source ///< Add comment. +); + + +/** pmSourceSextractType() + * + * This function converts the source classification into the closest available + * approximation to the Sextractor classification scheme. the correspondence is + * not yet defined (TBD) . + * + * XXX: Must code this. + * + */ +int pmSourceSextractType( + pmSource *source ///< Add comment. +); + +/** pmSourceFitModel_v5() + * + * Fit the requested model to the specified source. The starting guess for the + * model is given by the input source.model parameter values. The pixels of + * interest are specified by the source.pixelsand source.maskentries. This + * function calls psMinimizeLMChi2() on the image data. The function returns TRUE + * on success or FALSE on failure. + * + */ +bool pmSourceFitModel_v5( + pmSource *source, ///< The input pmSource + pmModel *model, ///< model to be fitted + const bool PSF ///< Treat model as PSF or FLT? +); + + +/** pmSourceFitModel_v7() + * + * Fit the requested model to the specified source. The starting guess for the + * model is given by the input source.model parameter values. The pixels of + * interest are specified by the source.pixelsand source.maskentries. This + * function calls psMinimizeLMChi2() on the image data. The function returns TRUE + * on success or FALSE on failure. + * + */ +bool pmSourceFitModel_v7( + pmSource *source, ///< The input pmSource + pmModel *model, ///< model to be fitted + const bool PSF ///< Treat model as PSF or FLT? +); + +#endif Index: /branches/eam_rel8_b2/psModules/src/objects/pmPSF.c =================================================================== --- /branches/eam_rel8_b2/psModules/src/objects/pmPSF.c (revision 5606) +++ /branches/eam_rel8_b2/psModules/src/objects/pmPSF.c (revision 5606) @@ -0,0 +1,194 @@ +/** @file pmPSF.c + * + * This file contains typedefs for the Point-Spread Function and prototypes + * for functions that calculate the PSF. + * + * @author EAM, IfA + * + * @version $Revision: 1.1.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-26 02:46:26 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +/*****************************************************************************/ +/* INCLUDE FILES */ +/*****************************************************************************/ + +#include +#include "psLibUtils.h" +#include "pmObjects.h" +#include "pmPSF.h" +#include "pmModelGroup.h" + +/*****************************************************************************/ +/* DEFINE STATEMENTS */ +/*****************************************************************************/ + +/*****************************************************************************/ +/* TYPE DEFINITIONS */ +/*****************************************************************************/ + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ + + + +/***************************************************************************** +pmPSFFree(psf): function to free a pmPSF structure + *****************************************************************************/ +static void pmPSFFree (pmPSF *psf) +{ + + if (psf == NULL) + return; + + psFree (psf->params); + return; +} + + + +/***************************************************************************** +pmPSFAlloc (type): allocate a pmPSF. + NOTE: a PSF always has 4 parameters fewer than the equivalent model. + This is because those 4 parameters are + X-center + Y-center + ???: Sky background value? + ???: Zo? + *****************************************************************************/ +pmPSF *pmPSFAlloc (pmModelType type) +{ + int Nparams; + + pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF)); + + psf->type = type; + psf->chisq = 0.0; + psf->ApResid = 0.0; + psf->dApResid = 0.0; + psf->skyBias = 0.0; + + Nparams = pmModelParameterCount (type); + if (!Nparams) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return(NULL); + } + + psf->params = psArrayAlloc(Nparams - 4); + for (int i = 0; i < psf->params->n; i++) { + // XXX EAM : make polynomial order user-specified? + if (1) { + psf->params->data[i] = psPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD); + } else { + psf->params->data[i] = psPolynomial2DAlloc(0, 0, PS_POLYNOMIAL_ORD); + } + } + + psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree); + return(psf); +} + + + +/***************************************************************************** +pmPSFFromModels (*psf, *models, *mask): build a PSF from a collection of +free-fitted models. The PSF ignores the first 4 (independent) model +parameters and constructs a polynomial fit to the remaining as a function of +image coordinate. + input: an array of pmModels, pre-allocated psf +Note: some of the array entries may be NULL (failed fits); ignore them. + *****************************************************************************/ +bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) +{ + + // construct the fit vectors from the collection of objects + // use the mask to ignore missing fits + psVector *x = psVectorAlloc (models->n, PS_TYPE_F64); + psVector *y = psVectorAlloc (models->n, PS_TYPE_F64); + psVector *z = psVectorAlloc (models->n, PS_TYPE_F64); + psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64); + + for (int i = 0; i < models->n; i++) { + pmModel *model = models->data[i]; + if (model == NULL) + continue; + + // XXX EAM : this is fragile: x and y must be stored consistently at 2,3 + // note that the data is provided in the F64 array + x->data.F64[i] = model->params->data.F32[2]; + y->data.F64[i] = model->params->data.F32[3]; + } + + // we are doing a robust fit. after each pass, we drop points which are + // more deviant than three sigma. + // mask is currently updated for each pass. this is inconsistent? + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + + for (int i = 0; i < psf->params->n; i++) { + for (int j = 0; j < models->n; j++) { + pmModel *model = models->data[j]; + if (model == NULL) + continue; + z->data.F64[j] = model->params->data.F32[i + 4]; + dz->data.F64[j] = 1; + // XXX EAM : need to use actual errors? + // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4 + } + + // XXX EAM : this is the expected API (cycle 7? cycle 8?) + psf->params->data[i] = psVectorClipFitPolynomial2D(psf->params->data[i], stats, mask, 0xff, z, dz, x, y); + + // XXX EAM : drop this when above is implemented... + // psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz); + + // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i); + // psPolynomial2DDump (psf->params->data[i]); + } + + psFree (stats); + psFree (x); + psFree (y); + psFree (z); + psFree (dz); + return (true); +} + + + +/***************************************************************************** +pmModelFromPSF (*modelFLT, *psf): use the model position parameters to +construct a realization of the PSF model at the object coordinates + *****************************************************************************/ +pmModel *pmModelFromPSF (pmModel *modelFLT, pmPSF *psf) +{ + + // need to define the relationship between the modelFLT and modelPSF ? + + // find function used to set the model parameters + pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type); + + // allocate a new pmModel to hold the PSF version + pmModel *modelPSF = pmModelAlloc (psf->type); + + // set model parameters for this source based on PSF information + modelFromPSFFunc (modelPSF, modelFLT, psf); + + return (modelPSF); +} Index: /branches/eam_rel8_b2/psModules/src/objects/pmPSFtry.c =================================================================== --- /branches/eam_rel8_b2/psModules/src/objects/pmPSFtry.c (revision 5606) +++ /branches/eam_rel8_b2/psModules/src/objects/pmPSFtry.c (revision 5606) @@ -0,0 +1,399 @@ +# include +# include "psLibUtils.h" +# include "pmObjects.h" +# include "psModulesUtils.h" +# include "pmPSF.h" +# include "pmPSFtry.h" +# include "pmModelGroup.h" + +// ******** pmPSFtry functions ************************************************** +// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely +// * fitted version of the model, the pmPSF fit to the fitted model parameters, +// * and the PSF fits to the source. It also includes the statistics from the +// * fits, both the individual sources, and the collection + +// free a pmPSFtry structure +static void pmPSFtryFree (pmPSFtry *test) +{ + + if (test == NULL) + return; + + psFree (test->psf); + psFree (test->sources); + psFree (test->modelFLT); + psFree (test->modelPSF); + psFree (test->metric); + psFree (test->fitMag); + psFree (test->mask); + return; +} + +// allocate a pmPSFtry based on the desired sources and the model (identified by name) +pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName) +{ + + pmModelType type; + + pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry)); + + // XXX probably need to increment ref counter + type = pmModelSetType (modelName); + test->psf = pmPSFAlloc (type); + test->sources = psMemCopy(sources); + test->modelFLT = psArrayAlloc (sources->n); + test->modelPSF = psArrayAlloc (sources->n); + test->metric = psVectorAlloc (sources->n, PS_TYPE_F64); + test->fitMag = psVectorAlloc (sources->n, PS_TYPE_F64); + test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); + + for (int i = 0; i < test->modelFLT->n; i++) { + test->mask->data.U8[i] = 0; + test->modelFLT->data[i] = NULL; + test->modelPSF->data[i] = NULL; + test->metric->data.F64[i] = 0; + test->fitMag->data.F64[i] = 0; + } + + psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); + return (test); +} + +// build a pmPSFtry for the given model: +// - fit each source with the free-floating model +// - construct the pmPSF from the collection of models +// - fit each source with the PSF-parameter models +// - measure the pmPSF quality metric (dApResid) + +// sources used in for pmPSFtry may be masked by the analysis +// mask values indicate the reason the source was rejected: + +pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS) +{ + bool status; + float obsMag; + float fitMag; + float x; + float y; + int Nflt = 0; + int Npsf = 0; + + pmPSFtry *try + = pmPSFtryAlloc (sources, modelName); + + // stage 1: fit an independent model (freeModel) to all sources + psTimerStart ("fit"); + for (int i = 0; i < try + ->sources->n; i++) { + + pmSource *source = try + ->sources->data[i]; + pmModel *model = pmSourceModelGuess (source, try + ->psf->type); + x = source->peak->x; + y = source->peak->y; + + // set temporary object mask and fit object + // fit model as FLT, not PSF + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); + status = pmSourceFitModel (source, model, false); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); + + // exclude the poor fits + if (!status) { + try + ->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL; + psFree (model); + continue; + } + try + ->modelFLT->data[i] = model; + Nflt ++; + } + psLogMsg ("psphot.psftry", 4, "fit flt: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); + psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n); + + // make this optional? + // DumpModelFits (try->modelFLT, "modelsFLT.dat"); + + // stage 2: construct a psf (pmPSF) from this collection of model fits + pmPSFFromModels (try + ->psf, try + ->modelFLT, try + ->mask); + + // stage 3: refit with fixed shape parameters + psTimerStart ("fit"); + for (int i = 0; i < try + ->sources->n; i++) { + // masked for: bad model fit, outlier in parameters + if (try + ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; + + pmSource *source = try + ->sources->data[i]; + pmModel *modelFLT = try + ->modelFLT->data[i]; + + // set shape for this model based on PSF + pmModel *modelPSF = pmModelFromPSF (modelFLT, try + ->psf); + x = source->peak->x; + y = source->peak->y; + + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); + status = pmSourceFitModel (source, modelPSF, true); + + // skip poor fits + if (!status) { + try + ->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL; + psFree (modelPSF); + goto next_source; + } + + // otherwise, save the resulting model + try + ->modelPSF->data[i] = modelPSF; + + // XXX : use a different aperture radius from the fit radius? + // XXX : use a different estimator for the local sky? + // XXX : pass 'source' as input? + if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) { + try + ->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; + goto next_source; + } + + try + ->metric->data.F64[i] = obsMag - fitMag; + try + ->fitMag->data.F64[i] = fitMag; + Npsf ++; + +next_source: + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); + + } + psLogMsg ("psphot.psftry", 4, "fit psf: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); + psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n); + + // make this optional + // DumpModelFits (try->modelPSF, "modelsPSF.dat"); + + // XXX this function wants aperture radius for pmSourcePhotometry + pmPSFtryMetric_Alt (try + , RADIUS); + psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n", + modelName, try + ->psf->ApResid, try + ->psf->dApResid, try + ->psf->skyBias); + + return (try + ); +} + + +bool pmPSFtryMetric (pmPSFtry *try + , float RADIUS) +{ + + float dBin; + int nKeep, nSkip; + + // the measured (aperture - fit) magnitudes (dA == try->metric) + // depend on both the true ap-fit (dAo) and the bias in the sky measurement: + // dA = dAo + dsky/flux + // where flux is the flux of the star + // we fit this trend to find the infinite flux aperture correction (dAo), + // the nominal sky bias (dsky), and the error on dAo + // the values of dA are contaminated by stars with close neighbors in the aperture + // we use an outlier rejection to avoid this bias + + // rflux = ten(0.4*fitMag); + psVector *rflux = psVectorAlloc (try + ->sources->n, PS_TYPE_F64); + for (int i = 0; i < try + ->sources->n; i++) { + if (try + ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; + rflux->data.F64[i] = pow(10.0, 0.4*try + ->fitMag->data.F64[i]); + } + + // find min and max of (1/flux): + psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); + psVectorStats (stats, rflux, NULL, try + ->mask, PSFTRY_MASK_ALL); + + // build binned versions of rflux, metric + dBin = (stats->max - stats->min) / 10.0; + psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64); + psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64); + psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8); + psFree (stats); + + psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin); + + // group data in daBin bins, measure lower 50% mean + for (int i = 0; i < daBin->n; i++) { + + psVector *tmp = psVectorAlloc (try + ->sources->n, PS_TYPE_F64); + tmp->n = 0; + + // accumulate data within bin range + for (int j = 0; j < try + ->sources->n; j++) { + // masked for: bad model fit, outlier in parameters + if (try + ->mask->data.U8[j] & PSFTRY_MASK_ALL) continue; + + // skip points with extreme dA values + if (fabs(try + ->metric->data.F64[j]) > 0.5) continue; + + // skip points outside of this bin + if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) + continue; + if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) + continue; + + tmp->data.F64[tmp->n] = try + ->metric->data.F64[j]; + tmp->n ++; + } + + // is this a valid point? + maskB->data.U8[i] = 0; + if (tmp->n < 2) { + maskB->data.U8[i] = 1; + psFree (tmp); + continue; + } + + // dA values are contaminated with low outliers + // measure statistics only on upper 50% of points + // this would be easier if we could sort in reverse: + + psVectorSort (tmp, tmp); + nKeep = 0.5*tmp->n; + nSkip = tmp->n - nKeep; + + psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64); + for (int j = 0; j < tmp2->n; j++) { + tmp2->data.F64[j] = tmp->data.F64[j + nSkip]; + } + + stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + psVectorStats (stats, tmp2, NULL, NULL, 0); + psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian); + + daBin->data.F64[i] = stats->sampleMedian; + + psFree (stats); + psFree (tmp); + psFree (tmp2); + } + + // linear clipped fit to rfBin, daBin + psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD); + psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin); + + psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin); + psVector *daResid = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit); + + stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV); + stats = psVectorStats (stats, daResid, NULL, maskB, 1); + + try + ->psf->ApResid = poly->coeff[0]; + try + ->psf->dApResid = stats->clippedStdev; + try + ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); + + psFree (rflux); + psFree (rfBin); + psFree (daBin); + psFree (maskB); + psFree (daBinFit); + psFree (daResid); + psFree (poly); + psFree (stats); + psFree (fitstat); + + return true; +} + +bool pmPSFtryMetric_Alt (pmPSFtry *try + , float RADIUS) +{ + + // the measured (aperture - fit) magnitudes (dA == try->metric) + // depend on both the true ap-fit (dAo) and the bias in the sky measurement: + // dA = dAo + dsky/flux + // where flux is the flux of the star + // we fit this trend to find the infinite flux aperture correction (dAo), + // the nominal sky bias (dsky), and the error on dAo + // the values of dA are contaminated by stars with close neighbors in the aperture + // we use an outlier rejection to avoid this bias + + // rflux = ten(0.4*fitMag); + psVector *rflux = psVectorAlloc (try + ->sources->n, PS_TYPE_F64); + for (int i = 0; i < try + ->sources->n; i++) { + if (try + ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; + rflux->data.F64[i] = pow(10.0, 0.4*try + ->fitMag->data.F64[i]); + } + + // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + + // XXX EAM + stats->min = 1.0; + stats->max = 3.0; + stats->clipIter = 3; + + // linear clipped fit to rfBin, daBin + psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD); + poly = psVectorClipFitPolynomial1D (poly, stats, try + ->mask, PSFTRY_MASK_ALL, try + ->metric, NULL, rflux); + fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); + + try + ->psf->ApResid = poly->coeff[0]; + try + ->psf->dApResid = stats->sampleStdev; + try + ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); + + FILE *f; + f = fopen ("apresid.dat", "w"); + if (f == NULL) + psAbort ("pmPSFtry", "can't open output file"); + + for (int i = 0; i < try + ->sources->n; i++) { + fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try + ->fitMag->data.F64[i], rflux->data.F64[i], try + ->metric->data.F64[i], try + ->mask->data.U8[i]); + } + fclose (f); + + psFree (rflux); + psFree (poly); + psFree (stats); + + // psFree (daFit); + // psFree (daResid); + + return true; +} Index: /branches/eam_rel8_b2/psModules/src/objects/pmPSFtry.h =================================================================== --- /branches/eam_rel8_b2/psModules/src/objects/pmPSFtry.h (revision 5606) +++ /branches/eam_rel8_b2/psModules/src/objects/pmPSFtry.h (revision 5606) @@ -0,0 +1,124 @@ +/** @file pmPSFtry.h + * + * This file contains code that allows the user to try to fit several + * PSF models to an image. + * + * @author EAM, IfA + * + * @version $Revision: 1.1.6.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-26 02:45:18 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_PSF_TRY_H +# define PM_PSF_TRY_H + + +/** + * + * This structure contains a pointer to the collection of sources which will + * be used to test the PSF model form. It lists the pmModelType type of model + * being tests, and contains an element to store the resulting psf + * representation. In addition, this structure carries the complete collection of + * FLT (floating parameter) and PSF (fixed parameter) model fits to each of the + * sources modelFLT and modelPSF. It also contains a mask which is set by the + * model fitting and psf fitting steps. For each model, the value of the quality + * metric is stored in the vector metric and the fitted instrumental magnitude is + * stored in fitMag. The quality metric for the PSF model is the aperture + * magnitude minus the fitted magnitude for each source. This collection of + * aperture residuals is examined in the analysis process, and a linear trend of + * the residual with the inverse object flux (ie, 100:4his structure contains a + * pointer to the collection of sources which will be used to test the PSF model + * form. It lists the pmModelType type of modmag) is fitted. The result of this + * fit is a measured sky bias (systematic error in the sky measured by the fits), + * an effective infinite-magnitude aperture correction (ApResid), and the scatter + * of the aperture correction for the ensemble of PSF stars (dApResid). The + * ultimate metric to intercompare multiple types of PSF models is the value of + * the aperture correction scatter. + * + * XXX: There are many more members in the SDRS then in the prototype code. + * I stuck with the prototype code. + * + * + */ +typedef struct +{ + pmPSF *psf; ///< Add comment. + psArray *sources; ///< pointers to the original sources + psArray *modelFLT; ///< model fits, floating parameters + psArray *modelPSF; ///< model fits, PSF parameters + psVector *mask; ///< Add comment. + psVector *metric; ///< Add comment. + psVector *fitMag; ///< Add comment. +} +pmPSFtry; + + +/** pmPSFtryMaskValues + * + * The following datatype defines the masks used by the pmPSFtry analysis to + * identify sources which should or should not be included in the analysis. + * + */ +enum { + PSFTRY_MASK_CLEAR = 0x00, ///< Add comment. + PSFTRY_MASK_OUTLIER = 0x01, ///< 1: outlier in psf polynomial fit (provided by psPolynomials) + PSFTRY_MASK_FLT_FAIL = 0x02, ///< 2: flt model failed to converge + PSFTRY_MASK_PSF_FAIL = 0x04, ///< 3: psf model failed to converge + PSFTRY_MASK_BAD_PHOT = 0x08, ///< 4: invalid source photometry + PSFTRY_MASK_ALL = 0x0f, ///< Add comment. +} pmPSFtryMaskValues; + + +/** pmPSFtryAlloc() + * + * Allocate a pmPSFtry data structure. + * + */ +pmPSFtry *pmPSFtryAlloc( + psArray *stars, ///< Add comment. + char *modelName ///< Add comment. +); + + +/** pmPSFtryModel() + * + * This function takes the input collection of sources and performs a complete + * analysis to determine a PSF model of the given type (specified by model name). + * The result is a pmPSFtry with the results of the analysis. + * + */ +pmPSFtry *pmPSFtryModel( + psArray *sources, ///< Add comment. + char *modelName, ///< Add comment. + float radius ///< Add comment. +); + + +/** pmPSFtryMetric() + * + * This function is used to measure the PSF model metric for the set of + * results contained in the pmPSFtry structure. + * + */ +bool pmPSFtryMetric( + pmPSFtry *try + , ///< Add comment. + float RADIUS ///< Add comment. +); + +/** pmPSFtryMetric_Alt() + * + * This function is used to measure the PSF model metric for the set of + * results contained in the pmPSFtry structure (alternative implementation). + * + */ +bool pmPSFtryMetric_Alt( + pmPSFtry *try + , ///< Add comment. + float RADIUS ///< Add comment. +); + +# endif