Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 25755) +++ trunk/psphot/src/psphot.h (revision 25852) @@ -182,5 +182,5 @@ bool psphotVisualShowFootprints (pmDetections *detections); bool psphotVisualShowMoments (psArray *sources); -bool psphotVisualPlotMoments (psMetadata *recipe, psArray *sources); +bool psphotVisualPlotMoments (psMetadata *recipe, psMetadata *analysis, psArray *sources); bool psphotVisualShowRoughClass (psArray *sources); bool psphotVisualShowPSFStars (psMetadata *recipe, pmPSF *psf, psArray *sources); @@ -193,5 +193,5 @@ bool psphotVisualShowResidualImage (pmReadout *readout); bool psphotVisualPlotApResid (psArray *sources, float mean, float error); -bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources); +bool psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources); bool psphotVisualShowPetrosians (psArray *sources); @@ -212,9 +212,9 @@ // bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef); // bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian); -// bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, -// psVector *refRadius, psVector *meanSB, -// psVector *petRatio, psVector *petRatioErr, psVector *fluxSum, -// float petRadius, float ratioForRadius, -// float petFlux, float radiusForFlux); +// bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, +// psVector *refRadius, psVector *meanSB, +// psVector *petRatio, psVector *petRatioErr, psVector *fluxSum, +// float petRadius, float ratioForRadius, +// float petFlux, float radiusForFlux); bool psphotImageQuality (psMetadata *recipe, psArray *sources); Index: trunk/psphot/src/psphotRoughClass.c =================================================================== --- trunk/psphot/src/psphotRoughClass.c (revision 25755) +++ trunk/psphot/src/psphotRoughClass.c (revision 25852) @@ -1,11 +1,11 @@ # include "psphotInternal.h" -# define CHECK_STATUS(S,MSG) { \ - if (!status) { \ - psError(PSPHOT_ERR_CONFIG, false, "missing PSF Clump entry: %s\n", MSG); \ - return false; \ - } } +# define CHECK_STATUS(S,MSG) { \ + if (!status) { \ + psError(PSPHOT_ERR_CONFIG, false, "missing PSF Clump entry: %s\n", MSG); \ + return false; \ + } } -bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *recipe, const bool havePSF); +bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *target, psMetadata *recipe, const bool havePSF); // 2006.02.02 : no leaks @@ -24,16 +24,16 @@ int nRegion = 0; for (int ix = 0; ix < NX; ix ++) { - for (int iy = 0; iy < NY; iy ++) { + for (int iy = 0; iy < NY; iy ++) { - psRegion *region = psRegionAlloc (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY); - if (!psphotRoughClassRegion (nRegion, region, sources, recipe, havePSF)) { - psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n", - region->x0, region->y0, region->x1, region->y1); - psFree (region); - continue; - } - psFree (region); - nRegion ++; - } + psRegion *region = psRegionAlloc (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY); + if (!psphotRoughClassRegion (nRegion, region, sources, readout->analysis, recipe, havePSF)) { + psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n", + region->x0, region->y0, region->x1, region->y1); + psFree (region); + continue; + } + psFree (region); + nRegion ++; + } } psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", nRegion); @@ -44,5 +44,5 @@ psLogMsg ("psphot.roughclass", PS_LOG_INFO, "rough classification: %f sec\n", psTimerMark ("psphot.rough")); - psphotVisualPlotMoments (recipe, sources); + psphotVisualPlotMoments (recipe, readout->analysis, sources); psphotVisualShowRoughClass (sources); // XXX better visualization: psphotVisualShowFlags (sources); @@ -51,5 +51,5 @@ } -bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *recipe, const bool havePSF) { +bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *target, psMetadata *recipe, const bool havePSF) { bool status; @@ -62,40 +62,40 @@ snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", nRegion); - psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + psMetadata *regionMD = psMetadataLookupPtr (&status, target, regionName); if (!regionMD) { - // allocate the region metadata folder and add this region to it. - regionMD = psMetadataAlloc(); - psMetadataAddMetadata (recipe, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); - psFree (regionMD); + // allocate the region metadata folder and add this region to it. + regionMD = psMetadataAlloc(); + psMetadataAddMetadata (target, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); + psFree (regionMD); } psMetadataAddPtr (regionMD, PS_LIST_TAIL, "REGION", PS_DATA_REGION | PS_META_REPLACE, "psf clump region", region); if (!havePSF) { - // determine the PSF parameters from the source moment values - // XXX why not save the psfClump as a PTR? - psfClump = pmSourcePSFClump (region, sources, recipe); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); + // determine the PSF parameters from the source moment values + // XXX why not save the psfClump as a PTR? + psfClump = pmSourcePSFClump (region, sources, target); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); } else { - // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y - psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); - if (!status) { - psLogMsg ("psphot", 4, "No PSF clump defined for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); - return false; - } - psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); psAssert (status, "missing PSF.CLUMP.Y"); - psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); psAssert (status, "missing PSF.CLUMP.DX"); - psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); psAssert (status, "missing PSF.CLUMP.DY"); + // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y + psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + if (!status) { + psLogMsg ("psphot", 4, "No PSF clump defined for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + return false; + } + psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); psAssert (status, "missing PSF.CLUMP.Y"); + psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); psAssert (status, "missing PSF.CLUMP.DX"); + psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); psAssert (status, "missing PSF.CLUMP.DY"); } if (psfClump.X < 0) { - psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourcePSFClump"); - return false; + psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourcePSFClump"); + return false; } if (!psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) { - psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); - return false; + psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + return false; } psLogMsg ("psphot", 3, "psf clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); @@ -104,6 +104,6 @@ // group into STAR, COSMIC, EXTENDED, SATURATED, etc. if (!pmSourceRoughClass (region, sources, recipe, psfClump, maskSat)) { - psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourceRoughClass"); - return false; + psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourceRoughClass"); + return false; } Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 25755) +++ trunk/psphot/src/psphotSourceSize.c (revision 25852) @@ -83,5 +83,5 @@ psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size")); - psphotVisualPlotSourceSize (recipe, sources); + psphotVisualPlotSourceSize (recipe, readout->analysis, sources); psphotVisualShowSourceSize (readout, sources); psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr); @@ -103,26 +103,26 @@ pmFootprint *footprint = peak->footprint; if (!footprint) { - // if we have not footprint, use the old code to mask by isophot - psphotMaskCosmicRayIsophot (source, maskVal, crMask); - return true; + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + return true; } if (!footprint->spans) { - // if we have no footprint, use the old code to mask by isophot - psphotMaskCosmicRayIsophot (source, maskVal, crMask); - return true; + // if we have no footprint, use the old code to mask by isophot + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + return true; } // mask all of the pixels covered by the spans of the footprint for (int j = 1; j < footprint->spans->n; j++) { - pmSpan *span1 = footprint->spans->data[j]; - - int iy = span1->y; - int xs = span1->x0; - int xe = span1->x1; - - for (int ix = xs; ix < xe; ix++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } + pmSpan *span1 = footprint->spans->data[j]; + + int iy = span1->y; + int xs = span1->x0; + int xe = span1->x1; + + for (int ix = xs; ix < xe; ix++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } } return true; @@ -147,14 +147,14 @@ // mark the pixels in this row to the left, then the right for (int ix = xo; ix >= 0; ix--) { - float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; - } + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } } for (int ix = xo + 1; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; - } + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } } @@ -162,33 +162,33 @@ // first go up: for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } } // next go down: for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } } return true; @@ -201,5 +201,5 @@ psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32); psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32); - + psImageMaskType maskVal = options->maskVal | options->markVal; @@ -208,6 +208,6 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + pmSource *source = sources->data[i]; + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; // replace object in image @@ -216,22 +216,22 @@ } - // clear the mask bit and set the circular mask pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); - - // XXX can we test if psfMag is set and calculate only if needed? - pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal - - // clear the mask bit - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); // re-subtract the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); - float apMag = -2.5*log10(source->moments->Sum); - float dMag = source->psfMag - apMag; - - psVectorAppend (Ap, 100, dMag); - psVectorAppend (ApErr, 100, source->errMag); + float apMag = -2.5*log10(source->moments->Sum); + float dMag = source->psfMag - apMag; + + psVectorAppend (Ap, 100, dMag); + psVectorAppend (ApErr, 100, source->errMag); } @@ -242,5 +242,5 @@ psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32); for (int i = 0; i < Ap->n; i++) { - dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian; + dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian; } @@ -266,29 +266,29 @@ psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip"); - int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS"); for (int i = 0; i < nRegions; i ++) { - snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); - psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); - psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion"); - - psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION"); - psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion"); - - // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y - psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); psAssert (status, "missing PSF.CLUMP.X"); - psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); psAssert (status, "missing PSF.CLUMP.Y"); - psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); psAssert (status, "missing PSF.CLUMP.DX"); - psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); psAssert (status, "missing PSF.CLUMP.DY"); - - if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) { - psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); - continue; - } - - if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) { - psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); - continue; - } - } + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); + psMetadata *regionMD = psMetadataLookupPtr (&status, readout->analysis, regionName); + psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion"); + + psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION"); + psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion"); + + // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y + psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); psAssert (status, "missing PSF.CLUMP.X"); + psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); psAssert (status, "missing PSF.CLUMP.Y"); + psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); psAssert (status, "missing PSF.CLUMP.DX"); + psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); psAssert (status, "missing PSF.CLUMP.DY"); + + if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) { + psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + continue; + } + + if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) { + psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + continue; + } + } return true; @@ -314,35 +314,35 @@ for (psS32 i = 0 ; i < sources->n ; i++) { - pmSource *source = (pmSource *) sources->data[i]; - - // psfClumps are found for image subregions: - // skip sources not in this region - if (source->peak->x < region->x0) continue; - if (source->peak->x >= region->x1) continue; - if (source->peak->y < region->y0) continue; - if (source->peak->y >= region->y1) continue; - - // skip source if it was already measured - if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { - psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); - continue; - } - - // source must have been subtracted - if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { - source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; - psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); - continue; - } - - // we are basically classifying by moments; use the default if not found - psAssert (source->moments, "why is this source missing moments?"); - if (source->mode & noMoments) { - Nskip ++; - continue; - } - - psF32 Mxx = source->moments->Mxx; - psF32 Myy = source->moments->Myy; + pmSource *source = (pmSource *) sources->data[i]; + + // psfClumps are found for image subregions: + // skip sources not in this region + if (source->peak->x < region->x0) continue; + if (source->peak->x >= region->x1) continue; + if (source->peak->y < region->y0) continue; + if (source->peak->y >= region->y1) continue; + + // skip source if it was already measured + if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { + psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); + continue; + } + + // source must have been subtracted + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { + source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; + psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); + continue; + } + + // we are basically classifying by moments; use the default if not found + psAssert (source->moments, "why is this source missing moments?"); + if (source->mode & noMoments) { + Nskip ++; + continue; + } + + psF32 Mxx = source->moments->Mxx; + psF32 Myy = source->moments->Myy; // replace object in image @@ -351,58 +351,58 @@ } - // clear the mask bit and set the circular mask pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); - - // XXX can we test if psfMag is set and calculate only if needed? - pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal - - // clear the mask bit - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); // re-subtract the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); - float apMag = -2.5*log10(source->moments->Sum); - float dMag = source->psfMag - apMag; - float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr); - - source->extNsigma = nSigma; - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; - - // Anything within this region is a probably PSF-like object. Saturated stars may land - // in this region, but are detected elsewhere on the basis of their peak value. - bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY); - if (isPSF) { - Npsf ++; - continue; - } - - // Defects may not always match CRs from peak curvature analysis - // Defects may also be marked as SATSTAR -- XXX deactivate this flag? - // XXX this rule is not great - if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) { - source->mode |= PM_SOURCE_MODE_DEFECT; - Ncr ++; - continue; - } - - // saturated star (determined in PSF fit). These may also be saturated galaxies, or - // just large saturated regions. - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - Nsat ++; - continue; - } - - // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)? - bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y)); - if (isEXT) { - source->mode |= PM_SOURCE_MODE_EXT_LIMIT; - Next ++; - continue; - } - - psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma); - Nmiss ++; + float apMag = -2.5*log10(source->moments->Sum); + float dMag = source->psfMag - apMag; + float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr); + + source->extNsigma = nSigma; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + + // Anything within this region is a probably PSF-like object. Saturated stars may land + // in this region, but are detected elsewhere on the basis of their peak value. + bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY); + if (isPSF) { + Npsf ++; + continue; + } + + // Defects may not always match CRs from peak curvature analysis + // Defects may also be marked as SATSTAR -- XXX deactivate this flag? + // XXX this rule is not great + if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) { + source->mode |= PM_SOURCE_MODE_DEFECT; + Ncr ++; + continue; + } + + // saturated star (determined in PSF fit). These may also be saturated galaxies, or + // just large saturated regions. + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + Nsat ++; + continue; + } + + // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)? + bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y)); + if (isEXT) { + source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + Next ++; + continue; + } + + psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma); + Nmiss ++; } @@ -417,5 +417,5 @@ // no longer used by psphot. float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, - psImageMaskType maskVal, const pmModel *model, float Ro) + psImageMaskType maskVal, const pmModel *model, float Ro) { psF32 *PAR = model->params->data.F32; // Model parameters @@ -432,6 +432,6 @@ float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0); if (Q < 0.0) { - // ellipse is imaginary - return NAN; + // ellipse is imaginary + return NAN; } @@ -441,44 +441,44 @@ for (int x = -radius; x <= radius; x++) { - // Polynomial coefficients - // XXX Should we be using the centre of the pixel as x or x+0.5? - float A = PS_SQR (1.0 / syy); - float B = x * sxy; - float C = PS_SQR (x / sxx) - Ro; - float T = PS_SQR(B) - 4*A*C; - if (T < 0.0) { - continue; - } - - // y position in source frame - float yP = (-B + sqrt (T)) / (2.0 * A); - float yM = (-B - sqrt (T)) / (2.0 * A); - - // Get the closest pixel positions (image frame) - int xPix = x + PAR[PM_PAR_XPOS] - image->col0 + 0.5; - int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5; - int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5; - - if (xPix < 0 || xPix >= image->numCols) { - continue; - } - - if (yPixM >= 0 && yPixM < image->numRows && - !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) { - float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]); - nSigma += dSigma; - nPts++; - } - - if (yPixM == yPixP) { - continue; - } - - if (yPixP >= 0 && yPixP < image->numRows && - !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) { - float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]); - nSigma += dSigma; - nPts++; - } + // Polynomial coefficients + // XXX Should we be using the centre of the pixel as x or x+0.5? + float A = PS_SQR (1.0 / syy); + float B = x * sxy; + float C = PS_SQR (x / sxx) - Ro; + float T = PS_SQR(B) - 4*A*C; + if (T < 0.0) { + continue; + } + + // y position in source frame + float yP = (-B + sqrt (T)) / (2.0 * A); + float yM = (-B - sqrt (T)) / (2.0 * A); + + // Get the closest pixel positions (image frame) + int xPix = x + PAR[PM_PAR_XPOS] - image->col0 + 0.5; + int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5; + int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5; + + if (xPix < 0 || xPix >= image->numCols) { + continue; + } + + if (yPixM >= 0 && yPixM < image->numRows && + !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]); + nSigma += dSigma; + nPts++; + } + + if (yPixM == yPixP) { + continue; + } + + if (yPixP >= 0 && yPixP < image->numRows && + !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]); + nSigma += dSigma; + nPts++; + } } nSigma /= nPts; @@ -491,128 +491,128 @@ // XXX use an internal flag to mark sources which have already been measured for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // skip source if it was already measured - if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { - psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); - continue; - } - - // source must have been subtracted - if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { - source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; - psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); - continue; - } - - psF32 **resid = source->pixels->data.F32; - psF32 **variance = source->variance->data.F32; - psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; - - // Integer position of peak - int xPeak = source->peak->xf - source->pixels->col0 + 0.5; - int yPeak = source->peak->yf - source->pixels->row0 + 0.5; - - // Skip sources which are too close to a boundary. These are mostly caught as DEFECT - if (xPeak < 1 || xPeak > source->pixels->numCols - 2 || - yPeak < 1 || yPeak > source->pixels->numRows - 2) { - psTrace("psphot", 7, "Not calculating crNsigma due to edge\n"); - continue; - } - - // Skip sources with masked pixels. These are mostly caught as DEFECT - bool keep = true; - for (int iy = -1; (iy <= +1) && keep; iy++) { - for (int ix = -1; (ix <= +1) && keep; ix++) { - if (mask[yPeak+iy][xPeak+ix] & options->maskVal) { - keep = false; - } - } - } - if (!keep) { - psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); - continue; - } - - // Compare the central pixel with those on either side, for the four possible lines through it. - - // Soften variances (add systematic error) - float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances - - // Across the middle: y = 0 - float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1]; - float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1]; - float nX = cX / sqrtf(dcX + softening); - - // Up the centre: x = 0 - float cY = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak+0] - resid[yPeak+1][xPeak+0]; - float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0]; - float nY = cY / sqrtf(dcY + softening); - - // Diagonal: x = y - float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; - float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1]; - float nL = cL / sqrtf(dcL + softening); - - // Diagonal: x = - y - float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; - float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1]; - float nR = cR / sqrtf(dcR + softening); - - // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) - // Ndof = 4 ? (four measurements, no free parameters) - // XXX this value is going to be biased low because of systematic errors. - // we need to calibrate it somehow - // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq); - - // not strictly accurate: overcounts the chisq contribution from the center pixel (by - // factor of 4); also biases a bit low if any pixels are masked - // XXX I am not sure I want to keep this value... - source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR); - - float fCR = 0.0; - int nCR = 0; - if (nX > 0.0) { - fCR += nX; - nCR ++; - } - if (nY > 0.0) { - fCR += nY; - nCR ++; - } - if (nL > 0.0) { - fCR += nL; - nCR ++; - } - if (nR > 0.0) { - fCR += nR; - nCR ++; - } - source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; - - if (!isfinite(source->crNsigma)) { - continue; - } - - // this source is thought to be a cosmic ray. flag the detection and mask the pixels - if (source->crNsigma > options->nSigmaCR) { - source->mode |= PM_SOURCE_MODE_CR_LIMIT; - // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask); - // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask); - } + pmSource *source = sources->data[i]; + + // skip source if it was already measured + if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { + psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); + continue; + } + + // source must have been subtracted + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { + source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; + psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); + continue; + } + + psF32 **resid = source->pixels->data.F32; + psF32 **variance = source->variance->data.F32; + psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + + // Integer position of peak + int xPeak = source->peak->xf - source->pixels->col0 + 0.5; + int yPeak = source->peak->yf - source->pixels->row0 + 0.5; + + // Skip sources which are too close to a boundary. These are mostly caught as DEFECT + if (xPeak < 1 || xPeak > source->pixels->numCols - 2 || + yPeak < 1 || yPeak > source->pixels->numRows - 2) { + psTrace("psphot", 7, "Not calculating crNsigma due to edge\n"); + continue; + } + + // Skip sources with masked pixels. These are mostly caught as DEFECT + bool keep = true; + for (int iy = -1; (iy <= +1) && keep; iy++) { + for (int ix = -1; (ix <= +1) && keep; ix++) { + if (mask[yPeak+iy][xPeak+ix] & options->maskVal) { + keep = false; + } + } + } + if (!keep) { + psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); + continue; + } + + // Compare the central pixel with those on either side, for the four possible lines through it. + + // Soften variances (add systematic error) + float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances + + // Across the middle: y = 0 + float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1]; + float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1]; + float nX = cX / sqrtf(dcX + softening); + + // Up the centre: x = 0 + float cY = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak+0] - resid[yPeak+1][xPeak+0]; + float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0]; + float nY = cY / sqrtf(dcY + softening); + + // Diagonal: x = y + float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; + float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1]; + float nL = cL / sqrtf(dcL + softening); + + // Diagonal: x = - y + float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; + float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1]; + float nR = cR / sqrtf(dcR + softening); + + // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) + // Ndof = 4 ? (four measurements, no free parameters) + // XXX this value is going to be biased low because of systematic errors. + // we need to calibrate it somehow + // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq); + + // not strictly accurate: overcounts the chisq contribution from the center pixel (by + // factor of 4); also biases a bit low if any pixels are masked + // XXX I am not sure I want to keep this value... + source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR); + + float fCR = 0.0; + int nCR = 0; + if (nX > 0.0) { + fCR += nX; + nCR ++; + } + if (nY > 0.0) { + fCR += nY; + nCR ++; + } + if (nL > 0.0) { + fCR += nL; + nCR ++; + } + if (nR > 0.0) { + fCR += nR; + nCR ++; + } + source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + + if (!isfinite(source->crNsigma)) { + continue; + } + + // this source is thought to be a cosmic ray. flag the detection and mask the pixels + if (source->crNsigma > options->nSigmaCR) { + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask); + // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask); + } } // now that we have masked pixels associated with CRs, we can grow the mask if (options->grow > 0) { - bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask - psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow); - psImageConvolveSetThreads(oldThreads); - if (!newMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask"); - return false; - } - psFree(readout->mask); - readout->mask = newMask; + bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask + psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow); + psImageConvolveSetThreads(oldThreads); + if (!newMask) { + psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask"); + return false; + } + psFree(readout->mask); + readout->mask = newMask; } return true; Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 25755) +++ trunk/psphot/src/psphotSourceStats.c (revision 25852) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources); +bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources); psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow) { @@ -68,8 +68,8 @@ if (setWindow) { - if (!psphotSetMomentsWindow(recipe, sources)) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!"); - return NULL; - } + if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!"); + return NULL; + } } @@ -242,5 +242,5 @@ // this function attempts to iteratively determine the best value for sigma of the moments weighting Gaussian -bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources) { +bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources) { bool status; @@ -259,44 +259,44 @@ for (int i = 0; i < 4; i++) { - // XXX move max source number to config - for (int j = 0; (j < sources->n) && (j < 400); j++) { - - pmSource *source = sources->data[j]; - psAssert (source->moments, "force moments to exist"); - source->moments->nPixels = 0; - - // skip faint sources for moments measurement - if (source->peak->SN < MIN_SN) { - continue; - } - - // measure basic source moments (no S/N clipping on input pixels) - status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0); - } - - // choose a grid scale that is a fixed fraction of the psf sigma^2 - psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i])); - psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i])); - psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i])); - - // determine the PSF parameters from the source moment values - pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe); - psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); - - psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); - psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000"); - if (!regionMD) { - regionMD = psMetadataAlloc(); - psMetadataAddMetadata (recipe, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD); - psFree (regionMD); - } - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); - psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); - - // psphotVisualPlotMoments (recipe, sources); - - Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i]; + // XXX move max source number to config + for (int j = 0; (j < sources->n) && (j < 400); j++) { + + pmSource *source = sources->data[j]; + psAssert (source->moments, "force moments to exist"); + source->moments->nPixels = 0; + + // skip faint sources for moments measurement + if (source->peak->SN < MIN_SN) { + continue; + } + + // measure basic source moments (no S/N clipping on input pixels) + status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0); + } + + // choose a grid scale that is a fixed fraction of the psf sigma^2 + psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i])); + psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i])); + psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i])); + + // determine the PSF parameters from the source moment values + pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe); + psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); + + psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, "PSF.CLUMP.REGION.000"); + if (!regionMD) { + regionMD = psMetadataAlloc(); + psMetadataAddMetadata (analysis, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD); + psFree (regionMD); + } + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); + + // psphotVisualPlotMoments (recipe, sources); + + Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i]; } @@ -307,6 +307,6 @@ float maxS = Sout[0]; for (int i = 0; i < 4; i++) { - minS = PS_MIN(Sout[i], minS); - maxS = PS_MAX(Sout[i], maxS); + minS = PS_MIN(Sout[i], minS); + maxS = PS_MAX(Sout[i], maxS); } if (minS > 0.65) Sigma = sigma[3]; @@ -314,10 +314,10 @@ for (int i = 0; i < 3; i++) { - if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue; - if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue; - Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]); + if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue; + if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue; + Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]); } psAssert (isfinite(Sigma), "did we miss a case?"); - + // choose a grid scale that is a fixed fraction of the psf sigma^2 psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma)); Index: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 25755) +++ trunk/psphot/src/psphotVisual.c (revision 25852) @@ -28,32 +28,32 @@ switch (channel) { case 1: - if (kapa1 == -1) { - kapa1 = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa1 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - } - } - return kapa1; + if (kapa1 == -1) { + kapa1 = KapaOpenNamedSocket ("kapa", "psphot:images"); + if (kapa1 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + } + } + return kapa1; case 2: - if (kapa2 == -1) { - kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - } - } - return kapa2; + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + } + } + return kapa2; case 3: - if (kapa3 == -1) { - kapa3 = KapaOpenNamedSocket ("kapa", "psphot:stamps"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - } - } - return kapa3; + if (kapa3 == -1) { + kapa3 = KapaOpenNamedSocket ("kapa", "psphot:stamps"); + if (kapa3 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + } + } + return kapa3; default: - psAbort ("unknown kapa channel"); + psAbort ("unknown kapa channel"); } psAbort ("unknown kapa channel"); @@ -315,5 +315,5 @@ // draw the top - // XXX need to allow top (and bottom) to have more than one span + // XXX need to allow top (and bottom) to have more than one span span = footprint->spans->data[0]; overlay[Noverlay].type = KII_OVERLAY_LINE; @@ -448,5 +448,5 @@ } -bool psphotVisualPlotMoments (psMetadata *recipe, psArray *sources) { +bool psphotVisualPlotMoments (psMetadata *recipe, psMetadata *analysis, psArray *sources) { bool status; @@ -469,29 +469,29 @@ float Ymin = 1000.0, Ymax = 0.0; { - int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); for (int n = 0; n < nRegions; n++) { char regionName[64]; snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); - psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); - - float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); - float X0 = psfX - 4.0*psfdX; - float X1 = psfX + 4.0*psfdX; - float Y0 = psfY - 4.0*psfdY; - float Y1 = psfY + 4.0*psfdY; - - if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } - if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } - if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } - if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } - } - } - Xmin = PS_MAX(Xmin, -0.1); - Ymin = PS_MAX(Ymin, -0.1); + float X0 = psfX - 4.0*psfdX; + float X1 = psfX + 4.0*psfdX; + float Y0 = psfY - 4.0*psfdY; + float Y1 = psfY + 4.0*psfdY; + + if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } + if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } + if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } + if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } + } + } + Xmin = PS_MAX(Xmin, -0.1); + Ymin = PS_MAX(Ymin, -0.1); // storage vectors for data to be plotted @@ -653,15 +653,15 @@ psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); - int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); - float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); + float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, analysis, "PSF_CLUMP_NSIGMA"); graphdata.color = KapaColorByName ("blue"); graphdata.style = 0; - graphdata.xmin = Xmin; - graphdata.ymin = Ymin; - graphdata.xmax = Xmax; - graphdata.ymax = Ymax; - KapaSetLimits (myKapa, &graphdata); + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); for (int n = 0; n < nRegions; n++) { @@ -669,5 +669,5 @@ char regionName[64]; snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); - psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); @@ -926,12 +926,12 @@ if (Xo == 0) { // place source alone on this row - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Yo += DY; @@ -943,12 +943,12 @@ Xo = 0; - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Xo = DX; @@ -957,11 +957,11 @@ } else { // extend this row - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Xo += DX; @@ -1020,5 +1020,5 @@ pmSource *source = sources->data[i]; - // only show "real" saturated stars (not defects) + // only show "real" saturated stars (not defects) if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;; if (source->mode & PM_SOURCE_MODE_DEFECT) continue;; @@ -1065,8 +1065,8 @@ pmSource *source = sources->data[i]; - // only show "real" saturated stars (not defects) + // only show "real" saturated stars (not defects) if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;; if (source->mode & PM_SOURCE_MODE_DEFECT) continue;; - nSAT ++; + nSAT ++; if (Xo + DX > NX) { @@ -1074,8 +1074,8 @@ if (Xo == 0) { // place source alone on this row - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); Yo += DY; @@ -1087,8 +1087,8 @@ Xo = 0; - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); Xo = DX; @@ -1097,8 +1097,8 @@ } else { // extend this row - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); Xo += DX; @@ -1406,11 +1406,11 @@ // if (source->type != type) continue; - if (mode) { - if (keep) { - if (!(source->mode & mode)) continue; - } else { - if (source->mode & mode) continue; - } - } + if (mode) { + if (keep) { + if (!(source->mode & mode)) continue; + } else { + if (source->mode & mode) continue; + } + } pmMoments *moments = source->moments; @@ -1463,5 +1463,5 @@ } -bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources) { +bool psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources) { bool status; @@ -1482,29 +1482,29 @@ float Ymin = 1000.0, Ymax = 0.0; { - int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); for (int n = 0; n < nRegions; n++) { char regionName[64]; snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); - psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); - - float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); - float X0 = psfX - 10.0*psfdX; - float X1 = psfX + 10.0*psfdX; - float Y0 = psfY - 10.0*psfdY; - float Y1 = psfY + 10.0*psfdY; - - if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } - if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } - if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } - if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } - } - } - Xmin = PS_MAX(Xmin, -0.1); - Ymin = PS_MAX(Ymin, -0.1); + float X0 = psfX - 10.0*psfdX; + float X1 = psfX + 10.0*psfdX; + float Y0 = psfY - 10.0*psfdY; + float Y1 = psfY + 10.0*psfdY; + + if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } + if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } + if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } + if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } + } + } + Xmin = PS_MAX(Xmin, -0.1); + Ymin = PS_MAX(Ymin, -0.1); // storage vectors for data to be plotted @@ -1544,42 +1544,42 @@ if (source->moments == NULL) continue; - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { - xCR->data.F32[nCR] = source->moments->Mxx; - yCR->data.F32[nCR] = source->moments->Myy; - mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum); - sCR->data.F32[nCR] = source->extNsigma; - nCR++; - } - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - xSAT->data.F32[nSAT] = source->moments->Mxx; - ySAT->data.F32[nSAT] = source->moments->Myy; - mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum); - sSAT->data.F32[nSAT] = source->extNsigma; - nSAT++; - } - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - xEXT->data.F32[nEXT] = source->moments->Mxx; - yEXT->data.F32[nEXT] = source->moments->Myy; - mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum); - sEXT->data.F32[nEXT] = source->extNsigma; - nEXT++; - continue; - } - if (source->mode & PM_SOURCE_MODE_DEFECT) { - xDEF->data.F32[nDEF] = source->moments->Mxx; - yDEF->data.F32[nDEF] = source->moments->Myy; - mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum); - sDEF->data.F32[nDEF] = source->extNsigma; - nDEF++; - continue; - } - if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) { - continue; - } - xPSF->data.F32[nPSF] = source->moments->Mxx; - yPSF->data.F32[nPSF] = source->moments->Myy; - mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum); - sPSF->data.F32[nPSF] = source->extNsigma; - nPSF++; + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { + xCR->data.F32[nCR] = source->moments->Mxx; + yCR->data.F32[nCR] = source->moments->Myy; + mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum); + sCR->data.F32[nCR] = source->extNsigma; + nCR++; + } + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + xSAT->data.F32[nSAT] = source->moments->Mxx; + ySAT->data.F32[nSAT] = source->moments->Myy; + mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum); + sSAT->data.F32[nSAT] = source->extNsigma; + nSAT++; + } + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + xEXT->data.F32[nEXT] = source->moments->Mxx; + yEXT->data.F32[nEXT] = source->moments->Myy; + mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum); + sEXT->data.F32[nEXT] = source->extNsigma; + nEXT++; + continue; + } + if (source->mode & PM_SOURCE_MODE_DEFECT) { + xDEF->data.F32[nDEF] = source->moments->Mxx; + yDEF->data.F32[nDEF] = source->moments->Myy; + mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum); + sDEF->data.F32[nDEF] = source->extNsigma; + nDEF++; + continue; + } + if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) { + continue; + } + xPSF->data.F32[nPSF] = source->moments->Mxx; + yPSF->data.F32[nPSF] = source->moments->Myy; + mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum); + sPSF->data.F32[nPSF] = source->extNsigma; + nPSF++; } xSAT->n = nSAT; @@ -1861,15 +1861,15 @@ psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); - int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); - float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); + float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, analysis, "PSF_CLUMP_NSIGMA"); graphdata.color = KapaColorByName ("blue"); graphdata.style = 0; - graphdata.xmin = Xmin; - graphdata.ymin = Ymin; - graphdata.xmax = Xmax; - graphdata.ymax = Ymax; - KapaSetLimits (myKapa, &graphdata); + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); for (int n = 0; n < nRegions; n++) { @@ -1877,5 +1877,5 @@ char regionName[64]; snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); - psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); @@ -2069,32 +2069,32 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - if (!source) continue; - if (!source->extpars) continue; - if (!source->extpars->profile) continue; - if (!source->extpars->petrosian_80) continue; - - pmSourceRadialProfile *profile = source->extpars->profile; - pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80; - - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = source->peak->xf; - overlay[Noverlay].y = source->peak->yf; - overlay[Noverlay].dx = 2.0*petrosian->radius; - overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major; - overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD; - overlay[Noverlay].text = NULL; - Noverlay ++; + pmSource *source = sources->data[i]; + + if (!source) continue; + if (!source->extpars) continue; + if (!source->extpars->profile) continue; + if (!source->extpars->petrosian_80) continue; + + pmSourceRadialProfile *profile = source->extpars->profile; + pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = source->peak->xf; + overlay[Noverlay].y = source->peak->yf; + overlay[Noverlay].dx = 2.0*petrosian->radius; + overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major; + overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = source->peak->xf; - overlay[Noverlay].y = source->peak->yf; - overlay[Noverlay].dx = 4.0*petrosian->radius; - overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major; - overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD; - overlay[Noverlay].text = NULL; - Noverlay ++; + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = source->peak->xf; + overlay[Noverlay].y = source->peak->yf; + overlay[Noverlay].dx = 4.0*petrosian->radius; + overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major; + overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); }