Index: /branches/eam_02_branch/psphot/src/Makefile.am =================================================================== --- /branches/eam_02_branch/psphot/src/Makefile.am (revision 12932) +++ /branches/eam_02_branch/psphot/src/Makefile.am (revision 12933) @@ -56,4 +56,5 @@ psphotDeblendSatstars.c \ psphotMosaicSubimage.c \ + psphotMakeResiduals.c \ psphotAddNoise.c Index: /branches/eam_02_branch/psphot/src/psphot.h =================================================================== --- /branches/eam_02_branch/psphot/src/psphot.h (revision 12932) +++ /branches/eam_02_branch/psphot/src/psphot.h (revision 12933) @@ -103,3 +103,5 @@ bool psphotDeblendSatstars (psArray *sources, psMetadata *recipe); +bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf); + #endif Index: /branches/eam_02_branch/psphot/src/psphotChoosePSF.c =================================================================== --- /branches/eam_02_branch/psphot/src/psphotChoosePSF.c (revision 12932) +++ /branches/eam_02_branch/psphot/src/psphotChoosePSF.c (revision 12933) @@ -155,4 +155,6 @@ try = models->data[bestN]; + // XXX calculate the PSF residual table here (and add to the pmPSF structure) + // measure and save the median value of dparams[PM_PAR_SXX],dparams[PM_PAR_SYY] // these are used by psphotEvalPSF to identify extended sources @@ -221,4 +223,7 @@ } + // build a PSF residual image + psphotMakeResiduals (try->sources, recipe, try->psf); + // XXX test dump of psf star data and psf-subtracted image if (psTraceGetLevel("psphot.psfstars") > 5) { Index: /branches/eam_02_branch/psphot/src/psphotMakeResiduals.c =================================================================== --- /branches/eam_02_branch/psphot/src/psphotMakeResiduals.c (revision 12933) +++ /branches/eam_02_branch/psphot/src/psphotMakeResiduals.c (revision 12933) @@ -0,0 +1,164 @@ +# include "psphotInternal.h" + +bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf) { + + bool status, isPSF; + double flux, dflux; + psU8 mflux; + + psTimerStart ("residuals"); + + int xBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.XBIN"); + PS_ASSERT (status, false); + + int yBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.YBIN"); + PS_ASSERT (status, false); + + // user parameters: + // size of aperture (determine from source images?) + // binning factor + + // select the subset of sources which are the PSFSTARs + // for each input source: + // - construct a residual image, renormalized + // - construct a renormalized weight image + // - construct a new mask image + + // construct the output residual table (Nx*DX,Ny*DY) + // for each output pixel: + // - construct a histogram of the values & weights (interpolate to the common pixel coordinate) + // - measure the robust median & sigma + // - reject (mask) input pixels which are outliers + // - re-measure the robust median & sigma + // - set output pixel, weight, and mask + + // determine the maximum image size from the input sources + int xSize = 0; + int ySize = 0; + + psVector *xC = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32); + + psArray *input = psArrayAllocEmpty (100); + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + + // which model to use? + pmModel *model = pmSourceGetModel (&isPSF, source); + if (model == NULL) continue; // model must be defined + + psImage *image = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + psImage *weight = psImageCopy (NULL, source->weight, PS_TYPE_F32); + psImage *mask = psImageCopy (NULL, source->mask, PS_TYPE_U8); + pmModelSub (image, mask, model, false, false); + psFree (mask); + + // re-normalize image and weight + float Io = model->params->data.F32[PM_PAR_I0]; + psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32)); + psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32)); + + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc( + PS_INTERPOLATE_BILINEAR, + image, weight, NULL, 0, 0.0, 0.0, 1, 0, 0.0); + psArrayAdd (input, 100, interp); + + // save the X,Y position for future reference + xC->data.F32[xC->n] = model->params->data.F32[PM_PAR_XPOS]; + yC->data.F32[yC->n] = model->params->data.F32[PM_PAR_YPOS]; + psVectorExtend (xC, 100, 1); + psVectorExtend (yC, 100, 1); + + xSize = PS_MAX (xSize, image->numCols); + ySize = PS_MAX (ySize, image->numRows); + + // free up the excess references + psFree (image); + psFree (weight); + psFree (interp); + } + pmResiduals *resid = pmResidualsAlloc (xSize, ySize, xBin, yBin); + + // x(resid) = (x(image) - Xo)*xBin + xCenter + + psVector *fluxes = psVectorAlloc (input->n, PS_TYPE_F32); + psVector *dfluxes = psVectorAlloc (input->n, PS_TYPE_F32); + psVector *fmasks = psVectorAlloc (input->n, PS_TYPE_U8); + + psStats *fluxClip = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psStats *fluxClipDef = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psStats *fluxStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + psStats *fluxStatsDef = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + + // build the residual image pixel-by-pixel + for (int oy = 0; oy < resid->image->numRows; oy++) { + fprintf (stderr, "."); + for (int ox = 0; ox < resid->image->numCols; ox++) { + + for (int i = 0; i < input->n; i++) { + + psImageInterpolateOptions *interp = input->data[i]; + + // fractional image position + float ix = (ox - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0; + float iy = (oy - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0; + + mflux = 0; + psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp); + fluxes->data.F32[i] = flux; + dfluxes->data.F32[i] = dflux; + fmasks->data.U8[i] = mflux; + // fprintf (stderr, "%f %f : %f %f (%d)\n", ix, iy, flux, dflux, fmasks->data.U8[i]); + } + + *fluxClip = *fluxClipDef; + psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff); + psErrorClear(); + + // mark input pixels which are more than N sigma from the median + for (int i = 0; i < fluxes->n; i++) { + float delta = fluxes->data.F32[i] - fluxClip->robustMedian; + float sigma = sqrt (dfluxes->data.F32[i]); + + float swing = fabs(delta) / sigma; + + // make this a user option + if (swing > 3.0) { + fmasks->data.U8[i] = 1; + } + } + + *fluxStats = *fluxStatsDef; + psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff); + psErrorClear(); + + resid->image->data.F32[oy][ox] = fluxStats->sampleMean; + resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev; + } + } + + psphotSaveImage (NULL, resid->image, "resid.im.fits"); + psphotSaveImage (NULL, resid->weight, "resid.wt.fits"); + psphotSaveImage (NULL, resid->mask, "resid.mk.fits"); + + psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals")); + + psFree (xC); + psFree (yC); + psFree (input); + + psFree (fluxes); + psFree (dfluxes); + psFree (fmasks); + + psFree (fluxStats); + psFree (fluxStatsDef); + psFree (fluxClip); + psFree (fluxClipDef); + + psf->residuals = resid; + return true; +} Index: /branches/eam_02_branch/psphot/src/psphotRadialPlot.c =================================================================== --- /branches/eam_02_branch/psphot/src/psphotRadialPlot.c (revision 12932) +++ /branches/eam_02_branch/psphot/src/psphotRadialPlot.c (revision 12933) @@ -8,5 +8,5 @@ static int nCount = 0; -bool psphotRadialPlot (int *kapa, char *filename, pmSource *source) { +bool psphotRadialPlot (int *kapa, const char *filename, pmSource *source) { Graphdata graphdata; @@ -95,7 +95,7 @@ sprintf (pagename, "%02d", nCount); if (nCount == 0) { - KiiPS (*kapa, false, KAPA_PS_NEWPLOT, filename, pagename); + KiiPS (*kapa, filename, false, KAPA_PS_NEWPLOT, pagename); } else { - KiiPS (*kapa, false, KAPA_PS_NEWPAGE, filename, pagename); + KiiPS (*kapa, filename, false, KAPA_PS_NEWPAGE, pagename); }