Index: /trunk/psModules/src/Makefile.am =================================================================== --- /trunk/psModules/src/Makefile.am (revision 3876) +++ /trunk/psModules/src/Makefile.am (revision 3877) @@ -1,10 +1,25 @@ lib_LIBRARIES = libpsmodule.a -libpsmodule_a_SOURCES = pmFlatField.c pmMaskBadPixels.c pmSubtractBias.c \ - pmSubtractBias.c pmNonLinear.c pmReadoutCombine.c pmSubtractSky.c pmObjects.c +libpsmodule_a_SOURCES = pmFlatField.c \ + pmImageCombine.c \ + pmImageSubtract.c \ + pmMaskBadPixels.c\ + pmNonLinear.c \ + pmObjects.c \ + pmReadoutCombine.c \ + pmSubtractBias.c \ + pmSubtractSky.c -libpsmodule_a_HEADERS = pmFlatFieldErrors.h pmMaskBadPixelsErrors.h \ - pmFlatField.h pmMaskBadPixels.h pmSubtractBias.h pmSubtractBias.h \ - pmNonLinear.h pmReadoutCombine.h pmSubtractSky.h pmObjects.h +libpsmodule_a_HEADERS = pmFlatFieldErrors.h \ + pmMaskBadPixelsErrors.h \ + pmFlatField.h \ + pmImageCombine.h \ + pmImageSubtract.h \ + pmMaskBadPixels.h \ + pmNonLinear.h \ + pmObjects.h \ + pmReadoutCombine.h \ + pmSubtractBias.h \ + pmSubtractSky.h EXTRA_DIST = psErrorCodes.dat Index: /trunk/psModules/src/pmImageCombine.c =================================================================== --- /trunk/psModules/src/pmImageCombine.c (revision 3877) +++ /trunk/psModules/src/pmImageCombine.c (revision 3877) @@ -0,0 +1,403 @@ +/** @file pmImageCombine.c + * + * This file will ... + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-05-10 23:48:19 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + * XXX: None of this has been tested. + * + * + * + * + * + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psConstants.h" +// XXX: Remove this, once psPixels.h is handled properly. +#include "psPixels.h" + + +/****************************************************************************** +pmCombineImages(combine, questionablePixels, images, errors, masks, maskVal, + pixels, numIter, sigmaClip, stats) + *****************************************************************************/ +psImage *pmCombineImages(psImage *combine, ///< Combined image (output) + psArray **questionablePixels, ///< Array of rejection masks + const psArray *images, ///< Array of input images + const psArray *errors, ///< Array of input error images + const psArray *masks, ///< Array of input masks + psU32 maskVal, ///< Mask value + const psPixels *pixels, ///< Pixels to combine + psS32 numIter, ///< Number of rejection iterations + psF32 sigmaClip, ///< Number of standard deviations at which to reject + const psStats *stats ///< Statistics to use in the combination + ) +{ + PS_PTR_CHECK_NULL(images, combine); + psU32 numImages = images->n; + + if (errors != NULL) { + if (images->n != errors->n) { + psError(PS_ERR_UNKNOWN, true, "images and errors args must have same length (%d != %d)\n", + images->n, errors->n); + return(combine); + } + } + if (masks != NULL) { + if (images->n != masks->n) { + psError(PS_ERR_UNKNOWN, true, "images and masks args must have same length (%d != %d)\n", + images->n, masks->n); + return(combine); + } + } + + psImage *tmpImg = (psImage *) images->data[0]; + psU32 numRows = tmpImg->numRows; + psU32 numCols = tmpImg->numCols; + + // + // Check that all images have the appropriate size and type. + // + for (psS32 i = 0 ; i < numImages ; i++) { + psImage *tmpDataImg = (psImage *) images->data[i]; + psImage *tmpErrorImg = (psImage *) errors->data[i]; + psImage *tmpMaskImg = (psImage *) masks->data[i]; + + PS_IMAGE_CHECK_NULL(tmpDataImg, combine); + PS_IMAGE_CHECK_NULL(tmpErrorImg, combine); + PS_IMAGE_CHECK_NULL(tmpMaskImg, combine); + + PS_IMAGE_CHECK_TYPE(tmpDataImg, PS_TYPE_F32, combine); + PS_IMAGE_CHECK_TYPE(tmpErrorImg, PS_TYPE_F32, combine); + PS_IMAGE_CHECK_TYPE(tmpMaskImg, PS_TYPE_U8, combine); + + if ((tmpDataImg->numRows != numRows) || (tmpDataImg->numCols != numCols)) { + psError(PS_ERR_UNKNOWN, true, "error image %d has size (%d, %d); should be (%d, %d).\n", + i, tmpDataImg->numRows, tmpDataImg->numCols, numRows, numCols); + } + + PS_IMAGE_CHECK_SIZE_EQUAL(tmpDataImg, tmpErrorImg, NULL); + PS_IMAGE_CHECK_SIZE_EQUAL(tmpDataImg, tmpMaskImg, NULL); + } + PS_PTR_CHECK_NULL(pixels, combine); + PS_PTR_CHECK_NULL(stats, combine); + + // Allocate an initialize the combined image, if necessary. + if (combine == NULL) { + combine = psImageAlloc(numCols, numRows, PS_TYPE_F32); + if (pixels != NULL) { + PS_IMAGE_SET_F32(combine, 0.0); + } + } + + // + // Allocate the necessary psVectors for the call to psVectorStats(). + // These vectors will be used whether we are combining a list of pixels, + // or every pixel in the input images. + // + psVector *pixelData = psVectorAlloc(numImages, PS_TYPE_F32); + + psVector *pixelMask = NULL; + if (masks != NULL) { + pixelMask = psVectorAlloc(numImages, PS_TYPE_F32); + PS_VECTOR_SET_U8(pixelMask, 0); + } + + psVector *pixelErrors = NULL; + if (errors != NULL) { + pixelErrors = psVectorAlloc(numImages, PS_TYPE_F32); + PS_VECTOR_SET_F32(pixelErrors, 1.0); + } + + if (pixels != NULL) { + // Only those specified pixels should be combined. + psStats *stdevStats = psStatsAlloc(PS_STAT_SAMPLE_STDEV); + for (psS32 p = 0 ; p < pixels->n ; p++) { + // Must initialize the mask to 0 for every pixel. + PS_VECTOR_SET_U8(pixelMask, 0); + psS32 col = (pixels->data[p]).x; + psS32 row = (pixels->data[p]).y; + + // + // Loop through each image, extract the pixel/mask/error data + // into psVectors. + // + for (psS32 im = 0 ; im < numImages ; im++) { + // Set the pixel data + pixelData->data.F32[im] = ((psImage *) images->data[im])->data.F32[row][col]; + + // Set the pixel mask data, if necessary + if (masks != NULL) { + pixelMask->data.F32[im] = ((psImage *) masks->data[im])->data.F32[row][col]; + } + + // Set the pixel error data, if necessary + if (errors != NULL) { + pixelErrors->data.F32[im] = ((psImage *) errors->data[im])->data.F32[row][col]; + } + } + + // + // Iterate on the pixels, rejecting outliers + // + for (psS32 iter = 0 ; iter < numIter ; iter++) { + // Combine all the pixels, using the specified stat. + stats = psVectorStats((psStats *) stats, pixelData, pixelErrors, pixelMask, maskVal); + psF64 combinedPixel; + psBool rc = p_psGetStatValue(stats, &combinedPixel); + if (rc != true) { + // XXX: Warning. + } + if (iter == 0) { + combine->data.F32[row][col] = (psF32) combinedPixel; + } + + // + // Reject all pixels that lie more that sigmaClip standard deviations from + // the combined pixel value. + // + for (psS32 im = 0 ; im < numImages ; im++) { + stats = psVectorStats(stdevStats, pixelData, pixelErrors, pixelMask, maskVal); + psF64 stdev; + psBool rc = p_psGetStatValue(stdevStats, &stdev); + if (rc != true) { + // XXX: Warning. + } + + if (fabs(pixelData->data.F32[im] - combinedPixel) > fabs(sigmaClip * stdev)) { + pixelMask->data.F32[im] = maskVal; + // + // XXX: Must Code This. + // XXX: Add this to the list of questionable pixels. + // + } + } + } + } + psFree(stdevStats); + } else { + // + // We get here if there is a NULL list of pixels to combine. + // Therefore, we combine all pixels in all images. + // + + // + // Loop over all pixels in all images, set the appropriate data, mask, + // error vectors, call psVectorStats(), and set the result in the combine + // image. + // + for (psS32 row = 0 ; row < numRows ; row++) { + for (psS32 col = 0 ; col < numCols ; col++) { + for (psS32 im = 0 ; im < numImages ; im++) { + // Set the pixel data + pixelData->data.F32[im] = ((psImage *) images->data[im])->data.F32[row][col]; + + // Set the pixel mask data, if necessary + if (masks != NULL) { + pixelMask->data.F32[im] = ((psImage *) masks->data[im])->data.F32[row][col]; + } + + // Set the pixel error data, if necessary + if (errors != NULL) { + pixelErrors->data.F32[im] = ((psImage *) errors->data[im])->data.F32[row][col]; + } + + } + // Combine all the pixels, using the specified stat. + stats = psVectorStats((psStats *) stats, pixelData, pixelErrors, pixelMask, maskVal); + psF64 tmpF64; + psBool rc = p_psGetStatValue(stats, &tmpF64); + combine->data.F32[row][col] = (psF32) tmpF64; + if (rc != true) { + // XXX: Warning. + } + } + } + } + + psFree(pixelData); + psFree(pixelMask); + psFree(pixelErrors); + + + return(combine); +} + + +/****************************************************************************** +XXX: Directly from Paul Price +XXX: Confirm that this gradient calculation is still what's required. + *****************************************************************************/ +static psF32 CalcGradient(psImage *image, + psS32 x, + psS32 y) +{ + int num = 0; + psVector *pixels = psVectorAlloc(8, PS_TYPE_F32); // Array of pixels + psVector *mask = psVectorAlloc(8, PS_TYPE_U8); // Corresponding mask + + // Get limits + int xMin = PS_MAX(x - 1, 0); + int xMax = PS_MIN(x + 1, image->numCols - 1); + int yMin = PS_MAX(y - 1, 0); + int yMax = PS_MIN(y + 1, image->numRows - 1); + for (int j = yMin; j <= yMax; j++) { + for (int i = xMin; i <= xMax; i++) { + if ((i != x) && (j != y)) { + pixels->data.F32[num] = image->data.F32[j][i]; + mask->data.U8[num] = 0; + num++; + } + } + } + pixels->n = num; + mask->n = num; + + // Get the median + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + (void)psVectorStats(stats, pixels, NULL, mask, 1); + float median = stats->sampleMedian; + psFree(stats); + psFree(pixels); + psFree(mask); + return median / image->data.F32[y][x]; +} + +/****************************************************************************** +pmRejectPixels(images, errors, inToOut, outToIn, rejThreshold, +gradLimit) + +XXX: Use static variables. +XXX: Optimization: we don't need to transform the entire mask image. +XXX: The inToOut and outToIn transforms are confusing. Verify that what I think + they mean syncs with PWP. + *****************************************************************************/ +psArray *pmRejectPixels(const psArray *images, ///< Array of input images + const psArray *errors, ///< The pixels which were rejected in the combination + const psArray *inToOut, ///< Transformation from input to output system + const psArray *outToIn, ///< Transformation from output to input system + psF32 rejThreshold, ///< Rejection threshold + psF32 gradLimit ///< Gradient limit + ) +{ + PS_PTR_CHECK_NULL(images, NULL); + PS_PTR_CHECK_NULL(errors, NULL); + PS_PTR_CHECK_NULL(inToOut, NULL); + PS_PTR_CHECK_NULL(outToIn, NULL); + // Ensure that the psArray parameters have an element for each image. + psS32 numImages = images->n; + PS_INT_CHECK_NON_EQUALS(numImages, errors->n, NULL); + PS_INT_CHECK_NON_EQUALS(numImages, inToOut->n, NULL); + PS_INT_CHECK_NON_EQUALS(numImages, outToIn->n, NULL); + // XXX: Loop through all images, ensure their sizes are equal? + + psS32 numCols = ((psImage *) images->data[0])->numCols; + psS32 numRows = ((psImage *) images->data[0])->numRows; + psRegion *myRegion = psRegionAlloc(0, numCols-1, 0, numRows-1); + psU32 maskVal = 1; + + psPlane *inCoords = psAlloc(sizeof(psPlane)); + psPlane *outCoords = psAlloc(sizeof(psPlane)); + + for (psS32 img = 0 ; img < numImages ; img++) { + // + // Extract data from psArrays. + // + psPixels *pixelList = (psPixels *) errors->data[img]; + psImage *currImage = (psImage *) images->data[img]; + myRegion->x0 = 0; + myRegion->x1 = currImage->numCols - 1; + myRegion->y0 = 0; + myRegion->y1 = currImage->numRows - 1; + psPlaneTransform *myInToOut = (psPlaneTransform *) inToOut->data[img]; + psPlaneTransform *myOutToIn = (psPlaneTransform *) outToIn->data[img]; + + // + // Create a psU8 mask image from the list of cosmic pixels. + // + psImage *maskImage = NULL; + maskImage = psPixelsToMask(maskImage, pixelList, myRegion, maskVal); + + // + // Transform that mask image into detector coordinate space + // + psImage *transformedImage = psImageTransform(NULL, maskImage, NULL, 0, myOutToIn, NULL, 0); + + // + // Loop over all cosmic pixels. Transform their coords to detector space. + // If the value of the transformed mask is larger than rejThreshold, then + // this might be a cosmic ray pixel. We then calculate the mean gradient + // in other images. + // + for (psS32 p = 0 ; p < pixelList->n ; p++) { + inCoords->x = 0.5 + (psF32) (pixelList->data)->x; + inCoords->y = 0.5 + (psF32) (pixelList->data)->y; + psPlaneTransformApply(outCoords, myInToOut, inCoords); + + psF32 maskVal = (psF32) psImagePixelInterpolate(transformedImage, outCoords->x, outCoords->y, + NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); + if (maskVal > rejThreshold) { + + // This is a possible cosmic array pixel. We must calculate the gradient + // at this location in all input images. + psF32 meanGrads = 0.0; + psS32 numGrads = 0; + + // + // Loop through all other images, calculate their mean gradient. + // + for (psS32 otherImg = 0 ; otherImg < numImages ; otherImg++) { + if (img != otherImg) { + // Map the outCoords to inCoords that for otherImg space. + psPlaneTransformApply(inCoords, (psPlaneTransform * )outToIn->data[otherImg], outCoords); + + psS32 xPix = (int)(inCoords->x + 0.5); + psS32 yPix = (int)(inCoords->y + 0.5); + if ((xPix >= 0) && (xPix <= ((psImage*)(images->data[otherImg]))->numCols - 1) && + (yPix >= 0) && (yPix <= ((psImage*)(images->data[otherImg]))->numRows - 1)) { + meanGrads += CalcGradient(images->data[otherImg], xPix, yPix); + numGrads++; + } + + } + } + if (numGrads > 0) { + meanGrads /= (psF32) numGrads; + } else { + // XXX: my idea. Verify with PWP: + meanGrads = 1.0 + gradLimit; + } + + // XXX: The SDRS and the prototype code differ significantly here: + // if (CalcGradient(inputs->data[i], pixelList->data.x, pixelList->data.y) < (gradLimit * meanGrads)) { + if (meanGrads < gradLimit) { + // + // XXX: Add this to the list of questionable pixels. + // XXX: Must Code This. + // + } + } + } + psFree(myInToOut); + psFree(myOutToIn); + psFree(maskImage); + psFree(transformedImage); + } + + psFree(myRegion); + psFree(inCoords); + psFree(outCoords); + return(NULL); +} +// This code is + Index: /trunk/psModules/src/pmObjects.c =================================================================== --- /trunk/psModules/src/pmObjects.c (revision 3876) +++ /trunk/psModules/src/pmObjects.c (revision 3877) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-04-19 23:44:54 $ + * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-05-10 23:48:19 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -20,5 +20,5 @@ /****************************************************************************** pmPeakAlloc(): Allocate the psPeak data structure and set appropriate members. - *****************************************************************************/ +*****************************************************************************/ psPeak *pmPeakAlloc(psS32 x, psS32 y, @@ -38,5 +38,5 @@ pmMomentsAlloc(): Allocate the psMoments structure and initialize the members to zero. - *****************************************************************************/ +*****************************************************************************/ psMoments *pmMomentsAlloc() { @@ -45,5 +45,5 @@ tmp->y = 0.0; tmp->Sx = 0.0; - tmp->Sy = 0.0; + tmp->Sx = 0.0; tmp->Sxy = 0.0; tmp->Sum = 0.0; @@ -64,5 +64,6 @@ pmModelAlloc(): Allocate the psModel structure, along with its parameters, and initialize the type member. Initialize the params to 0.0. - *****************************************************************************/ +XXX EAM: changing params and dparams to psVector +*****************************************************************************/ psModel *pmModelAlloc(psModelType type) { @@ -73,32 +74,26 @@ switch (type) { case PS_MODEL_GAUSS: - tmp->Nparams = 7; - tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); - tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); + tmp->params = psVectorAlloc(7, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(7, PS_TYPE_F32); break; case PS_MODEL_PGAUSS: - tmp->Nparams = 7; - tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); - tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); + tmp->params = psVectorAlloc(7, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(7, PS_TYPE_F32); break; case PS_MODEL_TWIST_GAUSS: - tmp->Nparams = 11; - tmp->params = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32)); - tmp->dparams = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32)); + tmp->params = psVectorAlloc(11, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(11, PS_TYPE_F32); break; case PS_MODEL_WAUSS: - tmp->Nparams = 9; - tmp->params = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32)); - tmp->dparams = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32)); + tmp->params = psVectorAlloc(9, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(9, PS_TYPE_F32); break; case PS_MODEL_SERSIC: - tmp->Nparams = 8; - tmp->params = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32)); - tmp->dparams = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32)); + tmp->params = psVectorAlloc(8, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(8, PS_TYPE_F32); break; case PS_MODEL_SERSIC_CORE: - tmp->Nparams = 12; - tmp->params = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32)); - tmp->dparams = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32)); + tmp->params = psVectorAlloc(12, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(12, PS_TYPE_F32); break; default: @@ -107,7 +102,7 @@ } - for (psS32 i = 0 ; i < tmp->Nparams ; i++) { - tmp->params[i] = 0.0; - tmp->dparams[i] = 0.0; + for (psS32 i = 0; i < tmp->params->n; i++) { + tmp->params->data.F32[i] = 0.0; + tmp->dparams->data.F32[i] = 0.0; } @@ -117,13 +112,14 @@ /****************************************************************************** -p_psSourceFree(tmp): a private function which frees the psSource struct. - *****************************************************************************/ -static void p_psSourceFree(psSource *tmpSrc) -{ - psFree(tmpSrc->peak); - psFree(tmpSrc->pixels); - psFree(tmpSrc->mask); - psFree(tmpSrc->moments); - psFree(tmpSrc->models); +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); } @@ -131,16 +127,17 @@ pmSourceAlloc(): Allocate the psSource structure and initialize its members to NULL. - *****************************************************************************/ +*****************************************************************************/ psSource *pmSourceAlloc() { - psSource *tmpSrc = (psSource *) psAlloc(sizeof(psSource)); - tmpSrc->peak = NULL; - tmpSrc->pixels = NULL; - tmpSrc->mask = NULL; - tmpSrc->moments = NULL; - tmpSrc->models = NULL; - psMemSetDeallocator(tmpSrc, (psFreeFcn) p_psSourceFree); - - return(tmpSrc); + psSource *tmp = (psSource *) psAlloc(sizeof(psSource)); + tmp->peak = NULL; + tmp->pixels = NULL; + tmp->mask = NULL; + tmp->moments = NULL; + tmp->models = NULL; + tmp->type = 0; + psMemSetDeallocator(tmp, (psFreeFcn) p_psSourceFree); + + return(tmp); } @@ -154,7 +151,6 @@ XXX: We currently step through the input vector twice; once to determine the size of the output vector, then to set the values of the output vector. -Depending upon actual use, this may need to be optimized. Use the function -which adds to a psVector. It's not clear which way is faster. - *****************************************************************************/ +Depending upon actual use, this may need to be optimized. +*****************************************************************************/ psVector *pmFindVectorPeaks(const psVector *vector, psF32 threshold) @@ -250,7 +246,5 @@ XXX: Is there a better way to do this? - -XXX: Use memcpy() on the data transfer. - *****************************************************************************/ +*****************************************************************************/ psVector *p_psGetRowVectorFromImage(psImage *image, psU32 row) @@ -267,202 +261,62 @@ /****************************************************************************** -p_psGetColVectorFromImage(): a private function which simply returns a -psVector containing the specified col of data from the psImage. - -XXX: Is there a better way to do this? - -XXX: Use memcpy() on the data transfer. - *****************************************************************************/ -psVector *p_psGetColVectorFromImage(psImage *image, - psU32 col) +MyListAddPeak(): A private function which allocates a psArray, if the list +argument is NULL, otherwise it adds the peak to that list. +XXX EAM : changed the output to psArray +XXX EAM : Switched row, col args +XXX EAM : NOTE: this was changed in the call, so the new code is consistent +*****************************************************************************/ +psArray *MyListAddPeak(psArray *list, + psS32 row, + psS32 col, + psF32 counts, + psPeakType type) +{ + psPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); + + if (list == NULL) { + list = psArrayAlloc(100); + list->n = 0; + } + psArrayAdd(list, 100, tmpPeak); + + return(list); +} + +/****************************************************************************** +pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage +above the given threshold. Returns a psArray containing location (x/y value) +of all peaks. + +XXX: I'm not convinced the peak type definition in the SDRS is mutually +exclusive. Some peaks can have multiple types. Edges for sure. Also, a +digonal line with the same value at each point will have a peak for every +point on that line. + +XXX: This does not work if image has either a single row, or a single column. + +XXX: In the output psArray elements, should we use the image row/column offsets? + Currently, we do not. +*****************************************************************************/ +psArray *pmFindImagePeaks(const psImage *image, + psF32 threshold) { PS_IMAGE_CHECK_NULL(image, NULL); PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL); - - psVector *tmpVector = psVectorAlloc(image->numRows, PS_TYPE_F32); - for (psU32 row = 0; row < image->numRows ; row++) { - tmpVector->data.F32[row] = image->data.F32[row][col]; - } - return(tmpVector); -} - -/****************************************************************************** -MyListAddPeak(): A private function which allocates a psList, if the list -argument is NULL, otherwise it adds the peak to that list. - *****************************************************************************/ -psList *MyListAddPeak(psList *list, - psS32 col, - psS32 row, - 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); -} - -/****************************************************************************** -pmFindImagePeaks(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. - *****************************************************************************/ -psList *pmFindImagePeaks(const psImage *image, - psF32 threshold) -{ - PS_IMAGE_CHECK_NULL(image, NULL); - PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL); - psPeakType myPeakClass = PM_PEAK_UNDEF; + 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; - - // - // Special case: a 1-by-1 image. - // - if ((image->numCols == 1) && (image->numRows == 1)) { - if (image->data.F32[row][col] > threshold) { - myPeakClass = PM_PEAK_EDGE | PM_PEAK_LONE; - list = MyListAddPeak(list, 0, 0, image->data.F32[row][col], myPeakClass); - return(list); - } else { - return(NULL); - } - } - - // - // Find peaks in row 0 only (single-row image). This is a special case since - // we can not test data at image->data.F32[row+1][...]) - // - if (image->numRows == 1) { - row = 0; - tmpRow = p_psGetRowVectorFromImage((psImage *) image, row); - psVector *row1 = pmFindVectorPeaks(tmpRow, threshold); - for (psU32 i = 0 ; i < row1->n ; i++ ) { - col = row1->data.U32[i]; - // - // Determine if pixel (0,0) is a peak. - // - if (col == 0) { - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] > image->data.F32[row][col+1])) { - myPeakClass = PM_PEAK_EDGE | PM_PEAK_LONE; - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); - } - } else if (col < (image->numCols - 1)) { - - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] >= image->data.F32[row][col-1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1])) { - myPeakClass = PM_PEAK_EDGE; - - if (image->data.F32[row][col] > image->data.F32[row][col-1]) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; - } - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); - } - - } else if (col == (image->numCols - 1)) { - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] >= image->data.F32[row][col-1])) { - myPeakClass = PM_PEAK_EDGE; - - if (image->data.F32[row][col] > image->data.F32[row][col-1]) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; - } - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); - } - } else { - psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); - } - - } - // - // We exit here since this is a single-row image. - // - // XXX: Why are we not getting memory leak errors? - // - // - // psFree(tmpRow); - // psFree(row1); - // - return(list); - } - - // - // Find peaks in col 0 only (single-col image). This is a special case since - // we can not test data at image->data.F32[...][col+1]) - // - if (image->numCols == 1) { - col = 0; - psVector *tmpCol = p_psGetColVectorFromImage((psImage *) image, col); - psVector *col1 = pmFindVectorPeaks(tmpCol, threshold); - for (psU32 i = 0 ; i < col1->n ; i++ ) { - row = col1->data.U32[i]; - // - // Determine if pixel (0,0) is a peak. - // - if (row == 0) { - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] > image->data.F32[row+1][col])) { - myPeakClass = PM_PEAK_EDGE | PM_PEAK_LONE; - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); - } - } else if (row < (image->numRows - 1)) { - - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] >= image->data.F32[row-1][col]) && - (image->data.F32[row][col] > image->data.F32[row+1][col])) { - myPeakClass = PM_PEAK_EDGE; - - if (image->data.F32[row][col] > image->data.F32[row-1][col]) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; - } - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); - } - - } else if (row == (image->numRows - 1)) { - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] >= image->data.F32[row-1][col])) { - myPeakClass = PM_PEAK_EDGE; - - if (image->data.F32[row][col] > image->data.F32[row-1][col]) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; - } - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); - } - } else { - psError(PS_ERR_UNKNOWN, true, "peak specified outside valid row range."); - } - - } - // - // We exit here since this is a single-column image. - // - // XXX: free tmpRow col1? - // - return(list); - } - - // - // Find peaks in row 0 only (multi-row image). + psArray *list = NULL; + + // + // Find peaks in row 0 only. // row = 0; tmpRow = p_psGetRowVectorFromImage((psImage *) image, row); psVector *row1 = pmFindVectorPeaks(tmpRow, threshold); + for (psU32 i = 0 ; i < row1->n ; i++ ) { col = row1->data.U32[i]; @@ -472,52 +326,41 @@ // if (col == 0) { - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] > image->data.F32[row][col+1]) && + if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) && (image->data.F32[row][col] > image->data.F32[row+1][col]) && (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { - myPeakClass = PM_PEAK_EDGE; - if (image->data.F32[row][col] > image->data.F32[row+1][col+1]) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; + if (image->data.F32[row][col] > threshold) { + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); } } else if (col < (image->numCols - 1)) { - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] >= image->data.F32[row][col-1]) && + if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && (image->data.F32[row][col] > image->data.F32[row][col+1]) && (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && (image->data.F32[row][col] > image->data.F32[row+1][col]) && (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { - - myPeakClass = PM_PEAK_EDGE; - if ( (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+1])) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; + if (image->data.F32[row][col] > threshold) { + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); } } else if (col == (image->numCols - 1)) { - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] >= image->data.F32[row][col-1]) && + if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && (image->data.F32[row][col] > image->data.F32[row+1][col]) && (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) { - myPeakClass = PM_PEAK_EDGE; - if ( (image->data.F32[row][col] > image->data.F32[row][col-1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col-1])) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; + if (image->data.F32[row][col] > threshold) { + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); } + } else { - psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); - } + psError(PS_ERR_UNKNOWN, true, "peak specified valid column range."); + } + } + + // + // Exit if this image has a single row. + // + if (image->numRows == 1) { + return(list); } @@ -531,30 +374,20 @@ // Step through all local peaks in this row. for (psU32 i = 0 ; i < row1->n ; i++ ) { + psPeakType myType = PM_PEAK_UNDEF; col = row1->data.U32[i]; if (col == 0) { // If col==0, then we can not read col-1 pixels - if ((image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && + if ((image->data.F32[row][col] > image->data.F32[row-1][col]) && (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && (image->data.F32[row][col] >= image->data.F32[row][col+1]) && (image->data.F32[row][col] >= image->data.F32[row+1][col]) && (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { - - myPeakClass = PM_PEAK_EDGE; - if ((image->data.F32[row][col] > image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col]) && - (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; - } - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); + myType = PM_PEAK_EDGE; + list = MyListAddPeak(list, row, col, image->data.F32[row][col], myType); } } else if (col < (image->numCols - 1)) { // This is an interior pixel - if ((image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && (image->data.F32[row][col] > image->data.F32[row-1][col]) && (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && @@ -564,22 +397,33 @@ (image->data.F32[row][col] >= image->data.F32[row+1][col]) && (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { - - if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col]) && - (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { - myPeakClass = PM_PEAK_LONE; - } else { - myPeakClass = PM_PEAK_FLAT; + if (image->data.F32[row][col] > threshold) { + if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] > image->data.F32[row][col+1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { + myType = PM_PEAK_LONE; + } + + if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) || + (image->data.F32[row][col] == image->data.F32[row-1][col]) || + (image->data.F32[row][col] == image->data.F32[row-1][col+1]) || + (image->data.F32[row][col] == image->data.F32[row][col-1]) || + (image->data.F32[row][col] == image->data.F32[row][col+1]) || + (image->data.F32[row][col] == image->data.F32[row+1][col-1]) || + (image->data.F32[row][col] == image->data.F32[row+1][col]) || + (image->data.F32[row][col] == image->data.F32[row+1][col+1])) { + myType = PM_PEAK_FLAT; + } + + list = MyListAddPeak(list, row, col, image->data.F32[row][col], myType); } - - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); } } else if (col == (image->numCols - 1)) { // If col==numCols - 1, then we can not read col+1 pixels - if ((image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && (image->data.F32[row][col] > image->data.F32[row-1][col]) && (image->data.F32[row][col] > image->data.F32[row][col-1]) && @@ -587,19 +431,11 @@ (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && (image->data.F32[row][col] >= image->data.F32[row+1][col])) { - - myPeakClass = PM_PEAK_EDGE; - if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col])) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; - } - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); + myType = PM_PEAK_EDGE; + list = MyListAddPeak(list, row, col, image->data.F32[row][col], myType); } } else { - psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); + psError(PS_ERR_UNKNOWN, true, "peak specified valid column range."); } + } } @@ -614,52 +450,29 @@ col = row1->data.U32[i]; if (col == 0) { - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && + if ( (image->data.F32[row][col] > image->data.F32[row-1][col]) && (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && (image->data.F32[row][col] > image->data.F32[row][col+1])) { - - myPeakClass = PM_PEAK_EDGE; - if (image->data.F32[row][col] > image->data.F32[row-1][col+1]) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; + if (image->data.F32[row][col] > threshold) { + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } - - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); } } else if (col < (image->numCols - 1)) { - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && (image->data.F32[row][col] > image->data.F32[row-1][col]) && (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && (image->data.F32[row][col] > image->data.F32[row][col-1]) && (image->data.F32[row][col] >= image->data.F32[row][col+1])) { - - myPeakClass = PM_PEAK_EDGE; - if ( (image->data.F32[row][col] > image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1])) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; + if (image->data.F32[row][col] > threshold) { + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } - - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); } } else if (col == (image->numCols - 1)) { - if ( (image->data.F32[row][col] > threshold) && - (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && (image->data.F32[row][col] > image->data.F32[row-1][col]) && (image->data.F32[row][col] > image->data.F32[row][col-1])) { - - myPeakClass = PM_PEAK_EDGE; - if (image->data.F32[row][col] > image->data.F32[row-1][col-1]) { - myPeakClass|= PM_PEAK_LONE; - } else { - myPeakClass|= PM_PEAK_FLAT; + if (image->data.F32[row][col] > threshold) { + list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); } - - list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); } } else { @@ -667,24 +480,34 @@ } } - return(list); } - -/****************************************************************************** -PS_REGION_CHECK(VALID, X, Y): this macro evaluates to TRUE if the coordinate -(X, Y) is within REGION, otherwise it evaluates to FALSE. - *****************************************************************************/ -#define PS_REGION_CHECK(VALID, X, Y) \ -(((X) >= (VALID)->x0) && \ - ((X) <= (VALID)->x1) && \ - ((Y) >= (VALID)->y0) && \ - ((Y) <= (VALID)->y1)) - - -/****************************************************************************** -psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psList that have +// 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); +} + + +/****************************************************************************** +psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have a peak value above the given maximum, or fall outside the valid region. - *****************************************************************************/ + +XXX: Should the sky value be used when comparing the maximum? + +XXX: warning message if valid is NULL? + +XXX: changed API to create a NEW output psArray (should change name as well) +*****************************************************************************/ psList *pmCullPeaks(psList *peaks, psF32 maxValue, @@ -692,16 +515,15 @@ { PS_PTR_CHECK_NULL(peaks, NULL); - if (valid == NULL) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: psRegion valid is NULL. Ignoring ...\n"); - } + // PS_PTR_CHECK_NULL(valid, NULL); psListElem *tmpListElem = (psListElem *) peaks->head; psS32 indexNum = 0; + // printf("pmCullPeaks(): list size is %d\n", peaks->size); while (tmpListElem != NULL) { psPeak *tmpPeak = (psPeak *) tmpListElem->data; if ((tmpPeak->counts > maxValue) || ((valid != NULL) && - (true == PS_REGION_CHECK(valid, tmpPeak->x, tmpPeak->y)))) { + (true == IsItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) { psListRemoveData(peaks, (psPtr) tmpPeak); } @@ -712,4 +534,28 @@ return(peaks); +} + +// XXX EAM: I changed this to return a new, subset array +// rather than alter the existing one +psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid) +{ + PS_PTR_CHECK_NULL(peaks, NULL); + + psArray *output = psArrayAlloc (200); + output->n = 0; + + psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n); + + for (int i = 0; i < peaks->n; i++) { + psPeak *tmpPeak = (psPeak *) peaks->data[i]; + if (tmpPeak->counts > maxValue) + continue; + if (valid != NULL) { + if (IsItInThisRegion(valid, tmpPeak->x, tmpPeak->y)) + continue; + } + psArrayAdd (output, 200, tmpPeak); + } + return(output); } @@ -726,4 +572,6 @@ psImageStats on that subImage+mask. +XXX: The subImage has width of 1+2*outerRadius. Verify with IfA. + XXX: Use static data structures for: subImage @@ -748,5 +596,5 @@ XXX: Don't use separate structs for the subimage and mask. Use the source-> members. - *****************************************************************************/ +*****************************************************************************/ psSource *pmSourceLocalSky(const psImage *image, const psPeak *peak, @@ -769,10 +617,13 @@ // 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; - psS32 SubImageStartRow = PS_MAX(SubImageCenterRow - outerRadiusS32, 0); - psS32 SubImageEndRow = PS_MIN(SubImageCenterRow + outerRadiusS32, image->numRows - 1); - psS32 SubImageStartCol = PS_MAX(SubImageCenterCol - outerRadiusS32, 0); - psS32 SubImageEndCol = PS_MIN(SubImageCenterCol + outerRadiusS32, image->numCols - 1); + + // 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. @@ -785,9 +636,12 @@ // 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); - } + // 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", @@ -805,4 +659,5 @@ return(NULL); } + # endif // @@ -830,4 +685,5 @@ // // 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++) { @@ -873,19 +729,41 @@ /****************************************************************************** -bool PS_RADIUS_CHECK(): a macro which evaluates to TRUE if the (x, y) point is -within the radius of the specified peak. - *****************************************************************************/ -#define PS_RADIUS_CHECK(PEAK, RADIUS, X, Y) \ -(PS_SQR((RADIUS)) >= ((psF32) (PS_SQR((X) - (PEAK)->x) + PS_SQR((Y) - (PEAK)->y)))) - - -/****************************************************************************** -bool PS_RADIUS_CHECK22(): a macro which evaluates to TRUE if the (x, y) point -is within the radius of the specified point. - *****************************************************************************/ -#define PS_RADIUS_CHECK2(X_CENTER, Y_CENTER, RADIUS, X, Y) \ -(PS_SQR((RADIUS)) >= ((psF32) (PS_SQR((X) - (X_CENTER)) + PS_SQR((Y) - (Y_CENTER))))) - - +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. +*****************************************************************************/ +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(): private function which simply determines if the (x, y) +point is within the radius of the specified peak. + +XXX: macro this for performance. +XXX: this is rather inefficient - at least compute and compare against radius^2 +*****************************************************************************/ +bool CheckRadius2(psF32 xCenter, + psF32 yCenter, + psF32 radius, + psF32 x, + psF32 y) +{ + /// XXX EAM should compare with hypot (x,y) for speed + if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { + return(true); + } + + return(false); +} /****************************************************************************** @@ -899,6 +777,8 @@ psSource->pixels +XXX: The peak calculations are done in image coords, not subImage coords. + XXX: mask values? - *****************************************************************************/ +*****************************************************************************/ psSource *pmSourceMoments(psSource *source, psF32 radius) @@ -909,9 +789,10 @@ 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) { - psError(PS_ERR_UNKNOWN, true, "Undefined source->psMoments is NULL"); - psFree(source); - return(NULL); + source->moments = pmMomentsAlloc(); } else { sky = source->moments->Sky; @@ -932,8 +813,14 @@ psF32 Y2 = 0.0; psF32 XY = 0.0; - // + psF32 x = 0; + psF32 y = 0; + // + // XXX why do I get different results for these two methods of finding Sx? + // XXX Sx, Sy would be better measured if we clip pixels close to sky + // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? // We loop through all pixels in this subimage (source->pixels), and for each // pixel that is not masked, AND within the radius of the peak pixel, we - // proceed with the moments calculation. + // proceed with the moments calculation. need to do two loops for a + // numerically stable result. first loop: get the sums. // for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -942,8 +829,8 @@ psS32 imgColCoord = col + source->pixels->col0; psS32 imgRowCoord = row + source->pixels->row0; - if (PS_RADIUS_CHECK(source->peak, - radius, - imgColCoord, - imgRowCoord)) { + if (CheckRadius(source->peak, + radius, + imgColCoord, + imgRowCoord)) { psF32 xDiff = (psF32) (imgColCoord - source->peak->x); psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); @@ -953,7 +840,8 @@ X1+= xDiff * pDiff; Y1+= yDiff * pDiff; + XY+= xDiff * yDiff * pDiff; + X2+= PS_SQR(xDiff) * pDiff; Y2+= PS_SQR(yDiff) * pDiff; - XY+= xDiff * yDiff * pDiff; if (source->pixels->data.F32[row][col] > peakPixel) { @@ -971,13 +859,84 @@ // 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)); + x = X1/Sum; + y = Y1/Sum; + source->moments->x = x + ((psF32) source->peak->x); + source->moments->y = y + ((psF32) source->peak->y); + source->moments->Sxy = XY/Sum; + source->moments->Sum = Sum; source->moments->Peak = peakPixel; source->moments->nPixels = numPixels; + // XXX EAM : these values can be negative, so we need to limit the range + source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); + source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); return(source); + + // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it + # if (0) + // + // second loop: get the difference sums + // + X2 = Y2 = 0; + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) { + psS32 imgColCoord = col + source->pixels->col0; + psS32 imgRowCoord = row + source->pixels->row0; + if (CheckRadius(source->peak, + radius, + imgColCoord, + imgRowCoord)) { + psF32 xDiff = (psF32) (imgColCoord - source->peak->x); + psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); + psF32 pDiff = source->pixels->data.F32[row][col] - sky; + + Sum+= pDiff; + X2+= PS_SQR(xDiff - x) * pDiff; + Y2+= PS_SQR(yDiff - y) * pDiff; + } + } + } + } + + // + // second moment X = sqrt (X2/Sum) + // + source->moments->Sx = (X2/Sum); + source->moments->Sy = (Y2/Sum); + return(source); + # endif +} + +// XXX EAM : I used +int pmComparePeakAscend (const void **a, const void **b) +{ + psPeak *A = *(psPeak **)a; + psPeak *B = *(psPeak **)b; + + psF32 diff; + + diff = A->counts - B->counts; + if (diff < FLT_EPSILON) + return (-1); + if (diff > FLT_EPSILON) + return (+1); + return (0); +} + +int pmComparePeakDescend (const void **a, const void **b) +{ + psPeak *A = *(psPeak **)a; + psPeak *B = *(psPeak **)b; + + psF32 diff; + + diff = A->counts - B->counts; + if (diff < FLT_EPSILON) + return (+1); + if (diff > FLT_EPSILON) + return (-1); + return (0); } @@ -990,80 +949,240 @@ XXX: The sigX and sigY stuff in the SDRS is unclear. - *****************************************************************************/ -#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); + +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_PTR_CHECK_NULL(sources, 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); + psArray *peaks = NULL; + psF32 clumpX = 0.0; + psF32 clumpDX = 0.0; + psF32 clumpY = 0.0; + psF32 clumpDY = 0.0; + + // find the sigmaX, sigmaY clump + { + psStats *stats = NULL; + psImage *splane = NULL; + int binX, binY; + + // construct a sigma-plane image + splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32); + + // place the sources in the sigma-plane image (ignore 0,0 values?) + for (psS32 i = 0 ; i < sources->n ; i++) + { + psSource *tmpSrc = (psSource *) sources->data[i]; + PS_PTR_CHECK_NULL(tmpSrc, false); // just skip this one? + PS_PTR_CHECK_NULL(tmpSrc->moments, false); // just skip this one? + + // Sx,Sy are limited at 0. a peak at 0,0 is artificial + if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { + continue; + } + + // for the moment, force splane dimensions to be 10x10 image pix + binX = tmpSrc->moments->Sx/SCALE; + if (binX < 0) + continue; + if (binX >= splane->numCols) + continue; + + binY = tmpSrc->moments->Sy/SCALE; + if (binY < 0) + continue; + if (binY >= splane->numRows) + continue; + + splane->data.F32[binY][binX] += 1.0; + } + + // find the peak in this image + stats = psStatsAlloc (PS_STAT_MAX); + stats = psImageStats (stats, splane, NULL, 0); + peaks = pmFindImagePeaks (splane, stats[0].max / 2); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); + + } + + // measure statistics on Sx, Sy if Sx, Sy within range of clump + { + psPeak *clump; + psF32 minSx, maxSx; + psF32 minSy, maxSy; + psVector *tmpSx = NULL; + psVector *tmpSy = NULL; + psStats *stats = NULL; + + // XXX EAM : this lets us takes the single highest peak + psArraySort (peaks, pmComparePeakDescend); + clump = peaks->data[0]; + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d\n", clump->x, clump->y); + + // define section window for clump + minSx = clump->x * SCALE - 0.2; + maxSx = clump->x * SCALE + 0.2; + minSy = clump->y * SCALE - 0.2; + maxSy = clump->y * SCALE + 0.2; + + tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); + tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); + tmpSx->n = 0; + tmpSy->n = 0; + + // XXX clip sources based on flux? + // create vectors with Sx, Sy values in window + for (psS32 i = 0 ; i < sources->n ; i++) + { + psSource *tmpSrc = (psSource *) sources->data[i]; + + if (tmpSrc->moments->Sx < minSx) + continue; + if (tmpSrc->moments->Sx > maxSx) + continue; + if (tmpSrc->moments->Sy < minSy) + continue; + if (tmpSrc->moments->Sy > maxSy) + continue; + tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; + tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; + tmpSx->n++; + tmpSy->n++; + if (tmpSx->n == tmpSx->nalloc) { + psVectorRealloc (tmpSx, tmpSx->nalloc + 100); + psVectorRealloc (tmpSy, tmpSy->nalloc + 100); + } + } + + // measures stats of Sx, Sy + stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + + stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); + clumpX = stats->clippedMean; + clumpDX = stats->clippedStdev; + + stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); + clumpY = stats->clippedMean; + clumpDY = stats->clippedStdev; + + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", clumpX, clumpY); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", clumpDX, clumpDY); + // these values should be pushed on the metadata somewhere + } + + int Nsat = 0; + int Ngal = 0; + int Nfaint = 0; + int Nstar = 0; + int Npsf = 0; + int Ncr = 0; + psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); + starsn->n = 0; + + // XXX allow clump size to be scaled relative to sigmas? + // make rough IDs based on clumpX,Y,DX,DY + for (psS32 i = 0 ; i < sources->n ; i++) { + + psSource *tmpSrc = (psSource *) sources->data[i]; + tmpSrc->peak->class = 0; - psF32 sigX = 0.0; - psF32 sigY = 0.0; - // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE. - psF32 clumpX = 0.0; - psF32 clumpDX = 0.0; - psF32 clumpY = 0.0; - psF32 clumpDY = 0.0; - + psF32 sigX = tmpSrc->moments->Sx; + psF32 sigY = tmpSrc->moments->Sy; + + // check return status value (do these exist?) + bool status; + psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE"); + psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN"); + psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE"); + + psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); + psF32 FAINT_SN_LIM = psMetadataLookupF32 (&status, metadata, "FAINT_SN_LIM"); + + // saturated object (star or single pixel not distinguished) if (tmpSrc->moments->Peak > SATURATE) { - tmpSrc->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; - } - // XXX: The SDRS is not real clear on how to calculate sigX, sigY. - if ((fabs(sigX - clumpX) < clumpDX) && - (fabs(sigY - clumpY) < clumpDY)) { - tmpSrc->peak->class|= PS_SOURCE_PSFSTAR; - } - - if ((sigX < (clumpX - clumpDX)) && - (sigY < (clumpY - clumpDY))) - tmpSrc->peak->class|= PS_SOURCE_DEFECT; - } - - if ((sigX > (clumpX + clumpDX)) && - (sigY > (clumpY + clumpDY))) { - tmpSrc->peak->class|= PS_SOURCE_GALAXY; - } - - if (tmpSrc->peak->class == 0) { - tmpSrc->peak->class|= PS_SOURCE_OTHER; - } - } + tmpSrc->type |= PS_SOURCE_SATURATED; + Nsat ++; + continue; + } + + // too small to be stellar + if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) { + tmpSrc->type |= PS_SOURCE_DEFECT; + Ncr ++; + continue; + } + + // possible galaxy + if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) { + tmpSrc->type |= PS_SOURCE_GALAXY; + Ngal ++; + continue; + } + + // the rest are probable stellar objects + psF32 S = tmpSrc->moments->Sum; + psF32 A = PS_PI * sigX * sigY; + psF32 B = tmpSrc->moments->Sky; + psF32 RT = PS_SQRT_F32(S + (A * B) + (A * PS_SQR(RDNOISE) / PS_SQRT_F32(GAIN))); + psF32 SN = (S * PS_SQRT_F32(GAIN) / RT); + + starsn->data.F32[starsn->n] = SN; + starsn->n ++; + Nstar ++; + + // faint star + if (SN < FAINT_SN_LIM) { + tmpSrc->type |= PS_SOURCE_FAINTSTAR; + Nfaint ++; + continue; + } + + // PSF star + if (SN > PSF_SN_LIM) { + tmpSrc->type |= PS_SOURCE_PSFSTAR; + Npsf ++; + continue; + } + + // random type of star + tmpSrc->type |= PS_SOURCE_OTHER; + } + + { + psStats *stats = NULL; + stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); + stats = psVectorStats (stats, starsn, NULL, NULL, 0); + psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); + } + + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nfaint: %3d\n", Nfaint); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal: %3d\n", Ngal); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); return(rc); } - - /****************************************************************************** pmSourceSetPixelCircle(source, image, radius) -XXX: Checking source->moments for NULL. Circle must be centered on the - centroid, not peak (from IfA 2005-04-06). - *****************************************************************************/ +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, const psImage *image, @@ -1073,6 +1192,6 @@ 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_PTR_CHECK_NULL(source->moments, false); + // PS_PTR_CHECK_NULL(source->peak, false); PS_FLOAT_COMPARE(0.0, radius, false); @@ -1086,15 +1205,19 @@ 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) { + // XXX EAM : for the circle to stay on the image + psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); + psS32 SubImageEndRow = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32); + psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); + psS32 SubImageEndCol = PS_MIN (image->numCols - 1, SubImageCenterCol + radiusS32); + + // XXX EAM : this should not be needed: we can never hit this error + # if (1) + + if (SubImageStartRow < 0) { + psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", + SubImageStartRow); + return(false); + } + if (SubImageEndRow >= image->numRows) { psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", SubImageEndRow); @@ -1106,13 +1229,16 @@ return(false); } - if (SubImageEndCol+1 >= image->numCols) { + if (SubImageEndCol >= image->numCols) { psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", SubImageEndCol); return(false); } + # endif // XXX: Must recycle image. + // XXX EAM: this message reflects a programming error we know about. + // i am setting it to a trace message which we can take out if (source->pixels != NULL) { - psLogMsg(__func__, PS_LOG_WARN, + psTrace (".psModule.pmObjects.pmSourceSetPixelCircle", 4, "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL. Freeing and reallocating.\n"); psFree(source->pixels); @@ -1121,6 +1247,6 @@ SubImageStartCol, SubImageStartRow, - SubImageEndCol+1, - SubImageEndRow+1); + SubImageEndCol, + SubImageEndRow); // XXX: Must recycle image. @@ -1128,22 +1254,22 @@ psFree(source->mask); } - source->mask = psImageAlloc(1 + 2 * radiusS32, - 1 + 2 * radiusS32, + source->mask = psImageAlloc(source->pixels->numCols, + source->pixels->numRows, PS_TYPE_F32); // // Loop through the subimage mask, initialize mask to 0 or 1. - // + // XXX EAM: valid pixels should have 0, not 1 for (psS32 row = 0 ; row < source->mask->numRows; row++) { for (psS32 col = 0 ; col < source->mask->numCols; col++) { - if (PS_RADIUS_CHECK2((psF32) radiusS32, - (psF32) radiusS32, - radius, - (psF32) col, - (psF32) row)) { + if (CheckRadius2((psF32) radiusS32, + (psF32) radiusS32, + radius, + (psF32) col, + (psF32) row)) { + source->mask->data.U8[row][col] = 0; + } else { source->mask->data.U8[row][col] = 1; - } else { - source->mask->data.U8[row][col] = 0; } } @@ -1151,5 +1277,4 @@ return(true); } - /****************************************************************************** @@ -1168,5 +1293,5 @@ image, not subImage coords. Therefore, the calls to the model evaluation functions will be in image, not subImage coords. Remember this. - *****************************************************************************/ +*****************************************************************************/ bool pmSourceModelGuess(psSource *source, const psImage *image, @@ -1184,48 +1309,35 @@ source->models = pmModelAlloc(model); + psVector *params = source->models->params; + switch (model) { case PS_MODEL_GAUSS: - 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; + 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: - 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; + 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_TWIST_GAUSS: - 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; - 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; + 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->models->params[7] = B2; @@ -1233,30 +1345,46 @@ 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: - 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; + 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->models->params[7] = Nexp; + //params->data.F32[7] = Nexp; return(true); case PS_MODEL_SERSIC_CORE: - 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; + 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? - //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; + // 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); @@ -1268,5 +1396,5 @@ /****************************************************************************** -evalModel(src, col, row): a private function which evaluates the +evalModel(source, level, row): a private function which evaluates the source->model function at the specified coords. The coords are subImage, not image coords. @@ -1274,4 +1402,6 @@ NOTE: The coords are in subImage source->pixel coords, not image coords. +XXX: reverse order of row,col args? + XXX: rename all coords in this file such that their name defines whether the coords is in subImage or image space. @@ -1285,8 +1415,8 @@ XXX: For a while, the first psVectorAlloc() was generating a seg fault during testing. Try to reproduce that and debug. - *****************************************************************************/ -psF32 evalModel(psSource *src, - psU32 col, - psU32 row) +*****************************************************************************/ +static psF32 evalModel(psSource *src, + psU32 row, + psU32 col) { PS_PTR_CHECK_NULL(src, false); @@ -1296,8 +1426,7 @@ // XXX: The following step will not be necessary if the models->params // member is a psVector. Suggest to IfA. - psVector *params = psVectorAlloc(7, PS_TYPE_F32); - for (psS32 i = 0 ; i < src->models->Nparams ; i++) { - params->data.F32[i] = src->models->params[i]; - } + + // XXX EAM: done: models->params is now a vector + psVector *params = src->models->params; // @@ -1332,6 +1461,4 @@ return(NAN); } - - psFree(params); psFree(x); return(tmpF); @@ -1339,18 +1466,20 @@ /****************************************************************************** -findValue(source, level, col, row, dir): a private function which determines +findValue(source, level, row, col, dir): a private function which determines the column coordinate of the model function which has the value "level". If dir equals 0, then you loop leftwards from the peak pixel, otherwise, rightwards. +XXX: reverse order of row,col args? + XXX: Input row/col are in image coords. XXX: The result is returned in image coords. - *****************************************************************************/ -psF32 findValue(psSource *source, - psF32 level, - psU32 col, - psU32 row, - psU32 dir) +*****************************************************************************/ +static psF32 findValue(psSource *source, + psF32 level, + psU32 row, + psU32 col, + psU32 dir) { // @@ -1370,5 +1499,5 @@ } - psF32 oldValue = evalModel(source, subCol, subRow); + psF32 oldValue = evalModel(source, subRow, subCol); if (oldValue == level) { return(((psF32) (subCol + source->pixels->col0))); @@ -1391,5 +1520,5 @@ while (subCol != lastColumn) { - psF32 newValue = evalModel(source, subCol, subRow); + psF32 newValue = evalModel(source, subRow, subCol); if (oldValue == level) { return((psF32) (subCol + source->pixels->col0)); @@ -1422,5 +1551,5 @@ XXX: What is momde? XXX: The top, bottom of the contour is not correctly determined. - *****************************************************************************/ +*****************************************************************************/ psArray *pmSourceContour(psSource *source, const psImage *image, @@ -1451,5 +1580,5 @@ // Starting at peak pixel, search leftwards for the column intercept. - psF32 leftIntercept = findValue(source, level, col, row, 0); + psF32 leftIntercept = findValue(source, level, row, col, 0); if (isnan(leftIntercept)) { psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); @@ -1462,5 +1591,5 @@ // Starting at peak pixel, search rightwards for the column intercept. - psF32 rightIntercept = findValue(source, level, col, row, 1); + psF32 rightIntercept = findValue(source, level, row, col, 1); if (isnan(rightIntercept)) { psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); @@ -1485,5 +1614,5 @@ // Starting at peak pixel, search leftwards for the column intercept. - psF32 leftIntercept = findValue(source, level, col, row, 0); + psF32 leftIntercept = findValue(source, level, row, col, 0); if (isnan(leftIntercept)) { psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); @@ -1495,5 +1624,5 @@ // Starting at peak pixel, search rightwards for the column intercept. - psF32 rightIntercept = findValue(source, level, col, row, 1); + psF32 rightIntercept = findValue(source, level, row, col, 1); if (isnan(rightIntercept)) { psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); @@ -1524,7 +1653,7 @@ 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 +// XXX EAM : these are better starting values, but should be available from metadata? +#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 20 +#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 /****************************************************************************** pmSourceFitModel(source, image): must create the appropiate arguments to the @@ -1533,5 +1662,5 @@ XXX: should there be a mask value? XXX: Probably should remove the "image" argument. - *****************************************************************************/ +*****************************************************************************/ bool pmSourceFitModel(psSource *source, const psImage *image) @@ -1545,7 +1674,11 @@ PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); psBool rc; + + // find the number of valid pixels + // XXX EAM : this loop and the loop below could just be one pass + // using the psArrayAdd and psVectorExtend functions psS32 count = 0; - for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { - for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { + 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++; @@ -1553,16 +1686,23 @@ } } + + // 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++) { + 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: - coord->data.F32[0] = (psF32) (i + source->pixels->row0); - coord->data.F32[1] = (psF32) (j + source->pixels->col0); + // XXX EAM: coord order is (x,y) == (col,row) + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); x->data[tmpCnt] = (psPtr *) coord; y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; + + // XXX EAM : this is approximate: need to apply the gain and rdnoise + yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j])); tmpCnt++; } @@ -1573,31 +1713,24 @@ PM_SOURCE_FIT_MODEL_TOLERANCE); - psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32); + psVector *params = source->models->params; switch (source->models->type) { case PS_MODEL_GAUSS: - - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec); + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D); break; case PS_MODEL_PGAUSS: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec); + 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, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec); + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D); break; case PS_MODEL_WAUSS: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec); + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D); break; case PS_MODEL_SERSIC: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec); + 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, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec); + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore); break; default: @@ -1605,9 +1738,13 @@ rc = false; } + // XXX EAM: we need to do something (give an error?) if rc is false + // XXX EAM: save the resulting chisq, nDOF, nIter + source->models->chisq = myMin->value; + source->models->nDOF = y->n - params->n; + source->models->nIter = myMin->iter; psFree(x); psFree(y); psFree(myMin); - psFree(params); return(rc); } @@ -1626,12 +1763,9 @@ PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); - psVector *params = psVectorAlloc(src->models->Nparams, PS_TYPE_F32); psVector *x = psVectorAlloc(2, PS_TYPE_F32); - for (psS32 i = 0 ; i < src->models->Nparams ; i++) { - params->data.F32[i] = src->models->params[i]; - } - - for (psS32 i = 0 ; i < src->pixels->numRows ; i++) { - for (psS32 j = 0 ; j < src->pixels->numCols ; j++) { + psVector *params = src->models->params; + + for (psS32 i = 0; i < src->pixels->numRows; i++) { + for (psS32 j = 0; j < src->pixels->numCols; j++) { psF32 pixelValue; // XXX: Should you be adding the pixels for the entire subImage, @@ -1667,5 +1801,4 @@ psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); psFree(x); - psFree(params); return(false); } @@ -1681,5 +1814,4 @@ } psFree(x); - psFree(params); return(true); } @@ -1716,10 +1848,14 @@ -/****************************************************************************** -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. - +/** + 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; @@ -1729,39 +1865,19 @@ 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. - -XXX: I added the following. Must conform with IfA. If deriv==NULL, then -we simply don't perform the derivative calculations. - -XXX: It is not clear whether the subImage coords, or the image coords should -be used when calling these functions. I don't think it really matters, as -long as we are consistent. - *****************************************************************************/ -psF32 pmMinLM_Gauss2D(psVector *deriv, +*****************************************************************************/ +psF64 pmMinLM_Gauss2D(psVector *deriv, psVector *params, psVector *x) { - 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]; + 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 f = params->data.F32[1]*r + params->data.F32[0]; + 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) { - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); - - psF32 q = params->data.F32[1]*r; deriv->data.F32[0] = +1.0; deriv->data.F32[1] = +r; @@ -1772,49 +1888,6 @@ deriv->data.F32[6] = -q*X*Y; } - return(f); } - -/****************************************************************************** -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]; - } - - psVector *tmpVec2 = (psVector *) x->data[i]; - tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow, - params, - tmpVec2); - - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - deriv->data.F32[i][j] = tmpRow->data.F32[j]; - } - } - - psFree(tmpRow); - return(tmpVec); -} - /****************************************************************************** @@ -1826,88 +1899,31 @@ params->data.F32[5] = sqrt(2) / SigmaY; params->data.F32[6] = Sxy; - *****************************************************************************/ -psF32 pmMinLM_PsuedoGauss2D(psVector *deriv, +*****************************************************************************/ +psF64 pmMinLM_PsuedoGauss2D(psVector *deriv, psVector *params, psVector *x) { - 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]; + 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]; + 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) { - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); - - // // 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[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); } - -/****************************************************************************** -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); -} - - - /****************************************************************************** @@ -1921,16 +1937,9 @@ params->data.F32[7] = B2; params->data.F32[8] = B3; - *****************************************************************************/ -psF32 pmMinLM_Wauss2D(psVector *deriv, +*****************************************************************************/ +psF64 pmMinLM_Wauss2D(psVector *deriv, psVector *params, psVector *x) { - 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]; @@ -1943,10 +1952,5 @@ if (deriv != NULL) { - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 9, NAN); - // // 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; @@ -1959,56 +1963,8 @@ 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); } - -/****************************************************************************** -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? @@ -2027,16 +1983,9 @@ params->data.F32[9] = SxyOuter; params->data.F32[10] = N; - *****************************************************************************/ -psF32 pmMinLM_TwistGauss2D(psVector *deriv, +*****************************************************************************/ +psF64 pmMinLM_TwistGauss2D(psVector *deriv, psVector *params, psVector *x) { - 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]; @@ -2049,11 +1998,7 @@ 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) { - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 11, NAN); - 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)); @@ -2063,16 +2008,12 @@ 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; @@ -2087,46 +2028,4 @@ return(f); } - -/****************************************************************************** -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); -} - - /****************************************************************************** @@ -2140,62 +2039,11 @@ params->data.F32[6] = Sxy; params->data.F32[7] = Nexp; - *****************************************************************************/ -psF32 pmMinLM_Sersic(psVector *deriv, +*****************************************************************************/ +psF64 pmMinLM_Sersic(psVector *deriv, psVector *params, psVector *x) { - 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); - - if (deriv != NULL) { - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 8, 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); } @@ -2214,72 +2062,10 @@ params->data.F32[10] = SxyOuter; params->data.F32[11] = Nexp; - *****************************************************************************/ -psF32 pmMinLM_SersicCore(psVector *deriv, +*****************************************************************************/ +psF64 pmMinLM_SersicCore(psVector *deriv, psVector *params, psVector *x) { - 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); - - if (deriv != NULL) { - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 12, 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); -} - - - -/****************************************************************************** - *****************************************************************************/ -psF32 pmMinLM_PsuedoSersic(psVector *deriv, - psVector *params, - psVector *x) -{ - return(0.0); -} Index: /trunk/psModules/src/pmObjects.h =================================================================== --- /trunk/psModules/src/pmObjects.h (revision 3876) +++ /trunk/psModules/src/pmObjects.h (revision 3877) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-04-14 02:42:41 $ + * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-05-10 23:48:19 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -24,8 +24,8 @@ typedef enum { - PM_PEAK_LONE = 0x00000001, ///< Isolated peak. - PM_PEAK_EDGE = 0x00000002, ///< Peak on edge. - PM_PEAK_FLAT = 0x00000004, ///< Peak has equal-value neighbors. - PM_PEAK_UNDEF = 0x00000008 ///< Undefined. + PM_PEAK_LONE, ///< Isolated peak. + PM_PEAK_EDGE, ///< Peak on edge. + PM_PEAK_FLAT, ///< Peak has equal-value neighbors. + PM_PEAK_UNDEF ///< Undefined. } psPeakType; @@ -54,5 +54,5 @@ typedef enum { - PS_MODEL_GAUSS, ///< Regular 2-D Gaussian + PS_MODEL_GAUSS, ///< Regular 2-D Gaussian PS_MODEL_PGAUSS, ///< Psuedo 2-D Gaussian PS_MODEL_TWIST_GAUSS, ///< 2-D Twisted Gaussian @@ -67,8 +67,9 @@ { psModelType type; ///< Model to be used. - psS32 Nparams; ///< Number of parameters. - psF32 *params; ///< Paramater values. - psF32 *dparams; ///< Parameter errors. + psVector *params; ///< Paramater values. + psVector *dparams; ///< Parameter errors. psF32 chisq; ///< Fit chi-squared. + psS32 nDOF; ///< number of degrees of freedom + psS32 nIter; ///< number of iterations to reach min } psModel; @@ -121,7 +122,7 @@ value) of all peaks. *****************************************************************************/ -psList *pmFindImagePeaks(const psImage *image, ///< The input image where peaks will be found (psF32) - psF32 threshold ///< Threshold above which to find a peak - ); +psArray *pmFindImagePeaks(const psImage *image, ///< The input image where peaks will be found (psF32) + psF32 threshold ///< Threshold above which to find a peak + ); /****************************************************************************** @@ -209,6 +210,7 @@ /****************************************************************************** XXX: Why only *x argument? - *****************************************************************************/ -psF32 pmMinLM_Gauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives +XXX EAM: psMinimizeLMChi2Func returns psF64, not psF32 + *****************************************************************************/ +psF64 pmMinLM_Gauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives psVector *params, ///< A psVector which holds the parameters of this function psVector *x ///< A psVector which holds the row/col coordinate @@ -218,5 +220,5 @@ XXX: Why only *x argument? *****************************************************************************/ -psF32 pmMinLM_PsuedoGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives +psF64 pmMinLM_PsuedoGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives psVector *params, ///< A psVector which holds the parameters of this function psVector *x ///< A psVector which holds the row/col coordinate @@ -225,5 +227,5 @@ /****************************************************************************** *****************************************************************************/ -psF32 pmMinLM_Wauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives +psF64 pmMinLM_Wauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives psVector *params, ///< A psVector which holds the parameters of this function psVector *x ///< A psVector which holds the row/col coordinate @@ -232,5 +234,5 @@ /****************************************************************************** *****************************************************************************/ -psF32 pmMinLM_TwistGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives +psF64 pmMinLM_TwistGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives psVector *params, ///< A psVector which holds the parameters of this function psVector *x ///< A psVector which holds the row/col coordinate @@ -239,5 +241,5 @@ /****************************************************************************** *****************************************************************************/ -psF32 pmMinLM_Sersic(psVector *deriv, ///< A possibly-NULL structure for the output derivatives +psF64 pmMinLM_Sersic(psVector *deriv, ///< A possibly-NULL structure for the output derivatives psVector *params, ///< A psVector which holds the parameters of this function psVector *x ///< A psVector which holds the row/col coordinate @@ -246,5 +248,5 @@ /****************************************************************************** *****************************************************************************/ -psF32 pmMinLM_SersicCore(psVector *deriv, ///< A possibly-NULL structure for the output derivatives +psF64 pmMinLM_SersicCore(psVector *deriv, ///< A possibly-NULL structure for the output derivatives psVector *params, ///< A psVector which holds the parameters of this function psVector *x ///< A psVector which holds the row/col coordinate @@ -253,5 +255,5 @@ /****************************************************************************** *****************************************************************************/ -psF32 pmMinLM_PsuedoSersic(psVector *deriv, ///< A possibly-NULL structure for the output derivatives +psF64 pmMinLM_PsuedoSersic(psVector *deriv, ///< A possibly-NULL structure for the output derivatives psVector *params, ///< A psVector which holds the parameters of this function psVector *x ///< A psVector which holds the row/col coordinate