Index: trunk/psphot/src/pmObjects_EAM.c =================================================================== --- trunk/psphot/src/pmObjects_EAM.c (revision 4949) +++ trunk/psphot/src/pmObjects_EAM.c (revision 4949) @@ -0,0 +1,2123 @@ +/** @file pmObjects.c + * + * This file will ... + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-06 06:43:59 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include "pslib.h" +#include "psConstants.h" +#include "pmObjects_EAM.h" + +/****************************************************************************** +pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members. +*****************************************************************************/ +pmPeak *pmPeakAlloc(psS32 x, + psS32 y, + psF32 counts, + pmPeakType class) +{ + pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak)); + tmp->x = x; + tmp->y = y; + tmp->counts = counts; + tmp->class = class; + + return(tmp); +} + +/****************************************************************************** +pmMomentsAlloc(): Allocate the pmMoments structure and initialize the members +to zero. +*****************************************************************************/ +pmMoments *pmMomentsAlloc() +{ + pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments)); + tmp->x = 0.0; + tmp->y = 0.0; + tmp->Sx = 0.0; + tmp->Sx = 0.0; + tmp->Sxy = 0.0; + tmp->Sum = 0.0; + tmp->Peak = 0.0; + tmp->Sky = 0.0; + tmp->nPixels = 0; + + return(tmp); +} + +static void modelFree(pmModel *tmp) +{ + psFree(tmp->params); + psFree(tmp->dparams); +} + +/****************************************************************************** +pmModelAlloc(): Allocate the pmModel structure, along with its parameters, +and initialize the type member. Initialize the params to 0.0. +XXX EAM: changing params and dparams to psVector +*****************************************************************************/ +pmModel *pmModelAlloc(pmModelType type) +{ + pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel)); + + tmp->type = type; + tmp->chisq = 0.0; + switch (type) { + case PS_MODEL_GAUSS: + tmp->params = psVectorAlloc(7, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(7, PS_TYPE_F32); + break; + case PS_MODEL_PGAUSS: + tmp->params = psVectorAlloc(7, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(7, PS_TYPE_F32); + break; + case PS_MODEL_TWIST_GAUSS: + tmp->params = psVectorAlloc(11, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(11, PS_TYPE_F32); + break; + case PS_MODEL_WAUSS: + tmp->params = psVectorAlloc(9, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(9, PS_TYPE_F32); + break; + case PS_MODEL_SERSIC: + tmp->params = psVectorAlloc(8, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(8, PS_TYPE_F32); + break; + case PS_MODEL_SERSIC_CORE: + tmp->params = psVectorAlloc(12, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(12, PS_TYPE_F32); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return(NULL); + } + + 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); + return(tmp); +} + +/****************************************************************************** +XXX: We don't free pixels and mask since that caused a memory error. +We might need to increase the reference counter and decrease it here. +*****************************************************************************/ +static void sourceFree(pmSource *tmp) +{ + psFree(tmp->peak); + // psFree(tmp->pixels); + // psFree(tmp->mask); + psFree(tmp->moments); + psFree(tmp->modelPSF); + psFree(tmp->modelFLT); +} + +/****************************************************************************** +pmSourceAlloc(): Allocate the pmSource structure and initialize its members +to NULL. +*****************************************************************************/ +pmSource *pmSourceAlloc() +{ + pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource)); + tmp->peak = NULL; + tmp->pixels = NULL; + tmp->mask = NULL; + tmp->moments = NULL; + tmp->modelPSF = NULL; + tmp->modelFLT = NULL; + tmp->type = 0; + psMemSetDeallocator(tmp, (psFreeFunc) sourceFree); + + 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) +{ + 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); + } + 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; + } + + 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? +*****************************************************************************/ +static psVector *getRowVectorFromImage(psImage *image, + psU32 row) +{ + 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]; + } + 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) +{ + pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); + + if (list == NULL) { + list = psArrayAlloc(100); + list->n = 0; + } + psArrayAdd(list, 100, tmpPeak); + + 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) +{ + 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."); + 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."); + } + } + + // + // Exit if this image has a single row. + // + if (image->numRows == 1) { + 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 valid column range."); + } + + } + } + + // + // 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 valid colum range."); + } + } + return(list); +} + +// XXX: Macro this. +static bool isItInThisRegion(const psRegion *valid, + psS32 x, + psS32 y) +{ + + if ((x >= valid->x0) && + (x <= valid->x1) && + (y >= valid->y0) && + (y <= valid->y1)) { + return(true); + } + + 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) +*****************************************************************************/ +psList *pmCullPeaks(psList *peaks, + psF32 maxValue, + const psRegion *valid) +{ + PS_ASSERT_PTR_NON_NULL(peaks, NULL); + // PS_ASSERT_PTR_NON_NULL(valid, 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) || + ((valid != NULL) && + (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) { + psListRemoveData(peaks, (psPtr) tmpPeak); + } + + indexNum++; + tmpListElem = tmpListElem->next; + } + + return(peaks); +} + +// XXX EAM: I changed this to return a new, subset array +// rather than alter the existing one +psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid) +{ + 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 (valid != NULL) { + if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)) + continue; + } + psArrayAdd (output, 200, tmpPeak); + } + 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. +*****************************************************************************/ +pmSource *pmSourceLocalSky(const psImage *image, + const pmPeak *peak, + psStatsOptions statsOptions, + psF32 innerRadius, + psF32 outerRadius) +{ + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); + PS_ASSERT_PTR_NON_NULL(peak, NULL); + PS_FLOAT_COMPARE(0.0, innerRadius, NULL); + PS_FLOAT_COMPARE(innerRadius, outerRadius, NULL); + psS32 innerRadiusS32 = (psS32) innerRadius; + psS32 outerRadiusS32 = (psS32) outerRadius; + + // + // 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). + // + // peak->x,y is guaranteed to be on image + psS32 SubImageCenterRow = peak->y; + psS32 SubImageCenterCol = peak->x; + + // XXX EAM : I added this code to stay on the image. So did George + psS32 SubImageStartRow = PS_MAX(0, SubImageCenterRow - outerRadiusS32); + psS32 SubImageEndRow = PS_MIN(image->numRows - 1, SubImageCenterRow + outerRadiusS32); + psS32 SubImageStartCol = PS_MAX(0, SubImageCenterCol - outerRadiusS32); + psS32 SubImageEndCol = PS_MIN(image->numCols - 1, SubImageCenterCol + outerRadiusS32); + // AnulusWidth == number of pixels width in the annulus. We add one since + // the pixels at the inner AND outher radius are included. + psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32); + // Example: assume an outer/inner radius of 20/10. Then the subimage + // should have width/length of 40. An 18-by-18 interior region will + // be masked. + // printf("pmSourceLocalSky(): innerRadiusS32 is %d\n", innerRadiusS32); + // printf("pmSourceLocalSky(): outerRadiusS32 is %d\n", outerRadiusS32); + // printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth); + + // XXX EAM : these tests should not be needed: we can never hit this error because of above + # if (1) + + if (SubImageStartRow < 0) { + psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", + SubImageStartRow); + return(NULL); + } + if (SubImageEndRow >= image->numRows) { + psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", + SubImageEndRow); + return(NULL); + } + if (SubImageStartCol < 0) { + psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n", + SubImageStartCol); + return(NULL); + } + if (SubImageEndCol >= image->numCols) { + psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", + SubImageEndCol); + return(NULL); + } + # endif + + // + // Grab a subimage of the original image of size (2 * outerRadius). + // + // XXX: Must fix for new psImageSubset + // psImage *subImage = psImageSubset((psImage *) image, + // SubImageStartCol, + // SubImageStartRow, + // SubImageEndCol, + // SubImageEndRow); + // printf("pmSourceLocalSky: subimage width/length is (%d, %d)\n", subImage->numCols, subImage->numRows); + psRegion tmpRegion = psRegionSet(SubImageStartCol, + SubImageEndCol, + SubImageStartRow, + SubImageEndRow); + psImage *subImage = psImageSubset((psImage *) image, tmpRegion); + + psImage *subImageMask = psImageAlloc(subImage->numCols, + subImage->numRows, + PS_TYPE_U8); + + // + // Loop through the subimage mask, initialize mask to 0. + // + for (psS32 row = 0 ; row < subImageMask->numRows; row++) { + for (psS32 col = 0 ; col < subImageMask->numCols; col++) { + subImageMask->data.U8[row][col] = 0; + } + } + + // + // Loop through the subimage, mask off pixels in the inner square. + // XXX this uses a static mask value of 1 + // + for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) { + for (psS32 col = AnulusWidth; col <= (subImageMask->numCols - AnulusWidth) - 1; col++) { + subImageMask->data.U8[row][col] = 1; + } + } + + + // for (psS32 row = 0 ; row < subImage->numRows; row++) { + // for (psS32 col = 0 ; col < subImage->numCols; col++) { + // printf("(%d) ", subImageMask->data.U8[row][col]); + // } + // printf("\n"); + // } + + // + // Allocate the myStats structure, then call psImageStats(), which will + // calculate the specified statistic. + // + psStats *myStats = psStatsAlloc(statsOptions); + myStats = psImageStats(myStats, subImage, subImageMask, 1); + + // + // Create the output mySource, and set appropriate members. + // + pmSource *mySource = pmSourceAlloc(); + mySource->peak = (pmPeak *) peak; + mySource->moments = pmMomentsAlloc(); + psF64 tmpF64; + p_psGetStatValue(myStats, &tmpF64); + mySource->moments->Sky = (psF32) tmpF64; + mySource->pixels = subImage; + mySource->mask = subImageMask; + + // + // Free things. XXX: This should be static memory. + // + psFree(myStats); + + return(mySource); +} + +/****************************************************************************** +bool checkRadius(*peak, radius, x, y): private function which simply +determines if the (x, y) point is within the radius of the specified peak. + +XXX: macro this for performance. +*****************************************************************************/ +static bool checkRadius(pmPeak *peak, + psF32 radius, + psS32 x, + psS32 y) +{ + if (PS_SQR(radius) >= (psF32) (PS_SQR(x - peak->x) + PS_SQR(y - peak->y))) { + return(true); + } + + return(false); +} + +/****************************************************************************** +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) +{ + /// XXX EAM should compare with hypot (x,y) for speed + if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { + return(true); + } + + 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 + +XXX: The peak calculations are done in image coords, not subImage coords. + +XXX: mask values? +*****************************************************************************/ +pmSource *pmSourceMoments(pmSource *source, + psF32 radius) +{ + PS_ASSERT_PTR_NON_NULL(source, NULL); + PS_ASSERT_PTR_NON_NULL(source->peak, NULL); + PS_ASSERT_PTR_NON_NULL(source->pixels, NULL); + PS_FLOAT_COMPARE(0.0, radius, NULL); + + // + // 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 Sum = 0.0; + psF32 peakPixel = -PS_MAX_F32; + psS32 numPixels = 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; + // + // 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. + // + 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] != 0)) { + psS32 imgColCoord = col + source->pixels->col0; + psS32 imgRowCoord = row + source->pixels->row0; + if (checkRadius(source->peak, + radius, + imgColCoord, + imgRowCoord)) { + psF32 xDiff = (psF32) (imgColCoord - source->peak->x); + psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); + psF32 pDiff = source->pixels->data.F32[row][col] - sky; + + Sum+= pDiff; + X1+= xDiff * pDiff; + Y1+= yDiff * pDiff; + XY+= xDiff * yDiff * pDiff; + + X2+= PS_SQR(xDiff) * pDiff; + Y2+= PS_SQR(yDiff) * pDiff; + + if (source->pixels->data.F32[row][col] > peakPixel) { + peakPixel = source->pixels->data.F32[row][col]; + } + 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; + source->moments->x = x + ((psF32) source->peak->x); + source->moments->y = y + ((psF32) source->peak->y); + + source->moments->Sxy = XY/Sum; + 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)); + return(source); + + // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it + # if (0) + // + // second loop: get the difference sums + // + X2 = Y2 = 0; + 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] != 0)) { + psS32 imgColCoord = col + source->pixels->col0; + psS32 imgRowCoord = row + source->pixels->row0; + if (checkRadius(source->peak, + radius, + imgColCoord, + imgRowCoord)) { + psF32 xDiff = (psF32) (imgColCoord - source->peak->x); + psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); + psF32 pDiff = source->pixels->data.F32[row][col] - sky; + + Sum+= pDiff; + X2+= PS_SQR(xDiff - x) * pDiff; + Y2+= PS_SQR(yDiff - y) * pDiff; + } + } + } + } + + // + // second moment X = sqrt (X2/Sum) + // + source->moments->Sx = (X2/Sum); + source->moments->Sy = (Y2/Sum); + return(source); + # endif +} + +// XXX EAM : I used +int pmComparePeakAscend (const void **a, const void **b) +{ + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 diff; + + diff = A->counts - B->counts; + if (diff < FLT_EPSILON) + return (-1); + if (diff > FLT_EPSILON) + return (+1); + return (0); +} + +int pmComparePeakDescend (const void **a, const void **b) +{ + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 diff; + + diff = A->counts - B->counts; + if (diff < FLT_EPSILON) + return (+1); + if (diff > FLT_EPSILON) + return (-1); + return (0); +} + +/****************************************************************************** +pmSourceRoughClass(source, metadata): make a guess at the source +classification. + +XXX: This is not useable code, as of the release date. There remains a fair +bit of coding to be completed. + +XXX: The sigX and sigY stuff in the SDRS is unclear. + +XXX: How can this function ever return FALSE? +*****************************************************************************/ + +# define NPIX 10 +# define SCALE 0.1 + +// XXX I am ignore memory freeing issues (EAM) +bool pmSourceRoughClass(psArray *sources, psMetadata *metadata) +{ + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(metadata, false); + psBool rc = true; + psArray *peaks = NULL; + psF32 clumpX = 0.0; + psF32 clumpDX = 0.0; + psF32 clumpY = 0.0; + psF32 clumpDY = 0.0; + + // find the sigmaX, sigmaY clump + { + psStats *stats = NULL; + psImage *splane = NULL; + int binX, binY; + + // construct a sigma-plane image + splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, 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]; + PS_ASSERT_PTR_NON_NULL(tmpSrc, false); // just skip this one? + PS_ASSERT_PTR_NON_NULL(tmpSrc->moments, false); // just skip this one? + + // 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); + + } + + // 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\n", clump->x, clump->y); + + // 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); + clumpX = stats->clippedMean; + clumpDX = stats->clippedStdev; + + stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); + clumpY = stats->clippedMean; + clumpDY = stats->clippedStdev; + + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", clumpX, clumpY); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", clumpDX, clumpDY); + // these values should be pushed on the metadata somewhere + } + + int Nsat = 0; + int Ngal = 0; + int Nfaint = 0; + int Nstar = 0; + int Npsf = 0; + int Ncr = 0; + psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); + starsn->n = 0; + + // 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; + + // check return status value (do these exist?) + bool status; + psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE"); + psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN"); + psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE"); + + psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); + psF32 FAINT_SN_LIM = psMetadataLookupF32 (&status, metadata, "FAINT_SN_LIM"); + + // saturated object (star or single pixel not distinguished) + if (tmpSrc->moments->Peak > SATURATE) { + tmpSrc->type |= PS_SOURCE_SATURATED; + Nsat ++; + continue; + } + + // too small to be stellar + if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) { + tmpSrc->type |= PS_SOURCE_DEFECT; + Ncr ++; + continue; + } + + // possible galaxy + if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) { + tmpSrc->type |= PS_SOURCE_GALAXY; + Ngal ++; + continue; + } + + // the rest are probable stellar objects + psF32 S = tmpSrc->moments->Sum; + psF32 A = M_PI * sigX * sigY; + psF32 B = tmpSrc->moments->Sky; + psF32 RT = sqrtf(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrtf(GAIN))); + psF32 SN = (S * sqrtf(GAIN) / RT); + + starsn->data.F32[starsn->n] = SN; + starsn->n ++; + Nstar ++; + + // faint star + if (SN < FAINT_SN_LIM) { + tmpSrc->type |= PS_SOURCE_FAINTSTAR; + Nfaint ++; + continue; + } + + // PSF star + if (SN > PSF_SN_LIM) { + tmpSrc->type |= PS_SOURCE_PSFSTAR; + Npsf ++; + continue; + } + + // random type of star + tmpSrc->type |= PS_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); + } + + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nfaint: %3d\n", Nfaint); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal: %3d\n", Ngal); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); + + return(rc); +} + +/****************************************************************************** +pmSourceSetPixelsCircle(source, image, radius) + +XXX: Why boolean output? + +XXX: Why are we checking source->moments for NULL? Should the circle be + centered on the centroid or the peak? + +XXX: The circle will have a diameter of (1+radius). This is different from + the pmSourceSetLocal() function. +*****************************************************************************/ +bool pmSourceSetPixelsCircle(pmSource *source, + const psImage *image, + psF32 radius) +{ + 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 + psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); + psS32 SubImageEndRow = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32); + psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); + psS32 SubImageEndCol = PS_MIN (image->numCols - 1, SubImageCenterCol + radiusS32); + + // XXX EAM : this should not be needed: we can never hit this error + # if (1) + + if (SubImageStartRow < 0) { + psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", + SubImageStartRow); + return(false); + } + if (SubImageEndRow >= image->numRows) { + psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", + SubImageEndRow); + return(false); + } + if (SubImageStartCol < 0) { + psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n", + SubImageStartCol); + return(false); + } + if (SubImageEndCol >= image->numCols) { + psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", + SubImageEndCol); + return(false); + } + # endif + + // 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); + } + // XXX: Must fix this. psImageSubset() has different parameters in latest CVS. + // source->pixels = psImageSubset((psImage *) image, + // 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_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; + } + } + } + return(true); +} + +/****************************************************************************** +pmSourceModelGuess(source, image, model): This function allocates a new +pmModel structure and stores it in the pmSource data structure specified in +the argument list. The model type is specified in the argument list. The +params array in that pmModel structure are allocated, and then set to the +appropriate values. This function returns true if everything was successful. + +XXX: Many of the initial parameters are set to 0.0 since I don't know what +the appropiate initial guesses are. + +XXX: The image argument is redundant. + +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. + +XXX: The source->models member used to be allocated here. Now I assume +->modelPSF should be allocated +*****************************************************************************/ +bool pmSourceModelGuess(pmSource *source, + const psImage *image, + pmModelType model) +{ + 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_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); + if (source->modelPSF != NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->modelPSF was non-NULL; calling psFree(source->modelPSF).\n"); + psFree(source->modelPSF); + } + if (!((model == PS_MODEL_GAUSS) || + (model == PS_MODEL_PGAUSS) || + (model == PS_MODEL_WAUSS) || + (model == PS_MODEL_TWIST_GAUSS) || + (model == PS_MODEL_SERSIC) || + (model == PS_MODEL_SERSIC_CORE))) { + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + return(false); + } + + source->modelPSF = pmModelAlloc(model); + + psVector *params = source->modelPSF->params; + + switch (model) { + case PS_MODEL_GAUSS: + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = sqrt(2.0) / source->moments->Sx; + params->data.F32[5] = sqrt(2.0) / source->moments->Sy; + params->data.F32[6] = source->moments->Sxy; + return(true); + + case PS_MODEL_PGAUSS: + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = sqrt(2.0) / source->moments->Sx; + params->data.F32[5] = sqrt(2.0) / source->moments->Sy; + params->data.F32[6] = source->moments->Sxy; + return(true); + + case PS_MODEL_WAUSS: + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = sqrt(2.0) / source->moments->Sx; + params->data.F32[5] = sqrt(2.0) / source->moments->Sy; + params->data.F32[6] = source->moments->Sxy; + // XXX: What are these? + // source->modelPSF->params[7] = B2; + // source->modelPSF->params[8] = B3; + return(true); + + // XXX EAM : I might drop this model (or rather, replace it) + case PS_MODEL_TWIST_GAUSS: + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + // XXX: What are these? + // params->data.F32[4] = SxInner; + // params->data.F32[5] = SyInner; + // params->data.F32[6] = SxyInner; + // params->data.F32[7] = SxOuter; + // params->data.F32[8] = SyOuter; + // params->data.F32[9] = SxyOuter; + // params->data.F32[10] = N; + return(true); + + case PS_MODEL_SERSIC: + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + params->data.F32[4] = sqrt(2.0) / source->moments->Sx; + params->data.F32[5] = sqrt(2.0) / source->moments->Sy; + params->data.F32[6] = source->moments->Sxy; + // XXX: What are these? + //params->data.F32[7] = Nexp; + return(true); + + case PS_MODEL_SERSIC_CORE: + params->data.F32[0] = source->moments->Sky; + params->data.F32[1] = source->peak->counts - source->moments->Sky; + params->data.F32[2] = source->moments->x; + params->data.F32[3] = source->moments->y; + // XXX: What are these? + // params->data.F32[4] SxInner; + // params->data.F32[5] SyInner; + // params->data.F32[6] SxyInner; + // params->data.F32[7] Zd; + // params->data.F32[8] SxOuter; + // params->data.F32[9] SyOuter; + // params->data.F32[10] = SxyOuter; + // params->data.F32[11] = Nexp; + return(true); + + default: + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + return(false); + } +} + +/****************************************************************************** +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. +*****************************************************************************/ +static psF32 evalModel(pmSource *src, + psU32 row, + psU32 col) +{ + PS_ASSERT_PTR_NON_NULL(src, false); + PS_ASSERT_PTR_NON_NULL(src->modelPSF, false); + PS_ASSERT_PTR_NON_NULL(src->modelPSF->params, false); + + // XXX: The following step will not be necessary if the modelPSF->params + // member is a psVector. Suggest to IfA. + + // XXX EAM: done: modelPSF->params is now a vector + psVector *params = src->modelPSF->params; + + // + // 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 + src->pixels->col0); + x->data.F32[1] = (psF32) (row + src->pixels->row0); + psF32 tmpF; + + switch (src->modelPSF->type) { + case PS_MODEL_GAUSS: + tmpF = pmMinLM_Gauss2D(NULL, params, x); + break; + case PS_MODEL_PGAUSS: + tmpF = pmMinLM_PsuedoGauss2D(NULL, params, x); + break; + case PS_MODEL_TWIST_GAUSS: + tmpF = pmMinLM_TwistGauss2D(NULL, params, x); + break; + case PS_MODEL_WAUSS: + tmpF = pmMinLM_Wauss2D(NULL, params, x); + break; + case PS_MODEL_SERSIC: + tmpF = pmMinLM_Sersic(NULL, params, x); + break; + case PS_MODEL_SERSIC_CORE: + tmpF = pmMinLM_SersicCore(NULL, params, x); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return(NAN); + } + psFree(x); + 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) +{ + // + // 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"); + return(NAN); + } + if (!((0 <= subRow) && (subRow < source->pixels->numRows))) { + psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range"); + return(NAN); + } + + psF32 oldValue = evalModel(source, subRow, subCol); + if (oldValue == level) { + 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 = evalModel(source, subRow, subCol); + if (oldValue == level) { + return((psF32) (subCol + source->pixels->col0)); + } + + if ((newValue <= level) && (level <= oldValue)) { + // This is simple linear interpolation. + return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); + } + + if ((oldValue <= level) && (level <= newValue)) { + // This is simple linear interpolation. + return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); + } + + subCol+=incr; + } + + 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. +*****************************************************************************/ +psArray *pmSourceContour(pmSource *source, + const psImage *image, + psF32 level, + pmContourType mode) +{ + 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->modelPSF, false); + + // + // 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); + 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); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + //printf("HERE 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); + 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); + 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; + return(tmpArray); +} + +#if 0 +static psVector *minLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x); +static psVector *minLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); +static psVector *minLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x); +static psVector *minLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); +static psVector *minLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x); +static psVector *minLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x); +#endif + +// XXX EAM : these are better starting values, but should be available from metadata? +#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 20 +#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 +/****************************************************************************** +pmSourceFitModel(source, image): 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: Probably should remove the "image" argument. +*****************************************************************************/ +bool pmSourceFitModel(pmSource *source, + const psImage *image) +{ + 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->modelPSF, false); + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); + psBool rc; + + // 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++; + } + } + } + + // 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]; + + // XXX EAM : this is approximate: need to apply the gain and rdnoise + yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j])); + tmpCnt++; + } + } + } + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + psVector *params = source->modelPSF->params; + + switch (source->modelPSF->type) { + case PS_MODEL_GAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D); + break; + case PS_MODEL_PGAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_PsuedoGauss2D); + break; + case PS_MODEL_TWIST_GAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D); + break; + case PS_MODEL_WAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D); + break; + case PS_MODEL_SERSIC: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Sersic); + break; + case PS_MODEL_SERSIC_CORE: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + rc = false; + } + // XXX EAM: we need to do something (give an error?) if rc is false + // XXX EAM: save the resulting chisq, nDOF, nIter + source->modelPSF->chisq = myMin->value; + source->modelPSF->nDOF = y->n - params->n; + source->modelPSF->nIter = myMin->iter; + + psFree(x); + psFree(y); + psFree(myMin); + return(rc); +} + +static bool sourceAddOrSubModel(psImage *image, + pmSource *src, + bool center, + psS32 flag) +{ + PS_ASSERT_PTR_NON_NULL(src, false); + PS_ASSERT_PTR_NON_NULL(src->moments, false); + PS_ASSERT_PTR_NON_NULL(src->peak, false); + PS_ASSERT_PTR_NON_NULL(src->pixels, false); + PS_ASSERT_PTR_NON_NULL(src->modelPSF, 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 = src->modelPSF->params; + + for (psS32 i = 0; i < src->pixels->numRows; i++) { + for (psS32 j = 0; j < src->pixels->numCols; j++) { + 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. + psS32 imageRow = i + src->pixels->row0; + psS32 imageCol = j + src->pixels->col0; + + x->data.F32[0] = (float) imageCol; + x->data.F32[1] = (float) imageRow; + switch (src->modelPSF->type) { + case PS_MODEL_GAUSS: + pixelValue = pmMinLM_Gauss2D(NULL, params, x); + break; + case PS_MODEL_PGAUSS: + pixelValue = pmMinLM_PsuedoGauss2D(NULL, params, x); + break; + case PS_MODEL_TWIST_GAUSS: + pixelValue = pmMinLM_TwistGauss2D(NULL, params, x); + break; + case PS_MODEL_WAUSS: + pixelValue = pmMinLM_Wauss2D(NULL, params, x); + break; + case PS_MODEL_SERSIC: + pixelValue = pmMinLM_Sersic(NULL, params, x); + break; + case PS_MODEL_SERSIC_CORE: + pixelValue = pmMinLM_SersicCore(NULL, params, x); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + psFree(x); + return(false); + } + 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[imageRow][imageCol]+= pixelValue; + } + } + psFree(x); + return(true); +} + + + +/****************************************************************************** + *****************************************************************************/ +bool pmSourceAddModel(psImage *image, + pmSource *src, + bool center) +{ + return(sourceAddOrSubModel(image, src, center, 0)); +} + +/****************************************************************************** + *****************************************************************************/ +bool pmSourceSubModel(psImage *image, + pmSource *src, + bool center) +{ + return(sourceAddOrSubModel(image, src, center, 1)); +} + + +// XXX: Put this is psConstants.h +#define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \ +if (VEC1->n != N) { \ + psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ + "psVector %s has size %d, should be %d.", \ + #VEC1, VEC1->n, N); \ + return(RVAL); \ +} + + +/** + all of these object representation functions have the same form : func(*deriv, *params, *x) + + the argument "x" contains a single "x,y" coordinate pair. The function computes the object + model, based on the parameters in "params" at the x,y point specified by *x, and returns the value. + The derivatives are also caculated and returned in the "deriv" argument. parameter error checking is + skipped because speed is most important. +**/ + +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; +*****************************************************************************/ +float pmMinLM_Gauss2D( + psVector *deriv, + const psVector *params, + const psVector *x) +{ + PS_ASSERT_VECTOR_NON_NULL(params, NAN); + PS_ASSERT_VECTOR_NON_NULL(x, NAN); + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[3]; + psF32 px = params->data.F32[4]*X; + psF32 py = params->data.F32[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; + psF32 r = exp(-z); + psF32 q = params->data.F32[1]*r; + psF32 f = q + params->data.F32[0]; + + if (deriv != NULL) { + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); + deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + } + return(f); +} + +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2) / SigmaX; + params->data.F32[5] = sqrt(2) / SigmaY; + params->data.F32[6] = Sxy; +*****************************************************************************/ +float pmMinLM_PsuedoGauss2D( + psVector *deriv, + const psVector *params, + const psVector *x) +{ + PS_ASSERT_VECTOR_NON_NULL(params, NAN); + PS_ASSERT_VECTOR_NON_NULL(x, NAN); + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[3]; + psF32 px = params->data.F32[4]*X; + psF32 py = params->data.F32[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; + psF32 t = 1 + z + 0.5*z*z; + psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */ + psF32 f = params->data.F32[1]*r + params->data.F32[0]; + + if (deriv != NULL) { + // note difference from a pure gaussian: q = params->data.F32[1]*r + psF32 q = params->data.F32[1]*r*r*t; + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); + deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + } + return(f); +} + +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = Sx; + params->data.F32[5] = Sy; + params->data.F32[6] = Sxy; + params->data.F32[7] = B2; + params->data.F32[8] = B3; +*****************************************************************************/ +float pmMinLM_Wauss2D( + psVector *deriv, + const psVector *params, + const psVector *x) +{ + PS_ASSERT_VECTOR_NON_NULL(params, NAN); + PS_ASSERT_VECTOR_NON_NULL(x, NAN); + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[2]; + psF32 px = params->data.F32[4]*X; + psF32 py = params->data.F32[5]*Y; + psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; + psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0); + psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */ + psF32 f = params->data.F32[1]*r + params->data.F32[0]; + + if (deriv != NULL) { + // note difference from gaussian: q = params->data.F32[1]*r + psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); + deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); + deriv->data.F32[4] = -2.0*q*px*X; + deriv->data.F32[5] = -2.0*q*py*Y; + deriv->data.F32[6] = -q*X*Y; + deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; + deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; + // The values of 100 dampen the swing of params->data.F32[7,8] */ + } + return(f); +} + +// XXX: What should these be? +#define FFACTOR 1.0 +#define FSCALE 1.0 +/****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = SxInner; + params->data.F32[5] = SyInner; + params->data.F32[6] = SxyInner; + params->data.F32[7] = SxOuter; + params->data.F32[8] = SyOuter; + params->data.F32[9] = SxyOuter; + params->data.F32[10] = N; +*****************************************************************************/ +float pmMinLM_TwistGauss2D( + psVector *deriv, + const psVector *params, + const psVector *x) +{ + PS_ASSERT_VECTOR_NON_NULL(params, NAN); + PS_ASSERT_VECTOR_NON_NULL(x, NAN); + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[3]; + psF32 px1 = params->data.F32[4]*X; + psF32 py1 = params->data.F32[5]*Y; + psF32 px2 = params->data.F32[7]*X; + psF32 py2 = params->data.F32[8]*Y; + psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y; + psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y; + psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10])); + + psF32 f = params->data.F32[5]*r + params->data.F32[6]; + + if (deriv != NULL) { + psF32 q1 = params->data.F32[5]*PS_SQR(r); + psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y); + deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X); + + // These fudge factors impede the growth of params->data.F32[4] beyond + // params->data.F32[7]. + psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); + psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; + deriv->data.F32[4] = -2.0*q1*px1*X*f2; + + // These fudge factors impede the growth of params->data.F32[5] beyond + // params->data.F32[8]. + f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); + f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; + deriv->data.F32[5] = -2.0*q1*py1*Y*f2; + deriv->data.F32[6] = -q1*X*Y; + deriv->data.F32[7] = -2.0*q2*px2*X; + deriv->data.F32[8] = -2.0*q2*py2*Y; + deriv->data.F32[9] = -q2*X*Y; + deriv->data.F32[10] = -q1*log(z2); + } + + return(f); +} + +/****************************************************************************** + float Sersic() + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = Sx; + params->data.F32[5] = Sy; + params->data.F32[6] = Sxy; + params->data.F32[7] = Nexp; +*****************************************************************************/ +float pmMinLM_Sersic( + psVector *deriv, + const psVector *params, + const psVector *x) +{ + PS_ASSERT_VECTOR_NON_NULL(params, NAN); + PS_ASSERT_VECTOR_NON_NULL(x, NAN); + psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); + return(0.0); +} + +/****************************************************************************** + float SersicBulge() + params->data.F32[0] So; + params->data.F32[1] Zo; + params->data.F32[2] Xo; + params->data.F32[3] Yo; + params->data.F32[4] SxInner; + params->data.F32[5] SyInner; + params->data.F32[6] SxyInner; + params->data.F32[7] Zd; + params->data.F32[8] SxOuter; + params->data.F32[9] SyOuter; + params->data.F32[10] = SxyOuter; + params->data.F32[11] = Nexp; +*****************************************************************************/ +float pmMinLM_SersicCore( + psVector *deriv, + const psVector *params, + const psVector *x) +{ + PS_ASSERT_VECTOR_NON_NULL(params, NAN); + PS_ASSERT_VECTOR_NON_NULL(x, NAN); + psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); + return(0.0); +} Index: trunk/psphot/src/pmObjects_EAM.h =================================================================== --- trunk/psphot/src/pmObjects_EAM.h (revision 4949) +++ trunk/psphot/src/pmObjects_EAM.h (revision 4949) @@ -0,0 +1,154 @@ + +typedef struct +{ + float X; + float dX; + float Y; + float dY; +} pmPSFClump; + +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; + int nPixels; ///< Number of pixels used. +} +pmMoments; + +typedef psS32 pmModelType; + +typedef struct +{ + pmModelType type; ///< Model to be used. + psVector *params; ///< Paramater values. + psVector *dparams; ///< Parameter errors. + float radius; ///< fit radius actually used + float chisq; ///< Fit chi-squared. + int nDOF; ///< number of degrees of freedom + int nIter; ///< number of iterations to reach min +} +pmModel; + +typedef enum { + PS_SOURCE_DEFECT, + PS_SOURCE_SATURATED, + + PS_SOURCE_SATSTAR, + PS_SOURCE_PSFSTAR, + PS_SOURCE_GOODSTAR, + + PS_SOURCE_POOR_FIT_PSF, + PS_SOURCE_FAIL_FIT_PSF, + PS_SOURCE_FAINTSTAR, + + PS_SOURCE_GALAXY, + PS_SOURCE_FAINT_GALAXY, + PS_SOURCE_DROP_GALAXY, + PS_SOURCE_FAIL_FIT_GAL, + PS_SOURCE_POOR_FIT_GAL, + + PS_SOURCE_OTHER, +} pmSourceType; + +typedef struct +{ + pmPeak *peak; ///< Description of peak pixel. + psImage *pixels; ///< Rectangular region including object pixels. + psImage *noise; ///< Mask which marks pixels associated with objects. + 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 Model fit (parameters and type). + pmSourceType type; ///< Best identification of object. +} +pmSource; + +// XXX EAM function to return pmModel +typedef psMinimizeLMChi2Func pmModelFunc; +typedef psF64 (*pmModelFlux)(const psVector *params); +typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); +typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf); +typedef psF64 (*pmModelRadius)(const psVector *params, double flux); + +typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max); +typedef bool (*pmModelFitStatusFunc)(pmModel *model); + +// XXX EAM structure to carry model options +typedef struct { + char *name; + int nParams; + pmModelFunc modelFunc; + pmModelFlux modelFlux; + pmModelRadius modelRadius; + pmModelLimits modelLimits; + pmModelGuessFunc modelGuessFunc; + pmModelFromPSFFunc modelFromPSFFunc; + pmModelFitStatusFunc modelFitStatusFunc; +} pmModelGroup; + +// XXX EAM : strucures to define elliptical shape parameters +typedef struct { + double major; + double minor; + double theta; +} EllipseAxes; + +typedef struct { + double x2; + double y2; + double xy; +} EllipseMoments; + +typedef struct { + double sx; + double sy; + double sxy; +} EllipseShape; + +EllipseAxes EllipseMomentsToAxes (EllipseMoments moments); +EllipseShape EllipseAxesToShape (EllipseAxes axes); +EllipseAxes EllipseShapeToAxes (EllipseShape shape); + +/****************************************************************************** +pmSourcePSFClump(pmArray *source, psMetaDeta *metadata): find the source PSF clump + *****************************************************************************/ +pmPSFClump pmSourcePSFClump(psArray *source, ///< The input psSource + psMetadata *metadata ///< Contains classification parameters + ); + +/****************************************************************************** +pmSourceRoughClass(pmArray *source, psMetaDeta *metadata): make a guess at the +source classification. + *****************************************************************************/ +bool pmSourceRoughClass(psArray *source, ///< The input psSource + psMetadata *metadata, ///< Contains classification parameters + pmPSFClump clump ///< Statistics about the PSF clump + ); + +pmModel *pmSourceModelGuess(pmSource *source, ///< The input psSource + pmModelType model ///< The type of model to be created. + ); + +// XXX EAM : added utility functions +pmModelType pmModelSetType (char *name); +pmModelFunc pmModelFunc_GetFunction (pmModelType type); +pmModelFlux pmModelFlux_GetFunction (pmModelType type); +pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type); +pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type); +pmModelRadius pmModelRadius_GetFunction (pmModelType type); +char *pmModelGetType (pmModelType type); +psS32 pmModelParameterCount (pmModelType type); +pmModelLimits pmModelLimits_GetFunction (pmModelType type); +pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type); +int pmCountSatPixels (psImage *image); + +pmMoments *pmMomentsAlloc(); +pmModel *pmModelAlloc(pmModelType type); +pmSource *pmSourceAlloc(); Index: trunk/psphot/src/pmSourceUtils.c =================================================================== --- trunk/psphot/src/pmSourceUtils.c (revision 4946) +++ trunk/psphot/src/pmSourceUtils.c (revision 4949) @@ -1,5 +1,5 @@ # include "psphot.h" -bool pmSourceDefinePixels(psSource *mySource, +bool pmSourceDefinePixels(pmSource *mySource, const psImageData *imdata, psF32 x, @@ -26,5 +26,5 @@ } -bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask) { +bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) { float obsSum = 0; @@ -32,5 +32,5 @@ float sky = model->params->data.F32[0]; - psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type); + pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); fitSum = modelFluxFunc (model->params); @@ -39,5 +39,4 @@ if (mask->data.U8[iy][ix]) continue; obsSum += image->data.F32[iy][ix] - sky; - // fitSum += psModelEval (model, image, ix, iy) - sky; } } Index: trunk/psphot/src/psModulesUtils.c =================================================================== --- trunk/psphot/src/psModulesUtils.c (revision 4946) +++ trunk/psphot/src/psModulesUtils.c (revision 4949) @@ -4,5 +4,5 @@ // this sets and clears bit 0x80 -bool pmSourceLocalSky_EAM (psSource *source, +bool pmSourceLocalSky_EAM (pmSource *source, psStatsOptions statsOptions, psF32 Radius) @@ -36,7 +36,7 @@ # define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 # define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 -bool pmSourceFitModel_EAM(psSource *source, - psModel *model, - const bool PSF) +bool pmSourceFitModel_EAM (pmSource *source, + pmModel *model, + const bool PSF) { PS_PTR_CHECK_NULL(source, false); @@ -55,6 +55,6 @@ // tests below could be conditions (!NULL) - psModelFunc modelFunc = psModelFunc_GetFunction (model->type); - psModelLimits modelLimits = psModelLimits_GetFunction (model->type); + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); psVector *params = model->params; @@ -177,11 +177,11 @@ /****************************************************************************** pmSourceMoments(source, radius): this function takes a subImage defined in the -psSource data structure, along with the peak location, and determines 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: - psSource - psSource->peak - psSource->pixels + pmSource + pmSource->peak + pmSource->pixels XXX: The peak calculations are done in image coords, not subImage coords. @@ -191,5 +191,5 @@ # define VALID_RADIUS(X,Y) (((R2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) -bool pmSourceMoments_EAM(psSource *source, +bool pmSourceMoments_EAM(pmSource *source, psF32 radius) { @@ -310,9 +310,9 @@ } -bool pmModelFitStatus (psModel *model) { +bool pmModelFitStatus (pmModel *model) { bool status; - psModelFitStatusFunc statusFunc = psModelFitStatusFunc_GetFunction (model->type); + pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type); status = statusFunc (model); Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 4946) +++ trunk/psphot/src/psphot.h (revision 4949) @@ -1,4 +1,6 @@ # include -# include +// # include +// my additions and modifications: +# include "pmObjects_EAM.h" # include # include @@ -7,5 +9,5 @@ # include -# define M_PI 3.14159265358979323846264338328 /* pi */ +// # define M_PI 3.14159265358979323846264338328 /* pi */ typedef struct { @@ -18,5 +20,5 @@ // data to test a given PSF model type typedef struct { - psModelType modelType; + pmModelType modelType; pmPSF *psf; psArray *sources; // pointers to the original sources @@ -61,6 +63,6 @@ void psphotOutput (psImageData *imdata, psMetadata *config, psArray *sources); -bool psphotMarkPSF (psSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE); -bool psphotSubtractPSF (psSource *source); +bool psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE); +bool psphotSubtractPSF (pmSource *source); int psphotSortBySN (const void **a, const void **b); @@ -69,10 +71,10 @@ // psf utilities -pmPSF *pmPSFAlloc (psModelType type); +pmPSF *pmPSFAlloc (pmModelType type); pmPSF_Test *pmPSF_TestAlloc (psArray *stars, char *modelName); pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float radius); bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask); -psModel *psModelFromPSF (psModel *model, pmPSF *psf); -bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask); +pmModel *pmModelFromPSF (pmModel *model, pmPSF *psf); +bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask); bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS); @@ -83,18 +85,18 @@ bool pmSourcesWriteCMF (psImageData *imdata, char *filename, psArray *sources); bool pmSourcesWriteSX (psImageData *imdata, char *filename, psArray *sources); -int pmSourcesDophotType (psSource *source); -bool psPeaksWriteText (psArray *sources, char *filename); -bool psMomentsWriteText (psArray *sources, char *filename); -bool psModelWritePSFs (psArray *sources, char *filename); -bool psModelWriteFLTs (psArray *sources, char *filename); -bool psModelWriteNULLs (psArray *sources, char *filename); +int pmSourcesDophotType (pmSource *source); +bool pmPeaksWriteText (psArray *sources, char *filename); +bool pmMomentsWriteText (psArray *sources, char *filename); +bool pmModelWritePSFs (psArray *sources, char *filename); +bool pmModelWriteFLTs (psArray *sources, char *filename); +bool pmModelWriteNULLs (psArray *sources, char *filename); // psModule extra utilities -bool pmSourceDefinePixels(psSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius); -bool pmSourceLocalSky_EAM (psSource *source, psStatsOptions statsOptions, psF32 Radius); -bool pmSourceFitModel_EAM(psSource *source, psModel *model, const bool PSF); -bool pmSourceMoments_EAM(psSource *source, psF32 radius); -bool pmModelFitStatus (psModel *model); -int pmSourceDophotType (psSource *source); +bool pmSourceDefinePixels(pmSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius); +bool pmSourceLocalSky_EAM (pmSource *source, psStatsOptions statsOptions, psF32 Radius); +bool pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF); +bool pmSourceMoments_EAM(pmSource *source, psF32 radius); +bool pmModelFitStatus (pmModel *model); +int pmSourceDophotType (pmSource *source); // minimize @@ -114,14 +116,14 @@ psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name); psVector *psVectorCreate (double lower, double upper, double delta, psElemType type); -void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue); -void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue); -void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue); -void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue); +//void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue); +//void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue); +//void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue); +//void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue); psVector *psGetRowVectorFromImage(psImage *image, psU32 row); // basic image functions bool psImageInit (psImage *image,...); -void psImageSmooth (psImage *image, float sigma, float Nsigma); -psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in); +//void psImageSmooth (psImage *image, float sigma, float Nsigma); +//psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in); psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius); int psphotSaveImage (psMetadata *header, psImage *image, char *filename); @@ -153,6 +155,2 @@ int fs_usage (); bool DumpImage (psImage *image, char *filename); - -// deprecated -// void pmSourceMaskRegion (psSource *source, psRegion *region); -//psVector *p_psGetRowVectorFromImage(psImage *image, psU32 row); Index: trunk/psphot/src/psphotApplyPSF.c =================================================================== --- trunk/psphot/src/psphotApplyPSF.c (revision 4946) +++ trunk/psphot/src/psphotApplyPSF.c (revision 4949) @@ -31,9 +31,9 @@ // this function specifies the radius at this the model hits the given flux - psModelRadius modelRadius = psModelRadius_GetFunction (psf->type); + pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type); for (int i = 0; i < sources->n; i++) { - psSource *source = sources->data[i]; + pmSource *source = sources->data[i]; // skip non-astronomical objects (very likely defects) @@ -43,8 +43,8 @@ // use the source moments, etc to guess basic model parameters - psModel *modelFLT = pmSourceModelGuess (source, psf->type); + pmModel *modelFLT = pmSourceModelGuess (source, psf->type); // set PSF parameters for this model - psModel *model = psModelFromPSF (modelFLT, psf); + pmModel *model = pmModelFromPSF (modelFLT, psf); x = model->params->data.F32[2]; y = model->params->data.F32[3]; @@ -66,7 +66,7 @@ // fit PSF model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, model->radius, OR, PSPHOT_MASK_KEEP); + psImageKeepCircle (source->mask, x, y, model->radius, "or", PSPHOT_MASK_KEEP); status = pmSourceFitModel (source, model, true); - psImageKeepCircle (source->mask, x, y, model->radius, AND, ~PSPHOT_MASK_KEEP); + psImageKeepCircle (source->mask, x, y, model->radius, "and", ~PSPHOT_MASK_KEEP); if (!status || (model->params->data.F32[1] < 0)) { psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter); Index: trunk/psphot/src/psphotChoosePSF.c =================================================================== --- trunk/psphot/src/psphotChoosePSF.c (revision 4946) +++ trunk/psphot/src/psphotChoosePSF.c (revision 4949) @@ -17,5 +17,5 @@ // select the candidate PSF stars (pointers to original sources) for (int i = 0; i < sources->n; i++) { - psSource *source = sources->data[i]; + pmSource *source = sources->data[i]; if (source->type != PS_SOURCE_PSFSTAR) continue; psArrayAdd (stars, 200, source); @@ -35,6 +35,5 @@ // set up an array to store the test results - psArray *tests = psArrayAlloc (list->size); - // XXX EAM : new value : psArray *tests = psArrayAlloc (list->n); + psArray *tests = psArrayAlloc (list->n); // test each model option listed in config @@ -68,5 +67,5 @@ // keep only the selected test: test = tests->data[bestN]; - modelName = psModelGetType (test->modelType); + modelName = pmModelGetType (test->modelType); psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, test->ApResid, test->dApResid); @@ -77,5 +76,5 @@ // set source->model based on best psf model fit for (int i = 0; i < test->sources->n; i++) { - psSource *source = test->sources->data[i]; + pmSource *source = test->sources->data[i]; // drop masked sources from PSFSTAR list if (test->mask->data.U8[i]) { Index: trunk/psphot/src/psphotFitGalaxies.c =================================================================== --- trunk/psphot/src/psphotFitGalaxies.c (revision 4946) +++ trunk/psphot/src/psphotFitGalaxies.c (revision 4949) @@ -21,11 +21,11 @@ char *modelName = psMetadataLookupPtr (&status, config, "GAL_MODEL"); - psModelType modelType = psModelSetType (modelName); - psModelRadius modelRadius = psModelRadius_GetFunction (modelType); + pmModelType modelType = pmModelSetType (modelName); + pmModelRadius modelRadius = pmModelRadius_GetFunction (modelType); psTimerStart ("psphot"); for (int i = 0; i < sources->n; i++) { - psSource *source = sources->data[i]; + pmSource *source = sources->data[i]; // sources which should not be fitted @@ -50,5 +50,5 @@ // recalculate the source moments using the larger galaxy moments radius - status = pmSourceMoments_EAM (source, GAL_MOMENTS_RAD); + status = pmSourceMoments (source, GAL_MOMENTS_RAD); if (!status) { source->type = PS_SOURCE_DROP_GALAXY; // better choice? @@ -57,5 +57,5 @@ // use the source moments, etc to guess basic model parameters - psModel *model = pmSourceModelGuess (source, modelType); + pmModel *model = pmSourceModelGuess (source, modelType); x = model->params->data.F32[2]; y = model->params->data.F32[3]; @@ -75,5 +75,5 @@ // fit FLT (not PSF) model (set/unset the pixel mask) psImageKeepCircle (source->mask, x, y, model->radius, OR, PSPHOT_MASK_KEEP); - status = pmSourceFitModel_EAM (source, model, false); + status = pmSourceFitModel (source, model, false); psImageKeepCircle (source->mask, x, y, model->radius, AND, ~PSPHOT_MASK_KEEP); if (!status) { Index: trunk/psphot/src/psphotMarkPSF.c =================================================================== --- trunk/psphot/src/psphotMarkPSF.c (revision 4946) +++ trunk/psphot/src/psphotMarkPSF.c (revision 4949) @@ -24,5 +24,5 @@ # define MIN_DS 0.01 -bool psphotMarkPSF (psSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE) +bool psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE) { int keep; Index: trunk/psphot/src/psphotOutput.c =================================================================== --- trunk/psphot/src/psphotOutput.c (revision 4946) +++ trunk/psphot/src/psphotOutput.c (revision 4949) @@ -52,14 +52,14 @@ sprintf (name, "%s.psf.dat", filename); - psModelWritePSFs (sources, name); + pmModelWritePSFs (sources, name); sprintf (name, "%s.flt.dat", filename); - psModelWriteFLTs (sources, name); + pmModelWriteFLTs (sources, name); sprintf (name, "%s.nul.dat", filename); - psModelWriteNULLs (sources, name); + pmModelWriteNULLs (sources, name); sprintf (name, "%s.mnt.dat", filename); - psMomentsWriteText (sources, name); + pmMomentsWriteText (sources, name); psFree (name); @@ -71,5 +71,5 @@ int i, type, status; - psModel *model; + pmModel *model; psF32 *PAR, *dPAR; float sky, dmag, apMag, fitMag; @@ -85,5 +85,5 @@ // write sources with models first for (i = 0; i < sources->n; i++) { - psSource *source = (psSource *) sources->data[i]; + pmSource *source = (pmSource *) sources->data[i]; // use the correct model @@ -162,5 +162,5 @@ int i, status, flags; - psModel *model; + pmModel *model; psF32 *PAR, *dPAR; float sky, dmag, apMag, fitMag; @@ -176,5 +176,5 @@ // write sources with models first for (i = 0; i < sources->n; i++) { - psSource *source = (psSource *) sources->data[i]; + pmSource *source = (pmSource *) sources->data[i]; flags = 0; @@ -260,5 +260,5 @@ int i, type, status; - psModel *model; + pmModel *model; psF32 *PAR, *dPAR; float sky, dmag, apMag, fitMag, lsky, Theta; @@ -282,5 +282,5 @@ // write sources with models first for (i = 0; i < sources->n; i++) { - psSource *source = (psSource *) sources->data[i]; + pmSource *source = (pmSource *) sources->data[i]; // use the correct model @@ -395,18 +395,18 @@ // write the moments to an output file -bool psMomentsWriteText (psArray *sources, char *filename) { +bool pmMomentsWriteText (psArray *sources, char *filename) { int i; FILE *f; - psSource *source; + pmSource *source; f = fopen (filename, "w"); if (f == NULL) { - psLogMsg ("psMomentsWriteText", 3, "can't open output file for moments%s\n", filename); - return false; - } - - for (i = 0; i < sources->n; i++) { - source = (psSource *) sources->data[i]; + psLogMsg ("pmMomentsWriteText", 3, "can't open output file for moments%s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + source = (pmSource *) sources->data[i]; if (source->moments == NULL) continue; fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %8.1f %7.1f %7.1f %7.1f %4d %2d\n", @@ -423,5 +423,5 @@ // dump the sources to an output file -bool psModelWritePSFs (psArray *sources, char *filename) { +bool pmModelWritePSFs (psArray *sources, char *filename) { double dP, flux; @@ -430,9 +430,9 @@ psVector *params; psVector *dparams; - psModel *model; + pmModel *model; f = fopen (filename, "w"); if (f == NULL) { - psLogMsg ("psModelWritePSFs", 3, "can't open output file for moments%s\n", filename); + psLogMsg ("pmModelWritePSFs", 3, "can't open output file for moments%s\n", filename); return false; } @@ -440,6 +440,6 @@ // write sources with models first for (i = 0; i < sources->n; i++) { - psSource *source = (psSource *) sources->data[i]; - model = (psModel *) source->modelPSF; + pmSource *source = (pmSource *) sources->data[i]; + model = (pmModel *) source->modelPSF; if (model == NULL) continue; @@ -459,5 +459,5 @@ dP = sqrt (dP); - psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type); + pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); flux = modelFluxFunc (params); @@ -489,5 +489,5 @@ // dump the sources to an output file -bool psModelWriteFLTs (psArray *sources, char *filename) { +bool pmModelWriteFLTs (psArray *sources, char *filename) { double dP; @@ -496,9 +496,9 @@ psVector *params; psVector *dparams; - psModel *model; + pmModel *model; f = fopen (filename, "w"); if (f == NULL) { - psLogMsg ("psModelWriteFLTs", 3, "can't open output file for moments%s\n", filename); + psLogMsg ("pmModelWriteFLTs", 3, "can't open output file for moments%s\n", filename); return false; } @@ -506,6 +506,6 @@ // write sources with models first for (i = 0; i < sources->n; i++) { - psSource *source = (psSource *) sources->data[i]; - model = (psModel *) source->modelFLT; + pmSource *source = (pmSource *) sources->data[i]; + model = (pmModel *) source->modelFLT; if (model == NULL) continue; if (source->type == PS_SOURCE_GALAXY) goto valid; @@ -551,5 +551,5 @@ // dump the sources to an output file -bool psModelWriteNULLs (psArray *sources, char *filename) +bool pmModelWriteNULLs (psArray *sources, char *filename) { @@ -563,9 +563,9 @@ } - psMoments *empty = pmMomentsAlloc (); + pmMoments *empty = pmMomentsAlloc (); // write sources with models first for (i = 0; i < sources->n; i++) { - psSource *source = (psSource *) sources->data[i]; + pmSource *source = (pmSource *) sources->data[i]; // skip these sources (in PSF or FLT) @@ -598,5 +598,5 @@ // translations between psphot object types and dophot object types -int pmSourceDophotType (psSource *source) { +int pmSourceDophotType (pmSource *source) { switch (source->type) { Index: trunk/psphot/src/psphotSetup.c =================================================================== --- trunk/psphot/src/psphotSetup.c (revision 4946) +++ trunk/psphot/src/psphotSetup.c (revision 4949) @@ -73,8 +73,8 @@ float YMIN = psMetadataLookupF32 (&status, config, "YMIN"); float YMAX = psMetadataLookupF32 (&status, config, "YMAX"); - psRegion *keep = psRegionAlloc (XMIN, XMAX, YMIN, YMAX); - keep = psRegionForImage (keep, image, keep); - psImageKeepRegion (mask, keep, OR, PSPHOT_MASK_INVALID); - psFree (keep); + psRegion keep = psRegionSet (XMIN, XMAX, YMIN, YMAX); + keep = psRegionForImage (image, &keep); + psImageKeepRegion (mask, &keep, OR, PSPHOT_MASK_INVALID); + // XXX EAM : psRegionForImage should not take psRegion *keep // mask the saturated pixels Index: trunk/psphot/src/psphotSortBySN.c =================================================================== --- trunk/psphot/src/psphotSortBySN.c (revision 4946) +++ trunk/psphot/src/psphotSortBySN.c (revision 4949) @@ -4,6 +4,6 @@ int psphotSortBySN (const void **a, const void **b) { - psSource *A = *(psSource **)a; - psSource *B = *(psSource **)b; + pmSource *A = *(pmSource **)a; + pmSource *B = *(pmSource **)b; psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN; Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 4946) +++ trunk/psphot/src/psphotSourceStats.c (revision 4949) @@ -19,6 +19,6 @@ // create a new source, add peak - psSource *source = pmSourceAlloc(); - source->peak = (psPeak *)psMemCopy(peaks->data[i]); + pmSource *source = pmSourceAlloc(); + source->peak = (pmPeak *)psMemCopy(peaks->data[i]); // allocate image, noise, mask arrays for each peak (square of radius OUTER) Index: trunk/psphot/src/psphotSubtractPSF.c =================================================================== --- trunk/psphot/src/psphotSubtractPSF.c (revision 4946) +++ trunk/psphot/src/psphotSubtractPSF.c (revision 4949) @@ -4,8 +4,8 @@ // need to control application of this with some thresholding -bool psphotSubtractPSF (psSource *source) +bool psphotSubtractPSF (pmSource *source) { float sky; - psModel *model; + pmModel *model; psImage *pixels; Index: trunk/psphot/src/pspsf.c =================================================================== --- trunk/psphot/src/pspsf.c (revision 4946) +++ trunk/psphot/src/pspsf.c (revision 4949) @@ -10,5 +10,5 @@ // a PSF always has 4 parameters fewer than the equivalent model -pmPSF *pmPSFAlloc (psModelType type) { +pmPSF *pmPSFAlloc (pmModelType type) { int Nparams; @@ -19,7 +19,7 @@ psf->chisq = 0.0; - Nparams = psModelParameterCount (type); + Nparams = pmModelParameterCount (type); if (!Nparams) { - psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); return(NULL); } @@ -54,5 +54,5 @@ // XXX probably need to increment ref counter - test->modelType = psModelSetType (modelName); + test->modelType = pmModelSetType (modelName); test->psf = pmPSFAlloc (test->modelType); test->sources = psMemCopy(sources); @@ -99,6 +99,6 @@ for (int i = 0; i < test->sources->n; i++) { - psSource *source = test->sources->data[i]; - psModel *model = pmSourceModelGuess (source, test->modelType); + pmSource *source = test->sources->data[i]; + pmModel *model = pmSourceModelGuess (source, test->modelType); x = source->peak->x; y = source->peak->y; @@ -134,9 +134,9 @@ if (test->mask->data.U8[i]) continue; - psSource *source = test->sources->data[i]; - psModel *modelFLT = test->modelFLT->data[i]; + pmSource *source = test->sources->data[i]; + pmModel *modelFLT = test->modelFLT->data[i]; // set shape for this model based on PSF - psModel *modelPSF = psModelFromPSF (modelFLT, test->psf); + pmModel *modelPSF = pmModelFromPSF (modelFLT, test->psf); x = source->peak->x; y = source->peak->y; @@ -185,5 +185,5 @@ } -// input: an array of psModels, pre-allocated psf +// input: an array of pmModels, pre-allocated psf // some of the array entries may be NULL, ignore them bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) { @@ -197,5 +197,5 @@ for (int i = 0; i < models->n; i++) { - psModel *model = models->data[i]; + pmModel *model = models->data[i]; if (model == NULL) continue; @@ -211,5 +211,5 @@ for (int i = 0; i < psf->params->n; i++) { for (int j = 0; j < models->n; j++) { - psModel *model = models->data[j]; + pmModel *model = models->data[j]; if (model == NULL) continue; z->data.F64[j] = model->params->data.F32[i + 4]; @@ -231,13 +231,13 @@ } -psModel *psModelFromPSF (psModel *modelFLT, pmPSF *psf) { +pmModel *pmModelFromPSF (pmModel *modelFLT, pmPSF *psf) { // need to define the relationship between the modelFLT and modelPSF ? // find function used to set the model parameters - psModelFromPSFFunc modelFromPSFFunc = psModelFromPSFFunc_GetFunction (psf->type); + pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type); // do we need a different model for the PSF vs FLT version? - psModel *modelPSF = psModelAlloc (psf->type); + pmModel *modelPSF = pmModelAlloc (psf->type); // set model parameters for this source based on PSF information