Index: /branches/eam_rel9_p0/psModules/src/objects/pmObjects.c =================================================================== --- /branches/eam_rel9_p0/psModules/src/objects/pmObjects.c (revision 6181) +++ /branches/eam_rel9_p0/psModules/src/objects/pmObjects.c (revision 6182) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.5.4.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-15 18:21:49 $ + * @version $Revision: 1.5.4.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-22 20:20:47 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -844,4 +844,6 @@ psF32 xPeak = source->peak->x; psF32 yPeak = source->peak->y; + psF32 xOff = source->pixels->col0 - source->peak->x; + psF32 yOff = source->pixels->row0 - source->peak->y; // XXX why do I get different results for these two methods of finding Sx? @@ -854,12 +856,19 @@ // XXX EAM : mask == 0 is valid + psF32 **vPix = source->pixels->data.F32; + psF32 **vWgt = source->weight->data.F32; + psU8 **vMsk = source->mask->data.U8; + 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])) { + if ((source->mask != NULL) && (vMsk[row][col])) { continue; } - psF32 xDiff = col + source->pixels->col0 - xPeak; - psF32 yDiff = row + source->pixels->row0 - yPeak; + // psF32 xDiff = col + source->pixels->col0 - xPeak; + // psF32 yDiff = row + source->pixels->row0 - yPeak; + + psF32 xDiff = col + xOff; + psF32 yDiff = row + yOff; // XXX EAM : calculate xDiff, yDiff up front; @@ -869,980 +878,976 @@ } - psF32 pDiff = source->pixels->data.F32[row][col] - sky; - psF32 wDiff = source->weight->data.F32[row][col]; + psF32 pDiff = vPix[row][col] - sky; + psF32 wDiff = vWgt[row][col]; // XXX EAM : check for valid S/N in pixel // XXX EAM : should this limit be user-defined? - if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) { + if (PS_SQR(pDiff) < wDiff) { + + Var += wDiff; + Sum += pDiff; + + psF32 xWght = xDiff * pDiff; + psF32 yWght = yDiff * pDiff; + + X1 += xWght; + Y1 += yWght; + XY += xDiff * yWght; + + X2 += xDiff * xWght; + Y2 += yDiff * yWght; + + peakPixel = PS_MAX (vPix[row][col], peakPixel); + numPixels++; + } + } + + // if we have less than (1/4) of the possible pixels, force a retry + // XXX EAM - the limit is a bit arbitrary. make it user defined? + if ((numPixels < 0.75*R2) || (Sum <= 0)) { + psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + 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", 3, + "large centroid swing; invalid peak %d, %d\n", + source->peak->x, source->peak->y); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + 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; + source->moments->SN = Sum / sqrt(Var); + 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); + + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); + } + + // XXX EAM : I used + int pmComparePeakAscend (const void **a, const void **b) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 diff; + + diff = A->counts - B->counts; + if (diff < FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); + return (-1); + } else if (diff > FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); + return (+1); + } + psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); + return (0); + } + + int pmComparePeakDescend (const void **a, const void **b) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 diff; + + diff = A->counts - B->counts; + if (diff < FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); + return (+1); + } else if (diff > FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); + return (-1); + } + psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); + return (0); + } + + /****************************************************************************** + pmSourcePSFClump(source, metadata): Find the likely PSF clump in the + sigma-x, sigma-y plane. return 0,0 clump in case of error. + *****************************************************************************/ + + // XXX EAM include a S/N cutoff in selecting the sources? + pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + # define NPIX 10 + # define SCALE 0.1 + + psArray *peaks = NULL; + pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; + pmPSFClump psfClump = emptyClump; + + PS_ASSERT_PTR_NON_NULL(sources, emptyClump); + PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); + + // find the sigmaX, sigmaY clump + { + psStats *stats = NULL; + psImage *splane = NULL; + int binX, binY; + bool status; + + psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX"); + if (!status) + SX_MAX = 10.0; + psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX"); + if (!status) + SY_MAX = 10.0; + + // construct a sigma-plane image + // psImageAlloc does zero the data + splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/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?) + for (psS32 i = 0 ; i < sources->n ; i++) + { + pmSource *tmpSrc = (pmSource *) sources->data[i]; + if (tmpSrc == NULL) { + continue; + } + if (tmpSrc->moments == NULL) { + continue; + } + + // Sx,Sy are limited at 0. a peak at 0,0 is artificial + if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { + continue; + } + + // for the moment, force splane dimensions to be 10x10 image pix + binX = tmpSrc->moments->Sx/SCALE; + if (binX < 0) + continue; + if (binX >= splane->numCols) + continue; + + binY = tmpSrc->moments->Sy/SCALE; + if (binY < 0) + continue; + if (binY >= splane->numRows) + continue; + + splane->data.F32[binY][binX] += 1.0; + } + + // find the peak in this image + stats = psStatsAlloc (PS_STAT_MAX); + stats = psImageStats (stats, splane, NULL, 0); + peaks = pmFindImagePeaks (splane, stats[0].max / 2); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); + + 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 + { + pmPeak *clump; + psF32 minSx, maxSx; + psF32 minSy, maxSy; + psVector *tmpSx = NULL; + psVector *tmpSy = NULL; + psStats *stats = NULL; + + // XXX EAM : this lets us takes the single highest peak + psArraySort (peaks, pmComparePeakDescend); + clump = peaks->data[0]; + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); + + // define section window for clump + minSx = clump->x * SCALE - 0.2; + maxSx = clump->x * SCALE + 0.2; + minSy = clump->y * SCALE - 0.2; + maxSy = clump->y * SCALE + 0.2; + + tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); + tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); + tmpSx->n = 0; + tmpSy->n = 0; + + // XXX clip sources based on flux? + // create vectors with Sx, Sy values in window + for (psS32 i = 0 ; i < sources->n ; i++) + { + pmSource *tmpSrc = (pmSource *) sources->data[i]; + + if (tmpSrc->moments->Sx < minSx) + continue; + if (tmpSrc->moments->Sx > maxSx) + continue; + if (tmpSrc->moments->Sy < minSy) + continue; + if (tmpSrc->moments->Sy > maxSy) + continue; + tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; + tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; + tmpSx->n++; + tmpSy->n++; + if (tmpSx->n == tmpSx->nalloc) { + psVectorRealloc (tmpSx, tmpSx->nalloc + 100); + psVectorRealloc (tmpSy, tmpSy->nalloc + 100); + } + } + + // measures stats of Sx, Sy + stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + + stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); + psfClump.X = stats->clippedMean; + psfClump.dX = stats->clippedStdev; + + stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); + 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 + + psFree (stats); + psFree (peaks); + psFree (tmpSx); + psFree (tmpSy); + } + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + 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? + + EAM: I moved S/N calculation to pmSourceMoments, using weight image + *****************************************************************************/ + + bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + psBool rc = true; + + int Nsat = 0; + int Next = 0; + int Nstar = 0; + int Npsf = 0; + int Ncr = 0; + int Nsatstar = 0; + // psRegion allArray = psRegionSet (0, 0, 0, 0); + psRegion inner; + + // report stats on S/N values for star-like objects + psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); + starsn->n = 0; + + // check return status value (do these exist?) + bool status; + psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); + + // XXX allow clump size to be scaled relative to sigmas? + // make rough IDs based on clumpX,Y,DX,DY + for (psS32 i = 0 ; i < sources->n ; i++) { + + pmSource *tmpSrc = (pmSource *) sources->data[i]; + + tmpSrc->peak->class = 0; + + psF32 sigX = tmpSrc->moments->Sx; + psF32 sigY = tmpSrc->moments->Sy; + + // XXX EAM : can we use the value of SATURATE if mask is NULL? + inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2); + int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED); + + // 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)); + big = true; + if ((Nsatpix > 1) && big) { + tmpSrc->type = PM_SOURCE_STAR; + tmpSrc->mode = PM_SOURCE_SATSTAR; + Nsatstar ++; continue; } - Var += wDiff; - 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++; - } - } - - // if we have less than (1/4) of the possible pixels, force a retry - // XXX EAM - the limit is a bit arbitrary. make it user defined? - if ((numPixels < 0.75*R2) || (Sum <= 0)) { - psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n"); - psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); - 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", 3, - "large centroid swing; invalid peak %d, %d\n", - source->peak->x, source->peak->y); - psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); - 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; - source->moments->SN = Sum / sqrt(Var); - 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); - - psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); - return(true); -} - -// XXX EAM : I used -int pmComparePeakAscend (const void **a, const void **b) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - pmPeak *A = *(pmPeak **)a; - pmPeak *B = *(pmPeak **)b; - - psF32 diff; - - diff = A->counts - B->counts; - if (diff < FLT_EPSILON) { - psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); - return (-1); - } else if (diff > FLT_EPSILON) { - psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); - return (+1); - } - psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); - return (0); -} - -int pmComparePeakDescend (const void **a, const void **b) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - pmPeak *A = *(pmPeak **)a; - pmPeak *B = *(pmPeak **)b; - - psF32 diff; - - diff = A->counts - B->counts; - if (diff < FLT_EPSILON) { - psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); - return (+1); - } else if (diff > FLT_EPSILON) { - psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); - return (-1); - } - psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); - return (0); -} - -/****************************************************************************** -pmSourcePSFClump(source, metadata): Find the likely PSF clump in the -sigma-x, sigma-y plane. return 0,0 clump in case of error. -*****************************************************************************/ - -// XXX EAM include a S/N cutoff in selecting the sources? -pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - - # define NPIX 10 - # define SCALE 0.1 - - psArray *peaks = NULL; - pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; - pmPSFClump psfClump = emptyClump; - - PS_ASSERT_PTR_NON_NULL(sources, emptyClump); - PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); - - // find the sigmaX, sigmaY clump + // saturated object (not a star, eg bleed trails, hot pixels) + if (Nsatpix > 1) { + tmpSrc->type = PM_SOURCE_SATURATED; + tmpSrc->mode = PM_SOURCE_DEFAULT; + Nsat ++; + continue; + } + + // 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 = PM_SOURCE_DEFECT; + tmpSrc->mode = PM_SOURCE_DEFAULT; + Ncr ++; + continue; + } + + // likely unsaturated extended source (too large to be stellar) + if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { + tmpSrc->type = PM_SOURCE_EXTENDED; + tmpSrc->mode = PM_SOURCE_DEFAULT; + Next ++; + continue; + } + + // the rest are probable stellar objects + starsn->data.F32[starsn->n] = tmpSrc->moments->SN; + starsn->n ++; + Nstar ++; + + // PSF star (within 1.5 sigma of clump center, S/N > limit) + psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); + if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) { + tmpSrc->type = PM_SOURCE_STAR; + tmpSrc->mode = PM_SOURCE_PSFSTAR; + Npsf ++; + continue; + } + + // random type of star + tmpSrc->type = PM_SOURCE_STAR; + tmpSrc->mode = PM_SOURCE_DEFAULT; + } + + { + psStats *stats = NULL; + stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); + stats = psVectorStats (stats, starsn, NULL, NULL, 0); + psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); + psFree (stats); + psFree (starsn); + } + + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Next: %3d\n", Next); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); + + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); + } + + /** pmSourceDefinePixels() + * + * Define psImage subarrays for the source located at coordinates x,y on the + * image set defined by readout. The pixels defined by this operation consist of + * a square window (of full width 2Radius+1) centered on the pixel which contains + * the given coordinate, in the frame of the readout. The window is defined to + * have limits which are valid within the boundary of the readout image, thus if + * the radius would fall outside the image pixels, the subimage is truncated to + * only consist of valid pixels. If readout->mask or readout->weight are not + * NULL, matching subimages are defined for those images as well. This function + * fails if no valid pixels can be defined (x or y less than Radius, for + * example). This function should be used to define a region of interest around a + * source, including both source and sky pixels. + * + * XXX: must code this. + * + */ + bool pmSourceDefinePixels( + pmSource *mySource, ///< Add comment. + pmReadout *readout, ///< Add comment. + psF32 x, ///< Add comment. + psF32 y, ///< Add comment. + psF32 Radius) ///< Add comment. { - psStats *stats = NULL; - psImage *splane = NULL; - int binX, binY; - bool status; - - psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX"); - if (!status) - SX_MAX = 10.0; - psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX"); - if (!status) - SY_MAX = 10.0; - - // construct a sigma-plane image - // psImageAlloc does zero the data - splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/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?) - for (psS32 i = 0 ; i < sources->n ; i++) - { - pmSource *tmpSrc = (pmSource *) sources->data[i]; - if (tmpSrc == NULL) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented. Returning FALSE.\n"); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(false); + } + + /****************************************************************************** + pmSourceSetPixelsCircle(source, image, radius) + + XXX: This was replaced by DefinePixels in SDRS. Remove it. + *****************************************************************************/ + bool pmSourceSetPixelsCircle(pmSource *source, + const psImage *image, + psF32 radius) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_FLOAT_COMPARE(0.0, radius, false); + + // + // 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). + // + psS32 radiusS32 = (psS32) radius; + psS32 SubImageCenterRow = source->peak->y; + 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, SubImageCenterRow + radiusS32 + 1); + psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); + psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); + + // XXX: Must recycle image. + // XXX EAM: this message reflects a programming error we know about. + // i am setting it to a trace message which we can take out + if (source->pixels != NULL) { + psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4, + "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL. Freeing and reallocating.\n"); + psFree(source->pixels); + } + source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, + SubImageStartRow, + SubImageEndCol, + SubImageEndRow)); + + // XXX: Must recycle image. + if (source->mask != NULL) { + psFree(source->mask); + } + source->mask = psImageAlloc(source->pixels->numCols, + source->pixels->numRows, + PS_TYPE_U8); // XXX EAM : type was F32 + + // + // Loop through the subimage mask, initialize mask to 0 or 1. + // XXX EAM: valid pixels should have 0, not 1 + for (psS32 row = 0 ; row < source->mask->numRows; row++) { + for (psS32 col = 0 ; col < source->mask->numCols; col++) { + + if (checkRadius2((psF32) radiusS32, + (psF32) radiusS32, + radius, + (psF32) col, + (psF32) row)) { + source->mask->data.U8[row][col] = 0; + } else { + source->mask->data.U8[row][col] = 1; + } + } + } + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); + } + + /****************************************************************************** + 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 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. + *****************************************************************************/ + pmModel *pmSourceModelGuess(pmSource *source, + pmModelType modelType) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + + pmModel *model = pmModelAlloc(modelType); + + pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); + modelGuessFunc(model, source); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(model); + } + + /****************************************************************************** + evalModel(source, level, row): a private function which evaluates the + source->modelPSF function at the specified coords. The coords are subImage, not + image coords. + + NOTE: The coords are in subImage source->pixel coords, not image coords. + + XXX: reverse order of row,col args? + + XXX: rename all coords in this file such that their name defines whether + the coords is in subImage or image space. + + XXX: This should probably be a public pmModules function. + + XXX: Use static vectors for x. + + XXX: Figure out if it's (row, col) or (col, row) for the model functions. + + XXX: For a while, the first psVectorAlloc() was generating a seg fault during + testing. Try to reproduce that and debug. + *****************************************************************************/ + + // XXX EAM : I have made this a public function + // XXX EAM : this now uses a pmModel as the input + // XXX EAM : it was using src->type to find the model, not model->type + psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(image, false); + PS_ASSERT_PTR_NON_NULL(model, false); + PS_ASSERT_PTR_NON_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 + image->col0); + x->data.F32[1] = (psF32) (row + image->row0); + psF32 tmpF; + pmModelFunc modelFunc; + + modelFunc = pmModelFunc_GetFunction (model->type); + tmpF = modelFunc (NULL, model->params, x); + psFree(x); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpF); + } + + /****************************************************************************** + pmSourceContour(src, img, level, mode): For an input subImage, and model, this + routine returns a psArray of coordinates that evaluate to the specified level. + + XXX: Probably should remove the "image" argument. + XXX: What type should the output coordinate vectors consist of? col,row? + XXX: Why a pmArray output? + XXX: doex x,y correspond with col,row or row/col? + XXX: What is mode? + XXX: The top, bottom of the contour is not correctly determined. + XXX EAM : this function is using the model for the contour, but it should + be using only the image counts + *****************************************************************************/ + psArray *pmSourceContour(pmSource *source, + const psImage *image, + psF32 level, + pmContourType mode) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(image, 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->modelEXT, false); + // XXX EAM : what is the purpose of modelPSF/modelEXT? + + // + // Allocate data for x/y pairs. + // + psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); + psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); + + // + // Start at the row with peak pixel, then decrement. + // + psS32 col = source->peak->x; + for (psS32 row = source->peak->y; row>= 0 ; row--) { + // XXX: yVec contain no real information. Do we really need it? + yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); + yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); + + // Starting at peak pixel, search leftwards for the column intercept. + psF32 leftIntercept = findValue(source, level, row, col, 0); + if (isnan(leftIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; + + // Starting at peak pixel, search rightwards for the column intercept. + + psF32 rightIntercept = findValue(source, level, row, col, 1); + if (isnan(rightIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); + xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; + + // Set starting column for next row + col = (psS32) ((leftIntercept + rightIntercept) / 2.0); + } + // + // Start at the row (+1) with peak pixel, then increment. + // + col = source->peak->x; + for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) { + // XXX: yVec contain no real information. Do we really need it? + yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); + yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); + + // Starting at peak pixel, search leftwards for the column intercept. + psF32 leftIntercept = findValue(source, level, row, col, 0); + if (isnan(leftIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; + + // Starting at peak pixel, search rightwards for the column intercept. + psF32 rightIntercept = findValue(source, level, row, col, 1); + if (isnan(rightIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; + + // Set starting column for next row + col = (psS32) ((leftIntercept + rightIntercept) / 2.0); + } + + // + // Allocate an array for result, store coord vectors there. + // + psArray *tmpArray = psArrayAlloc(2); + tmpArray->data[0] = (psPtr *) yVec; + tmpArray->data[1] = (psPtr *) xVec; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpArray); + } + + // save a static values so they may be set externally + static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; + static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; + + bool pmSourceFitModelInit (float nIter, float tol) { + + PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter; + PM_SOURCE_FIT_MODEL_TOLERANCE = tol; + } + + bool pmSourceFitModel (pmSource *source, + pmModel *model, + const bool PSF) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + 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->weight, false); + + // XXX EAM : is it necessary for the mask & weight to exist? the + // tests below could be conditions (!NULL) + + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + psVector *params = model->params; + psVector *dparams = model->dparams; + psVector *paramMask = NULL; + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + + int nParams = PSF ? 4 : params->n; + + // maximum number of valid pixels + psS32 nPix = source->pixels->numRows * source->pixels->numCols; + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(nPix); + psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); + + nPix = 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]) { + continue; + } + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + // Convert i/j to image space: + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[nPix] = (psPtr *) coord; + y->data.F32[nPix] = source->pixels->data.F32[i][j]; + + // psMinimizeLMChi2 takes wt = 1/dY^2 + if (source->weight->data.F32[i][j] == 0) { + continue; + } + yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + nPix++; + } + } + if (nPix < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + model->status = PM_MODEL_BADARGS; + psFree (x); + psFree (y); + psFree (yErr); + return(false); + } + x->n = nPix; + y->n = nPix; + yErr->n = nPix; + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + // PSF model only fits first 4 parameters, EXT 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; + } + } + + // Set the parameter range checks + 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(myMin, covar, params, paramMask, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if ((paramMask != NULL) && paramMask->data.U8[i]) continue; - } - if (tmpSrc->moments == NULL) { - continue; - } - - // Sx,Sy are limited at 0. a peak at 0,0 is artificial - if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { - continue; - } - - // for the moment, force splane dimensions to be 10x10 image pix - binX = tmpSrc->moments->Sx/SCALE; - if (binX < 0) - continue; - if (binX >= splane->numCols) - continue; - - binY = tmpSrc->moments->Sy/SCALE; - if (binY < 0) - continue; - if (binY >= splane->numRows) - continue; - - splane->data.F32[binY][binX] += 1.0; - } - - // find the peak in this image - stats = psStatsAlloc (PS_STAT_MAX); - stats = psImageStats (stats, splane, NULL, 0); - peaks = pmFindImagePeaks (splane, stats[0].max / 2); - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); - - 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 - { - pmPeak *clump; - psF32 minSx, maxSx; - psF32 minSy, maxSy; - psVector *tmpSx = NULL; - psVector *tmpSy = NULL; - psStats *stats = NULL; - - // XXX EAM : this lets us takes the single highest peak - psArraySort (peaks, pmComparePeakDescend); - clump = peaks->data[0]; - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); - - // define section window for clump - minSx = clump->x * SCALE - 0.2; - maxSx = clump->x * SCALE + 0.2; - minSy = clump->y * SCALE - 0.2; - maxSy = clump->y * SCALE + 0.2; - - tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); - tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); - tmpSx->n = 0; - tmpSy->n = 0; - - // XXX clip sources based on flux? - // create vectors with Sx, Sy values in window - for (psS32 i = 0 ; i < sources->n ; i++) - { - pmSource *tmpSrc = (pmSource *) sources->data[i]; - - if (tmpSrc->moments->Sx < minSx) - continue; - if (tmpSrc->moments->Sx > maxSx) - continue; - if (tmpSrc->moments->Sy < minSy) - continue; - if (tmpSrc->moments->Sy > maxSy) - continue; - tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; - tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; - tmpSx->n++; - tmpSy->n++; - if (tmpSx->n == tmpSx->nalloc) { - psVectorRealloc (tmpSx, tmpSx->nalloc + 100); - psVectorRealloc (tmpSy, tmpSy->nalloc + 100); - } - } - - // measures stats of Sx, Sy - stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - - stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); - psfClump.X = stats->clippedMean; - psfClump.dX = stats->clippedStdev; - - stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); - 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 - - psFree (stats); - psFree (peaks); - psFree (tmpSx); - psFree (tmpSy); - } - - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - 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? - -EAM: I moved S/N calculation to pmSourceMoments, using weight image -*****************************************************************************/ - -bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - - psBool rc = true; - - int Nsat = 0; - int Next = 0; - int Nstar = 0; - int Npsf = 0; - int Ncr = 0; - int Nsatstar = 0; - // psRegion allArray = psRegionSet (0, 0, 0, 0); - psRegion inner; - - // report stats on S/N values for star-like objects - psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); - starsn->n = 0; - - // check return status value (do these exist?) - bool status; - psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); - - // XXX allow clump size to be scaled relative to sigmas? - // make rough IDs based on clumpX,Y,DX,DY - for (psS32 i = 0 ; i < sources->n ; i++) { - - pmSource *tmpSrc = (pmSource *) sources->data[i]; - - tmpSrc->peak->class = 0; - - psF32 sigX = tmpSrc->moments->Sx; - psF32 sigY = tmpSrc->moments->Sy; - - // XXX EAM : can we use the value of SATURATE if mask is NULL? - inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2); - int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED); - - // 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)); - big = true; - if ((Nsatpix > 1) && big) { - tmpSrc->type = PM_SOURCE_STAR; - tmpSrc->mode = PM_SOURCE_SATSTAR; - Nsatstar ++; - continue; - } - - // saturated object (not a star, eg bleed trails, hot pixels) - if (Nsatpix > 1) { - tmpSrc->type = PM_SOURCE_SATURATED; - tmpSrc->mode = PM_SOURCE_DEFAULT; - Nsat ++; - continue; - } - - // 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 = PM_SOURCE_DEFECT; - tmpSrc->mode = PM_SOURCE_DEFAULT; - Ncr ++; - continue; - } - - // likely unsaturated extended source (too large to be stellar) - if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { - tmpSrc->type = PM_SOURCE_EXTENDED; - tmpSrc->mode = PM_SOURCE_DEFAULT; - Next ++; - continue; - } - - // the rest are probable stellar objects - starsn->data.F32[starsn->n] = tmpSrc->moments->SN; - starsn->n ++; - Nstar ++; - - // PSF star (within 1.5 sigma of clump center, S/N > limit) - psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); - if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) { - tmpSrc->type = PM_SOURCE_STAR; - tmpSrc->mode = PM_SOURCE_PSFSTAR; - Npsf ++; - continue; - } - - // random type of star - tmpSrc->type = PM_SOURCE_STAR; - tmpSrc->mode = PM_SOURCE_DEFAULT; - } - - { - psStats *stats = NULL; - stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); - stats = psVectorStats (stats, starsn, NULL, NULL, 0); - psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); - psFree (stats); - psFree (starsn); - } - - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Next: %3d\n", Next); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); - - psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); - return(rc); -} - -/** pmSourceDefinePixels() - * - * Define psImage subarrays for the source located at coordinates x,y on the - * image set defined by readout. The pixels defined by this operation consist of - * a square window (of full width 2Radius+1) centered on the pixel which contains - * the given coordinate, in the frame of the readout. The window is defined to - * have limits which are valid within the boundary of the readout image, thus if - * the radius would fall outside the image pixels, the subimage is truncated to - * only consist of valid pixels. If readout->mask or readout->weight are not - * NULL, matching subimages are defined for those images as well. This function - * fails if no valid pixels can be defined (x or y less than Radius, for - * example). This function should be used to define a region of interest around a - * source, including both source and sky pixels. - * - * XXX: must code this. - * - */ -bool pmSourceDefinePixels( - pmSource *mySource, ///< Add comment. - pmReadout *readout, ///< Add comment. - psF32 x, ///< Add comment. - psF32 y, ///< Add comment. - psF32 Radius) ///< Add comment. -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented. Returning FALSE.\n"); - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(false); -} - -/****************************************************************************** -pmSourceSetPixelsCircle(source, image, radius) - -XXX: This was replaced by DefinePixels in SDRS. Remove it. -*****************************************************************************/ -bool pmSourceSetPixelsCircle(pmSource *source, - const psImage *image, - psF32 radius) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_FLOAT_COMPARE(0.0, radius, false); - - // - // 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). - // - psS32 radiusS32 = (psS32) radius; - psS32 SubImageCenterRow = source->peak->y; - 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, SubImageCenterRow + radiusS32 + 1); - psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); - psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); - - // XXX: Must recycle image. - // XXX EAM: this message reflects a programming error we know about. - // i am setting it to a trace message which we can take out - if (source->pixels != NULL) { - psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4, - "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL. Freeing and reallocating.\n"); - psFree(source->pixels); - } - source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, - SubImageStartRow, - SubImageEndCol, - SubImageEndRow)); - - // XXX: Must recycle image. - if (source->mask != NULL) { - psFree(source->mask); - } - source->mask = psImageAlloc(source->pixels->numCols, - source->pixels->numRows, - PS_TYPE_U8); // XXX EAM : type was F32 - - // - // Loop through the subimage mask, initialize mask to 0 or 1. - // XXX EAM: valid pixels should have 0, not 1 - for (psS32 row = 0 ; row < source->mask->numRows; row++) { - for (psS32 col = 0 ; col < source->mask->numCols; col++) { - - if (checkRadius2((psF32) radiusS32, - (psF32) radiusS32, - radius, - (psF32) col, - (psF32) row)) { - source->mask->data.U8[row][col] = 0; - } else { - source->mask->data.U8[row][col] = 1; - } - } - } - psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); - return(true); -} - -/****************************************************************************** -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 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. -*****************************************************************************/ -pmModel *pmSourceModelGuess(pmSource *source, - pmModelType modelType) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - - pmModel *model = pmModelAlloc(modelType); - - pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); - modelGuessFunc(model, source); - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(model); -} - -/****************************************************************************** -evalModel(source, level, row): a private function which evaluates the -source->modelPSF function at the specified coords. The coords are subImage, not -image coords. - -NOTE: The coords are in subImage source->pixel coords, not image coords. - -XXX: reverse order of row,col args? - -XXX: rename all coords in this file such that their name defines whether -the coords is in subImage or image space. - -XXX: This should probably be a public pmModules function. - -XXX: Use static vectors for x. - -XXX: Figure out if it's (row, col) or (col, row) for the model functions. - -XXX: For a while, the first psVectorAlloc() was generating a seg fault during -testing. Try to reproduce that and debug. -*****************************************************************************/ - -// XXX EAM : I have made this a public function -// XXX EAM : this now uses a pmModel as the input -// XXX EAM : it was using src->type to find the model, not model->type -psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(image, false); - PS_ASSERT_PTR_NON_NULL(model, false); - PS_ASSERT_PTR_NON_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 + image->col0); - x->data.F32[1] = (psF32) (row + image->row0); - psF32 tmpF; - pmModelFunc modelFunc; - - modelFunc = pmModelFunc_GetFunction (model->type); - tmpF = modelFunc (NULL, model->params, x); - psFree(x); - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmpF); -} - -/****************************************************************************** -pmSourceContour(src, img, level, mode): For an input subImage, and model, this -routine returns a psArray of coordinates that evaluate to the specified level. - -XXX: Probably should remove the "image" argument. -XXX: What type should the output coordinate vectors consist of? col,row? -XXX: Why a pmArray output? -XXX: doex x,y correspond with col,row or row/col? -XXX: What is mode? -XXX: The top, bottom of the contour is not correctly determined. -XXX EAM : this function is using the model for the contour, but it should - be using only the image counts -*****************************************************************************/ -psArray *pmSourceContour(pmSource *source, - const psImage *image, - psF32 level, - pmContourType mode) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(image, 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->modelEXT, false); - // XXX EAM : what is the purpose of modelPSF/modelEXT? - - // - // Allocate data for x/y pairs. - // - psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); - psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); - - // - // Start at the row with peak pixel, then decrement. - // - psS32 col = source->peak->x; - for (psS32 row = source->peak->y; row>= 0 ; row--) { - // XXX: yVec contain no real information. Do we really need it? - yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); - yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); - - // Starting at peak pixel, search leftwards for the column intercept. - psF32 leftIntercept = findValue(source, level, row, col, 0); - if (isnan(leftIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; - - // Starting at peak pixel, search rightwards for the column intercept. - - psF32 rightIntercept = findValue(source, level, row, col, 1); - if (isnan(rightIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); - xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; - - // Set starting column for next row - col = (psS32) ((leftIntercept + rightIntercept) / 2.0); - } - // - // Start at the row (+1) with peak pixel, then increment. - // - col = source->peak->x; - for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) { - // XXX: yVec contain no real information. Do we really need it? - yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); - yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); - - // Starting at peak pixel, search leftwards for the column intercept. - psF32 leftIntercept = findValue(source, level, row, col, 0); - if (isnan(leftIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; - - // Starting at peak pixel, search rightwards for the column intercept. - psF32 rightIntercept = findValue(source, level, row, col, 1); - if (isnan(rightIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; - - // Set starting column for next row - col = (psS32) ((leftIntercept + rightIntercept) / 2.0); - } - - // - // Allocate an array for result, store coord vectors there. - // - psArray *tmpArray = psArrayAlloc(2); - tmpArray->data[0] = (psPtr *) yVec; - tmpArray->data[1] = (psPtr *) xVec; - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmpArray); -} - -// XXX EAM : these are better starting values, but should be available from metadata? -#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 -#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 -bool pmSourceFitModel (pmSource *source, - pmModel *model, - const bool PSF) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - 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->weight, false); - - // XXX EAM : is it necessary for the mask & weight to exist? the - // tests below could be conditions (!NULL) - - psBool fitStatus = true; - psBool onPic = true; - psBool rc = true; - - psVector *params = model->params; - psVector *dparams = model->dparams; - psVector *paramMask = NULL; - - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - - int nParams = PSF ? 4 : params->n; - - // maximum number of valid pixels - psS32 nPix = source->pixels->numRows * source->pixels->numCols; - - // construct the coordinate and value entries - psArray *x = psArrayAlloc(nPix); - psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); - psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); - - nPix = 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]) { - continue; - } - psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - - // Convert i/j to image space: - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); - x->data[nPix] = (psPtr *) coord; - y->data.F32[nPix] = source->pixels->data.F32[i][j]; - - // psMinimizeLMChi2 takes wt = 1/dY^2 - if (source->weight->data.F32[i][j] == 0) { - continue; - } - yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; - nPix++; - } - } - if (nPix < nParams + 1) { - psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); - psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); - model->status = PM_MODEL_BADARGS; - psFree (x); - psFree (y); - psFree (yErr); - return(false); - } - x->n = nPix; - y->n = nPix; - yErr->n = nPix; - - psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, - PM_SOURCE_FIT_MODEL_TOLERANCE); - - // PSF model only fits first 4 parameters, EXT 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; - } - } - - // Set the parameter range checks - 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(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]); - } - - // save the resulting chisq, nDOF, nIter - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nDOF = y->n - nParams; - - // 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]; - } - } - - // set the model success or failure status - if (!fitStatus) { - model->status = PM_MODEL_NONCONVERGE; - } else { - model->status = PM_MODEL_SUCCESS; - } - - // 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); - if (!onPic) { - model->status = PM_MODEL_OFFIMAGE; - } - - source->mode |= PM_SOURCE_FITTED; - - psFree(x); - psFree(y); - psFree(yErr); - psFree(myMin); - psFree(covar); - psFree(paramMask); - - rc = (onPic && fitStatus); - psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); - return(rc); -} - -bool p_pmSourceAddOrSubModel(psImage *image, - psImage *mask, - pmModel *model, - bool center, - bool sky, - bool add - ) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - - PS_ASSERT_PTR_NON_NULL(model, false); - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); - - psVector *x = psVectorAlloc(2, PS_TYPE_F32); - psVector *params = model->params; - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - psS32 imageCol; - psS32 imageRow; - psF32 skyValue = params->data.F32[0]; - psF32 pixelValue; - - 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; - - // XXX: Should you be adding the pixels for the entire subImage, - // or a radius of pixels around it? - - // Convert i/j to imace coord space: - // XXX: Make sure you have col/row order correct. - // 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; - - // set the appropriate pixel value for this coordinate - if (sky) { - pixelValue = modelFunc (NULL, params, x); - } else { - pixelValue = modelFunc (NULL, params, x) - skyValue; - } - - - // add or subtract the value - if (add - ) { - image->data.F32[i][j] += pixelValue; - } - else { - image->data.F32[i][j] -= pixelValue; - } - } - } - psFree(x); - psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); - return(true); -} - - - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceAddModel(psImage *image, - psImage *mask, - pmModel *model, - bool center, - bool sky) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true); - psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); - return(rc); -} - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceSubModel(psImage *image, - psImage *mask, - pmModel *model, - bool center, - bool sky) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false); - psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); - return(rc); -} - -bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) -{ - - float obsSum = 0; - float fitSum = 0; - float sky = model->params->data.F32[0]; - - pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); - fitSum = modelFluxFunc (model->params); - - for (int ix = 0; ix < image->numCols; ix++) { - for (int iy = 0; iy < image->numRows; iy++) { - if (mask->data.U8[iy][ix]) - continue; - obsSum += image->data.F32[iy][ix] - sky; - } - } - if (obsSum <= 0) - return false; - if (fitSum <= 0) - return false; - - *fitMag = -2.5*log10(fitSum); - *obsMag = -2.5*log10(obsSum); - return (true); -} - + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // save the resulting chisq, nDOF, nIter + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // 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]; + } + } + + // set the model success or failure status + if (!fitStatus) { + model->status = PM_MODEL_NONCONVERGE; + } else { + model->status = PM_MODEL_SUCCESS; + } + + // 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); + if (!onPic) { + model->status = PM_MODEL_OFFIMAGE; + } + + source->mode |= PM_SOURCE_FITTED; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(paramMask); + + rc = (onPic && fitStatus); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); + } + + bool p_pmSourceAddOrSubModel(psImage *image, + psImage *mask, + pmModel *model, + bool center, + bool sky, + bool add + ) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + PS_ASSERT_PTR_NON_NULL(model, false); + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); + + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + psVector *params = model->params; + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + psS32 imageCol; + psS32 imageRow; + psF32 skyValue = params->data.F32[0]; + psF32 pixelValue; + + 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; + + // XXX: Should you be adding the pixels for the entire subImage, + // or a radius of pixels around it? + + // Convert i/j to imace coord space: + // XXX: Make sure you have col/row order correct. + // 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; + + // set the appropriate pixel value for this coordinate + if (sky) { + pixelValue = modelFunc (NULL, params, x); + } else { + pixelValue = modelFunc (NULL, params, x) - skyValue; + } + + + // add or subtract the value + if (add + ) { + image->data.F32[i][j] += pixelValue; + } + else { + image->data.F32[i][j] -= pixelValue; + } + } + } + psFree(x); + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); + } + + + + /****************************************************************************** + *****************************************************************************/ + bool pmSourceAddModel(psImage *image, + psImage *mask, + pmModel *model, + bool center, + bool sky) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); + } + + /****************************************************************************** + *****************************************************************************/ + bool pmSourceSubModel(psImage *image, + psImage *mask, + pmModel *model, + bool center, + bool sky) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); + } + + bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) { + + float obsSum = 0; + float fitSum = 0; + float sky = model->params->data.F32[0]; + + pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); + fitSum = modelFluxFunc (model->params); + + for (int ix = 0; ix < image->numCols; ix++) { + for (int iy = 0; iy < image->numRows; iy++) { + if (mask->data.U8[iy][ix]) + continue; + obsSum += image->data.F32[iy][ix] - sky; + } + } + if (obsSum <= 0) + return false; + if (fitSum <= 0) + return false; + + *fitMag = -2.5*log10(fitSum); + *obsMag = -2.5*log10(obsSum); + return (true); + } + Index: /branches/eam_rel9_p0/psModules/src/objects/pmObjects.h =================================================================== --- /branches/eam_rel9_p0/psModules/src/objects/pmObjects.h (revision 6181) +++ /branches/eam_rel9_p0/psModules/src/objects/pmObjects.h (revision 6182) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.4.4.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-20 02:38:28 $ + * @version $Revision: 1.4.4.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-22 20:20:47 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -475,4 +475,9 @@ +bool pmSourceFitModelInit( + float nIter; ///< max number of allowed iterations + float tol ///< convergence criterion +); + /** pmSourceFitModel() *