Index: /trunk/psphot/doc/notes.txt =================================================================== --- /trunk/psphot/doc/notes.txt (revision 5047) +++ /trunk/psphot/doc/notes.txt (revision 5048) @@ -50,4 +50,6 @@ - fixed 'center' option in pmSourceAdd/Sub + - pmPSF params should be F64 + lingering concerns: - Polynomials and related functions are STILL messed up. this is a Index: /trunk/psphot/src/pmObjects_EAM.c =================================================================== --- /trunk/psphot/src/pmObjects_EAM.c (revision 5047) +++ /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); Index: /trunk/psphot/src/pmObjects_EAM.h =================================================================== --- /trunk/psphot/src/pmObjects_EAM.h (revision 5047) +++ /trunk/psphot/src/pmObjects_EAM.h (revision 5048) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-09 18:31:29 $ + * @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 @@ -210,11 +210,9 @@ *****************************************************************************/ -pmSource *pmSourceLocalSky( - const psImage *image, ///< The input image (float) - const pmPeak *peak, ///< The peak for which the pmSource struct is created. +bool pmSourceLocalSky( + pmSource *source, ///< The input image (float) psStatsOptions statsOptions, ///< The statistic used in calculating the background sky - float innerRadius, ///< The inner radius of the suqare annulus for calculating sky - float outerRadius ///< The outer radius of the suqare annulus for calculating sky -); + float Radius ///< The inner radius of the square annulus to exclude + ); /****************************************************************************** @@ -298,3 +296,15 @@ ); +bool pmSourceFitModel_v5( + pmSource *source, ///< The input pmSource + pmModel *model, ///< model to be fitted + const bool PSF ///< Treat model as PSF or FLT? +); + +bool pmSourceFitModel_v7( + pmSource *source, ///< The input pmSource + pmModel *model, ///< model to be fitted + const bool PSF ///< Treat model as PSF or FLT? +); + #endif Index: /trunk/psphot/src/psLibUtils.c =================================================================== --- /trunk/psphot/src/psLibUtils.c (revision 5047) +++ /trunk/psphot/src/psLibUtils.c (revision 5048) @@ -147,6 +147,7 @@ } +# if (0) // alternate implementation of this function from pmObjects.c -// XXX EAM : move this to pmObjects_EAM.c +// now defined in psLib SDRS as psImageRow psVector *psGetRowVectorFromImage(psImage *image, psU32 row) @@ -159,4 +160,5 @@ return(tmpVector); } +# endif // XXX EAM : this is now in psLib Index: /trunk/psphot/src/psLibUtils.h =================================================================== --- /trunk/psphot/src/psLibUtils.h (revision 5047) +++ /trunk/psphot/src/psLibUtils.h (revision 5048) @@ -39,5 +39,5 @@ int psArgumentRemove (int N, int *argc, char **argv); // added to SDRS psVector *psVectorCreate (double lower, double upper, double delta, psElemType type); // added to SDRS -psVector *psGetRowVectorFromImage(psImage *image, psU32 row); // added to SDRS (as psImageRow) +// psVector *psGetRowVectorFromImage(psImage *image, psU32 row); // added to SDRS (as psImageRow) int psImageCountPixelMask (psImage *mask, psU8 value); // added to SDRS @@ -56,4 +56,5 @@ // polynomial functions -- all in SDRS, not done for v.15 +// XXX the pslib 0.7.0 polynomials are *still* nTerm, not nOrder!!! psPolynomial1D *Polynomial1DAlloc_EAM(psPolynomialType type, psS32 nOrder); void Polynomial1DDump_EAM(psPolynomial1D *poly); Index: /trunk/psphot/src/psMinimize.c =================================================================== --- /trunk/psphot/src/psMinimize.c (revision 5047) +++ /trunk/psphot/src/psMinimize.c (revision 5048) @@ -95,5 +95,5 @@ // iterate until the tolerance is reached, or give up - while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { + while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) { // set a new guess for Alpha, Beta, Params Index: /trunk/psphot/src/psModulesUtils.c =================================================================== --- /trunk/psphot/src/psModulesUtils.c (revision 5047) +++ /trunk/psphot/src/psModulesUtils.c (revision 5048) @@ -2,4 +2,5 @@ // extract config informatin from config data or from image header, if specified +// XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) { @@ -48,315 +49,4 @@ } -// XXX EAM : these are my alternate implementations of psModule functions - -// this sets and clears bit 0x80 -bool pmSourceLocalSky_EAM (pmSource *source, - psStatsOptions statsOptions, - psF32 Radius) -{ - - psImage *image = source->pixels; - psImage *mask = source->mask; - pmPeak *peak = source->peak; - psRegion srcRegion; - - // XXX EAM : psRegionXXX funcs need 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, image, mask, 0xff); - psImageMaskRegion (mask, &srcRegion, "AND", ~PSPHOT_MASK_MARKED); - - psF64 tmpF64; - p_psGetStatValue(myStats, &tmpF64); - psFree (myStats); - - if (isnan(tmpF64)) return (false); - source->moments = pmMomentsAlloc(); - source->moments->Sky = (psF32) tmpF64; - return (true); -} - -# define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 -# define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 -bool pmSourceFitModel_EAM (pmSource *source, - pmModel *model, - const bool PSF) -{ - // PS_PTR_CHECK_NULL(source, false); - // PS_PTR_CHECK_NULL(source->moments, false); - // PS_PTR_CHECK_NULL(source->peak, false); - // PS_PTR_CHECK_NULL(source->pixels, false); - // PS_PTR_CHECK_NULL(source->mask, false); - // PS_PTR_CHECK_NULL(source->noise, false); - psBool fitStatus = true; - psBool onPic = true; - psBool rc = true; - psF32 Ro, ymodel; - - - // XXX EAM : is it necessary for the mask & noise to exist? the - // tests below could be conditions (!NULL) - - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); - - psVector *params = model->params; - psVector *dparams = model->dparams; - psVector *paramMask = NULL; - - psVector *beta_lim = NULL; - psVector *params_min = NULL; - psVector *params_max = NULL; - - 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]; - - 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)); - yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * Ro); - // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then - // the minimization function calculates sq() - 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; - } - } - - 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); - return(rc); -} - -/****************************************************************************** -pmSourceMoments(source, radius): this function takes a subImage defined in the -pmSource data structure, along with the peak location, and determines the -various moments associated with that peak. - -Requires the following to have been created: - pmSource - pmSource->peak - pmSource->pixels - -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 -*****************************************************************************/ -# define VALID_RADIUS(X,Y) (((R2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) - -bool pmSourceMoments_EAM(pmSource *source, - psF32 radius) -{ - // 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); - - // - // XXX: Verify the setting for sky if source->moments == NULL. - // - psF32 sky = 0.0; - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } else { - sky = source->moments->Sky; - } - - // - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // X2 = SUM (x - xc)^2 * (z - sky) - // XY = SUM (x - xc)*(y - yc)*(z - sky) - // - psF32 peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; - psF32 Sum = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; - psF32 X2 = 0.0; - psF32 Y2 = 0.0; - psF32 XY = 0.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 - // 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. need to do two loops for a - // numerically stable result. first loop: get the sums. - // XXX EAM : mask == 0 is valid - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - for (psS32 col = 0; col < source->pixels->numCols ; col++) { - if ((source->mask != NULL) && (source->mask->data.U8[row][col])) continue; - - psF32 xDiff = col + source->pixels->col0 - xPeak; - psF32 yDiff = row + source->pixels->row0 - yPeak; - - if (!VALID_RADIUS(xDiff, yDiff)) continue; - - psF32 pDiff = source->pixels->data.F32[row][col] - sky; - - // XXX EAM : check for valid S/N in pixel - 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, - "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); - - // - // first moment X = X1/Sum + xc - // second moment X = sqrt (X2/Sum - (X1/Sum)^2) - // Sxy = XY / Sum - // - 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 + xPeak; - source->moments->y = y + yPeak; - - source->moments->Sxy = XY/Sum - x*y; - 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)); - - 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); -} - bool pmModelFitStatus (pmModel *model) { Index: /trunk/psphot/src/psModulesUtils.h =================================================================== --- /trunk/psphot/src/psModulesUtils.h (revision 5047) +++ /trunk/psphot/src/psModulesUtils.h (revision 5048) @@ -10,7 +10,5 @@ // psModule extra utilities -bool pmSourceLocalSky_EAM (pmSource *source, psStatsOptions statsOptions, psF32 Radius); -bool pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF); -bool pmSourceMoments_EAM(pmSource *source, psF32 radius); +// bool pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF); bool pmModelFitStatus (pmModel *model); int pmSourceDophotType (pmSource *source); Index: /trunk/psphot/src/psphotApplyPSF.c =================================================================== --- /trunk/psphot/src/psphotApplyPSF.c (revision 5047) +++ /trunk/psphot/src/psphotApplyPSF.c (revision 5048) @@ -39,6 +39,6 @@ // skip non-astronomical objects (very likely defects) // XXX EAM : should we try these anyway? - if (source->type == PS_SOURCE_DEFECT) continue; - if (source->type == PS_SOURCE_SATURATED) continue; + if (source->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; // use the source moments, etc to guess basic model parameters @@ -65,5 +65,5 @@ } - if (i > 66) psTraceSetLevel (".psLib.dataManip.psMinimizeLMChi2", 5); + // if (i > 66) psTraceSetLevel (".psLib.dataManip.psMinimizeLMChi2", 5); // fit PSF model (set/unset the pixel mask) @@ -73,5 +73,5 @@ if (!status || (model->params->data.F32[1] < 0)) { psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter); - source->type = PS_SOURCE_FAIL_FIT_PSF; // better choice? + source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? psFree (model); continue; Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 5047) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 5048) @@ -18,5 +18,5 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - if (source->type != PS_SOURCE_PSFSTAR) continue; + if (source->type != PM_SOURCE_PSFSTAR) continue; psArrayAdd (stars, 200, source); } Index: /trunk/psphot/src/psphotFitGalaxies.c =================================================================== --- /trunk/psphot/src/psphotFitGalaxies.c (revision 5047) +++ /trunk/psphot/src/psphotFitGalaxies.c (revision 5048) @@ -31,28 +31,26 @@ // sources which should not be fitted // skip all valid stars - if (source->type == PS_SOURCE_PSFSTAR) continue; - if (source->type == PS_SOURCE_SATSTAR) continue; - if (source->type == PS_SOURCE_GOODSTAR) continue; + if (source->type == PM_SOURCE_PSFSTAR) continue; + if (source->type == PM_SOURCE_SATSTAR) continue; + if (source->type == PM_SOURCE_GOODSTAR) continue; // skip all likely defects - if (source->type == PS_SOURCE_DEFECT) continue; - if (source->type == PS_SOURCE_SATURATED) continue; + if (source->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; // skip poorly fitted stars - if (source->type == PS_SOURCE_FAINTSTAR) continue; - if (source->type == PS_SOURCE_POOR_FIT_PSF) continue; + if (source->type == PM_SOURCE_FAINTSTAR) continue; + if (source->type == PM_SOURCE_POOR_FIT_PSF) continue; // XXX EAM when do we pick these up again? if (source->moments->SN < GAL_MIN_SN) { - source->type = PS_SOURCE_FAINT_GALAXY; // better choice? + source->type = PM_SOURCE_FAINT_GALAXY; // better choice? continue; } // recalculate the source moments using the larger galaxy moments radius - // XXX EAM : 0.5.0 code used _EAM version - // status = pmSourceMoments (source, GAL_MOMENTS_RAD); - status = pmSourceMoments_EAM (source, GAL_MOMENTS_RAD); + status = pmSourceMoments (source, GAL_MOMENTS_RAD); if (!status) { - source->type = PS_SOURCE_DROP_GALAXY; // better choice? + source->type = PM_SOURCE_DROP_GALAXY; // better choice? continue; } @@ -84,5 +82,5 @@ // if the fit fails, we need to change the classification psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); - source->type = PS_SOURCE_FAIL_FIT_GAL; // better choice? + source->type = PM_SOURCE_FAIL_FIT_GAL; // better choice? source->modelFLT = model; Nfail ++; @@ -95,5 +93,5 @@ // if the fit fails, we need to change the classification psLogMsg ("psphot", 3, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); - source->type = PS_SOURCE_POOR_FIT_GAL; // better choice? + source->type = PM_SOURCE_POOR_FIT_GAL; // better choice? source->modelFLT = model; Nfail ++; @@ -102,5 +100,5 @@ // accept the model - source->type = PS_SOURCE_GALAXY; + source->type = PM_SOURCE_GALAXY; source->modelFLT = model; Niter += model[0].nIter; Index: /trunk/psphot/src/psphotMarkPSF.c =================================================================== --- /trunk/psphot/src/psphotMarkPSF.c (revision 5047) +++ /trunk/psphot/src/psphotMarkPSF.c (revision 5048) @@ -18,5 +18,5 @@ // we also reject objects with S/N too low or ChiSquare to high -// any object which meets the criterion is marked as PS_SOURCE_BRIGHTSTAR +// any object which meets the criterion is marked as PM_SOURCE_BRIGHTSTAR // floor for DS value @@ -35,13 +35,13 @@ // remember: fit does not use saturated pixels (masked) if (source->modelPSF->params->data.F32[1] >= SATURATE) { - if (source->type == PS_SOURCE_PSFSTAR) { + if (source->type == PM_SOURCE_PSFSTAR) { psLogMsg ("psphot", 3, "PSFSTAR marked SATSTAR\n"); } - source->type = PS_SOURCE_SATSTAR; + source->type = PM_SOURCE_SATSTAR; return (true); } - if (source->type == PS_SOURCE_SATSTAR) { + if (source->type == PM_SOURCE_SATSTAR) { psLogMsg ("psphot", 4, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n"); - source->type = PS_SOURCE_GOODSTAR; + source->type = PM_SOURCE_GOODSTAR; } @@ -57,5 +57,5 @@ nSy = dSY / hypot (MIN_DS, 1 / (SY * SN)); - // assign PS_SOURCE_GOODSTAR to bright objects within PSF region of dparams[] + // assign PM_SOURCE_GOODSTAR to bright objects within PSF region of dparams[] keep = TRUE; keep &= (fabs(nSx) < shapeNsigma); @@ -64,10 +64,10 @@ keep &= (Chi < maxChi); if (keep) { - if (source->type == PS_SOURCE_PSFSTAR) return (true); - source->type = PS_SOURCE_GOODSTAR; + if (source->type == PM_SOURCE_PSFSTAR) return (true); + source->type = PM_SOURCE_GOODSTAR; return (true); } - if (source->type == PS_SOURCE_PSFSTAR) { + if (source->type == PM_SOURCE_PSFSTAR) { psLogMsg ("psphot", 3, "PSFSTAR demoted based on fit quality\n"); } @@ -75,5 +75,5 @@ // object appears to be small, suspected defect if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) { - source->type = PS_SOURCE_DEFECT; + source->type = PM_SOURCE_DEFECT; return (false); } @@ -81,5 +81,5 @@ // object appears to be large, suspected galaxy if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) { - source->type = PS_SOURCE_GALAXY; + source->type = PM_SOURCE_GALAXY; return (false); } @@ -87,10 +87,10 @@ // object appears to be extremely faint: what is this? if (SN < minSN) { - source->type = PS_SOURCE_FAINTSTAR; + source->type = PM_SOURCE_FAINTSTAR; return (false); } // these are pooly fitted, probable stars near other stars? - source->type = PS_SOURCE_POOR_FIT_PSF; + source->type = PM_SOURCE_POOR_FIT_PSF; return (false); } Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 5047) +++ /trunk/psphot/src/psphotOutput.c (revision 5048) @@ -91,22 +91,22 @@ // use PSF model with stars - case PS_SOURCE_PSFSTAR: - case PS_SOURCE_SATSTAR: - case PS_SOURCE_GOODSTAR: + case PM_SOURCE_PSFSTAR: + case PM_SOURCE_SATSTAR: + case PM_SOURCE_GOODSTAR: model = source->modelPSF; break; // use FLT model with galaxies - case PS_SOURCE_GALAXY: + case PM_SOURCE_GALAXY: model = source->modelFLT; break; // skip defects and poorly fitted stars or galaxies - case PS_SOURCE_DEFECT: - case PS_SOURCE_SATURATED: - case PS_SOURCE_FAINTSTAR: - case PS_SOURCE_POOR_FIT_PSF: - case PS_SOURCE_POOR_FIT_GAL: - case PS_SOURCE_FAIL_FIT_GAL: + case PM_SOURCE_DEFECT: + case PM_SOURCE_SATURATED: + case PM_SOURCE_FAINTSTAR: + case PM_SOURCE_POOR_FIT_PSF: + case PM_SOURCE_POOR_FIT_GAL: + case PM_SOURCE_FAIL_FIT_GAL: default: model = NULL; @@ -183,26 +183,26 @@ // use PSF model with stars - case PS_SOURCE_SATSTAR: + case PM_SOURCE_SATSTAR: flags = 4; model = source->modelPSF; break; - case PS_SOURCE_PSFSTAR: - case PS_SOURCE_GOODSTAR: + case PM_SOURCE_PSFSTAR: + case PM_SOURCE_GOODSTAR: model = source->modelPSF; break; // use FLT model with galaxies - case PS_SOURCE_GALAXY: + case PM_SOURCE_GALAXY: model = source->modelFLT; break; // skip defects and poorly fitted stars or galaxies - case PS_SOURCE_DEFECT: - case PS_SOURCE_SATURATED: - case PS_SOURCE_FAINTSTAR: - case PS_SOURCE_POOR_FIT_PSF: - case PS_SOURCE_POOR_FIT_GAL: - case PS_SOURCE_FAIL_FIT_GAL: + case PM_SOURCE_DEFECT: + case PM_SOURCE_SATURATED: + case PM_SOURCE_FAINTSTAR: + case PM_SOURCE_POOR_FIT_PSF: + case PM_SOURCE_POOR_FIT_GAL: + case PM_SOURCE_FAIL_FIT_GAL: default: model = NULL; @@ -288,22 +288,22 @@ // use PSF model with stars - case PS_SOURCE_PSFSTAR: - case PS_SOURCE_SATSTAR: - case PS_SOURCE_GOODSTAR: + case PM_SOURCE_PSFSTAR: + case PM_SOURCE_SATSTAR: + case PM_SOURCE_GOODSTAR: model = source->modelPSF; break; // use FLT model with galaxies - case PS_SOURCE_GALAXY: + case PM_SOURCE_GALAXY: model = source->modelFLT; break; // skip defects and poorly fitted stars or galaxies - case PS_SOURCE_DEFECT: - case PS_SOURCE_SATURATED: - case PS_SOURCE_FAINTSTAR: - case PS_SOURCE_POOR_FIT_PSF: - case PS_SOURCE_POOR_FIT_GAL: - case PS_SOURCE_FAIL_FIT_GAL: + case PM_SOURCE_DEFECT: + case PM_SOURCE_SATURATED: + case PM_SOURCE_FAINTSTAR: + case PM_SOURCE_POOR_FIT_PSF: + case PM_SOURCE_POOR_FIT_GAL: + case PM_SOURCE_FAIL_FIT_GAL: default: model = NULL; @@ -445,7 +445,7 @@ // valid source types for this function - if (source->type == PS_SOURCE_SATSTAR) goto valid; - if (source->type == PS_SOURCE_PSFSTAR) goto valid; - if (source->type == PS_SOURCE_GOODSTAR) goto valid; + if (source->type == PM_SOURCE_SATSTAR) goto valid; + if (source->type == PM_SOURCE_PSFSTAR) goto valid; + if (source->type == PM_SOURCE_GOODSTAR) goto valid; continue; @@ -509,5 +509,5 @@ model = (pmModel *) source->modelFLT; if (model == NULL) continue; - if (source->type == PS_SOURCE_GALAXY) goto valid; + if (source->type == PM_SOURCE_GALAXY) goto valid; continue; @@ -570,8 +570,8 @@ // skip these sources (in PSF or FLT) - if (source->type == PS_SOURCE_GALAXY) continue; - if (source->type == PS_SOURCE_PSFSTAR) continue; - if (source->type == PS_SOURCE_SATSTAR) continue; - if (source->type == PS_SOURCE_GOODSTAR) continue; + if (source->type == PM_SOURCE_GALAXY) continue; + if (source->type == PM_SOURCE_PSFSTAR) continue; + if (source->type == PM_SOURCE_SATSTAR) continue; + if (source->type == PM_SOURCE_GOODSTAR) continue; if (source->moments == NULL) { @@ -602,30 +602,30 @@ switch (source->type) { - case PS_SOURCE_DEFECT: - case PS_SOURCE_SATURATED: + case PM_SOURCE_DEFECT: + case PM_SOURCE_SATURATED: return (8); - case PS_SOURCE_SATSTAR: + case PM_SOURCE_SATSTAR: return (10); - case PS_SOURCE_PSFSTAR: - case PS_SOURCE_GOODSTAR: + case PM_SOURCE_PSFSTAR: + case PM_SOURCE_GOODSTAR: return (1); - case PS_SOURCE_POOR_FIT_PSF: + case PM_SOURCE_POOR_FIT_PSF: return (7); - case PS_SOURCE_FAIL_FIT_PSF: - case PS_SOURCE_FAINTSTAR: + case PM_SOURCE_FAIL_FIT_PSF: + case PM_SOURCE_FAINTSTAR: return (4); - case PS_SOURCE_GALAXY: - case PS_SOURCE_FAINT_GALAXY: - case PS_SOURCE_DROP_GALAXY: - case PS_SOURCE_FAIL_FIT_GAL: - case PS_SOURCE_POOR_FIT_GAL: + case PM_SOURCE_GALAXY: + case PM_SOURCE_FAINT_GALAXY: + case PM_SOURCE_DROP_GALAXY: + case PM_SOURCE_FAIL_FIT_GAL: + case PM_SOURCE_POOR_FIT_GAL: return (2); - case PS_SOURCE_OTHER: + case PM_SOURCE_OTHER: default: return (0); Index: /trunk/psphot/src/psphotSourceStats.c =================================================================== --- /trunk/psphot/src/psphotSourceStats.c (revision 5047) +++ /trunk/psphot/src/psphotSourceStats.c (revision 5048) @@ -28,5 +28,5 @@ // measure a local sky value // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken - status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER); + status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER); if (!status) { psFree (source); @@ -36,5 +36,5 @@ // measure basic source moments // XXX EAM : choose between these two versions - status = pmSourceMoments_EAM (source, RADIUS); + status = pmSourceMoments (source, RADIUS); if (!status) { psFree (source); Index: /trunk/psphot/src/psphotSubtractPSF.c =================================================================== --- /trunk/psphot/src/psphotSubtractPSF.c (revision 5047) +++ /trunk/psphot/src/psphotSubtractPSF.c (revision 5048) @@ -11,7 +11,7 @@ // only subtract successful fits - if (source->type == PS_SOURCE_SATSTAR) goto valid; - if (source->type == PS_SOURCE_PSFSTAR) goto valid; - if (source->type == PS_SOURCE_GOODSTAR) goto valid; + if (source->type == PM_SOURCE_SATSTAR) goto valid; + if (source->type == PM_SOURCE_PSFSTAR) goto valid; + if (source->type == PM_SOURCE_GOODSTAR) goto valid; return (false);