Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 36555) +++ trunk/psphot/src/Makefile.am (revision 36623) @@ -208,4 +208,5 @@ psphotMakeGrowthCurve.c \ psphotMagnitudes.c \ + psphotLensing.c \ psphotSetMaskBits.c \ psphotSkyReplace.c \ Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 36555) +++ trunk/psphot/src/psphot.h (revision 36623) @@ -144,4 +144,8 @@ bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf); bool psphotMagnitudes_Threaded (psThreadJob *job); + +bool psphotLensing (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotLensingReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotLensingPSFtrendsReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); bool psphotEfficiency (pmConfig *config, const pmFPAview *view, const char *filerule); Index: trunk/psphot/src/psphotApResid.c =================================================================== --- trunk/psphot/src/psphotApResid.c (revision 36555) +++ trunk/psphot/src/psphotApResid.c (revision 36623) @@ -214,5 +214,5 @@ // XXX make this user-configurable? - if (source->psfMagErr > 0.01) continue; + if (source->psfMagErr > 0.03) continue; // aperture residual for this source @@ -267,10 +267,11 @@ // set the max order (0 = constant) which the number of psf stars can support: - // rule of thumb: require 3 stars per 'cell' (order+1)^2 + // we require only 3 stars for n = 0, increase stars / cell for higher order int MaxOrderForStars = 0; - if (Npsf >= 12) MaxOrderForStars = 1; // 4 cells - if (Npsf >= 27) MaxOrderForStars = 2; // 9 cells - if (Npsf >= 48) MaxOrderForStars = 3; // 16 cells - if (Npsf > 75) MaxOrderForStars = 4; // 25 cells + if (Npsf >= 16) MaxOrderForStars = 1; // 4 cells + if (Npsf >= 54) MaxOrderForStars = 2; // 9 cells + if (Npsf >= 128) MaxOrderForStars = 3; // 16 cells + if (Npsf >= 300) MaxOrderForStars = 4; // 25 cells + if (Npsf > 576) MaxOrderForStars = 5; // 36 cells pmTrend2DMode mode = PM_TREND_MAP; Index: trunk/psphot/src/psphotChoosePSF.c =================================================================== --- trunk/psphot/src/psphotChoosePSF.c (revision 36555) +++ trunk/psphot/src/psphotChoosePSF.c (revision 36623) @@ -399,5 +399,5 @@ char *modelName = pmModelClassGetName (psf->type); psLogMsg ("psphot.pspsf", PS_LOG_WARN, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf")); - psLogMsg ("psphot.pspsf", PS_LOG_INFO, "psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid); + psLogMsg ("psphot.pspsf", PS_LOG_INFO, "psf model %s, ApResid: %f +/- %f (%d x %d model)\n", modelName, psf->ApResid, psf->dApResid, psf->trendNx, psf->trendNy); psFree (options); Index: trunk/psphot/src/psphotLensing.c =================================================================== --- trunk/psphot/src/psphotLensing.c (revision 36623) +++ trunk/psphot/src/psphotLensing.c (revision 36623) @@ -0,0 +1,319 @@ +# include "psphotInternal.h" + +// calculate lensing parameters (we have the moments already) +bool psphotLensing (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Lensing Parameters ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + if (!psphotLensingReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure lensing parameters for %s entry %d", filerule, i); + return false; + } + if (!psphotLensingPSFtrendsReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure lensing parameters for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +// this block measures the lensing parameters for all objects +bool psphotLensingReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) +{ + bool status; + pmSource *source; + + psTimerStart ("psphot.lensing"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // if we have measured the window, we will be saving the modified version of these recipe values on readout->analysis + float SIGMA = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_GAUSS_SIGMA"); + if (!status) { + SIGMA = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA"); + } + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping lensing"); + return true; + } + + for (int i = 0; i < sources->n; i++) { + source = sources->data[i]; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + + pmMoments *moments = source->moments; + + if (!source->lensingOBJ) { + source->lensingOBJ = pmSourceLensingAlloc (); + } + + // XXX big objects (eg, saturated stars) use 3*SIGMA. I'm ignoring this for now + pmSourceLensingShearFromMoments (source->lensingOBJ, moments, SIGMA); + pmSourceLensingSmearFromMoments (source->lensingOBJ, moments, SIGMA); + } + + psLogMsg ("psphot.lensing", PS_LOG_DETAIL, "calculate lensing parameters for %d objects\n", (int) sources->n); + + return true; +} + +// generate image maps for the Xij,ei elements based only on good stars +bool psphotLensingPSFtrendsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) +{ + int Nfail = 0; + bool status; + pmSource *source; + + psTimerStart ("psphot.lensing"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping lensing psf trends"); + return true; + } + + // gather the stats to assess the aperture residuals + int Npsf = 0; + psVector *psfX11sm = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *psfX12sm = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *psfX22sm = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *psf_e1sm = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *psf_e2sm = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *psfX11sh = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *psfX12sh = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *psfX22sh = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *psf_e1sh = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *psf_e2sh = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *xPos = psVectorAllocEmpty (300, PS_TYPE_F32); + psVector *yPos = psVectorAllocEmpty (300, PS_TYPE_F32); + + for (int i = 0; i < sources->n; i++) { + source = sources->data[i]; + + // only use good stars: + if (source->type != PM_SOURCE_TYPE_STAR) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + if (source->mode & PM_SOURCE_MODE_BLEND) continue; + if (source->mode & PM_SOURCE_MODE_FAIL) continue; + if (source->mode & PM_SOURCE_MODE_POOR) continue; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) continue; + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + if (source->mode2 & PM_SOURCE_MODE2_MATCHED) continue; // XXX ??? + + // use psf - ap to limit to good stars? + // XXX need to be sure if this is before or after ap corrections are done + if (!isfinite(source->apMag)) continue; + if (!isfinite(source->psfMag)) continue; + + // XXX make this user-configurable? + if (source->psfMagErr > 0.03) continue; + + // aperture residual for this source + float dap = source->apMag - source->psfMag; + + // good stars have |dap| < 0.2 + if (fabs(dap) > 0.2) { + Nfail ++; + psTrace ("psphot", 3, "fail : bad dap: %f", dap); + continue; + } + + pmModel *model = source->modelPSF; + + psVectorAppend (xPos, model->params->data.F32[PM_PAR_XPOS]); + psVectorAppend (yPos, model->params->data.F32[PM_PAR_YPOS]); + psVectorAppend (psfX11sm, source->lensingOBJ->smear->X11); + psVectorAppend (psfX12sm, source->lensingOBJ->smear->X12); + psVectorAppend (psfX22sm, source->lensingOBJ->smear->X22); + psVectorAppend (psf_e1sm, source->lensingOBJ->smear->e1); + psVectorAppend (psf_e2sm, source->lensingOBJ->smear->e2); + psVectorAppend (psfX11sh, source->lensingOBJ->shear->X11); + psVectorAppend (psfX12sh, source->lensingOBJ->shear->X12); + psVectorAppend (psfX22sh, source->lensingOBJ->shear->X22); + psVectorAppend (psf_e1sh, source->lensingOBJ->shear->e1); + psVectorAppend (psf_e2sh, source->lensingOBJ->shear->e2); + Npsf ++; + } + + // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different + // definition for 'order' (order_MAP = order_POLY + 1). in addition, we have a + // user-specified MAX order, which we should respect, regardless of the mode + + // set the max order (0 = constant) which the number of psf stars can support: + // we require only 3 stars for n = 0, increase stars / cell for higher order + int LensingTrendOrder = 0; + if (Npsf >= 16) LensingTrendOrder = 1; // 4 cells + if (Npsf >= 54) LensingTrendOrder = 2; // 9 cells + if (Npsf >= 128) LensingTrendOrder = 3; // 16 cells + if (Npsf >= 300) LensingTrendOrder = 4; // 25 cells + if (Npsf > 576) LensingTrendOrder = 5; // 36 cells + + pmTrend2DMode mode = PM_TREND_MAP; + if (mode == PM_TREND_MAP) { + LensingTrendOrder ++; + } + + // Nx,Ny are the superpixel size; but I need to choose them to allow a max of LensingTrendOrder values + // I'm not sure this is really sensible, but... + int NX = readout->image->numCols; + int NY = readout->image->numRows; + int Nx, Ny; + if (NX > NY) { + Nx = LensingTrendOrder; + Ny = PS_MAX (1, (int)(LensingTrendOrder * (NY / (float)(NX)) + 0.5)); + } else { + Ny = LensingTrendOrder; + Nx = PS_MAX (1, (int)(LensingTrendOrder * (NX / (float)(NY)) + 0.5)); + } + + psLogMsg ("psphot.lensing", PS_LOG_DETAIL, "generate lensing maps using %d objects (%d x %d grid)\n", Npsf, Nx, Ny); + + // XXX allow user to set this optionally? + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + + // measure Trend2D for the current spatial scale + pmTrend2D *trendX11sm = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + pmTrend2D *trendX12sm = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + pmTrend2D *trendX22sm = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + pmTrend2D *trend_e1sm = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + pmTrend2D *trend_e2sm = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + + pmTrend2D *trendX11sh = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + pmTrend2D *trendX12sh = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + pmTrend2D *trendX22sh = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + pmTrend2D *trend_e1sh = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + pmTrend2D *trend_e2sh = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + + bool goodFit = false; + if (!pmTrend2DFit (&goodFit, trendX11sm, NULL, 0xff, xPos, yPos, psfX11sm, NULL)) { psWarning ("failed to measure X11 smear trend"); } + if (!pmTrend2DFit (&goodFit, trendX12sm, NULL, 0xff, xPos, yPos, psfX12sm, NULL)) { psWarning ("failed to measure X12 smear trend"); } + if (!pmTrend2DFit (&goodFit, trendX22sm, NULL, 0xff, xPos, yPos, psfX22sm, NULL)) { psWarning ("failed to measure X22 smear trend"); } + if (!pmTrend2DFit (&goodFit, trend_e1sm, NULL, 0xff, xPos, yPos, psf_e1sm, NULL)) { psWarning ("failed to measure e1 smear trend"); } + if (!pmTrend2DFit (&goodFit, trend_e2sm, NULL, 0xff, xPos, yPos, psf_e2sm, NULL)) { psWarning ("failed to measure e2 smear trend"); } + + if (!pmTrend2DFit (&goodFit, trendX11sh, NULL, 0xff, xPos, yPos, psfX11sh, NULL)) { psWarning ("failed to measure X11 shear trend"); } + if (!pmTrend2DFit (&goodFit, trendX12sh, NULL, 0xff, xPos, yPos, psfX12sh, NULL)) { psWarning ("failed to measure X12 shear trend"); } + if (!pmTrend2DFit (&goodFit, trendX22sh, NULL, 0xff, xPos, yPos, psfX22sh, NULL)) { psWarning ("failed to measure X22 shear trend"); } + if (!pmTrend2DFit (&goodFit, trend_e1sh, NULL, 0xff, xPos, yPos, psf_e1sh, NULL)) { psWarning ("failed to measure e1 shear trend"); } + if (!pmTrend2DFit (&goodFit, trend_e2sh, NULL, 0xff, xPos, yPos, psf_e2sh, NULL)) { psWarning ("failed to measure e2 shear trend"); } + + // evaluate the PSF trend maps at the location of all sources + for (int i = 0; i < sources->n; i++) { + source = sources->data[i]; + pmPeak *peak = source->peak; + if (!peak) continue; + + float xPos = peak->xf; + float yPos = peak->yf; + + if (!source->lensingPSF) { + source->lensingPSF = pmSourceLensingAlloc (); + } + pmSourceLensing *lensing = source->lensingPSF; + + if (!lensing->shear) { + lensing->shear = pmLensingParsAlloc(); + } + pmLensingPars *shear = lensing->shear; + + shear->X11 = pmTrend2DEval (trendX11sh, xPos, yPos); + shear->X12 = pmTrend2DEval (trendX12sh, xPos, yPos); + shear->X22 = pmTrend2DEval (trendX22sh, xPos, yPos); + shear->e1 = pmTrend2DEval (trend_e1sh, xPos, yPos); + shear->e2 = pmTrend2DEval (trend_e2sh, xPos, yPos); + + if (!lensing->smear) { + lensing->smear = pmLensingParsAlloc(); + } + pmLensingPars *smear = lensing->smear; + + smear->X11 = pmTrend2DEval (trendX11sm, xPos, yPos); + smear->X12 = pmTrend2DEval (trendX12sm, xPos, yPos); + smear->X22 = pmTrend2DEval (trendX22sm, xPos, yPos); + smear->e1 = pmTrend2DEval (trend_e1sm, xPos, yPos); + smear->e2 = pmTrend2DEval (trend_e2sm, xPos, yPos); + } + + psFree (trendX11sm); + psFree (trendX12sm); + psFree (trendX22sm); + psFree (trend_e1sm); + psFree (trend_e2sm); + + psFree (trendX11sh); + psFree (trendX12sh); + psFree (trendX22sh); + psFree (trend_e1sh); + psFree (trend_e2sh); + + psFree (psfX11sm); + psFree (psfX12sm); + psFree (psfX22sm); + psFree (psf_e1sm); + psFree (psf_e2sm); + psFree (psfX11sh); + psFree (psfX12sh); + psFree (psfX22sh); + psFree (psf_e1sh); + psFree (psf_e2sh); + psFree (xPos ); + psFree (yPos ); + + psFree (stats); + + psLogMsg ("psphot.lensing", PS_LOG_DETAIL, "calculate lensing parameters for %d objects: %f sec\n", (int) sources->n, psTimerMark ("psphot.lensing")); + + return true; +} + Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 36555) +++ trunk/psphot/src/psphotReadout.c (revision 36623) @@ -340,7 +340,14 @@ // calculate source magnitudes if (!psphotMagnitudes(config, view, filerule)) { // pass 1 (detections->allSources) - psErrorStackPrint(stderr, "Unable to do magnitudes."); + psErrorStackPrint(stderr, "Unable to do magnitudes."); psErrorClear(); } + + // calculate lensing parameters + if (!psphotLensing(config, view, filerule)) { + psErrorStackPrint(stderr, "Unable to do lensing parameters."); + psErrorClear(); + } + if (!psphotEfficiency(config, view, filerule)) { // pass 1 psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); @@ -353,12 +360,12 @@ // replace background in residual image if (!psphotSkyReplace (config, view, filerule)) { // pass 1 - psErrorStackPrint(stderr, "Unable to replace sky"); - psErrorClear(); + psErrorStackPrint(stderr, "Unable to replace sky"); + psErrorClear(); } // drop the references to the image pixels held by each source if (!psphotSourceFreePixels (config, view, filerule)) { // pass 1 - psErrorStackPrint(stderr, "Unable to free source pixels"); - psErrorClear(); + psErrorStackPrint(stderr, "Unable to free source pixels"); + psErrorClear(); }