Index: /trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- /trunk/psModules/src/objects/pmPSF_IO.c (revision 19905) +++ /trunk/psModules/src/objects/pmPSF_IO.c (revision 19906) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-07-17 22:38:15 $ + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-06 13:05:13 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -319,13 +319,27 @@ pmPSFClump psfClump; - psfClump.X = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X"); assert (status); - psfClump.Y = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y"); assert (status); - psfClump.dX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX"); assert (status); - psfClump.dY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY"); assert (status); - - psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLX", 0, "psf clump center", psfClump.X); - psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLY", 0, "psf clump center", psfClump.Y); - psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDX", 0, "psf clump size", psfClump.dX); - psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDY", 0, "psf clump size", psfClump.dY); + // we now save clump parameters for each region : need to save all of those + int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + psMetadataAddS32 (header, PS_LIST_TAIL, "PSF_CLN", PS_META_REPLACE, "number of psf clump regions", nRegions); + for (int i = 0; i < nRegions; i++) { + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); + psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + + psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); assert (status); + psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); assert (status); + psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); assert (status); + psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); assert (status); + + char key[16]; + snprintf (key, 9, "CLX_%03d", i); + psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.X); + snprintf (key, 9, "CLY_%03d", i); + psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.Y); + snprintf (key, 9, "CLDX_%03d", i); + psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dX); + snprintf (key, 9, "CLDY_%03d", i); + psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dY); + } // save the dimensions of each parameter @@ -648,16 +662,65 @@ } - // psf clump data - pmPSFClump psfClump; - - psfClump.X = psMetadataLookupF32 (&status, header, "PSF_CLX" ); assert(status); - psfClump.Y = psMetadataLookupF32 (&status, header, "PSF_CLY" ); assert(status); - psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX"); assert(status); - psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY"); assert(status); - - psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.X" , PS_META_REPLACE, "psf clump center", psfClump.X); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.Y" , PS_META_REPLACE, "psf clump center", psfClump.Y); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump size", psfClump.dX); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump size", psfClump.dY); + // read the psf clump data for each region + int nRegions = psMetadataLookupS32 (&status, header, "PSF_CLN"); + if (!status) { + // read old-style psf clump data + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.000"); + psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); + + psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + if (!regionMD) { + regionMD = psMetadataAlloc(); + psMetadataAddMetadata (recipe, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); + psFree (regionMD); + } + + // psf clump data + pmPSFClump psfClump; + + psfClump.X = psMetadataLookupF32 (&status, header, "PSF_CLX" ); assert(status); + psfClump.Y = psMetadataLookupF32 (&status, header, "PSF_CLY" ); assert(status); + psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX"); assert(status); + psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY"); assert(status); + + 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 { + psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", nRegions); + + for (int i = 0; i < nRegions; i++) { + char key[10]; + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); + + psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + if (!regionMD) { + regionMD = psMetadataAlloc(); + psMetadataAddMetadata (recipe, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); + psFree (regionMD); + } + + // psf clump data + pmPSFClump psfClump; + + snprintf (key, 9, "CLX_%03d", i); + psfClump.X = psMetadataLookupF32 (&status, header, key); assert(status); + snprintf (key, 9, "CLY_%03d", i); + psfClump.Y = psMetadataLookupF32 (&status, header, key); assert(status); + snprintf (key, 9, "CLDX_%03d", i); + psfClump.dX = psMetadataLookupF32 (&status, header, key); assert(status); + snprintf (key, 9, "CLDY_%03d", i); + psfClump.dY = psMetadataLookupF32 (&status, header, key); assert(status); + + 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); + } + } options->poissonErrorsPhotLMM = psMetadataLookupBool (&status, header, "ERR_LMM"); Index: /trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtry.c (revision 19905) +++ /trunk/psModules/src/objects/pmPSFtry.c (revision 19906) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.61 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-09-11 00:38:23 $ + * @version $Revision: 1.62 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-06 13:05:13 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -150,5 +150,5 @@ } psLogMsg ("psphot.psftry", 4, "fit ext: %f sec for %d of %ld sources\n", psTimerMark ("fit"), Next, sources->n); - psTrace ("psphot.psftry", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n); + psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n); if (Next == 0) { @@ -215,5 +215,5 @@ psfTry->metricErr->data.F32[i] = source->errMag; - psTrace ("psphot.psftry", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n); + psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n); Npsf ++; } @@ -221,5 +221,5 @@ psLogMsg ("psphot.psftry", 4, "fit psf: %f sec for %d of %ld sources\n", psTimerMark ("fit"), Npsf, sources->n); - psTrace ("psphot.psftry", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n); + psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n); if (Npsf == 0) { @@ -568,5 +568,5 @@ // check the fit residuals and increase Nx,Ny until the error is minimized // pmPSFParamTrend increases the number along the longer of x or y - for (int i = 1; i < PS_MAX (psf->trendNx, psf->trendNy); i++) { + for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) { psVectorInit (mask, 0); @@ -578,5 +578,7 @@ // use the bright end scatter from the constant fit to set the scale for the // versions with 2D variation. potentially scale by poisson errors... - if (i == 1) { + // if (i == 1) { + // XXX let's drop this for the moment: relies on valid result from errorFloor + if (0) { // if (psf->poissonErrorsParams) do something else.. dz = psVectorAlloc (sources->n, PS_TYPE_F32); @@ -629,12 +631,30 @@ for (int i = 0; i < psf->psfTrendStats->clipIter; i++) { + psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); + psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); + + pmTrend2D *trend = NULL; + float mean, stdev; + // XXX we are using the same stats structure on each pass: do we need to re-init it? bool status = true; - status &= pmTrend2DFit (psf->params->data[PM_PAR_E0], srcMask, 0xff, x, y, e0, NULL); - psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n); - status &= pmTrend2DFit (psf->params->data[PM_PAR_E1], srcMask, 0xff, x, y, e1, NULL); - psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n); - status &= pmTrend2DFit (psf->params->data[PM_PAR_E2], srcMask, 0xff, x, y, e2, NULL); - psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n); + + trend = psf->params->data[PM_PAR_E0]; + status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n); + + trend = psf->params->data[PM_PAR_E1]; + status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n); + + trend = psf->params->data[PM_PAR_E2]; + status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n); if (!status) { @@ -676,8 +696,12 @@ if (psf->fieldNx > psf->fieldNy) { Nx = scale; - Ny = PS_MAX (1, Nx * (psf->fieldNy / psf->fieldNx)); + float AR = psf->fieldNy / (float) psf->fieldNx; + Ny = (int) (Nx * AR + 0.5); + Ny = PS_MAX (1, Ny); } else { Ny = scale; - Nx = PS_MAX (1, Ny * (psf->fieldNx / psf->fieldNy)); + float AR = psf->fieldNx / (float) psf->fieldNy; + Nx = (int) (Ny * AR + 0.5); + Nx = PS_MAX (1, Nx); } @@ -716,10 +740,38 @@ for (int i = 0; i < nIter; i++) { // XXX we are using the same stats structure on each pass: do we need to re-init it? - pmTrend2DFit (psf->params->data[PM_PAR_E0], mask, 0xff, x, y, e0obs, dz); - psTrace ("psModules.objects", 5, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0obs->n); - pmTrend2DFit (psf->params->data[PM_PAR_E1], mask, 0xff, x, y, e1obs, dz); - psTrace ("psModules.objects", 5, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1obs->n); - pmTrend2DFit (psf->params->data[PM_PAR_E2], mask, 0xff, x, y, e2obs, dz); - psTrace ("psModules.objects", 5, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2obs->n); + // pmTrend2DFit (psf->params->data[PM_PAR_E0], mask, 0xff, x, y, e0obs, dz); + // psTrace ("psModules.objects", 5, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0obs->n); + // pmTrend2DFit (psf->params->data[PM_PAR_E1], mask, 0xff, x, y, e1obs, dz); + // psTrace ("psModules.objects", 5, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1obs->n); + // pmTrend2DFit (psf->params->data[PM_PAR_E2], mask, 0xff, x, y, e2obs, dz); + // psTrace ("psModules.objects", 5, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2obs->n); + + psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); + psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); + + pmTrend2D *trend = NULL; + float mean, stdev; + + // XXX we are using the same stats structure on each pass: do we need to re-init it? + bool status = true; + + trend = psf->params->data[PM_PAR_E0]; + status &= pmTrend2DFit (trend, mask, 0xff, x, y, e0obs, dz); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs->n); + + trend = psf->params->data[PM_PAR_E1]; + status &= pmTrend2DFit (trend, mask, 0xff, x, y, e1obs, dz); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs->n); + + trend = psf->params->data[PM_PAR_E2]; + status &= pmTrend2DFit (trend, mask, 0xff, x, y, e2obs, dz); + mean = psStatsGetValue (trend->stats, meanOption); + stdev = psStatsGetValue (trend->stats, stdevOption); + psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n); + } psf->psfTrendStats->clipIter = nIter; // restore default setting @@ -734,6 +786,33 @@ psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs, "-", (void *) e2fit); +// XXX this code determines the formal error on the map, not the scatter +# if (0) // measure errors on the map pixels pmTrend2D *trend; + psStatsOptions meanOption; + psStatsOptions stdevOption; + float mapErrorSum = 0.0; + float mapError = 0.0; + + trend = psf->params->data[PM_PAR_E0]; + meanOption = psStatsMeanOption (trend->stats->options); + stdevOption = psStatsStdevOption (trend->stats->options); + mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption)); + mapErrorSum += mapError; + psTrace ("psModules.objects", 4, "E0 error: %f\n", sqrt(mapError)); + + trend = psf->params->data[PM_PAR_E1]; + meanOption = psStatsMeanOption (trend->stats->options); + stdevOption = psStatsStdevOption (trend->stats->options); + mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption)); + mapErrorSum += mapError; + psTrace ("psModules.objects", 4, "E1 error: %f\n", sqrt(mapError)); + + trend = psf->params->data[PM_PAR_E2]; + meanOption = psStatsMeanOption (trend->stats->options); + stdevOption = psStatsStdevOption (trend->stats->options); + mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption)); + mapErrorSum += mapError; + psTrace ("psModules.objects", 4, "E2 error: %f\n", sqrt(mapError)); trend = psf->params->data[PM_PAR_E0]; @@ -769,4 +848,5 @@ psTrace ("psModules.objects", PS_LOG_INFO, "total map error: %f\n", sqrt(mapErrorSum)); +# endif // measure systematic errorFloor & systematic / photon scale factor @@ -776,5 +856,6 @@ psStatsStdevOption(psf->psfTrendStats->options)); - *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum); + // XXX Bogus: *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum); + *errorTotal = *errorFloor; psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup); @@ -815,4 +896,5 @@ for (int i = 0; i < nBin; i++) { int j; + int nValid = 0; for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) { int N = index->data.U32[n]; @@ -822,5 +904,8 @@ mkSubset->data.U8[j] = mask->data.U8[N]; - } + if (!mask->data.U8[N]) nValid ++; + } + if (nValid < 3) continue; + dE0subset->n = j; dE1subset->n = j; @@ -849,4 +934,6 @@ } } + *errorFloor = min; + psFree (dE0subset); psFree (dE1subset); Index: /trunk/psModules/src/objects/pmSource.c =================================================================== --- /trunk/psModules/src/objects/pmSource.c (revision 19905) +++ /trunk/psModules/src/objects/pmSource.c (revision 19906) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.57 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-10-03 20:59:16 $ + * @version $Revision: 1.58 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-06 13:05:13 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -274,5 +274,6 @@ // XXX this function should probably accept the values, not a recipe. wrap with a // psphot-specific function which applies the recipe values -pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *recipe) +// only apply selection to sources within specified region +pmPSFClump pmSourcePSFClump(psRegion *region, psArray *sources, psMetadata *recipe) { psTrace("psModules.objects", 5, "---- begin ----\n"); @@ -329,4 +330,9 @@ continue; + if (tmpSrc->peak->x < region->x0) continue; + if (tmpSrc->peak->x > region->x1) continue; + if (tmpSrc->peak->y < region->y0) continue; + if (tmpSrc->peak->y > region->y1) continue; + // Sx,Sy are limited at 0. a peak at 0,0 is artificial if (fabs(tmpSrc->moments->Mxx) < 0.05) @@ -427,4 +433,9 @@ continue; + if (tmpSrc->peak->x < region->x0) continue; + if (tmpSrc->peak->x > region->x1) continue; + if (tmpSrc->peak->y < region->y0) continue; + if (tmpSrc->peak->y > region->y1) continue; + if (tmpSrc->moments->Mxx < minSx) continue; @@ -475,5 +486,5 @@ *****************************************************************************/ -bool pmSourceRoughClass(psArray *sources, psMetadata *recipe, pmPSFClump clump, psMaskType maskSat) +bool pmSourceRoughClass(psRegion *region, psArray *sources, psMetadata *recipe, pmPSFClump clump, psMaskType maskSat) { psTrace("psModules.objects", 5, "---- begin ----"); @@ -509,4 +520,9 @@ pmSource *source = (pmSource *) sources->data[i]; + + 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; source->peak->type = 0; Index: /trunk/psModules/src/objects/pmSource.h =================================================================== --- /trunk/psModules/src/objects/pmSource.h (revision 19905) +++ /trunk/psModules/src/objects/pmSource.h (revision 19906) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-04-08 18:35:38 $ + * @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-10-06 13:05:13 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -186,4 +186,5 @@ */ pmPSFClump pmSourcePSFClump( + psRegion *region, ///< restrict measurement to specified region psArray *source, ///< The input pmSource psMetadata *metadata ///< Contains classification parameters @@ -201,5 +202,6 @@ */ bool pmSourceRoughClass( - psArray *source, ///< The input pmSource + psRegion *region, ///< restrict measurement to specified region + psArray *sources, ///< The input pmSources psMetadata *metadata, ///< Contains classification parameters pmPSFClump clump, ///< Statistics about the PSF clump