Index: trunk/psphot/src/pmObjects_EAM.c =================================================================== --- trunk/psphot/src/pmObjects_EAM.c (revision 4977) +++ trunk/psphot/src/pmObjects_EAM.c (revision 5048) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-08 04:32:40 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-14 01:35:20 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -603,131 +603,31 @@ members. *****************************************************************************/ -pmSource *pmSourceLocalSky(const psImage *image, - const pmPeak *peak, - psStatsOptions statsOptions, - psF32 innerRadius, - psF32 outerRadius) -{ - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - PS_ASSERT_PTR_NON_NULL(peak, NULL); - PS_FLOAT_COMPARE(0.0, innerRadius, NULL); - PS_FLOAT_COMPARE(innerRadius, outerRadius, NULL); - psS32 innerRadiusS32 = (psS32) innerRadius; - psS32 outerRadiusS32 = (psS32) outerRadius; - - // - // We define variables for code readability. - // - // XXX: Since the peak->xy coords are in image, not subImage coords, - // these variables should be renamed for clarity (imageCenterRow, etc). - // - // peak->x,y is guaranteed to be on image - psS32 SubImageCenterRow = peak->y; - psS32 SubImageCenterCol = peak->x; - - // XXX EAM : I added this code to stay on the image. So did George - // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) - // XXX EAM : dropped the Annulus width (not needed) - // XXX EAM : dropped off-boundary tests (not needed) - 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); - // - // Grab a subimage of the original image of size (2 * outerRadius). - // - // XXX EAM : merged psImageSubset & psRegionSet - // XXX EAM : cast for image is needed because of const (above) - psImage *subImage = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, - SubImageEndCol, - SubImageStartRow, - SubImageEndRow)); - - psImage *subImageMask = psImageAlloc(subImage->numCols, - subImage->numRows, - PS_TYPE_U8); - // XXX EAM : for consistency, mask needs col0,row0 set to match image - // XXX EAM : CONFLICT between psLib and code requirement. FIX PSLIB!!! - // subImageMask->col0 = subImage->col0; - // subImageMask->row0 = subImage->row0; - psAbort ("pmObjects", "must fix this error before psphot will work!!!\n"); - - // - // Loop through the subimage mask & initialize PSPHOT_MASK_INVALID - // - for (psS32 row = 0 ; row < subImageMask->numRows; row++) { - for (psS32 col = 0 ; col < subImageMask->numCols; col++) { - subImageMask->data.U8[row][col] = PSPHOT_MASK_INVALID; - } - } - - // - // Loop through the subimage, mask off pixels in the inner square. - // 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 - // XXX EAM : this masks the pixels in the center region - // XXX EAM : should we be using psImageMaskRegion ??? - 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] = PSPHOT_MASK_CLEAR; - } - } - - // XXX EAM : make this trace information?? - // for (psS32 row = 0 ; row < subImage->numRows; row++) { - // for (psS32 col = 0 ; col < subImage->numCols; col++) { - // printf("(%d) ", subImageMask->data.U8[row][col]); - // } - // printf("\n"); - // } - - // - // Allocate the myStats structure, then call psImageStats(), which will - // calculate the specified statistic. - // +bool pmSourceLocalSky (pmSource *source, + psStatsOptions statsOptions, + psF32 Radius) +{ + + psImage *image = source->pixels; + psImage *mask = source->mask; + pmPeak *peak = source->peak; + psRegion srcRegion; + + // XXX EAM : psRegionXXX funcs should take value not ptr + srcRegion = psRegionForSquare (peak->x, peak->y, Radius); + srcRegion = psRegionForImage (mask, &srcRegion); + + psImageMaskRegion (mask, &srcRegion, "OR", PSPHOT_MASK_MARKED); psStats *myStats = psStatsAlloc(statsOptions); - myStats = psImageStats(myStats, subImage, subImageMask, 1); - - // - // Create the output mySource, and set appropriate members. - // - pmSource *mySource = pmSourceAlloc(); - mySource->peak = (pmPeak *) peak; - mySource->moments = pmMomentsAlloc(); + myStats = psImageStats(myStats, image, mask, 0xff); + psImageMaskRegion (mask, &srcRegion, "AND", ~PSPHOT_MASK_MARKED); + psF64 tmpF64; p_psGetStatValue(myStats, &tmpF64); - mySource->moments->Sky = (psF32) tmpF64; - mySource->pixels = subImage; - mySource->mask = subImageMask; - - // - // Free things. XXX: This should be static memory. - // - psFree(myStats); - - return(mySource); -} - -/****************************************************************************** -bool checkRadius(*peak, radius, x, y): private function which simply -determines if the (x, y) point is within the radius of the specified peak. - -XXX: macro this for performance. -*****************************************************************************/ -static bool checkRadius(pmPeak *peak, - psF32 radius, - psS32 x, - psS32 y) -{ - if (PS_SQR(radius) >= (psF32) (PS_SQR(x - peak->x) + PS_SQR(y - peak->y))) { - return(true); - } - - return(false); + psFree (myStats); + + if (isnan(tmpF64)) return (false); + source->moments = pmMomentsAlloc(); + source->moments->Sky = (psF32) tmpF64; + return (true); } @@ -762,15 +662,21 @@ pmSource->peak pmSource->pixels + pmSource->noise + pmSource->mask XXX: The peak calculations are done in image coords, not subImage coords. -XXX: mask values? -*****************************************************************************/ +XXX EAM : this version clips input pixels on S/N +XXX EAM : this version returns false for several reasons +*****************************************************************************/ +# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) + bool pmSourceMoments(pmSource *source, psF32 radius) { - PS_ASSERT_PTR_NON_NULL(source, NULL); - PS_ASSERT_PTR_NON_NULL(source->peak, NULL); - PS_ASSERT_PTR_NON_NULL(source->pixels, NULL); + // PS_PTR_CHECK_NULL(source, NULL); + // PS_PTR_CHECK_NULL(source->peak, NULL); + // PS_PTR_CHECK_NULL(source->pixels, NULL); + // PS_PTR_CHECK_NULL(source->noise, NULL); PS_FLOAT_COMPARE(0.0, radius, NULL); @@ -791,7 +697,7 @@ // XY = SUM (x - xc)*(y - yc)*(z - sky) // - psF32 Sum = 0.0; psF32 peakPixel = -PS_MAX_F32; psS32 numPixels = 0; + psF32 Sum = 0.0; psF32 X1 = 0.0; psF32 Y1 = 0.0; @@ -799,7 +705,11 @@ psF32 Y2 = 0.0; psF32 XY = 0.0; - psF32 x = 0; - psF32 y = 0; - // + psF32 x = 0; + psF32 y = 0; + psF32 R2 = PS_SQR(radius); + + psF32 xPeak = source->peak->x; + psF32 yPeak = source->peak->y; + // XXX why do I get different results for these two methods of finding Sx? // XXX Sx, Sy would be better measured if we clip pixels close to sky @@ -809,38 +719,39 @@ // proceed with the moments calculation. need to do two loops for a // numerically stable result. first loop: get the sums. - // + // XXX EAM : mask == 0 is valid + for (psS32 row = 0; row < source->pixels->numRows ; row++) { for (psS32 col = 0; col < source->pixels->numCols ; col++) { - // 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++; - } + + psF32 xDiff = col + source->pixels->col0 - xPeak; + psF32 yDiff = row + source->pixels->row0 - yPeak; + + // XXX EAM : calculate xDiff, yDiff up front; + // radius is just a function of (xDiff, yDiff) + if (!VALID_RADIUS(xDiff, yDiff, R2)) continue; + + psF32 pDiff = source->pixels->data.F32[row][col] - sky; + + // XXX EAM : check for valid S/N in pixel + // XXX EAM : should this limit be user-defined? + if (pDiff / sqrt(source->noise->data.F32[row][col]) < 1) continue; + + Sum += pDiff; + X1 += xDiff * pDiff; + Y1 += yDiff * pDiff; + XY += xDiff * yDiff * pDiff; + + X2 += PS_SQR(xDiff) * pDiff; + Y2 += PS_SQR(yDiff) * pDiff; + + peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel); + numPixels++; } } // XXX EAM - the limit is a bit arbitrary. make it user defined? if ((numPixels < 3) || (Sum <= 0)) { - psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n"); - return (false); + psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n"); + return (false); } @@ -857,13 +768,14 @@ 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); - + psTrace (".psModules.pmSourceMoments", 5, + "large centroid swing; invalid peak %d, %d\n", + source->peak->x, source->peak->y); + return (false); + } + + source->moments->x = x + xPeak; + source->moments->y = y + yPeak; + + // XXX EAM : Sxy needs to have x*y subtracted source->moments->Sxy = XY/Sum - x*y; source->moments->Sum = Sum; @@ -919,4 +831,5 @@ *****************************************************************************/ +// XXX EAM include a S/N cutoff in selecting the sources? pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) { @@ -1123,5 +1036,5 @@ bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); if ((Nsatpix > 1) && big) { - tmpSrc->type |= PS_SOURCE_SATSTAR; + tmpSrc->type |= PM_SOURCE_SATSTAR; Nsatstar ++; continue; @@ -1130,5 +1043,5 @@ // saturated object (not a star, eg bleed trails, hot pixels) if (Nsatpix > 1) { - tmpSrc->type |= PS_SOURCE_SATURATED; + tmpSrc->type |= PM_SOURCE_SATURATED; Nsat ++; continue; @@ -1139,5 +1052,5 @@ // only set candidate defects if if ((sigX < 0.05) || (sigY < 0.05)) { - tmpSrc->type |= PS_SOURCE_DEFECT; + tmpSrc->type |= PM_SOURCE_DEFECT; Ncr ++; continue; @@ -1146,5 +1059,5 @@ // 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; + tmpSrc->type |= PM_SOURCE_GALAXY; Ngal ++; continue; @@ -1159,5 +1072,5 @@ 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; + tmpSrc->type |= PM_SOURCE_PSFSTAR; Npsf ++; continue; @@ -1165,5 +1078,5 @@ // random type of star - tmpSrc->type |= PS_SOURCE_OTHER; + tmpSrc->type |= PM_SOURCE_OTHER; } @@ -1197,4 +1110,6 @@ XXX: The circle will have a diameter of (1+radius). This is different from the pmSourceSetLocal() function. + +XXX: this should probably use the psImageMaskCircle Function for consistency *****************************************************************************/ bool pmSourceSetPixelsCircle(pmSource *source, @@ -1537,7 +1452,7 @@ XXX EAM : fit the specified model (not necessarily the one in source) *****************************************************************************/ -bool pmSourceFitModel(pmSource *source, - pmModel *model, - const bool PSF) +bool pmSourceFitModel_v5(pmSource *source, + pmModel *model, + const bool PSF) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -1554,9 +1469,9 @@ // tests below could be conditions (!NULL) - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - psVector *params = model->params; psVector *dparams = model->dparams; psVector *paramMask = NULL; + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); int nParams = PSF ? params->n - 4 : params->n; @@ -1615,5 +1530,6 @@ } - psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); + // XXX EAM : covar must be F64? + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); @@ -1621,5 +1537,5 @@ for (int i = 0; i < dparams->n; i++) { if ((paramMask != NULL) && paramMask->data.U8[i]) continue; - dparams->data.F32[i] = sqrt(covar->data.F32[i][i]); + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); } @@ -1640,9 +1556,9 @@ // XXX EAM get the Gauss-Newton distance for fixed model parameters if (paramMask != NULL) { - psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32); + 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.F32[i]; + dparams->data.F32[i] = delta->data.F64[i]; } psFree (delta); @@ -1655,4 +1571,152 @@ psFree(covar); psFree(paramMask); + + rc = (onPic && fitStatus); + return(rc); +} + +// XXX EAM : new version with parameter range limits and noise enhancement +bool pmSourceFitModel (pmSource *source, + pmModel *model, + const bool PSF) +{ + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->noise, false); + + // XXX EAM : is it necessary for the mask & noise to exist? the + // tests below could be conditions (!NULL) + + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + psF32 Ro, ymodel; + + psVector *params = model->params; + psVector *dparams = model->dparams; + psVector *paramMask = NULL; + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + + // XXX EAM : I need to use the sky value to constrain the noise model + int nParams = PSF ? params->n - 4 : params->n; + psF32 So = params->data.F32[0]; + + // find the number of valid pixels + // XXX EAM : this loop and the loop below could just be one pass + // using the psArrayAdd and psVectorExtend functions + psS32 count = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j] == 0) { + count++; + } + } + } + if (count < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + return(false); + } + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(count); + psVector *y = psVectorAlloc(count, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); + psS32 tmpCnt = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j] == 0) { + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + // XXX: Convert i/j to image space: + // XXX EAM: coord order is (x,y) == (col,row) + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[tmpCnt] = (psPtr *) coord; + y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; + + // compare observed flux to model flux to adjust weight + ymodel = modelFunc (NULL, model->params, coord); + + // this test enhances the noise based on deviation from the model flux + Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So)); + + // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then + // the minimization function calculates sq(). + // should psMinimizeLMChi2 take dY^2? + yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * Ro); + tmpCnt++; + } + } + } + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + // PSF model only fits first 4 parameters, FLT model fits all + if (PSF) { + paramMask = psVectorAlloc (params->n, PS_TYPE_U8); + for (int i = 0; i < 4; i++) { + paramMask->data.U8[i] = 0; + } + for (int i = 4; i < paramMask->n; i++) { + paramMask->data.U8[i] = 1; + } + } + + // XXX EAM : I've added three types of parameter range checks + // XXX EAM : this requires my new psMinimization functions + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); + psVector *beta_lim = NULL; + psVector *params_min = NULL; + psVector *params_max = NULL; + + // XXX EAM : in this implementation, I pass in the limits with the covar matrix. + // in the SDRS, I define a new psMinimization which will take these in + psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); + modelLimits (&beta_lim, ¶ms_min, ¶ms_max); + for (int i = 0; i < params->n; i++) { + covar->data.F64[0][i] = beta_lim->data.F32[i]; + covar->data.F64[1][i] = params_min->data.F32[i]; + covar->data.F64[2][i] = params_max->data.F32[i]; + } + + psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2_EAM(myMin, covar, params, paramMask, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if ((paramMask != NULL) && paramMask->data.U8[i]) continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // XXX EAM: we need to do something (give an error?) if rc is false + // XXX EAM: psMinimizeLMChi2 does not check convergence + + // XXX models can go insane: reject these + onPic &= (params->data.F32[2] >= source->pixels->col0); + onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (params->data.F32[3] >= source->pixels->row0); + onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + + // XXX EAM: save the resulting chisq, nDOF, nIter + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // XXX EAM get the Gauss-Newton distance for fixed model parameters + if (paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if (!paramMask->data.U8[i]) continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + } + + psFree(paramMask); + psFree(x); + psFree(y); + psFree(myMin); rc = (onPic && fitStatus);