Index: /branches/eam_rel9_p0/psModules/src/objects/pmObjects.c =================================================================== --- /branches/eam_rel9_p0/psModules/src/objects/pmObjects.c (revision 6182) +++ /branches/eam_rel9_p0/psModules/src/objects/pmObjects.c (revision 6183) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.5.4.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-22 20:20:47 $ + * @version $Revision: 1.5.4.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-22 21:12:19 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -884,274 +884,279 @@ // XXX EAM : should this limit be user-defined? 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); - } - - /****************************************************************************** + continue; + } + + 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 +*****************************************************************************/ + +// 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++) { - 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 + 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++) { - 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); - } - - /****************************************************************************** + 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. @@ -1162,224 +1167,226 @@ 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 +*****************************************************************************/ + +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; - } - - { - 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. + 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; + } + { - 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); - } - - /****************************************************************************** + 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); - } - - /****************************************************************************** +*****************************************************************************/ +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. @@ -1390,20 +1397,21 @@ 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); - } - - /****************************************************************************** +*****************************************************************************/ +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 @@ -1425,31 +1433,32 @@ 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); - } - - /****************************************************************************** +*****************************************************************************/ + +// 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. @@ -1463,391 +1472,399 @@ 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); +*****************************************************************************/ +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; + return true; +} + +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 != NULL) && paramMask->data.U8[i]) + if (!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] = 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 6182) +++ /branches/eam_rel9_p0/psModules/src/objects/pmObjects.h (revision 6183) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.4.4.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-22 20:20:47 $ + * @version $Revision: 1.4.4.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-22 21:12:20 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -476,6 +476,6 @@ bool pmSourceFitModelInit( - float nIter; ///< max number of allowed iterations - float tol ///< convergence criterion + float nIter, ///< max number of allowed iterations + float tol ///< convergence criterion );