Index: trunk/psphot/src/pmObjects_EAM.c =================================================================== --- trunk/psphot/src/pmObjects_EAM.c (revision 4949) +++ trunk/psphot/src/pmObjects_EAM.c (revision 4950) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-06 06:43:59 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-06 08:05:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -15,4 +15,5 @@ #include #include +#include #include "pslib.h" #include "psConstants.h" @@ -66,4 +67,5 @@ and initialize the type member. Initialize the params to 0.0. XXX EAM: changing params and dparams to psVector +XXX EAM: simplifying code with psModelParameterCount *****************************************************************************/ pmModel *pmModelAlloc(pmModelType type) @@ -73,33 +75,13 @@ tmp->type = type; tmp->chisq = 0.0; - switch (type) { - case PS_MODEL_GAUSS: - tmp->params = psVectorAlloc(7, PS_TYPE_F32); - tmp->dparams = psVectorAlloc(7, PS_TYPE_F32); - break; - case PS_MODEL_PGAUSS: - tmp->params = psVectorAlloc(7, PS_TYPE_F32); - tmp->dparams = psVectorAlloc(7, PS_TYPE_F32); - break; - case PS_MODEL_TWIST_GAUSS: - tmp->params = psVectorAlloc(11, PS_TYPE_F32); - tmp->dparams = psVectorAlloc(11, PS_TYPE_F32); - break; - case PS_MODEL_WAUSS: - tmp->params = psVectorAlloc(9, PS_TYPE_F32); - tmp->dparams = psVectorAlloc(9, PS_TYPE_F32); - break; - case PS_MODEL_SERSIC: - tmp->params = psVectorAlloc(8, PS_TYPE_F32); - tmp->dparams = psVectorAlloc(8, PS_TYPE_F32); - break; - case PS_MODEL_SERSIC_CORE: - tmp->params = psVectorAlloc(12, PS_TYPE_F32); - tmp->dparams = psVectorAlloc(12, PS_TYPE_F32); - break; - default: - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + tmp->nIter = 0; + psS32 Nparams = psModelParameterCount (type); + if (Nparams == 0) { + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); return(NULL); } + + tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32); for (psS32 i = 0; i < tmp->params->n; i++) { @@ -113,12 +95,12 @@ /****************************************************************************** -XXX: We don't free pixels and mask since that caused a memory error. -We might need to increase the reference counter and decrease it here. +XXX EAM : we can now free these pixels - memory ref is incremented now *****************************************************************************/ static void sourceFree(pmSource *tmp) { psFree(tmp->peak); - // psFree(tmp->pixels); - // psFree(tmp->mask); + psFree(tmp->pixels); + psFree(tmp->noise); + psFree(tmp->mask); psFree(tmp->moments); psFree(tmp->modelPSF); @@ -135,4 +117,5 @@ tmp->peak = NULL; tmp->pixels = NULL; + tmp->noise = NULL; tmp->mask = NULL; tmp->moments = NULL; @@ -248,5 +231,6 @@ psVector containing the specified row of data from the psImage. -XXX: Is there a better way to do this? +XXX: Is there a better way to do this? +XXX EAM: does this really need to alloc a new vector??? *****************************************************************************/ static psVector *getRowVectorFromImage(psImage *image, @@ -283,4 +267,6 @@ } psArrayAdd(list, 100, tmpPeak); + psFree (tmpPeak); + // XXX EAM : is this free appropriate? (does psArrayAdd increment memory counter?) return(list); @@ -364,4 +350,6 @@ } } + psFree (tmpRow); + psFree (row1); // @@ -446,4 +434,6 @@ } + psFree (tmpRow); + psFree (row1); } @@ -487,4 +477,6 @@ } } + psFree (tmpRow); + psFree (row1); return(list); } @@ -496,9 +488,12 @@ { + // XXX EAM this function should check the status of valid + if (valid == NULL) return (true); + if ((x >= valid->x0) && - (x <= valid->x1) && - (y >= valid->y0) && - (y <= valid->y1)) { - return(true); + (x <= valid->x1) && + (y >= valid->y0) && + (y <= valid->y1)) { + return(true); } @@ -629,11 +624,12 @@ // 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); + // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) + psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - outerRadiusS32); + psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + outerRadiusS32 + 1); + psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - outerRadiusS32); + psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + outerRadiusS32 + 1); // AnulusWidth == number of pixels width in the annulus. We add one since // the pixels at the inner AND outher radius are included. - psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32); + // XXX EAM : not used : psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32); // Example: assume an outer/inner radius of 20/10. Then the subimage // should have width/length of 40. An 18-by-18 interior region will @@ -651,5 +647,5 @@ return(NULL); } - if (SubImageEndRow >= image->numRows) { + if (SubImageEndRow > image->numRows) { psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", SubImageEndRow); @@ -661,5 +657,5 @@ return(NULL); } - if (SubImageEndCol >= image->numCols) { + if (SubImageEndCol > image->numCols) { psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", SubImageEndCol); @@ -683,15 +679,20 @@ SubImageEndRow); psImage *subImage = psImageSubset((psImage *) image, tmpRegion); + // XXX EAM : can merge these: psImageSubset (image, psRegionSet ()); psImage *subImageMask = psImageAlloc(subImage->numCols, subImage->numRows, PS_TYPE_U8); - - // - // Loop through the subimage mask, initialize mask to 0. + // XXX EAM : for consistency, mask needs col0,row0 set to match image + subImageMask->col0 = subImage->col0; + subImageMask->row0 = subImage->row0; + + // + // Loop through the subimage mask, initialize mask to 1 (invalid pixel) + // XXX EAM : use PSPHOT_MASK_INVALID? // for (psS32 row = 0 ; row < subImageMask->numRows; row++) { for (psS32 col = 0 ; col < subImageMask->numCols; col++) { - subImageMask->data.U8[row][col] = 0; + subImageMask->data.U8[row][col] = 1; } } @@ -699,9 +700,13 @@ // // Loop through the subimage, mask off pixels in the inner square. - // XXX this uses a static mask value of 1 - // - for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) { - for (psS32 col = AnulusWidth; col <= (subImageMask->numCols - AnulusWidth) - 1; col++) { - subImageMask->data.U8[row][col] = 1; + // XXX this uses a static mask value of 0 + // XXX EAM : this was wrong: it masked the wrong pixels at the edge of the image + psS32 StartRow = PS_MAX (0, SubImageCenterRow - subImageMask->row0 - innerRadiusS32); + psS32 EndRow = PS_MIN (subImageMask->numRows, SubImageCenterRow - subImageMask->row0 + innerRadiusS32 + 1); + psS32 StartCol = PS_MAX (0, SubImageCenterCol - subImageMask->col0 - innerRadiusS32); + psS32 EndCol = PS_MIN (subImageMask->numCols, SubImageCenterCol - subImageMask->col0 + innerRadiusS32 + 1); + for (psS32 row = StartRow; row < EndRow; row++) { + for (psS32 col = StartCol; col < EndCol; col++) { + subImageMask->data.U8[row][col] = 0; } } @@ -795,6 +800,6 @@ XXX: mask values? *****************************************************************************/ -pmSource *pmSourceMoments(pmSource *source, - psF32 radius) +bool pmSourceMoments(psSource *source, + psF32 radius) { PS_ASSERT_PTR_NON_NULL(source, NULL); @@ -840,31 +845,40 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { for (psS32 col = 0; col < source->pixels->numCols ; col++) { - if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) { - psS32 imgColCoord = col + source->pixels->col0; - psS32 imgRowCoord = row + source->pixels->row0; - if (checkRadius(source->peak, - radius, - imgColCoord, - imgRowCoord)) { - psF32 xDiff = (psF32) (imgColCoord - source->peak->x); - psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); - psF32 pDiff = source->pixels->data.F32[row][col] - sky; - - Sum+= pDiff; - X1+= xDiff * pDiff; - Y1+= yDiff * pDiff; - XY+= xDiff * yDiff * pDiff; - - X2+= PS_SQR(xDiff) * pDiff; - Y2+= PS_SQR(yDiff) * pDiff; - - if (source->pixels->data.F32[row][col] > peakPixel) { - peakPixel = source->pixels->data.F32[row][col]; - } - numPixels++; - } + // XXX EAM : mask == 0 is valid + if ((source->mask != NULL) && (source->mask->data.U8[row][col])) continue; + psS32 imgColCoord = col + source->pixels->col0; + psS32 imgRowCoord = row + source->pixels->row0; + if (CheckRadius(source->peak, + radius, + imgColCoord, + imgRowCoord)) { + psF32 xDiff = (psF32) (imgColCoord - source->peak->x); + psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); + psF32 pDiff = source->pixels->data.F32[row][col] - sky; + + Sum+= pDiff; + X1+= xDiff * pDiff; + Y1+= yDiff * pDiff; + XY+= xDiff * yDiff * pDiff; + + X2+= PS_SQR(xDiff) * pDiff; + Y2+= PS_SQR(yDiff) * pDiff; + + if (source->pixels->data.F32[row][col] > peakPixel) { + peakPixel = source->pixels->data.F32[row][col]; + } + numPixels++; } } } + // XXX EAM - the limit is a bit arbitrary. make it user defined? + if ((numPixels < 3) || (Sum <= 0)) { + psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n"); + return (false); + } + + psTrace (".psModules.pmSourceMoments", 5, + "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", + sky, Sum, X1, Y1, X2, Y2, XY, numPixels); // @@ -875,8 +889,15 @@ x = X1/Sum; y = Y1/Sum; + if ((fabs(x) > radius) || (fabs(y) > radius)) { + psTrace (".psModules.pmSourceMoments", 5, + "large centroid swing; invalid peak %d, %d\n", + source->peak->x, source->peak->y); + return (false); + } + source->moments->x = x + ((psF32) source->peak->x); source->moments->y = y + ((psF32) source->peak->y); - source->moments->Sxy = XY/Sum; + source->moments->Sxy = XY/Sum - x*y; source->moments->Sum = Sum; source->moments->Peak = peakPixel; @@ -886,5 +907,11 @@ 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); + + psTrace (".psModules.pmSourceMoments", 4, + "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", + sky, Sum, source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, source->moments->Sxy); + + return(true); // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it @@ -920,5 +947,5 @@ source->moments->Sx = (X2/Sum); source->moments->Sy = (Y2/Sum); - return(source); + return(true); # endif } @@ -956,29 +983,20 @@ /****************************************************************************** -pmSourceRoughClass(source, metadata): make a guess at the source -classification. - -XXX: This is not useable code, as of the release date. There remains a fair -bit of coding to be completed. - -XXX: The sigX and sigY stuff in the SDRS is unclear. - -XXX: How can this function ever return FALSE? -*****************************************************************************/ - -# define NPIX 10 -# define SCALE 0.1 - -// XXX I am ignore memory freeing issues (EAM) -bool pmSourceRoughClass(psArray *sources, psMetadata *metadata) -{ - PS_ASSERT_PTR_NON_NULL(sources, false); - PS_ASSERT_PTR_NON_NULL(metadata, false); - psBool rc = true; +pmSourcePSFClump(source, metadata): Find the likely PSF clump in the +sigma-x, sigma-y plane. return 0,0 clump in case of error. +*****************************************************************************/ + +pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) +{ + + # define NPIX 10 + # define SCALE 0.1 + psArray *peaks = NULL; - psF32 clumpX = 0.0; - psF32 clumpDX = 0.0; - psF32 clumpY = 0.0; - psF32 clumpDY = 0.0; + pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; + pmPSFClump psfClump = emptyClump; + + // PS_PTR_CHECK_NULL(sources, emptyClump); + // PS_PTR_CHECK_NULL(metadata, emptyClump); // find the sigmaX, sigmaY clump @@ -989,5 +1007,9 @@ // construct a sigma-plane image + // psImageAlloc does zero the data splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32); + for (int i = 0; i < splane->numRows; i++) { + memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); + } // place the sources in the sigma-plane image (ignore 0,0 values?) @@ -995,6 +1017,10 @@ { pmSource *tmpSrc = (pmSource *) sources->data[i]; - PS_ASSERT_PTR_NON_NULL(tmpSrc, false); // just skip this one? - PS_ASSERT_PTR_NON_NULL(tmpSrc->moments, false); // just skip this one? + if (tmpSrc == NULL) { + continue; + } + if (tmpSrc->moments == NULL) { + continue; + } // Sx,Sy are limited at 0. a peak at 0,0 is artificial @@ -1025,5 +1051,11 @@ psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); - } + psFree (splane); + psFree (stats); + + } + // XXX EAM : possible errors: + // 1) no peak in splane + // 2) no significant peak in splane // measure statistics on Sx, Sy if Sx, Sy within range of clump @@ -1039,5 +1071,5 @@ psArraySort (peaks, pmComparePeakDescend); clump = peaks->data[0]; - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d\n", clump->x, clump->y); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); // define section window for clump @@ -1080,24 +1112,54 @@ stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); - clumpX = stats->clippedMean; - clumpDX = stats->clippedStdev; + psfClump.X = stats->clippedMean; + psfClump.dX = stats->clippedStdev; stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); - clumpY = stats->clippedMean; - clumpDY = stats->clippedStdev; - - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", clumpX, clumpY); - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", clumpDX, clumpDY); + psfClump.Y = stats->clippedMean; + psfClump.dY = stats->clippedStdev; + + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); // these values should be pushed on the metadata somewhere - } - - int Nsat = 0; - int Ngal = 0; - int Nfaint = 0; - int Nstar = 0; - int Npsf = 0; - int Ncr = 0; + + psFree (stats); + psFree (peaks); + psFree (tmpSx); + psFree (tmpSy); + } + + return (psfClump); +} + +/****************************************************************************** +pmSourceRoughClass(source, metadata): make a guess at the source +classification. + +XXX: push the clump info into the metadata? + +XXX: How can this function ever return FALSE? +*****************************************************************************/ + +bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) +{ + + psBool rc = true; + + int Nsat = 0; + int Ngal = 0; + int Nstar = 0; + int Npsf = 0; + int Ncr = 0; + int Nsatstar = 0; + psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); starsn->n = 0; + + // check return status value (do these exist?) + bool status; + psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE"); + psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN"); + psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); + // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE"); // XXX allow clump size to be scaled relative to sigmas? @@ -1112,15 +1174,26 @@ 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) { + // calculate and save signal-to-noise estimates + psF32 S = tmpSrc->moments->Sum; + psF32 A = 4 * M_PI * sigX * sigY; + psF32 B = tmpSrc->moments->Sky; + psF32 RT = PS_SQRT_F32(S + (A * B) + (A * PS_SQR(RDNOISE) / PS_SQRT_F32(GAIN))); + psF32 SN = (S * PS_SQRT_F32(GAIN) / RT); + tmpSrc->moments->SN = SN; + + // XXX EAM : can we use the value of SATIRATE if mask is NULL? + int Nsatpix = pmCountSatPixels (tmpSrc->mask); + + // saturated star (size consistent with PSF or larger) + // Nsigma should be user-configured parameter + bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); + if ((Nsatpix > 1) && big) { + tmpSrc->type |= PS_SOURCE_SATSTAR; + Nsatstar ++; + continue; + } + + // saturated object (not a star, eg bleed trails, hot pixels) + if (Nsatpix > 1) { tmpSrc->type |= PS_SOURCE_SATURATED; Nsat ++; @@ -1128,6 +1201,8 @@ } - // too small to be stellar - if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) { + // likely defect (too small to be stellar) (push out to 3 sigma) + // low S/N objects which are small are probably stellar + // only set candidate defects if + if ((sigX < 0.05) || (sigY < 0.05)) { tmpSrc->type |= PS_SOURCE_DEFECT; Ncr ++; @@ -1135,6 +1210,6 @@ } - // possible galaxy - if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) { + // likely unsaturated galaxy (too large to be stellar) + if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { tmpSrc->type |= PS_SOURCE_GALAXY; Ngal ++; @@ -1143,23 +1218,11 @@ // the rest are probable stellar objects - psF32 S = tmpSrc->moments->Sum; - psF32 A = M_PI * sigX * sigY; - psF32 B = tmpSrc->moments->Sky; - psF32 RT = sqrtf(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrtf(GAIN))); - psF32 SN = (S * sqrtf(GAIN) / RT); - starsn->data.F32[starsn->n] = SN; starsn->n ++; Nstar ++; - // faint star - if (SN < FAINT_SN_LIM) { - tmpSrc->type |= PS_SOURCE_FAINTSTAR; - Nfaint ++; - continue; - } - - // PSF star - if (SN > PSF_SN_LIM) { + // PSF star (within 1.5 sigma of clump center + psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); + if ((SN > PSF_SN_LIM) && (radius < 1.5)) { tmpSrc->type |= PS_SOURCE_PSFSTAR; Npsf ++; @@ -1176,14 +1239,30 @@ stats = psVectorStats (stats, starsn, NULL, NULL, 0); psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); - } - - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nfaint: %3d\n", Nfaint); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal: %3d\n", Ngal); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); + psFree (stats); + psFree (starsn); + } + + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal: %3d\n", Ngal); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); return(rc); +} + +int pmCountSatPixels (psImage *image) +{ + int Nsatpix = 0; + + for (int i = 0; i < image->numRows; i++) { + for (int j = 0; j < image->numCols; j++) { + if (image->data.U8[i][j] & 0x02) { + Nsatpix ++; + } + } + } + return (Nsatpix); } @@ -1220,8 +1299,9 @@ psS32 SubImageCenterCol = source->peak->x; // XXX EAM : for the circle to stay on the image + // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); - psS32 SubImageEndRow = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32); + psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1); psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); - psS32 SubImageEndCol = PS_MIN (image->numCols - 1, SubImageCenterCol + radiusS32); + psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); // XXX EAM : this should not be needed: we can never hit this error @@ -1233,5 +1313,5 @@ return(false); } - if (SubImageEndRow >= image->numRows) { + if (SubImageEndRow > image->numRows) { psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", SubImageEndRow); @@ -1243,5 +1323,5 @@ return(false); } - if (SubImageEndCol >= image->numCols) { + if (SubImageEndCol > image->numCols) { psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", SubImageEndCol); @@ -1264,4 +1344,8 @@ // SubImageEndCol, // SubImageEndRow); + source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, + SubImageStartRow, + SubImageEndCol, + SubImageEndRow)); // XXX: Must recycle image. @@ -1271,5 +1355,5 @@ source->mask = psImageAlloc(source->pixels->numCols, source->pixels->numRows, - PS_TYPE_F32); + PS_TYPE_U8); // XXX EAM : type was F32 // @@ -1294,131 +1378,28 @@ /****************************************************************************** -pmSourceModelGuess(source, image, model): This function allocates a new -pmModel structure and stores it in the pmSource data structure specified in -the argument list. The model type is specified in the argument list. The -params array in that pmModel structure are allocated, and then set to the -appropriate values. This function returns true if everything was successful. +pmSourceModelGuess(source, model): This function allocates a new +pmModel structure based on the given modelType specified in the argument list. +The corresponding pmModelGuess function is returned, and used to +supply the values of the params array in the pmModel structure. XXX: Many of the initial parameters are set to 0.0 since I don't know what the appropiate initial guesses are. - -XXX: The image argument is redundant. XXX: Many parameters are based on the src->moments structure, which is in image, not subImage coords. Therefore, the calls to the model evaluation functions will be in image, not subImage coords. Remember this. - -XXX: The source->models member used to be allocated here. Now I assume -->modelPSF should be allocated -*****************************************************************************/ -bool pmSourceModelGuess(pmSource *source, - const psImage *image, - pmModelType model) -{ - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); - if (source->modelPSF != NULL) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->modelPSF was non-NULL; calling psFree(source->modelPSF).\n"); - psFree(source->modelPSF); - } - if (!((model == PS_MODEL_GAUSS) || - (model == PS_MODEL_PGAUSS) || - (model == PS_MODEL_WAUSS) || - (model == PS_MODEL_TWIST_GAUSS) || - (model == PS_MODEL_SERSIC) || - (model == PS_MODEL_SERSIC_CORE))) { - psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); - return(false); - } - - source->modelPSF = pmModelAlloc(model); - - psVector *params = source->modelPSF->params; - - switch (model) { - case PS_MODEL_GAUSS: - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = sqrt(2.0) / source->moments->Sx; - params->data.F32[5] = sqrt(2.0) / source->moments->Sy; - params->data.F32[6] = source->moments->Sxy; - return(true); - - case PS_MODEL_PGAUSS: - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = sqrt(2.0) / source->moments->Sx; - params->data.F32[5] = sqrt(2.0) / source->moments->Sy; - params->data.F32[6] = source->moments->Sxy; - return(true); - - case PS_MODEL_WAUSS: - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = sqrt(2.0) / source->moments->Sx; - params->data.F32[5] = sqrt(2.0) / source->moments->Sy; - params->data.F32[6] = source->moments->Sxy; - // XXX: What are these? - // source->modelPSF->params[7] = B2; - // source->modelPSF->params[8] = B3; - return(true); - - // XXX EAM : I might drop this model (or rather, replace it) - case PS_MODEL_TWIST_GAUSS: - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - // XXX: What are these? - // params->data.F32[4] = SxInner; - // params->data.F32[5] = SyInner; - // params->data.F32[6] = SxyInner; - // params->data.F32[7] = SxOuter; - // params->data.F32[8] = SyOuter; - // params->data.F32[9] = SxyOuter; - // params->data.F32[10] = N; - return(true); - - case PS_MODEL_SERSIC: - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = sqrt(2.0) / source->moments->Sx; - params->data.F32[5] = sqrt(2.0) / source->moments->Sy; - params->data.F32[6] = source->moments->Sxy; - // XXX: What are these? - //params->data.F32[7] = Nexp; - return(true); - - case PS_MODEL_SERSIC_CORE: - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - // XXX: What are these? - // params->data.F32[4] SxInner; - // params->data.F32[5] SyInner; - // params->data.F32[6] SxyInner; - // params->data.F32[7] Zd; - // params->data.F32[8] SxOuter; - // params->data.F32[9] SyOuter; - // params->data.F32[10] = SxyOuter; - // params->data.F32[11] = Nexp; - return(true); - - default: - psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); - return(false); - } +*****************************************************************************/ +psModel *pmSourceModelGuess(psSource *source, + psModelType modelType) +{ + PS_PTR_CHECK_NULL(source, false); + PS_PTR_CHECK_NULL(source->moments, false); + PS_PTR_CHECK_NULL(source->peak, false); + + psModel *model = psModelAlloc(modelType); + + psModelGuessFunc modelGuessFunc = psModelGuessFunc_GetFunction (modelType); + modelGuessFunc (model, source); + return(model); } @@ -1452,41 +1433,24 @@ PS_ASSERT_PTR_NON_NULL(src->modelPSF->params, false); - // XXX: The following step will not be necessary if the modelPSF->params - // member is a psVector. Suggest to IfA. - - // XXX EAM: done: modelPSF->params is now a vector - psVector *params = src->modelPSF->params; - - // +XXX EAM : I've made this a public function +XXX EAM : this now uses a psModel as the input +XXX EAM : it was using src->type to find the model, not model->type +*****************************************************************************/ +psF32 psModelEval(psModel *model, psImage *image, psS32 col, psS32 row) +{ + PS_PTR_CHECK_NULL(image, false); + PS_PTR_CHECK_NULL(model, false); + PS_PTR_CHECK_NULL(model->params, false); + // Allocate the x coordinate structure and convert row/col to image space. // psVector *x = psVectorAlloc(2, PS_TYPE_F32); - x->data.F32[0] = (psF32) (col + src->pixels->col0); - x->data.F32[1] = (psF32) (row + src->pixels->row0); + x->data.F32[0] = (psF32) (col + image->col0); + x->data.F32[1] = (psF32) (row + image->row0); psF32 tmpF; - - switch (src->modelPSF->type) { - case PS_MODEL_GAUSS: - tmpF = pmMinLM_Gauss2D(NULL, params, x); - break; - case PS_MODEL_PGAUSS: - tmpF = pmMinLM_PsuedoGauss2D(NULL, params, x); - break; - case PS_MODEL_TWIST_GAUSS: - tmpF = pmMinLM_TwistGauss2D(NULL, params, x); - break; - case PS_MODEL_WAUSS: - tmpF = pmMinLM_Wauss2D(NULL, params, x); - break; - case PS_MODEL_SERSIC: - tmpF = pmMinLM_Sersic(NULL, params, x); - break; - case PS_MODEL_SERSIC_CORE: - tmpF = pmMinLM_SersicCore(NULL, params, x); - break; - default: - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return(NAN); - } + psModelFunc modelFunc; + + modelFunc = psModelFunc_GetFunction (model->type); + tmpF = modelFunc (NULL, model->params, x); psFree(x); return(tmpF); @@ -1527,5 +1491,7 @@ } - psF32 oldValue = evalModel(source, subRow, subCol); + // XXX EAM : i changed this to match psModelEval above, but see + // XXX EAM the note below in pmSourceContour + psF32 oldValue = psModelEval(source->modelFLT, source->pixels, subCol, subRow); if (oldValue == level) { return(((psF32) (subCol + source->pixels->col0))); @@ -1548,5 +1514,5 @@ while (subCol != lastColumn) { - psF32 newValue = evalModel(source, subRow, subCol); + psF32 newValue = psModelEval(source->modelFLT, source->pixels, subCol, subRow); if (oldValue == level) { return((psF32) (subCol + source->pixels->col0)); @@ -1579,4 +1545,6 @@ XXX: What is mode? XXX: The top, bottom of the contour is not correctly determined. +XXX EAM : this functions is using the model for the contour, but it should + be using only the image counts *****************************************************************************/ psArray *pmSourceContour(pmSource *source, @@ -1590,5 +1558,6 @@ PS_ASSERT_PTR_NON_NULL(source->peak, false); PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->modelPSF, false); + PS_ASSERT_PTR_NON_NULL(source->modelFLT, false); + // XXX EAM : what is the purpose of modelPSF/modelFLT? // @@ -1678,25 +1647,17 @@ } -#if 0 -static psVector *minLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x); -static psVector *minLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); -static psVector *minLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x); -static psVector *minLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); -static psVector *minLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x); -static psVector *minLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x); -#endif - // XXX EAM : these are better starting values, but should be available from metadata? -#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 20 +#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 /****************************************************************************** -pmSourceFitModel(source, image): must create the appropiate arguments to the +pmSourceFitModel(source, model): must create the appropiate arguments to the LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. XXX: should there be a mask value? -XXX: Probably should remove the "image" argument. +XXX EAM : fit the specified model (not necessarily the one in source) *****************************************************************************/ bool pmSourceFitModel(pmSource *source, - const psImage *image) + pmModel *model, + const bool PSF) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -1704,8 +1665,20 @@ PS_ASSERT_PTR_NON_NULL(source->peak, false); PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->modelPSF, false); - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); - psBool rc; + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->noise, false); + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + // XXX EAM : is it necessary for the mask & noise to exist? the + // tests below could be conditions (!NULL) + + psModelFunc modelFunc = psModelFunc_GetFunction (model->type); + + psVector *params = model->params; + psVector *dparams = model->dparams; + psVector *paramMask = NULL; + + int nParams = PSF ? params->n - 4 : params->n; // find the number of valid pixels @@ -1719,4 +1692,8 @@ } } + } + if (count < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + return(false); } @@ -1736,7 +1713,7 @@ x->data[tmpCnt] = (psPtr *) coord; y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; - - // XXX EAM : this is approximate: need to apply the gain and rdnoise - yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j])); + yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j]); + // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then + // the minimization function calculates sq() tmpCnt++; } @@ -1747,59 +1724,81 @@ PM_SOURCE_FIT_MODEL_TOLERANCE); - psVector *params = source->modelPSF->params; - - switch (source->modelPSF->type) { - case PS_MODEL_GAUSS: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D); - break; - case PS_MODEL_PGAUSS: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_PsuedoGauss2D); - break; - case PS_MODEL_TWIST_GAUSS: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D); - break; - case PS_MODEL_WAUSS: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D); - break; - case PS_MODEL_SERSIC: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Sersic); - break; - case PS_MODEL_SERSIC_CORE: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore); - break; - default: - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - rc = false; - } + // PSF model only fits first 4 parameters, FLT model fits all + if (PSF) { + paramMask = psVectorAlloc (params->n, PS_TYPE_U8); + for (int i = 0; i < 4; i++) { + paramMask->data.U8[i] = 0; + } + for (int i = 4; i < paramMask->n; i++) { + paramMask->data.U8[i] = 1; + } + } + + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); + + psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if ((paramMask != NULL) && paramMask->data.U8[i]) continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + // XXX EAM: we need to do something (give an error?) if rc is false + // XXX EAM: psMinimizeLMChi2 does not check convergence + + // XXX models can go insane: reject these + onPic &= (params->data.F32[2] >= source->pixels->col0); + onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (params->data.F32[3] >= source->pixels->row0); + onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + // XXX EAM: save the resulting chisq, nDOF, nIter - source->modelPSF->chisq = myMin->value; - source->modelPSF->nDOF = y->n - params->n; - source->modelPSF->nIter = myMin->iter; + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // XXX EAM get the Gauss-Newton distance for fixed model parameters + if (paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if (!paramMask->data.U8[i]) continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + psFree (delta); + } psFree(x); psFree(y); + psFree(yErr); psFree(myMin); + psFree(covar); + psFree(paramMask); + + rc = (onPic && fitStatus); return(rc); } -static bool sourceAddOrSubModel(psImage *image, - pmSource *src, - bool center, - psS32 flag) -{ - PS_ASSERT_PTR_NON_NULL(src, false); - PS_ASSERT_PTR_NON_NULL(src->moments, false); - PS_ASSERT_PTR_NON_NULL(src->peak, false); - PS_ASSERT_PTR_NON_NULL(src->pixels, false); - PS_ASSERT_PTR_NON_NULL(src->modelPSF, false); - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); +bool p_pmSourceAddOrSubModel(psImage *image, + psImage *mask, + psModel *model, + bool center, + psS32 flag) +{ + + // XXX EAM : convert to ASSERTS + // PS_PTR_CHECK_NULL(model, false); + // PS_IMAGE_CHECK_NULL(image, false); + // PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); psVector *x = psVectorAlloc(2, PS_TYPE_F32); - psVector *params = src->modelPSF->params; - - for (psS32 i = 0; i < src->pixels->numRows; i++) { - for (psS32 j = 0; j < src->pixels->numCols; j++) { + psVector *params = model->params; + psModelFunc modelFunc = psModelFunc_GetFunction (model->type); + psS32 imageCol; + psS32 imageRow; + + for (psS32 i = 0; i < image->numRows; i++) { + for (psS32 j = 0; j < image->numCols; j++) { + if ((mask != NULL) && mask->data.U8[i][j]) continue; psF32 pixelValue; // XXX: Should you be adding the pixels for the entire subImage, @@ -1808,33 +1807,19 @@ // Convert i/j to imace coord space: // XXX: Make sure you have col/row order correct. - psS32 imageRow = i + src->pixels->row0; - psS32 imageCol = j + src->pixels->col0; + // XXX EAM : 'center' option changes this + // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2] + if (center) { + imageCol = j - 0.5*image->numCols + model->params->data.F32[2]; + imageRow = i - 0.5*image->numRows + model->params->data.F32[3]; + } else { + imageCol = j + image->col0; + imageRow = i + image->row0; + } x->data.F32[0] = (float) imageCol; x->data.F32[1] = (float) imageRow; - switch (src->modelPSF->type) { - case PS_MODEL_GAUSS: - pixelValue = pmMinLM_Gauss2D(NULL, params, x); - break; - case PS_MODEL_PGAUSS: - pixelValue = pmMinLM_PsuedoGauss2D(NULL, params, x); - break; - case PS_MODEL_TWIST_GAUSS: - pixelValue = pmMinLM_TwistGauss2D(NULL, params, x); - break; - case PS_MODEL_WAUSS: - pixelValue = pmMinLM_Wauss2D(NULL, params, x); - break; - case PS_MODEL_SERSIC: - pixelValue = pmMinLM_Sersic(NULL, params, x); - break; - case PS_MODEL_SERSIC_CORE: - pixelValue = pmMinLM_SersicCore(NULL, params, x); - break; - default: - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - psFree(x); - return(false); - } + pixelValue = modelFunc (NULL, params, x); + // fprintf (stderr, "%f %f %d %d %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue); + if (flag == 1) { pixelValue = -pixelValue; @@ -1844,5 +1829,5 @@ // how to use the boolean "center" flag. - image->data.F32[imageRow][imageCol]+= pixelValue; + image->data.F32[i][j]+= pixelValue; } } @@ -1856,8 +1841,9 @@ *****************************************************************************/ bool pmSourceAddModel(psImage *image, - pmSource *src, + psImage *mask, + psModel *model, bool center) { - return(sourceAddOrSubModel(image, src, center, 0)); + return(p_pmSourceAddOrSubModel(image, mask, model, center, 0)); } @@ -1865,8 +1851,9 @@ *****************************************************************************/ bool pmSourceSubModel(psImage *image, - pmSource *src, + psImage *mask, + psModel *model, bool center) { - return(sourceAddOrSubModel(image, src, center, 1)); + return(p_pmSourceAddOrSubModel(image, mask, model, center, 1)); } @@ -1880,244 +1867,2 @@ return(RVAL); \ } - - -/** - all of these object representation functions have the same form : func(*deriv, *params, *x) - - the argument "x" contains a single "x,y" coordinate pair. The function computes the object - model, based on the parameters in "params" at the x,y point specified by *x, and returns the value. - The derivatives are also caculated and returned in the "deriv" argument. parameter error checking is - skipped because speed is most important. -**/ - -/****************************************************************************** - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = sqrt(2.0) / SigmaX; - params->data.F32[5] = sqrt(2.0) / SigmaY; - params->data.F32[6] = Sxy; -*****************************************************************************/ -float pmMinLM_Gauss2D( - psVector *deriv, - const psVector *params, - const psVector *x) -{ - PS_ASSERT_VECTOR_NON_NULL(params, NAN); - PS_ASSERT_VECTOR_NON_NULL(x, NAN); - psF32 X = x->data.F32[0] - params->data.F32[2]; - psF32 Y = x->data.F32[1] - params->data.F32[3]; - psF32 px = params->data.F32[4]*X; - psF32 py = params->data.F32[5]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; - psF32 r = exp(-z); - psF32 q = params->data.F32[1]*r; - psF32 f = q + params->data.F32[0]; - - if (deriv != NULL) { - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); - deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); - deriv->data.F32[4] = -2.0*q*px*X; - deriv->data.F32[5] = -2.0*q*py*Y; - deriv->data.F32[6] = -q*X*Y; - } - return(f); -} - -/****************************************************************************** - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = sqrt(2) / SigmaX; - params->data.F32[5] = sqrt(2) / SigmaY; - params->data.F32[6] = Sxy; -*****************************************************************************/ -float pmMinLM_PsuedoGauss2D( - psVector *deriv, - const psVector *params, - const psVector *x) -{ - PS_ASSERT_VECTOR_NON_NULL(params, NAN); - PS_ASSERT_VECTOR_NON_NULL(x, NAN); - psF32 X = x->data.F32[0] - params->data.F32[2]; - psF32 Y = x->data.F32[1] - params->data.F32[3]; - psF32 px = params->data.F32[4]*X; - psF32 py = params->data.F32[5]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; - psF32 t = 1 + z + 0.5*z*z; - psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */ - psF32 f = params->data.F32[1]*r + params->data.F32[0]; - - if (deriv != NULL) { - // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 q = params->data.F32[1]*r*r*t; - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); - deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); - deriv->data.F32[4] = -2.0*q*px*X; - deriv->data.F32[5] = -2.0*q*py*Y; - deriv->data.F32[6] = -q*X*Y; - } - return(f); -} - -/****************************************************************************** - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = Sx; - params->data.F32[5] = Sy; - params->data.F32[6] = Sxy; - params->data.F32[7] = B2; - params->data.F32[8] = B3; -*****************************************************************************/ -float pmMinLM_Wauss2D( - psVector *deriv, - const psVector *params, - const psVector *x) -{ - PS_ASSERT_VECTOR_NON_NULL(params, NAN); - PS_ASSERT_VECTOR_NON_NULL(x, NAN); - psF32 X = x->data.F32[0] - params->data.F32[2]; - psF32 Y = x->data.F32[1] - params->data.F32[2]; - psF32 px = params->data.F32[4]*X; - psF32 py = params->data.F32[5]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; - psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0); - psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */ - psF32 f = params->data.F32[1]*r + params->data.F32[0]; - - if (deriv != NULL) { - // note difference from gaussian: q = params->data.F32[1]*r - psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); - deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); - deriv->data.F32[4] = -2.0*q*px*X; - deriv->data.F32[5] = -2.0*q*py*Y; - deriv->data.F32[6] = -q*X*Y; - deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; - deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; - // The values of 100 dampen the swing of params->data.F32[7,8] */ - } - return(f); -} - -// XXX: What should these be? -#define FFACTOR 1.0 -#define FSCALE 1.0 -/****************************************************************************** - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = SxInner; - params->data.F32[5] = SyInner; - params->data.F32[6] = SxyInner; - params->data.F32[7] = SxOuter; - params->data.F32[8] = SyOuter; - params->data.F32[9] = SxyOuter; - params->data.F32[10] = N; -*****************************************************************************/ -float pmMinLM_TwistGauss2D( - psVector *deriv, - const psVector *params, - const psVector *x) -{ - PS_ASSERT_VECTOR_NON_NULL(params, NAN); - PS_ASSERT_VECTOR_NON_NULL(x, NAN); - psF32 X = x->data.F32[0] - params->data.F32[2]; - psF32 Y = x->data.F32[1] - params->data.F32[3]; - psF32 px1 = params->data.F32[4]*X; - psF32 py1 = params->data.F32[5]*Y; - psF32 px2 = params->data.F32[7]*X; - psF32 py2 = params->data.F32[8]*Y; - psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y; - psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y; - psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10])); - - psF32 f = params->data.F32[5]*r + params->data.F32[6]; - - if (deriv != NULL) { - psF32 q1 = params->data.F32[5]*PS_SQR(r); - psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y); - deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X); - - // These fudge factors impede the growth of params->data.F32[4] beyond - // params->data.F32[7]. - psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); - psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; - deriv->data.F32[4] = -2.0*q1*px1*X*f2; - - // These fudge factors impede the growth of params->data.F32[5] beyond - // params->data.F32[8]. - f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); - f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; - deriv->data.F32[5] = -2.0*q1*py1*Y*f2; - deriv->data.F32[6] = -q1*X*Y; - deriv->data.F32[7] = -2.0*q2*px2*X; - deriv->data.F32[8] = -2.0*q2*py2*Y; - deriv->data.F32[9] = -q2*X*Y; - deriv->data.F32[10] = -q1*log(z2); - } - - return(f); -} - -/****************************************************************************** - float Sersic() - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = Sx; - params->data.F32[5] = Sy; - params->data.F32[6] = Sxy; - params->data.F32[7] = Nexp; -*****************************************************************************/ -float pmMinLM_Sersic( - psVector *deriv, - const psVector *params, - const psVector *x) -{ - PS_ASSERT_VECTOR_NON_NULL(params, NAN); - PS_ASSERT_VECTOR_NON_NULL(x, NAN); - psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); - return(0.0); -} - -/****************************************************************************** - float SersicBulge() - params->data.F32[0] So; - params->data.F32[1] Zo; - params->data.F32[2] Xo; - params->data.F32[3] Yo; - params->data.F32[4] SxInner; - params->data.F32[5] SyInner; - params->data.F32[6] SxyInner; - params->data.F32[7] Zd; - params->data.F32[8] SxOuter; - params->data.F32[9] SyOuter; - params->data.F32[10] = SxyOuter; - params->data.F32[11] = Nexp; -*****************************************************************************/ -float pmMinLM_SersicCore( - psVector *deriv, - const psVector *params, - const psVector *x) -{ - PS_ASSERT_VECTOR_NON_NULL(params, NAN); - PS_ASSERT_VECTOR_NON_NULL(x, NAN); - psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); - return(0.0); -}