Index: trunk/psModules/src/objects/pmObjects.c =================================================================== --- trunk/psModules/src/objects/pmObjects.c (revision 6511) +++ trunk/psModules/src/objects/pmObjects.c (revision 6873) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-04 01:01:33 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -19,2012 +19,3 @@ #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__); -} - -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__); -} - -/****************************************************************************** -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); - tmpVector->n = tmpVector->nalloc; - 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); -} - - -/****************************************************************************** -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); -} - -// 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); -} - -/****************************************************************************** -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); -} - -/****************************************************************************** -pmModelAlloc(): Allocate the pmModel structure, along with its parameters, -and initialize the type member. Initialize the params to 0.0. -XXX EAM: simplifying code with pmModelParameterCount -*****************************************************************************/ -pmModel *pmModelAlloc(pmModelType type) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel)); - - tmp->type = type; - tmp->chisq = 0.0; - tmp->nIter = 0; - tmp->radius = 0; - tmp->status = PM_MODEL_UNTRIED; - - psS32 Nparams = pmModelParameterCount(type); - if (Nparams == 0) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return(NULL); - } - - tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32); - tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32); - tmp->params->n = tmp->params->nalloc; - tmp->dparams->n = tmp->dparams->nalloc; - - 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 -*****************************************************************************/ - -/****************************************************************************** -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->n = 1; - 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); - tmpVector->n = tmpVector->nalloc; - count = 0; - - // - // Determine if first pixel is a peak - // - if ((vector->data.F32[0] > vector->data.F32[1]) && - (vector->data.F32[0] > threshold)) { - tmpVector->data.U32[count++] = 0; - } - - // - // Determine if interior pixels are peaks - // - for (psU32 i = 1; i < (n-1) ; i++) { - if ((vector->data.F32[i] > vector->data.F32[i-1]) && - (vector->data.F32[i] > vector->data.F32[i+1]) && - (vector->data.F32[i] > threshold)) { - tmpVector->data.U32[count++] = i; - } - } - - // - // Determine if last pixel is a peak - // - if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && - (vector->data.F32[n-1] > threshold)) { - tmpVector->data.U32[count++] = n-1; - } - - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmpVector); -} - - -/****************************************************************************** -pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage -above the given threshold. Returns a psArray containing location (x/y value) -of all peaks. - -XXX: I'm not convinced the peak type definition in the SDRS is mutually -exclusive. Some peaks can have multiple types. Edges for sure. Also, a -digonal line with the same value at each point will have a peak for every -point on that line. - -XXX: This does not work if image has either a single row, or a single column. - -XXX: In the output psArray elements, should we use the image row/column offsets? - Currently, we do not. - -XXX: 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); -} - - -/****************************************************************************** -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); -} - -/****************************************************************************** -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; - bool status; - - psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX"); - if (!status) - SX_MAX = 10.0; - psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX"); - if (!status) - SY_MAX = 10.0; - - // construct a sigma-plane image - // psImageAlloc does zero the data - splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32); - for (int i = 0; i < splane->numRows; i++) - { - memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); - } - - // place the sources in the sigma-plane image (ignore 0,0 values?) - for (psS32 i = 0 ; i < sources->n ; i++) - { - pmSource *tmpSrc = (pmSource *) sources->data[i]; - if (tmpSrc == NULL) { - continue; - } - if (tmpSrc->moments == NULL) { - continue; - } - - // Sx,Sy are limited at 0. a peak at 0,0 is artificial - if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { - continue; - } - - // for the moment, force splane dimensions to be 10x10 image pix - binX = tmpSrc->moments->Sx/SCALE; - if (binX < 0) - continue; - if (binX >= splane->numCols) - continue; - - binY = tmpSrc->moments->Sy/SCALE; - if (binY < 0) - continue; - if (binY >= splane->numRows) - continue; - - splane->data.F32[binY][binX] += 1.0; - } - - // find the peak in this image - stats = psStatsAlloc (PS_STAT_MAX); - stats = psImageStats (stats, splane, NULL, 0); - peaks = pmFindImagePeaks (splane, stats[0].max / 2); - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); - - psFree (splane); - psFree (stats); - - } - // XXX EAM : possible errors: - // 1) no peak in splane - // 2) no significant peak in splane - - // measure statistics on Sx, Sy if Sx, Sy within range of clump - { - pmPeak *clump; - psF32 minSx, maxSx; - psF32 minSy, maxSy; - psVector *tmpSx = NULL; - psVector *tmpSy = NULL; - psStats *stats = NULL; - - // XXX EAM : this lets us takes the single highest peak - psArraySort (peaks, pmComparePeakDescend); - clump = peaks->data[0]; - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); - - // define section window for clump - minSx = clump->x * SCALE - 0.2; - maxSx = clump->x * SCALE + 0.2; - minSy = clump->y * SCALE - 0.2; - maxSy = clump->y * SCALE + 0.2; - - tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); - tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); - tmpSx->n = 0; - tmpSy->n = 0; - - // XXX clip sources based on flux? - // create vectors with Sx, Sy values in window - for (psS32 i = 0 ; i < sources->n ; i++) - { - pmSource *tmpSrc = (pmSource *) sources->data[i]; - - if (tmpSrc->moments->Sx < minSx) - continue; - if (tmpSrc->moments->Sx > maxSx) - continue; - if (tmpSrc->moments->Sy < minSy) - continue; - if (tmpSrc->moments->Sy > maxSy) - continue; - tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; - tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; - tmpSx->n++; - tmpSy->n++; - if (tmpSx->n == tmpSx->nalloc) { - psVectorRealloc (tmpSx, tmpSx->nalloc + 100); - psVectorRealloc (tmpSy, tmpSy->nalloc + 100); - } - } - - // measures stats of Sx, Sy - stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - - stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); - psfClump.X = stats->clippedMean; - psfClump.dX = stats->clippedStdev; - - stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); - psfClump.Y = stats->clippedMean; - psfClump.dY = stats->clippedStdev; - - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); - // these values should be pushed on the metadata somewhere - - psFree (stats); - psFree (peaks); - psFree (tmpSx); - psFree (tmpSy); - } - - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return (psfClump); -} - -/****************************************************************************** -pmSourceRoughClass(source, metadata): make a guess at the source -classification. - -XXX: push the clump info into the metadata? - -XXX: How can this function ever return FALSE? - -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_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(radius, 0.0, 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->n = x->nalloc; - 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); -} - -/****************************************************************************** -pmSourceContour(src, img, level, mode): For an input subImage, and model, this -routine returns a psArray of coordinates that evaluate to the specified level. - -XXX: Probably should remove the "image" argument. -XXX: What type should the output coordinate vectors consist of? col,row? -XXX: Why a pmArray output? -XXX: doex x,y correspond with col,row or row/col? -XXX: What is mode? -XXX: The top, bottom of the contour is not correctly determined. -XXX EAM : this function is using the model for the contour, but it should - be using only the image counts -*****************************************************************************/ -psArray *pmSourceContour(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); - xVec->n = xVec->nalloc; - yVec->n = yVec->nalloc; - // - // 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->n = 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__); - model->status = PM_MODEL_BADARGS; - 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); - x->n = x->nalloc; - y->n = y->nalloc; - yErr->n = yErr->nalloc; - 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); - coord->n = 2; - // 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); - paramMask->n = paramMask->nalloc; - 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"); - - psMinConstrain *constrain = psMinConstrainAlloc(); - constrain->paramMask = paramMask; - fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, - x, y, yErr, modelFunc); - psFree(constrain); - for (int i = 0; i < dparams->n; i++) { - if ((paramMask != NULL) && paramMask->data.U8[i]) - continue; - dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); - } - - // save the resulting chisq, nDOF, nIter - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nDOF = y->n - nParams; - - // get the Gauss-Newton distance for fixed model parameters - if (paramMask != NULL) { - psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); - delta->n = delta->nalloc; - psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc); - for (int i = 0; i < dparams->n; i++) { - if (!paramMask->data.U8[i]) - continue; - dparams->data.F32[i] = delta->data.F64[i]; - } - psFree (delta); - } - - // set the model success or failure status - if (!fitStatus) { - model->status = PM_MODEL_NONCONVERGE; - } else { - model->status = PM_MODEL_SUCCESS; - } - - // models can go insane: reject these - onPic &= (params->data.F32[2] >= source->pixels->col0); - onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); - onPic &= (params->data.F32[3] >= source->pixels->row0); - onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); - if (!onPic) { - model->status = PM_MODEL_OFFIMAGE; - } - - 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__); - model->status = PM_MODEL_BADARGS; - 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); - x->n = x->nalloc; - y->n = y->nalloc; - yErr->n = yErr->nalloc; - 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); - coord->n = 2; - // 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 enhances 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_EAM 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); - paramMask->n = paramMask->nalloc; - 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"); - psMinConstrain *constrain = psMinConstrainAlloc(); - constrain->paramMask = paramMask; - fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, - x, y, yErr, modelFunc); - psFree(constrain); - 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 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); - delta->n = delta->nalloc; - 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]; - } - } - - // set the model success or failure status - if (!fitStatus) { - model->status = PM_MODEL_NONCONVERGE; - } else { - model->status = PM_MODEL_SUCCESS; - } - - // 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); - if (!onPic) { - model->status = PM_MODEL_OFFIMAGE; - } - - 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); -} - -bool p_pmSourceAddOrSubModel(psImage *image, - psImage *mask, - pmModel *model, - bool center, - bool sky, - bool add - ) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - - PS_ASSERT_PTR_NON_NULL(model, false); - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); - - psVector *x = psVectorAlloc(2, PS_TYPE_F32); - x->n = 2; - psVector *params = model->params; - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - psS32 imageCol; - psS32 imageRow; - psF32 skyValue = params->data.F32[0]; - psF32 pixelValue; - - for (psS32 i = 0; i < image->numRows; i++) { - for (psS32 j = 0; j < image->numCols; j++) { - if ((mask != NULL) && mask->data.U8[i][j]) - continue; - - // 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; - - // set the appropriate pixel value for this coordinate - if (sky) { - pixelValue = modelFunc (NULL, params, x); - } else { - pixelValue = modelFunc (NULL, params, x) - skyValue; - } - - - // add or subtract the value - if (add - ) { - image->data.F32[i][j] += pixelValue; - } - else { - image->data.F32[i][j] -= pixelValue; - } - } - } - psFree(x); - psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); - return(true); -} - - - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceAddModel(psImage *image, - psImage *mask, - pmModel *model, - bool center, - bool sky) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true); - psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); - return(rc); -} - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceSubModel(psImage *image, - psImage *mask, - pmModel *model, - bool center, - bool sky) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false); - psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); - return(rc); -} - -bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) -{ - - float obsSum = 0; - float fitSum = 0; - float sky = model->params->data.F32[0]; - - pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); - fitSum = modelFluxFunc (model->params); - - for (int ix = 0; ix < image->numCols; ix++) { - for (int iy = 0; iy < image->numRows; iy++) { - if (mask->data.U8[iy][ix]) - continue; - obsSum += image->data.F32[iy][ix] - sky; - } - } - if (obsSum <= 0) - return false; - if (fitSum <= 0) - return false; - - *fitMag = -2.5*log10(fitSum); - *obsMag = -2.5*log10(obsSum); - return (true); -} -