Index: trunk/psModules/src/objects/pmObjects.c =================================================================== --- trunk/psModules/src/objects/pmObjects.c (revision 5170) +++ trunk/psModules/src/objects/pmObjects.c (revision 5255) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-28 20:43:52 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-10-10 19:53:40 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -13,10 +13,10 @@ */ -#include -#include +#include +#include +#include #include "pslib.h" -#include "psConstants.h" #include "pmObjects.h" - +#include "pmModelGroup.h" /****************************************************************************** pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members. @@ -27,4 +27,5 @@ pmPeakType class) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak)); tmp->x = x; @@ -33,4 +34,5 @@ tmp->class = class; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return(tmp); } @@ -42,9 +44,10 @@ pmMoments *pmMomentsAlloc() { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments)); tmp->x = 0.0; tmp->y = 0.0; tmp->Sx = 0.0; - tmp->Sx = 0.0; + tmp->Sy = 0.0; tmp->Sxy = 0.0; tmp->Sum = 0.0; @@ -53,4 +56,5 @@ tmp->nPixels = 0; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return(tmp); } @@ -58,6 +62,8 @@ static void modelFree(pmModel *tmp) { + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); psFree(tmp->params); psFree(tmp->dparams); + psTrace(__func__, 4, "---- %s() end ----\n", __func__); } @@ -65,41 +71,22 @@ 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 +XXX EAM: simplifying code with pmModelParameterCount *****************************************************************************/ pmModel *pmModelAlloc(pmModelType type) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel)); tmp->type = type; tmp->chisq = 0.0; - 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: + tmp->nIter = 0; + psS32 Nparams = pmModelParameterCount(type); + if (Nparams == 0) { psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); return(NULL); } + + tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32); for (psS32 i = 0; i < tmp->params->n; i++) { @@ -109,19 +96,22 @@ psMemSetDeallocator(tmp, (psFreeFunc) modelFree); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); 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. +XXX EAM : we can now free these pixels - memory ref is incremented now *****************************************************************************/ static void sourceFree(pmSource *tmp) { + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); psFree(tmp->peak); - // psFree(tmp->pixels); - // psFree(tmp->mask); + psFree(tmp->pixels); + psFree(tmp->weight); + psFree(tmp->mask); psFree(tmp->moments); psFree(tmp->modelPSF); psFree(tmp->modelFLT); + psTrace(__func__, 4, "---- %s() end ----\n", __func__); } @@ -132,7 +122,9 @@ pmSource *pmSourceAlloc() { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource)); tmp->peak = NULL; tmp->pixels = NULL; + tmp->weight = NULL; tmp->mask = NULL; tmp->moments = NULL; @@ -142,4 +134,5 @@ psMemSetDeallocator(tmp, (psFreeFunc) sourceFree); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return(tmp); } @@ -159,4 +152,5 @@ psF32 threshold) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); PS_ASSERT_VECTOR_NON_NULL(vector, NULL); PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL); @@ -177,4 +171,5 @@ tmpVector = psVectorAlloc(0, PS_TYPE_U32); } + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return(tmpVector); } @@ -241,4 +236,5 @@ } + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return(tmpVector); } @@ -248,9 +244,11 @@ psVector containing the specified row of data from the psImage. -XXX: Is there a better way to do this? +XXX: Is there a better way to do this? +XXX EAM: does this really need to alloc a new vector??? *****************************************************************************/ static psVector *getRowVectorFromImage(psImage *image, psU32 row) { + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); PS_ASSERT_IMAGE_NON_NULL(image, NULL); PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); @@ -260,4 +258,5 @@ tmpVector->data.F32[col] = image->data.F32[row][col]; } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); return(tmpVector); } @@ -276,4 +275,5 @@ pmPeakType type) { + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); @@ -283,5 +283,8 @@ } psArrayAdd(list, 100, tmpPeak); - + psFree (tmpPeak); + // XXX EAM : is this free appropriate? (does psArrayAdd increment memory counter?) + + psTrace(__func__, 4, "---- %s() end ----\n", __func__); return(list); } @@ -308,8 +311,10 @@ psF32 threshold) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); PS_ASSERT_IMAGE_NON_NULL(image, NULL); PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); if ((image->numRows == 1) || (image->numCols == 1)) { psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns."); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); return(NULL); } @@ -364,4 +369,6 @@ } } + psFree (tmpRow); + psFree (row1); // @@ -369,4 +376,5 @@ // if (image->numRows == 1) { + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return(list); } @@ -442,8 +450,10 @@ } } else { - psError(PS_ERR_UNKNOWN, true, "peak specified valid column range."); + psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); } } + psFree (tmpRow); + psFree (row1); } @@ -484,7 +494,10 @@ } } else { - psError(PS_ERR_UNKNOWN, true, "peak specified valid colum range."); - } - } + psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); + } + } + psFree (tmpRow); + psFree (row1); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return(list); } @@ -495,12 +508,13 @@ psS32 y) { - + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); if ((x >= valid.x0) && (x <= valid.x1) && (y >= valid.y0) && (y <= valid.y1)) { + psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); return(true); } - + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); return(false); } @@ -516,4 +530,8 @@ XXX: changed API to create a NEW output psArray (should change name as well) + +XXX: Do we free the psList elements of those culled peaks? + +XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset? *****************************************************************************/ psList *pmCullPeaks(psList *peaks, @@ -521,6 +539,6 @@ const psRegion valid) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); PS_ASSERT_PTR_NON_NULL(peaks, NULL); - // PS_ASSERT_PTR_NON_NULL(valid, NULL); psListElem *tmpListElem = (psListElem *) peaks->head; @@ -539,4 +557,5 @@ } + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return(peaks); } @@ -544,6 +563,11 @@ // 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) -{ +// XXX: Fix the *valid pointer. +psArray *pmPeaksSubset( + psArray *peaks, + psF32 maxValue, + const psRegion valid) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); PS_ASSERT_PTR_NON_NULL(peaks, NULL); @@ -561,4 +585,5 @@ psArrayAdd (output, 200, tmpPeak); } + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return(output); } @@ -601,158 +626,51 @@ 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. - // + + + + + + + + + +bool pmSourceLocalSky( + pmSource *source, + psStatsOptions statsOptions, + psF32 Radius) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_IMAGE_NON_NULL(source->pixels, false); + PS_ASSERT_IMAGE_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_INT_POSITIVE(Radius, false); + PS_ASSERT_INT_NONNEGATIVE(Radius, false); + + psImage *image = source->pixels; + psImage *mask = source->mask; + pmPeak *peak = source->peak; + psRegion srcRegion; + + srcRegion = psRegionForSquare(peak->x, peak->y, Radius); + srcRegion = psRegionForImage(mask, srcRegion); + + psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED); psStats *myStats = psStatsAlloc(statsOptions); - myStats = psImageStats(myStats, subImage, subImageMask, 1); - - // - // Create the output mySource, and set appropriate members. - // - pmSource *mySource = pmSourceAlloc(); - mySource->peak = (pmPeak *) peak; - mySource->moments = pmMomentsAlloc(); + myStats = psImageStats(myStats, image, mask, 0xff); + psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED); + 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); + if (isnan(tmpF64)) { + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return(false); + } + source->moments = pmMomentsAlloc(); + source->moments->Sky = (psF32) tmpF64; + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return (true); } @@ -770,4 +688,5 @@ psF32 y) { + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); /// XXX EAM should compare with hypot (x,y) for speed if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { @@ -775,4 +694,5 @@ } + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); return(false); } @@ -787,16 +707,23 @@ pmSource->peak pmSource->pixels + pmSource->weight + pmSource->mask 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 EAM : this version clips input pixels on S/N +XXX EAM : this version returns false for several reasons +*****************************************************************************/ +# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) + +bool pmSourceMoments(pmSource *source, + psF32 radius) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); // @@ -816,7 +743,7 @@ // XY = SUM (x - xc)*(y - yc)*(z - sky) // - psF32 Sum = 0.0; psF32 peakPixel = -PS_MAX_F32; psS32 numPixels = 0; + psF32 Sum = 0.0; psF32 X1 = 0.0; psF32 Y1 = 0.0; @@ -824,7 +751,11 @@ psF32 Y2 = 0.0; psF32 XY = 0.0; - psF32 x = 0; - psF32 y = 0; - // + psF32 x = 0; + psF32 y = 0; + psF32 R2 = PS_SQR(radius); + + psF32 xPeak = source->peak->x; + psF32 yPeak = source->peak->y; + // XXX why do I get different results for these two methods of finding Sx? // XXX Sx, Sy would be better measured if we clip pixels close to sky @@ -834,34 +765,48 @@ // proceed with the moments calculation. need to do two loops for a // numerically stable result. first loop: get the sums. - // + // XXX EAM : mask == 0 is valid + for (psS32 row = 0; row < source->pixels->numRows ; row++) { for (psS32 col = 0; col < source->pixels->numCols ; col++) { - if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 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++; - } - } - } - } + if ((source->mask != NULL) && (source->mask->data.U8[row][col])) + continue; + + psF32 xDiff = col + source->pixels->col0 - xPeak; + psF32 yDiff = row + source->pixels->row0 - yPeak; + + // XXX EAM : calculate xDiff, yDiff up front; + // radius is just a function of (xDiff, yDiff) + if (!VALID_RADIUS(xDiff, yDiff, R2)) + continue; + + psF32 pDiff = source->pixels->data.F32[row][col] - sky; + + // XXX EAM : check for valid S/N in pixel + // XXX EAM : should this limit be user-defined? + if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) + continue; + + Sum += pDiff; + X1 += xDiff * pDiff; + Y1 += yDiff * pDiff; + XY += xDiff * yDiff * pDiff; + + X2 += PS_SQR(xDiff) * pDiff; + Y2 += PS_SQR(yDiff) * pDiff; + + peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel); + numPixels++; + } + } + // XXX EAM - the limit is a bit arbitrary. make it user defined? + if ((numPixels < 3) || (Sum <= 0)) { + psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return (false); + } + + psTrace (".psModules.pmSourceMoments", 5, + "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", + sky, Sum, X1, Y1, X2, Y2, XY, numPixels); // @@ -872,8 +817,17 @@ 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; + if ((fabs(x) > radius) || (fabs(y) > radius)) { + psTrace (".psModules.pmSourceMoments", 5, + "large centroid swing; invalid peak %d, %d\n", + source->peak->x, source->peak->y); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return (false); + } + + source->moments->x = x + xPeak; + source->moments->y = y + yPeak; + + // XXX EAM : Sxy needs to have x*y subtracted + source->moments->Sxy = XY/Sum - x*y; source->moments->Sum = Sum; source->moments->Peak = peakPixel; @@ -883,40 +837,12 @@ 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 + + psTrace (".psModules.pmSourceMoments", 4, + "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", + sky, Sum, source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, source->moments->Sxy); + + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); } @@ -924,4 +850,5 @@ int pmComparePeakAscend (const void **a, const void **b) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); pmPeak *A = *(pmPeak **)a; pmPeak *B = *(pmPeak **)b; @@ -930,8 +857,12 @@ diff = A->counts - B->counts; - if (diff < FLT_EPSILON) + if (diff < FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); return (-1); - if (diff > FLT_EPSILON) + } else if (diff > FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); return (+1); + } + psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); return (0); } @@ -939,4 +870,5 @@ int pmComparePeakDescend (const void **a, const void **b) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); pmPeak *A = *(pmPeak **)a; pmPeak *B = *(pmPeak **)b; @@ -945,37 +877,34 @@ diff = A->counts - B->counts; - if (diff < FLT_EPSILON) + if (diff < FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); return (+1); - if (diff > FLT_EPSILON) + } else if (diff > FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); return (-1); + } + psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); return (0); } /****************************************************************************** -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; +pmSourcePSFClump(source, metadata): Find the likely PSF clump in the +sigma-x, sigma-y plane. return 0,0 clump in case of error. +*****************************************************************************/ + +// XXX EAM include a S/N cutoff in selecting the sources? +pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + # define NPIX 10 + # define SCALE 0.1 + psArray *peaks = NULL; - psF32 clumpX = 0.0; - psF32 clumpDX = 0.0; - psF32 clumpY = 0.0; - psF32 clumpDY = 0.0; + pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; + pmPSFClump psfClump = emptyClump; + + PS_ASSERT_PTR_NON_NULL(sources, emptyClump); + PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); // find the sigmaX, sigmaY clump @@ -986,5 +915,10 @@ // construct a sigma-plane image + // psImageAlloc does zero the data splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32); + for (int i = 0; i < splane->numRows; i++) + { + memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); + } // place the sources in the sigma-plane image (ignore 0,0 values?) @@ -992,6 +926,10 @@ { 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? + if (tmpSrc == NULL) { + continue; + } + if (tmpSrc->moments == NULL) { + continue; + } // Sx,Sy are limited at 0. a peak at 0,0 is artificial @@ -1022,5 +960,11 @@ psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); - } + psFree (splane); + psFree (stats); + + } + // XXX EAM : possible errors: + // 1) no peak in splane + // 2) no significant peak in splane // measure statistics on Sx, Sy if Sx, Sy within range of clump @@ -1036,5 +980,5 @@ psArraySort (peaks, pmComparePeakDescend); clump = peaks->data[0]; - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d\n", clump->x, clump->y); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); // define section window for clump @@ -1077,24 +1021,58 @@ stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); - clumpX = stats->clippedMean; - clumpDX = stats->clippedStdev; + psfClump.X = stats->clippedMean; + psfClump.dX = 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); + psfClump.Y = stats->clippedMean; + psfClump.dY = stats->clippedStdev; + + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); // these values should be pushed on the metadata somewhere - } - - int Nsat = 0; - int Ngal = 0; - int Nfaint = 0; - int Nstar = 0; - int Npsf = 0; - int Ncr = 0; + + psFree (stats); + psFree (peaks); + psFree (tmpSx); + psFree (tmpSy); + } + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return (psfClump); +} + +/****************************************************************************** +pmSourceRoughClass(source, metadata): make a guess at the source +classification. + +XXX: push the clump info into the metadata? + +XXX: How can this function ever return FALSE? + +XXX EAM : add the saturated mask value to metadata +*****************************************************************************/ + +bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + psBool rc = true; + + int Nsat = 0; + int Ngal = 0; + int Nstar = 0; + int Npsf = 0; + int Ncr = 0; + int Nsatstar = 0; + psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); starsn->n = 0; + + // check return status value (do these exist?) + bool status; + psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE"); + psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN"); + psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); + // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE"); // XXX allow clump size to be scaled relative to sigmas? @@ -1109,30 +1087,48 @@ 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; + // calculate and save signal-to-noise estimates + psF32 S = tmpSrc->moments->Sum; + psF32 A = 4 * M_PI * sigX * sigY; + psF32 B = tmpSrc->moments->Sky; + psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN))); + psF32 SN = (S * sqrt(GAIN) / RT); + tmpSrc->moments->SN = SN; + + // XXX EAM : can we use the value of SATURATE if mask is NULL? + // + // XXX: Must verify this region (the region argument was added to psImageCountPixelMask() + // after EAM wrote this code. + // + psRegion tmpRegion = psRegionSet(0, tmpSrc->mask->numCols, 0, tmpSrc->mask->numRows); + int Nsatpix = psImageCountPixelMask(tmpSrc->mask, tmpRegion, PSPHOT_MASK_SATURATED); + + // saturated star (size consistent with PSF or larger) + // Nsigma should be user-configured parameter + bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); + if ((Nsatpix > 1) && big) { + tmpSrc->type |= PM_SOURCE_SATSTAR; + Nsatstar ++; + continue; + } + + // saturated object (not a star, eg bleed trails, hot pixels) + if (Nsatpix > 1) { + tmpSrc->type |= PM_SOURCE_SATURATED; Nsat ++; continue; } - // too small to be stellar - if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) { - tmpSrc->type |= PS_SOURCE_DEFECT; + // likely defect (too small to be stellar) (push out to 3 sigma) + // low S/N objects which are small are probably stellar + // only set candidate defects if + if ((sigX < 0.05) || (sigY < 0.05)) { + tmpSrc->type |= PM_SOURCE_DEFECT; Ncr ++; continue; } - // possible galaxy - if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) { - tmpSrc->type |= PS_SOURCE_GALAXY; + // likely unsaturated galaxy (too large to be stellar) + if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { + tmpSrc->type |= PM_SOURCE_GALAXY; Ngal ++; continue; @@ -1140,24 +1136,12 @@ // 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; + // PSF star (within 1.5 sigma of clump center + psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); + if ((SN > PSF_SN_LIM) && (radius < 1.5)) { + tmpSrc->type |= PM_SOURCE_PSFSTAR; Npsf ++; continue; @@ -1165,5 +1149,5 @@ // random type of star - tmpSrc->type |= PS_SOURCE_OTHER; + tmpSrc->type |= PM_SOURCE_OTHER; } @@ -1173,26 +1157,53 @@ 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); - + psFree (stats); + psFree (starsn); + } + + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal: %3d\n", Ngal); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); + + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); return(rc); } +/** pmSourceDefinePixels() + * + * Define psImage subarrays for the source located at coordinates x,y on the + * image set defined by readout. The pixels defined by this operation consist of + * a square window (of full width 2Radius+1) centered on the pixel which contains + * the given coordinate, in the frame of the readout. The window is defined to + * have limits which are valid within the boundary of the readout image, thus if + * the radius would fall outside the image pixels, the subimage is truncated to + * only consist of valid pixels. If readout->mask or readout->weight are not + * NULL, matching subimages are defined for those images as well. This function + * fails if no valid pixels can be defined (x or y less than Radius, for + * example). This function should be used to define a region of interest around a + * source, including both source and sky pixels. + * + * XXX: must code this. + * + */ +bool pmSourceDefinePixels( + pmSource *mySource, ///< Add comment. + pmReadout *readout, ///< Add comment. + psF32 x, ///< Add comment. + psF32 y, ///< Add comment. + psF32 Radius) ///< Add comment. +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented. Returning FALSE.\n"); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(false); +} + /****************************************************************************** pmSourceSetPixelsCircle(source, image, radius) -XXX: 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. +XXX: This was replaced by DefinePixels in SDRS. Remove it. *****************************************************************************/ bool pmSourceSetPixelsCircle(pmSource *source, @@ -1200,9 +1211,10 @@ psF32 radius) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); PS_ASSERT_IMAGE_NON_NULL(image, false); PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); PS_ASSERT_PTR_NON_NULL(source, false); PS_ASSERT_PTR_NON_NULL(source->moments, false); - // PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); PS_FLOAT_COMPARE(0.0, radius, false); @@ -1217,33 +1229,9 @@ psS32 SubImageCenterCol = source->peak->x; // XXX EAM : for the circle to stay on the image + // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); - psS32 SubImageEndRow = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32); + psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1); 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 + psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); // XXX: Must recycle image. @@ -1255,10 +1243,8 @@ psFree(source->pixels); } - // XXX: Must fix this. psImageSubset() has different parameters in latest CVS. - // source->pixels = psImageSubset((psImage *) image, - // SubImageStartCol, - // SubImageStartRow, - // SubImageEndCol, - // SubImageEndRow); + source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, + SubImageStartRow, + SubImageEndCol, + SubImageEndRow)); // XXX: Must recycle image. @@ -1268,5 +1254,5 @@ source->mask = psImageAlloc(source->pixels->numCols, source->pixels->numRows, - PS_TYPE_F32); + PS_TYPE_U8); // XXX EAM : type was F32 // @@ -1287,135 +1273,31 @@ } } + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 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. +pmSourceModelGuess(source, model): This function allocates a new +pmModel structure based on the given modelType specified in the argument list. +The corresponding pmModelGuess function is returned, and used to +supply the values of the params array in the pmModel structure. XXX: Many parameters are based on the src->moments structure, which is in image, not subImage coords. Therefore, the calls to the model evaluation functions will be in image, not subImage coords. Remember this. - -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); +*****************************************************************************/ +pmModel *pmSourceModelGuess(pmSource *source, + pmModelType modelType) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); PS_ASSERT_PTR_NON_NULL(source->moments, false); PS_ASSERT_PTR_NON_NULL(source->peak, false); - 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); - } + + pmModel *model = pmModelAlloc(modelType); + + pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); + modelGuessFunc(model, source); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(model); } @@ -1441,50 +1323,27 @@ 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; - - // + +// XXX EAM : I have made this a public function +// XXX EAM : this now uses a pmModel as the input +// XXX EAM : it was using src->type to find the model, not model->type +psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(image, false); + PS_ASSERT_PTR_NON_NULL(model, false); + PS_ASSERT_PTR_NON_NULL(model->params, false); + // Allocate the x coordinate structure and convert row/col to image space. // psVector *x = psVectorAlloc(2, PS_TYPE_F32); - x->data.F32[0] = (psF32) (col + src->pixels->col0); - x->data.F32[1] = (psF32) (row + src->pixels->row0); + x->data.F32[0] = (psF32) (col + image->col0); + x->data.F32[1] = (psF32) (row + image->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); - } + pmModelFunc modelFunc; + + modelFunc = pmModelFunc_GetFunction (model->type); + tmpF = modelFunc (NULL, model->params, x); psFree(x); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return(tmpF); } @@ -1508,4 +1367,5 @@ psU32 dir) { + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); // // Convert coords to subImage space. @@ -1517,13 +1377,18 @@ if (!((0 <= subCol) && (subCol < source->pixels->numCols))) { psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); return(NAN); } if (!((0 <= subRow) && (subRow < source->pixels->numRows))) { + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range"); return(NAN); } - psF32 oldValue = evalModel(source, subRow, subCol); + // XXX EAM : i changed this to match pmModelEval above, but see + // XXX EAM the note below in pmSourceContour + psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); if (oldValue == level) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); return(((psF32) (subCol + source->pixels->col0))); } @@ -1545,6 +1410,7 @@ while (subCol != lastColumn) { - psF32 newValue = evalModel(source, subRow, subCol); + psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); if (oldValue == level) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); return((psF32) (subCol + source->pixels->col0)); } @@ -1552,4 +1418,5 @@ if ((newValue <= level) && (level <= oldValue)) { // This is simple linear interpolation. + psTrace(__func__, 4, "---- %s() end ----\n", __func__); return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); } @@ -1557,4 +1424,5 @@ if ((oldValue <= level) && (level <= newValue)) { // This is simple linear interpolation. + psTrace(__func__, 4, "---- %s() end ----\n", __func__); return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); } @@ -1563,4 +1431,5 @@ } + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); return(NAN); } @@ -1576,4 +1445,6 @@ XXX: What is mode? XXX: The top, bottom of the contour is not correctly determined. +XXX EAM : this functions is using the model for the contour, but it should + be using only the image counts *****************************************************************************/ psArray *pmSourceContour(pmSource *source, @@ -1582,4 +1453,5 @@ pmContourType mode) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); PS_ASSERT_PTR_NON_NULL(source, false); PS_ASSERT_PTR_NON_NULL(image, false); @@ -1587,5 +1459,6 @@ 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_PTR_NON_NULL(source->modelFLT, false); + // XXX EAM : what is the purpose of modelPSF/modelFLT? // @@ -1610,4 +1483,5 @@ psFree(xVec); psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); return(NULL); //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); @@ -1622,8 +1496,9 @@ psFree(xVec); psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); return(NULL); //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); } - //printf("The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); + psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; @@ -1646,4 +1521,5 @@ psFree(xVec); psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); return(NULL); //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); @@ -1657,4 +1533,5 @@ psFree(xVec); psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); return(NULL); //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); @@ -1672,37 +1549,43 @@ tmpArray->data[0] = (psPtr *) yVec; tmpArray->data[1] = (psPtr *) xVec; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); 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_NUM_ITERATIONS 15 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 /****************************************************************************** -pmSourceFitModel(source, image): must create the appropiate arguments to the +pmSourceFitModel(source, model): must create the appropiate arguments to the LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. XXX: should there be a mask value? -XXX: Probably should remove the "image" argument. -*****************************************************************************/ -bool pmSourceFitModel(pmSource *source, - const psImage *image) -{ +XXX EAM : fit the specified model (not necessarily the one in source) +*****************************************************************************/ +bool pmSourceFitModel_v5(pmSource *source, + pmModel *model, + const bool PSF) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); PS_ASSERT_PTR_NON_NULL(source, false); PS_ASSERT_PTR_NON_NULL(source->moments, false); PS_ASSERT_PTR_NON_NULL(source->peak, false); PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->modelPSF, false); - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); - psBool rc; + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->weight, false); + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + // XXX EAM : is it necessary for the mask & weight to exist? the + // tests below could be conditions (!NULL) + + psVector *params = model->params; + psVector *dparams = model->dparams; + psVector *paramMask = NULL; + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + + int nParams = PSF ? params->n - 4 : params->n; // find the number of valid pixels @@ -1716,4 +1599,9 @@ } } + } + if (count < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return(false); } @@ -1733,7 +1621,7 @@ 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])); + yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]); + // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then + // the minimization function calculates sq() tmpCnt++; } @@ -1744,59 +1632,241 @@ 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; - } + // PSF model only fits first 4 parameters, FLT model fits all + if (PSF) { + paramMask = psVectorAlloc (params->n, PS_TYPE_U8); + for (int i = 0; i < 4; i++) { + paramMask->data.U8[i] = 0; + } + for (int i = 4; i < paramMask->n; i++) { + paramMask->data.U8[i] = 1; + } + } + + // XXX EAM : covar must be F64? + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); + + psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if ((paramMask != NULL) && paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + // XXX EAM: we need to do something (give an error?) if rc is false + // XXX EAM: psMinimizeLMChi2 does not check convergence + + // XXX models can go insane: reject these + onPic &= (params->data.F32[2] >= source->pixels->col0); + onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (params->data.F32[3] >= source->pixels->row0); + onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + // XXX EAM: save the resulting chisq, nDOF, nIter - source->modelPSF->chisq = myMin->value; - source->modelPSF->nDOF = y->n - params->n; - source->modelPSF->nIter = myMin->iter; - + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // XXX EAM get the Gauss-Newton distance for fixed model parameters + if (paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if (!paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + psFree (delta); + } + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(paramMask); + + rc = (onPic && fitStatus); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +// XXX EAM : new version with parameter range limits and weight enhancement +bool pmSourceFitModel (pmSource *source, + pmModel *model, + const bool PSF) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->weight, false); + + // XXX EAM : is it necessary for the mask & weight to exist? the + // tests below could be conditions (!NULL) + + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + psF32 Ro, ymodel; + + psVector *params = model->params; + psVector *dparams = model->dparams; + psVector *paramMask = NULL; + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + + // XXX EAM : I need to use the sky value to constrain the weight model + int nParams = PSF ? params->n - 4 : params->n; + psF32 So = params->data.F32[0]; + + // find the number of valid pixels + // XXX EAM : this loop and the loop below could just be one pass + // using the psArrayAdd and psVectorExtend functions + psS32 count = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j] == 0) { + count++; + } + } + } + if (count < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return(false); + } + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(count); + psVector *y = psVectorAlloc(count, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); + psS32 tmpCnt = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j] == 0) { + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + // XXX: Convert i/j to image space: + // XXX EAM: coord order is (x,y) == (col,row) + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[tmpCnt] = (psPtr *) coord; + y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; + + // compare observed flux to model flux to adjust weight + ymodel = modelFunc (NULL, model->params, coord); + + // this test enhances the weight based on deviation from the model flux + Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So)); + + // psMinimizeLMChi2_EAM takes wt = 1/dY^2 + if (source->weight->data.F32[i][j] == 0) { + yErr->data.F32[tmpCnt] = 0.0; + } else { + yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro); + } + tmpCnt++; + } + } + } + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + // PSF model only fits first 4 parameters, FLT model fits all + if (PSF) { + paramMask = psVectorAlloc (params->n, PS_TYPE_U8); + for (int i = 0; i < 4; i++) { + paramMask->data.U8[i] = 0; + } + for (int i = 4; i < paramMask->n; i++) { + paramMask->data.U8[i] = 1; + } + } + + // XXX EAM : I've added three types of parameter range checks + // XXX EAM : this requires my new psMinimization functions + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); + psVector *beta_lim = NULL; + psVector *params_min = NULL; + psVector *params_max = NULL; + + // XXX EAM : in this implementation, I pass in the limits with the covar matrix. + // in the SDRS, I define a new psMinimization which will take these in + psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); + modelLimits (&beta_lim, ¶ms_min, ¶ms_max); + for (int i = 0; i < params->n; i++) { + covar->data.F64[0][i] = beta_lim->data.F32[i]; + covar->data.F64[1][i] = params_min->data.F32[i]; + covar->data.F64[2][i] = params_max->data.F32[i]; + } + + psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if ((paramMask != NULL) && paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // XXX EAM: we need to do something (give an error?) if rc is false + // XXX EAM: psMinimizeLMChi2 does not check convergence + + // XXX models can go insane: reject these + onPic &= (params->data.F32[2] >= source->pixels->col0); + onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (params->data.F32[3] >= source->pixels->row0); + onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + + // XXX EAM: save the resulting chisq, nDOF, nIter + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // XXX EAM get the Gauss-Newton distance for fixed model parameters + if (paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if (!paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + } + + psFree(paramMask); psFree(x); psFree(y); psFree(myMin); + + rc = (onPic && fitStatus); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); return(rc); } -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); +bool p_pmSourceAddOrSubModel(psImage *image, + psImage *mask, + pmModel *model, + bool center, + psS32 flag) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + PS_ASSERT_PTR_NON_NULL(model, false); PS_ASSERT_IMAGE_NON_NULL(image, false); PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); psVector *x = psVectorAlloc(2, PS_TYPE_F32); - psVector *params = src->modelPSF->params; - - for (psS32 i = 0; i < src->pixels->numRows; i++) { - for (psS32 j = 0; j < src->pixels->numCols; j++) { + psVector *params = model->params; + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + psS32 imageCol; + psS32 imageRow; + + for (psS32 i = 0; i < image->numRows; i++) { + for (psS32 j = 0; j < image->numCols; j++) { + if ((mask != NULL) && mask->data.U8[i][j]) + continue; psF32 pixelValue; // XXX: Should you be adding the pixels for the entire subImage, @@ -1805,33 +1875,19 @@ // 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; + // XXX EAM : 'center' option changes this + // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2] + if (center) { + imageCol = j - 0.5*image->numCols + model->params->data.F32[2]; + imageRow = i - 0.5*image->numRows + model->params->data.F32[3]; + } else { + imageCol = j + image->col0; + imageRow = i + image->row0; + } x->data.F32[0] = (float) imageCol; x->data.F32[1] = (float) imageRow; - 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); - } + pixelValue = modelFunc (NULL, params, x); + // fprintf (stderr, "%f %f %d %d %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue); + if (flag == 1) { pixelValue = -pixelValue; @@ -1841,8 +1897,9 @@ // how to use the boolean "center" flag. - image->data.F32[imageRow][imageCol]+= pixelValue; + image->data.F32[i][j]+= pixelValue; } } psFree(x); + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); return(true); } @@ -1853,8 +1910,12 @@ *****************************************************************************/ bool pmSourceAddModel(psImage *image, - pmSource *src, + psImage *mask, + pmModel *model, bool center) { - return(sourceAddOrSubModel(image, src, center, 0)); + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); } @@ -1862,259 +1923,13 @@ *****************************************************************************/ bool pmSourceSubModel(psImage *image, - pmSource *src, + psImage *mask, + pmModel *model, 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); -} + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +