Index: /trunk/psModules/src/pmObjects.c =================================================================== --- /trunk/psModules/src/pmObjects.c (revision 3497) +++ /trunk/psModules/src/pmObjects.c (revision 3498) @@ -5,9 +5,9 @@ * @author GLG, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-02-15 23:59:05 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-24 22:36:52 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii -xd * + * */ @@ -17,4 +17,81 @@ #include "psConstants.h" #include "pmObjects.h" +/****************************************************************************** + *****************************************************************************/ + +psPeak *pmPeakAlloc(psS32 x, + psS32 y, + psF32 counts, + psPeakType class) +{ + psPeak *tmp = (psPeak *) psAlloc(sizeof(psPeak)); + tmp->x = x; + tmp->y = y; + tmp->counts = counts; + tmp->class = class; + + return(tmp); +} + +psMoments *pmMomentsAlloc() +{ + psMoments *tmp = (psMoments *) psAlloc(sizeof(psMoments)); + 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 p_psModelFree(psModel *tmp) +{ + psFree(tmp->params); + psFree(tmp->dparams); +} + +psModel *pmModelAlloc(psModelType type) +{ + psModel *tmp = (psModel *) psAlloc(sizeof(psModel)); + + tmp->type = type; + tmp->Nparams = 0; + tmp->params = NULL; + tmp->dparams = NULL; + tmp->chisq = 0.0; + + p_psMemSetDeallocator(tmp, (psFreeFcn) p_psModelFree); + 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 p_psSourceFree(psSource *tmp) +{ + psFree(tmp->peak); + // psFree(tmp->pixels); + // psFree(tmp->mask); + psFree(tmp->moments); + psFree(tmp->models); +} + +psSource *pmSourceAlloc() +{ + psSource *tmp = (psSource *) psAlloc(sizeof(psSource)); + tmp->peak = NULL; + tmp->pixels = NULL; + tmp->mask = NULL; + tmp->moments = NULL; + tmp->models = NULL; + p_psMemSetDeallocator(tmp, (psFreeFcn) p_psSourceFree); + + return(tmp); +} /****************************************************************************** @@ -23,5 +100,5 @@ the location (x value) of all peaks. -XXX: What types should be supported? Only F32 us implemented. +XXX: What types should be supported? Only F32 is implemented. XXX: We currently step through the input vector twice; once to determine the @@ -33,8 +110,27 @@ { PS_VECTOR_CHECK_NULL(vector, NULL); + PS_VECTOR_CHECK_EMPTY(vector, NULL); PS_VECTOR_CHECK_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)) { @@ -42,4 +138,7 @@ } + // + // Determine if interior pixels are peaks + // for (psU32 i = 1; i < n-1 ; i++) { if ((vector->data.F32[i] > vector->data.F32[i-1]) && @@ -49,4 +148,8 @@ } } + + // + // Determine if last pixel is a peak + // if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && (vector->data.F32[n-1] > threshold)) { @@ -54,10 +157,22 @@ } + // + // 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]) && @@ -67,4 +182,8 @@ } } + + // + // Determine if last pixel is a peak + // if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && (vector->data.F32[n-1] > threshold)) { @@ -91,8 +210,36 @@ } +// XXX: Switch row, col args? +psList *MyListAddPeak(psList *list, + psS32 row, + psS32 col, + psF32 counts, + psPeakType type) +{ + psPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); + + if (list == NULL) { + list = psListAlloc(tmpPeak); + } else { + psListAdd(list, PS_LIST_HEAD, tmpPeak); + } + + return(list); +} + /****************************************************************************** pmFindImagePeeks(image, threshold): Find all local peaks in the given psImage above the given threshold. Returns a psList 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 psList elements, should we use the image row/column offsets? + Currently, we do not. *****************************************************************************/ psList *pmFindImagePeeks(const psImage *image, @@ -101,8 +248,15 @@ PS_IMAGE_CHECK_NULL(image, NULL); PS_IMAGE_CHECK_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."); + } psVector *tmpRow = NULL; psU32 col = 0; psU32 row = 0; - + psList *list = NULL; + + // + // Find peaks in row 0 only. + // row = 0; tmpRow = p_psGetRowVectorFromImage((psImage *) image, row); @@ -110,4 +264,8 @@ 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]) && @@ -115,5 +273,5 @@ (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { if (image->data.F32[row][col] > threshold) { - // Add peak at location (row, col) + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } } @@ -125,5 +283,5 @@ (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { if (image->data.F32[row][col] > threshold) { - // Add peak at location (row, col) + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } } @@ -134,13 +292,22 @@ (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) { if (image->data.F32[row][col] > threshold) { - // Add peak at location (row, col) + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } } } else { - printf("XXX: ERROR!\n"); - } - } - + psError(PS_ERR_UNKNOWN, true, "peak specified valid colum 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 = p_psGetRowVectorFromImage((psImage *) image, 0); @@ -160,5 +327,28 @@ (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { if (image->data.F32[row][col] > threshold) { - // Add peak at location (row, col) + + psPeakType myType = PM_PEAK_UNDEF; + 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); } } @@ -166,4 +356,7 @@ } + // + // Find peaks in the last row only. + // row = image->numRows - 1; tmpRow = p_psGetRowVectorFromImage((psImage *) image, row); @@ -176,5 +369,5 @@ (image->data.F32[row][col] > image->data.F32[row][col+1])) { if (image->data.F32[row][col] > threshold) { - // Add peak at location (row, col) + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } } @@ -186,5 +379,5 @@ (image->data.F32[row][col] >= image->data.F32[row][col+1])) { if (image->data.F32[row][col] > threshold) { - // Add peak at location (row, col) + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } } @@ -195,14 +388,31 @@ (image->data.F32[row][col] > image->data.F32[row][col-1])) { if (image->data.F32[row][col] > threshold) { - // Add peak at location (row, col) + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } } } else { - printf("XXX: ERROR!\n"); - } - } - - return(NULL); -} + psError(PS_ERR_UNKNOWN, true, "peak specified valid colum range."); + } + } + + return(list); +} + +// XXX: Macro this. +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); +} + /****************************************************************************** @@ -210,5 +420,7 @@ a peak value above the given maximum, or fall outside the valid region. -XXX: Do we free the psList elements of those culled peaks? +XXX: Should the sky value be used when comparing the maximum? + +XXX: warning message if valid is NULL? *****************************************************************************/ psList *pmCullPeeks(psList *peaks, @@ -216,41 +428,372 @@ const psRegion *valid) { - return(NULL); -} - -/****************************************************************************** -psSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): + PS_PTR_CHECK_NULL(peaks, NULL); + // PS_PTR_CHECK_NULL(valid, NULL); + + psListElem *tmpListElem = (psListElem *) peaks->head; + psS32 indexNum = 0; + + // printf("pmCullPeeks(): list size is %d\n", peaks->size); + while (tmpListElem != NULL) { + psPeak *tmpPeak = (psPeak *) 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); +} + +/****************************************************************************** +psSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this +routine creates a new psSource data structure and sets the following members: + ->psPeak + ->psMoments->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 psSource->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. *****************************************************************************/ psSource *pmSourceLocalSky(const psImage *image, const psPeak *peak, + psStatsOptions statsOptions, psF32 innerRadius, psF32 outerRadius) { - return(NULL); -} - -/****************************************************************************** + PS_IMAGE_CHECK_NULL(image, NULL); + PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL); + PS_PTR_CHECK_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. + // + psS32 SubImageCenterRow = peak->y; + psS32 SubImageCenterCol = peak->x; + psS32 SubImageStartRow = SubImageCenterRow - outerRadiusS32; + psS32 SubImageEndRow = SubImageCenterRow + outerRadiusS32; + psS32 SubImageStartCol = SubImageCenterCol - outerRadiusS32; + psS32 SubImageEndCol = 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); + + 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); + } + + // + // Grab a subimage of the original image of size (2 * outerRadius). + // + psImage *subImage = psImageSubset((psImage *) image, + SubImageStartCol, + SubImageStartRow, + SubImageEndCol, + SubImageEndRow); + // printf("pmSourceLocalSky: subimage width/length is (%d, %d)\n", subImage->numCols, subImage->numRows); + 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. + // + 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. + // + psSource *mySource = pmSourceAlloc(); + mySource->peak = (psPeak *) 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. + +XXX: should arguments be (y, x) order? + *****************************************************************************/ +bool CheckRadius(psPeak *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(psF32 xCenter, + psF32 yCenter, + psF32 radius, + psF32 x, + psF32 y) +{ + if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { + return(true); + } + + return(false); +} + +/****************************************************************************** +pmSourceMoments(source, radius) + +Requires the following to have been created: + psSource + psSource->peak + psSource->pixels + +XXX: mask values? *****************************************************************************/ psSource *pmSourceMoments(psSource *source, - const psImage *image, psF32 radius) { - return(NULL); -} - -/****************************************************************************** -pmSourceRoughClass(source, saturate, SNlim, valid): make a guessat the source + PS_PTR_CHECK_NULL(source, NULL); + PS_PTR_CHECK_NULL(source->peak, NULL); + PS_PTR_CHECK_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; + // + // 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. + // + 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; + X2+= PS_SQR(xDiff) * pDiff; + Y2+= PS_SQR(yDiff) * pDiff; + XY+= xDiff * 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 + // + source->moments->x = X1/Sum + ((psF32) source->peak->x); + source->moments->y = Y1/Sum + ((psF32) source->peak->y); + source->moments->Sx = sqrt(X2/Sum - PS_SQR(X1/Sum)); + source->moments->Sy = sqrt(Y2/Sum - PS_SQR(Y1/Sum)); + source->moments->Sxy = XY/Sum; + source->moments->Peak = peakPixel; + source->moments->nPixels = numPixels; + + return(source); +} + +/****************************************************************************** +pmSourceRoughClass(source, saturate, SNlim, valid): make a guess at the source classification. - *****************************************************************************/ + psSource *pmSourceRoughClass(psSource *source, psF32 saturate, float SNlim, const psRegion *valid) -{ - return(NULL); -} +XXX: Waiting for sample code from IfA. + +XXX: Most code this. + *****************************************************************************/ +#define SATURATE 0.0 +#define FAINT_SN_LIM 0.0 +#define PSF_SN_LIM 0.0 +#define SATURATE 0.0 +#define SATURATE 0.0 + +bool pmSourceRoughClass(psArray *source, + psMetadata *metadata) +{ + PS_PTR_CHECK_NULL(source, false); + PS_PTR_CHECK_NULL(metadata, false); + psBool rc = true; + + for (psS32 i = 0 ; i < source->n ; i++) { + psSource *tmpSrc = (psSource *) source->data[i]; + PS_PTR_CHECK_NULL(tmpSrc->moments, false); + + if (tmpSrc->moments->Peak > SATURATE) { + tmpSrc->peak->class = PS_SOURCE_SATURATED; + } else { + // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE. + psF32 gain = 0.0; + psF32 readNoise = 0.0; + psF32 S = tmpSrc->moments->Sum; + psF32 A = PS_PI * tmpSrc->moments->Sx * tmpSrc->moments->Sy; + psF32 B = tmpSrc->moments->Sky; + psF32 SN = (PS_SQRT_F32(gain) * S) / + PS_SQRT_F32(S + (A * B) + ((A * readNoise * readNoise) / PS_SQRT_F32(gain))); + if (SN < FAINT_SN_LIM) { + tmpSrc->peak->class = PS_SOURCE_FAINTSTAR; + } + if (SN < PSF_SN_LIM) { + tmpSrc->peak->class = PS_SOURCE_FAINTSTAR; + } + } + } + + return(rc); +} + + /****************************************************************************** pmSourceSetPixelCircle(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 pmSourceSetPixelCircle(psSource *source, @@ -258,4 +801,89 @@ psF32 radius) { + PS_IMAGE_CHECK_NULL(image, false); + PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); + PS_PTR_CHECK_NULL(source, false); + // PS_PTR_CHECK_NULL(source->moments, false); + PS_PTR_CHECK_NULL(source->peak, false); + PS_FLOAT_COMPARE(0.0, radius, false); + + // + // We define variables for code readability. + // + psS32 radiusS32 = (psS32) radius; + psS32 SubImageCenterRow = source->peak->y; + psS32 SubImageCenterCol = source->peak->x; + psS32 SubImageStartRow = SubImageCenterRow - radiusS32; + psS32 SubImageEndRow = SubImageCenterRow + radiusS32; + psS32 SubImageStartCol = SubImageCenterCol - radiusS32; + psS32 SubImageEndCol = SubImageCenterCol + radiusS32; + + if (SubImageStartRow < 0) { + psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", + SubImageStartRow); + return(false); + } + if (SubImageEndRow+1 >= 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+1 >= image->numCols) { + psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", + SubImageEndCol); + return(false); + } + + // XXX: Must recycle image. + if (source->pixels != NULL) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL. Freeing and reallocating.\n"); + psFree(source->pixels); + } + source->pixels = psImageSubset((psImage *) image, + SubImageStartCol, + SubImageStartRow, + SubImageEndCol+1, + SubImageEndRow+1); + + // XXX: Must recycle image. + if (source->mask != NULL) { + psFree(source->mask); + } + source->mask = psImageAlloc(1 + 2 * radiusS32, 1 + 2 * radiusS32, PS_TYPE_F32); + + // + // Loop through the subimage mask, initialize mask to 0 or 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] = 1; + } else { + source->mask->data.U8[row][col] = 1; + } + } + } + + /* + for (psS32 row = SubImageCenterRow - radiusS32; row <= SubImageCenterRow + radiusS32; row++) { + for (psS32 col = SubImageCenterCol - radiusS32; col <= SubImageCenterCol + radiusS32; col++) { + if (CheckRadius(source->peak, radius, (psF32) col, (psF32) row)) { + source->mask->data.U8[row-SubImageCenterRow][col-SubImageCenterCol] = 1; + } + } + } + */ + return(true); } @@ -263,13 +891,170 @@ /****************************************************************************** +pmSourceModelGuess(source, image, model): This function allocates a new +psModel structure and store it in the psSource data structure specified in the +argument list. The model type is specified in the argument list. The params +array in that psModel 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. *****************************************************************************/ bool pmSourceModelGuess(psSource *source, - const psImage *image) -{ - return(true); -} - - -/****************************************************************************** + const psImage *image, + psModelType model) +{ + PS_PTR_CHECK_NULL(source, false); + PS_PTR_CHECK_NULL(source->moments, false); + PS_PTR_CHECK_NULL(source->peak, false); + PS_IMAGE_CHECK_NULL(image, false); + PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); + if (source->models != NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->models was non-NULL; calling psFree(source->models).\n"); + psFree(source->models); + } + source->models = pmModelAlloc(PS_MODEL_UNDEFINED); + + switch (model) { + case PS_MODEL_GAUSS: + source->models->type = PS_MODEL_GAUSS; + source->models->Nparams = 7; + source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); + for (psS32 i = 0 ; i < 7 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + source->models->params[4] = sqrt(2.0) / source->moments->Sx; + source->models->params[5] = sqrt(2.0) / source->moments->Sy; + source->models->params[6] = source->moments->Sxy; + source->models->chisq = 0.0; + return(true); + case PS_MODEL_PGAUSS: + source->models->type = PS_MODEL_PGAUSS; + source->models->Nparams = 7; + source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); + for (psS32 i = 0 ; i < 7 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + source->models->params[4] = sqrt(2.0) / source->moments->Sx; + source->models->params[5] = sqrt(2.0) / source->moments->Sy; + source->models->params[6] = source->moments->Sxy; + source->models->chisq = 0.0; + return(true); + case PS_MODEL_TWIST_GAUSS: + source->models->type = PS_MODEL_TWIST_GAUSS; + source->models->Nparams = 11; + source->models->params = (psF32 *) psAlloc(11 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(11 * sizeof(psF32)); + for (psS32 i = 0 ; i < 11 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + // XXX: What are these? + // source->models->params[4] = SxInner; + // source->models->params[5] = SyInner; + // source->models->params[6] = SxyInner; + // source->models->params[7] = SxOuter; + // source->models->params[8] = SyOuter; + // source->models->params[9] = SxyOuter; + // source->models->params[10] = N; + + source->models->chisq = 0.0; + return(true); + case PS_MODEL_WAUSS: + + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + source->models->params[4] = sqrt(2.0) / source->moments->Sx; + source->models->params[5] = sqrt(2.0) / source->moments->Sy; + source->models->params[6] = source->moments->Sxy; + // XXX: What are these? + // source->models->params[7] = B2; + // source->models->params[8] = B3; + + source->models->type = PS_MODEL_WAUSS; + source->models->Nparams = 9; + source->models->params = (psF32 *) psAlloc(9 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(9 * sizeof(psF32)); + for (psS32 i = 0 ; i < 9 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + source->models->chisq = 0.0; + return(true); + case PS_MODEL_SERSIC: + source->models->type = PS_MODEL_SERSIC; + source->models->Nparams = 8; + source->models->params = (psF32 *) psAlloc(8 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(8 * sizeof(psF32)); + for (psS32 i = 0 ; i < 8 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + source->models->params[4] = sqrt(2.0) / source->moments->Sx; + source->models->params[5] = sqrt(2.0) / source->moments->Sy; + source->models->params[6] = source->moments->Sxy; + // XXX: What are these? + //source->models->params[7] = Nexp; + + source->models->chisq = 0.0; + return(true); + case PS_MODEL_SERSIC_CORE: + source->models->type = PS_MODEL_SERSIC_CORE; + source->models->Nparams = 12; + source->models->params = (psF32 *) psAlloc(12 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(12 * sizeof(psF32)); + for (psS32 i = 0 ; i < 12 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + // XXX: What are these? + //source->models->params[4] SxInner; + //source->models->params[5] SyInner; + //source->models->params[6] SxyInner; + //source->models->params[7] Zd; + //source->models->params[8] SxOuter; + //source->models->params[9] SyOuter; + //source->models->params[10] = SxyOuter; + //source->models->params[11] = Nexp; + + source->models->chisq = 0.0; + return(true); + default: + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + return(false); + } +} + + +/****************************************************************************** + *****************************************************************************/ psArray *pmSourceContour(psSource *source, @@ -281,11 +1066,169 @@ } -/****************************************************************************** +psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x); +psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); +psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x); +psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); +psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x); +psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x); + +//XXX: What should these values be? +#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100 +#define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0 +/****************************************************************************** +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? *****************************************************************************/ bool pmSourceFitModel(psSource *source, const psImage *image) { + PS_PTR_CHECK_NULL(source, false); + PS_PTR_CHECK_NULL(source->moments, false); + PS_PTR_CHECK_NULL(source->peak, false); + PS_PTR_CHECK_NULL(source->pixels, false); + PS_PTR_CHECK_NULL(source->models, false); + PS_IMAGE_CHECK_NULL(image, false); + PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); + + psBool rc; + 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++; + } + } + } + psArray *x = psArrayAlloc(count); + psVector *y = psVectorAlloc(count, PS_TYPE_F32); + 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: Should we use the subimage offsets here, or does it not matter? + coord->data.F32[0] = (psF32) (i); + coord->data.F32[1] = (psF32) (j); + x->data[count] = (psPtr *) coord; + y->data.F32[count] = source->pixels->data.F32[i][j]; + } + } + } + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32); + + switch (source->models->type) { + case PS_MODEL_GAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec); + break; + case PS_MODEL_PGAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec); + break; + case PS_MODEL_TWIST_GAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec); + break; + case PS_MODEL_WAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec); + break; + case PS_MODEL_SERSIC: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec); + break; + case PS_MODEL_SERSIC_CORE: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + rc = false; + } + + psFree(x); + psFree(y); + psFree(myMin); + psFree(params); + return(rc); +} + +bool p_pmSourceAddOrSubModel(psImage *image, + psSource *source, + bool center, + psS32 flag) +{ + PS_PTR_CHECK_NULL(source, false); + PS_PTR_CHECK_NULL(source->moments, false); + PS_PTR_CHECK_NULL(source->peak, false); + PS_PTR_CHECK_NULL(source->pixels, false); + PS_PTR_CHECK_NULL(source->models, false); + PS_IMAGE_CHECK_NULL(image, false); + PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); + + // XXX: make static, document. + psVector *deriv = psVectorAlloc(100, PS_TYPE_F32); + psVector *params = psVectorAlloc(100, PS_TYPE_F32); + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + for (psS32 i = 0 ; i < source->models->Nparams ; i++) { + params->data.F32[i] = source->models->params[i]; + } + + for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { + for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { + psF32 pixelValue; + // XXX: Should you use offsets here? + // XXX: Make sure you have col/row order correct. + x->data.F32[0] = (float) j; + x->data.F32[1] = (float) i; + switch (source->models->type) { + case PS_MODEL_GAUSS: + pixelValue = pmMinLM_Gauss2D(deriv, params, x); + break; + case PS_MODEL_PGAUSS: + pixelValue = pmMinLM_PsuedoGauss2D(deriv, params, x); + break; + case PS_MODEL_TWIST_GAUSS: + pixelValue = pmMinLM_TwistGauss2D(deriv, params, x); + break; + case PS_MODEL_WAUSS: + pixelValue = pmMinLM_Wauss2D(deriv, params, x); + break; + case PS_MODEL_SERSIC: + pixelValue = pmMinLM_Sersic(deriv, params, x); + break; + case PS_MODEL_SERSIC_CORE: + pixelValue = pmMinLM_SersicCore(deriv, params, x); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + psFree(x); + psFree(deriv); + psFree(params); + 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. + psS32 imageRow = 0; + psS32 imageCol = 0; + image->data.F32[imageRow][imageCol]+= pixelValue; + } + } + psFree(x); + psFree(deriv); + psFree(params); return(true); } + + /****************************************************************************** @@ -295,18 +1238,32 @@ bool center) { - return(true); -} - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceSubModel(psSource *source) -{ - return(true); -} - -/****************************************************************************** -XXX: Why only *x argument? - -May x[0] is x and x[1] is y? + return(p_pmSourceAddOrSubModel(image, source, center, 0)); +} + +/****************************************************************************** + *****************************************************************************/ +bool pmSourceSubModel(psImage *image, + psSource *source, + bool center) +{ + return(p_pmSourceAddOrSubModel(image, source, 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); \ +} + + +/****************************************************************************** +pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x, +y" coordinate pair. This function computes the gaussian, specified by the +parameters in "params" at that x,y point and returns the value. The +derivatives are also caculated and returned in the "deriv" argument. params->data.F32[0] = So; @@ -317,4 +1274,7 @@ params->data.F32[5] = sqrt(2.0) / SigmaY; params->data.F32[6] = Sxy; + +XXX: Consider getting rid of the parameter checks since this might consume +a significant fraction of this function CPU time. *****************************************************************************/ psF32 pmMinLM_Gauss2D(psVector *deriv, @@ -322,4 +1282,14 @@ psVector *x) { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 7, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + psF32 X = x->data.F32[0] - params->data.F32[2]; psF32 Y = x->data.F32[1] - params->data.F32[3]; @@ -343,4 +1313,45 @@ /****************************************************************************** +p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above +function in a form that is usable in the LM minimization routines. + *****************************************************************************/ +psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, + psVector *params, + psArray *x) +{ + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); +} + + +/****************************************************************************** params->data.F32[0] = So; params->data.F32[1] = Zo; @@ -355,4 +1366,14 @@ psVector *x) { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 7, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + psF32 X = x->data.F32[0] - params->data.F32[2]; psF32 Y = x->data.F32[1] - params->data.F32[3]; @@ -382,4 +1403,47 @@ /****************************************************************************** +p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the +above function in a form that is usable in the LM minimization routines. + *****************************************************************************/ +psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, + psVector *params, + psArray *x) +{ + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); +} + + + + +/****************************************************************************** params->data.F32[0] = So; params->data.F32[1] = Zo; @@ -396,4 +1460,14 @@ psVector *x) { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 9, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 9, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + psF32 X = x->data.F32[0] - params->data.F32[2]; psF32 Y = x->data.F32[1] - params->data.F32[2]; @@ -424,4 +1498,49 @@ return(f); } + +/****************************************************************************** +p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above +function in a form that is usable in the LM minimization routines. + *****************************************************************************/ +psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, + psVector *params, + psArray *x) +{ + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); +} + + + + + // XXX: What should these be? @@ -445,4 +1564,14 @@ psVector *x) { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 11, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 11, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + psF32 X = x->data.F32[0] - params->data.F32[2]; psF32 Y = x->data.F32[1] - params->data.F32[3]; @@ -489,4 +1618,46 @@ /****************************************************************************** +p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above +function in a form that is usable in the LM minimization routines. + *****************************************************************************/ +psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, + psVector *params, + psArray *x) +{ + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); +} + + + +/****************************************************************************** float Sersic() params->data.F32[0] = So; @@ -503,5 +1674,55 @@ psVector *x) { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 8, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 8, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + + psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); return(0.0); +} +/****************************************************************************** +p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above +function in a form that is usable in the LM minimization routines. + *****************************************************************************/ +psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, + psVector *params, + psArray *x) +{ + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); } @@ -525,6 +1746,58 @@ psVector *x) { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 12, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 12, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + + psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); return(0.0); } +/****************************************************************************** +p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above +function in a form that is usable in the LM minimization routines. + *****************************************************************************/ +psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, + psVector *params, + psArray *x) +{ + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); +} + + /****************************************************************************** Index: /trunk/psModules/src/pmObjects.h =================================================================== --- /trunk/psModules/src/pmObjects.h (revision 3497) +++ /trunk/psModules/src/pmObjects.h (revision 3498) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-01-25 02:45:43 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-24 22:36:52 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -28,5 +28,6 @@ PM_PEAK_LONE, // Isolated peak. PM_PEAK_EDGE, // Peak on edge. - PM_PEAK_FLAT // Peak has equal-value neighbors. + PM_PEAK_FLAT, // Peak has equal-value neighbors. + PM_PEAK_UNDEF // Undefined. } psPeakType; @@ -60,5 +61,6 @@ PS_MODEL_WAUSS, PS_MODEL_SERSIC, - PS_MODEL_SERSIC_CORE + PS_MODEL_SERSIC_CORE, + PS_MODEL_UNDEFINED } psModelType; @@ -75,10 +77,12 @@ // XXX: What is this enum? typedef enum { - PS_SOURCE_TYPE_1, - PS_SOURCE_TYPE_2, - PS_SOURCE_TYPE_3, - PS_SOURCE_TYPE_4, - PS_SOURCE_TYPE_5, - PS_SOURCE_TYPE_6 + PS_SOURCE_PSFSTAR, + PS_SOURCE_GALAXY, + PS_SOURCE_DEFECT, + PS_SOURCE_SATURATED, + PS_SOURCE_SATSTAR, + PS_SOURCE_FAINTSTAR, + PS_SOURCE_BRIGHTSTAR, + PS_SOURCE_OTHER } psSourceType; @@ -94,4 +98,11 @@ psSource; +psPeak *pmPeakAlloc(psS32 x, + psS32 y, + psF32 counts, + psPeakType class); +psMoments *pmMomentsAlloc(); +psModel *pmModelAlloc(psModelType type); +psSource *pmSourceAlloc(); /****************************************************************************** @@ -131,4 +142,5 @@ psSource *pmSourceLocalSky(const psImage *image, const psPeak *peak, + psStatsOptions statsOptions, psF32 innerRadius, psF32 outerRadius); @@ -137,17 +149,12 @@ *****************************************************************************/ psSource *pmSourceMoments(psSource *source, - const psImage *image, psF32 radius); /****************************************************************************** -pmSourceRoughClass(source, saturate, SNlim, valid): make a guessat the source -classification. - *****************************************************************************/ -psSource *pmSourceRoughClass(psSource *source, - psF32 saturate, - float SNlim, - const psRegion *valid); - - +pmSourceRoughClass(pmArray *source, psMetaDeta *metadata): make a guess at the +source classification. + *****************************************************************************/ +bool pmSourceRoughClass(psArray *source, + psMetadata *metadata); /****************************************************************************** pmSourceSetPixelCircle(source, image, radius) @@ -161,5 +168,6 @@ *****************************************************************************/ bool pmSourceModelGuess(psSource *source, - const psImage *image); + const psImage *image, + psModelType model); /****************************************************************************** @@ -183,5 +191,7 @@ /****************************************************************************** *****************************************************************************/ -bool pmSourceSubModel(psSource *source); +bool pmSourceSubModel(psImage *image, + psSource *source, + bool center); /******************************************************************************