Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 30621) +++ trunk/psModules/src/objects/Makefile.am (revision 31153) @@ -32,4 +32,5 @@ pmSourceFitSet.c \ pmSourcePhotometry.c \ + pmSourceOutputs.c \ pmSourceIO.c \ pmSourceIO_RAW.c \ @@ -102,4 +103,5 @@ pmSourceFitSet.h \ pmSourcePhotometry.h \ + pmSourceOutputs.h \ pmSourceIO.h \ pmSourcePlots.h \ Index: trunk/psModules/src/objects/models/pmModel_DEV.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_DEV.c (revision 30621) +++ trunk/psModules/src/objects/models/pmModel_DEV.c (revision 31153) @@ -217,44 +217,29 @@ bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { - pmMoments *moments = source->moments; - pmPeak *peak = source->peak; - psF32 *PAR = model->params->data.F32; - - psEllipseMoments emoments; - emoments.x2 = moments->Mxx; - emoments.xy = moments->Mxy; - emoments.y2 = moments->Myy; - - // force the axis ratio to be < 20.0 - psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); - - if (!isfinite(axes.major)) return false; - if (!isfinite(axes.minor)) return false; - if (!isfinite(axes.theta)) return false; - - psEllipseShape shape = psEllipseAxesToShape (axes); - - if (!isfinite(shape.sx)) return false; - if (!isfinite(shape.sy)) return false; - if (!isfinite(shape.sxy)) return false; - - // the other parameters depend on the guess for PAR_7 + psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 + PAR[PM_PAR_SKY] = 0.0; + + // set the shape parameters + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + return false; + } + + // the normalization is modified by the slope float index = 0.5 / ALPHA; float bn = 1.9992*index - 0.3271; - // float fR = 1.0 / (sqrt(2.0) * pow (bn, index)); float Io = exp(0.5*bn); - float Sxx = PS_MAX(0.5, shape.sx); - float Syy = PS_MAX(0.5, shape.sy); - - PAR[PM_PAR_SKY] = 0.0; - PAR[PM_PAR_I0] = peak->flux / Io; - PAR[PM_PAR_XPOS] = peak->xf; - PAR[PM_PAR_YPOS] = peak->yf; - // PAR[PM_PAR_SXX] = Sxx * fR; - // PAR[PM_PAR_SYY] = Syy * fR; - PAR[PM_PAR_SXX] = Sxx; - PAR[PM_PAR_SYY] = Syy; - PAR[PM_PAR_SXY] = shape.sxy; + // set the model normalization + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + PAR[PM_PAR_I0] /= Io; + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } return(true); @@ -322,5 +307,7 @@ psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA); - psAssert (isfinite(radius), "fix this code: z should not be nan"); + psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]); + return (radius); } Index: trunk/psModules/src/objects/models/pmModel_EXP.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_EXP.c (revision 30621) +++ trunk/psModules/src/objects/models/pmModel_EXP.c (revision 31153) @@ -210,33 +210,24 @@ bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { - pmMoments *moments = source->moments; - pmPeak *peak = source->peak; - psF32 *PAR = model->params->data.F32; - - psEllipseMoments emoments; - emoments.x2 = moments->Mxx; - emoments.xy = moments->Mxy; - emoments.y2 = moments->Myy; - - // force the axis ratio to be < 20.0 - psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); - - if (!isfinite(axes.major)) return false; - if (!isfinite(axes.minor)) return false; - if (!isfinite(axes.theta)) return false; - - psEllipseShape shape = psEllipseAxesToShape (axes); - - if (!isfinite(shape.sx)) return false; - if (!isfinite(shape.sy)) return false; - if (!isfinite(shape.sxy)) return false; - + psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 PAR[PM_PAR_SKY] = 0.0; - PAR[PM_PAR_I0] = peak->flux; - PAR[PM_PAR_XPOS] = peak->xf; - PAR[PM_PAR_YPOS] = peak->yf; - PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); - PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); - PAR[PM_PAR_SXY] = shape.sxy; + + // set the shape parameters + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + return false; + } + + // set the model normalization + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + return(true); } @@ -303,5 +294,6 @@ psF64 radius = axes.major * sqrt (2.0) * zn; - psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]); + psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]); return (radius); } Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 30621) +++ trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 31153) @@ -193,24 +193,24 @@ bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { - pmMoments *moments = source->moments; - pmPeak *peak = source->peak; - psF32 *PAR = model->params->data.F32; - - psEllipseMoments emoments; - emoments.x2 = moments->Mxx; - emoments.y2 = moments->Myy; - emoments.xy = moments->Mxy; - - // force the axis ratio to be < 20.0 - psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); - psEllipseShape shape = psEllipseAxesToShape (axes); - + psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 PAR[PM_PAR_SKY] = 0.0; - PAR[PM_PAR_I0] = peak->flux; - PAR[PM_PAR_XPOS] = peak->xf; - PAR[PM_PAR_YPOS] = peak->yf; - PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); - PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); - PAR[PM_PAR_SXY] = shape.sxy; + + // set the shape parameters + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + return false; + } + + // set the model normalization + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + return(true); } @@ -258,5 +258,6 @@ psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux)); - psAssert (isfinite(radius), "fix this code: radius should not be nan for %f", PAR[PM_PAR_I0]); + psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]); return (radius); Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 30621) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 31153) @@ -194,23 +194,24 @@ bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { - pmMoments *moments = source->moments; - pmPeak *peak = source->peak; - psF32 *PAR = model->params->data.F32; - - psEllipseMoments emoments; - emoments.x2 = moments->Mxx; - emoments.xy = moments->Mxy; - emoments.y2 = moments->Myy; - - psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); - psEllipseShape shape = psEllipseAxesToShape (axes); - + psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 PAR[PM_PAR_SKY] = 0.0; - PAR[PM_PAR_I0] = peak->flux; - PAR[PM_PAR_XPOS] = peak->xf; - PAR[PM_PAR_YPOS] = peak->yf; - PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); - PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); - PAR[PM_PAR_SXY] = shape.sxy; + + // set the shape parameters + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + return false; + } + + // set the model normalization + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + return(true); } @@ -285,7 +286,11 @@ // choose a z value guaranteed to be beyond our limit float z0 = pow((1.0 / limit), (1.0 / 3.0)); - psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_I0]); + psAssert (isfinite(z0), "fix this code: z0 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]); + float z1 = (1.0 / limit); - psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_I0]); + psAssert (isfinite(z1), "fix this code: z1 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]); + z1 = PS_MAX (z0, z1); z0 = 0.0; Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 30621) +++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 31153) @@ -213,33 +213,25 @@ bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { - pmMoments *moments = source->moments; - pmPeak *peak = source->peak; - psF32 *PAR = model->params->data.F32; - - psEllipseMoments emoments; - emoments.x2 = moments->Mxx; - emoments.xy = moments->Mxy; - emoments.y2 = moments->Myy; - - // force the axis ratio to be < 20.0 - psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); - - if (!isfinite(axes.major)) return false; - if (!isfinite(axes.minor)) return false; - if (!isfinite(axes.theta)) return false; - - psEllipseShape shape = psEllipseAxesToShape (axes); - - if (!isfinite(shape.sx)) return false; - if (!isfinite(shape.sy)) return false; - if (!isfinite(shape.sxy)) return false; - + psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 PAR[PM_PAR_SKY] = 0.0; - PAR[PM_PAR_I0] = peak->flux; - PAR[PM_PAR_XPOS] = peak->xf; - PAR[PM_PAR_YPOS] = peak->yf; - PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); - PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); - PAR[PM_PAR_SXY] = shape.sxy; + + // set the shape parameters + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + return false; + } + + // set the model normalization + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + + // extra parameter PAR[PM_PAR_7] = 0.5; @@ -314,5 +306,6 @@ float z0 = 0.0; float z1 = pow((1.0 / limit), (1.0 / ALPHA)); - psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]); + psAssert (isfinite(z1), "fix this code: z1 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]); if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0; Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 30621) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 31153) @@ -214,33 +214,25 @@ bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { - pmMoments *moments = source->moments; - pmPeak *peak = source->peak; - psF32 *PAR = model->params->data.F32; - - psEllipseMoments emoments; - emoments.x2 = moments->Mxx; - emoments.xy = moments->Mxy; - emoments.y2 = moments->Myy; - - // force the axis ratio to be < 20.0 - psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); - - if (!isfinite(axes.major)) return false; - if (!isfinite(axes.minor)) return false; - if (!isfinite(axes.theta)) return false; - - psEllipseShape shape = psEllipseAxesToShape (axes); - - if (!isfinite(shape.sx)) return false; - if (!isfinite(shape.sy)) return false; - if (!isfinite(shape.sxy)) return false; - + psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 PAR[PM_PAR_SKY] = 0.0; - PAR[PM_PAR_I0] = peak->flux; - PAR[PM_PAR_XPOS] = peak->xf; - PAR[PM_PAR_YPOS] = peak->yf; - PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); - PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); - PAR[PM_PAR_SXY] = shape.sxy; + + // set the shape parameters + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + return false; + } + + // set the model normalization + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + + // extra parameter PAR[PM_PAR_7] = 1.0; @@ -315,5 +307,6 @@ float z0 = 0.0; float z1 = pow((1.0 / limit), (1.0 / ALPHA)); - psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]); + psAssert (isfinite(z1), "fix this code: z1 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]); if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0; Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 30621) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 31153) @@ -203,33 +203,25 @@ bool PM_MODEL_GUESS (pmModel *model, pmSource *source) { - pmMoments *moments = source->moments; - pmPeak *peak = source->peak; - psF32 *PAR = model->params->data.F32; - - psEllipseMoments emoments; - emoments.x2 = moments->Mxx; - emoments.xy = moments->Mxy; - emoments.y2 = moments->Myy; - - // force the axis ratio to be < 20.0 - psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); - - if (!isfinite(axes.major)) return false; - if (!isfinite(axes.minor)) return false; - if (!isfinite(axes.theta)) return false; - - psEllipseShape shape = psEllipseAxesToShape (axes); - - if (!isfinite(shape.sx)) return false; - if (!isfinite(shape.sy)) return false; - if (!isfinite(shape.sxy)) return false; - + psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 PAR[PM_PAR_SKY] = 0.0; - PAR[PM_PAR_I0] = peak->flux; - PAR[PM_PAR_XPOS] = peak->xf; - PAR[PM_PAR_YPOS] = peak->yf; - PAR[PM_PAR_SXX] = PS_MAX(0.5, shape.sx); - PAR[PM_PAR_SYY] = PS_MAX(0.5, shape.sy); - PAR[PM_PAR_SXY] = shape.sxy; + + // set the shape parameters + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + return false; + } + + // set the model normalization + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + + // extra parameter PAR[PM_PAR_7] = 1.5; @@ -305,7 +297,11 @@ // choose a z value guaranteed to be beyond our limit float z0 = pow((1.0 / limit), (1.0 / PAR[PM_PAR_7])); - psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_7]); + psAssert (isfinite(z0), "fix this code: z0 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], PAR[PM_PAR_7]); + float z1 = (1.0 / limit); - psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]); + psAssert (isfinite(z1), "fix this code: z1 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], PAR[PM_PAR_7]); + z1 = PS_MAX (z0, z1); z0 = 0.0; Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 30621) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 31153) @@ -222,6 +222,8 @@ { pmMoments *moments = source->moments; - pmPeak *peak = source->peak; psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 + PAR[PM_PAR_SKY] = 0.0; // the other parameters depend on the guess for PAR_7 @@ -236,4 +238,5 @@ float Zero = 1.16 - 0.615 * PAR[PM_PAR_7]; + // Sersic shape is a bit special psEllipseMoments emoments; emoments.x2 = moments->Mxx; @@ -273,11 +276,18 @@ float Syy = PS_MAX(0.5, shape.sy); - PAR[PM_PAR_SKY] = 0.0; - PAR[PM_PAR_I0] = peak->flux / Io; - PAR[PM_PAR_XPOS] = peak->xf; - PAR[PM_PAR_YPOS] = peak->yf; PAR[PM_PAR_SXX] = Sxx; PAR[PM_PAR_SYY] = Syy; PAR[PM_PAR_SXY] = shape.sxy; + + // set the model normalization (adjust for Sersic best guess) + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + PAR[PM_PAR_I0] /= Io; + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } return(true); @@ -346,7 +356,6 @@ psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]); - // fprintf (stderr, "sersic model %f %f, n %f, radius: %f, zn: %f, f/Io: %f, major: %f\n", PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], PAR[PM_PAR_7], radius, zn, flux/PAR[PM_PAR_I0], axes.major); - - psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]); + psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], PAR[PM_PAR_7]); return (radius); } Index: trunk/psModules/src/objects/pmFootprint.c =================================================================== --- trunk/psModules/src/objects/pmFootprint.c (revision 30621) +++ trunk/psModules/src/objects/pmFootprint.c (revision 31153) @@ -98,7 +98,12 @@ } +// XXX not actually used anywhere pmFootprint *pmFootprintNormalize(pmFootprint *fp) { if (fp != NULL && !fp->normalized) { - fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN); + if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE) { + fp->peaks = psArraySort(fp->peaks, pmPeaksSortBySmoothFluxDescend); + } else { + fp->peaks = psArraySort(fp->peaks, pmPeaksSortByRawFluxDescend); + } fp->normalized = true; } Index: trunk/psModules/src/objects/pmFootprint.h =================================================================== --- trunk/psModules/src/objects/pmFootprint.h (revision 30621) +++ trunk/psModules/src/objects/pmFootprint.h (revision 31153) @@ -10,4 +10,10 @@ #ifndef PM_FOOTPRINT_H #define PM_FOOTPRINT_H + +// We need to choose up front if the culling algorithm uses the raw or smoothed image. +// depending on which we choose, we should produce sorted peaks based on peak->rawFlux or +// peak->smoothFlux + +# define PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE 1 typedef struct { @@ -84,5 +90,7 @@ const float nsigma_delta, // how many sigma above local background a peak needs to be to survive const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance - const float min_threshold // minimum permitted coll height + const float min_threshold, // minimum permitted coll height + const float max_threshold,// maximum permitted coll height + const bool isWeightVar // the input weight may be variance (sigma^2) or S/N (1/sigma) ); Index: trunk/psModules/src/objects/pmFootprintAssignPeaks.c =================================================================== --- trunk/psModules/src/objects/pmFootprintAssignPeaks.c (revision 30621) +++ trunk/psModules/src/objects/pmFootprintAssignPeaks.c (revision 31153) @@ -46,10 +46,10 @@ */ psImage *ids = pmSetFootprintArrayIDs(footprints, true); - assert (ids != NULL); - assert (ids->type.type == PS_TYPE_S32); - const int row0 = ids->row0; - const int col0 = ids->col0; - const int numRows = ids->numRows; - const int numCols = ids->numCols; + if (ids) { assert (ids->type.type == PS_TYPE_S32); } + + const int row0 = ids ? ids->row0 : 0; + const int col0 = ids ? ids->col0 : 0; + const int numRows = ids ? ids->numRows : -1; + const int numCols = ids ? ids->numCols : -1; for (int i = 0; i < peaks->n; i++) { @@ -58,6 +58,7 @@ const int y = peak->y - row0; - assert (x >= 0 && x < numCols && y >= 0 && y < numRows); - int id = ids->data.S32[y][x - col0]; + if (ids) { assert (x >= 0 && x < numCols && y >= 0 && y < numRows);} + int id = ids ? ids->data.S32[y][x - col0] : 0; + // XXX I think the '[x - col0]' above is just wrong (should be [x], but never gets triggerd. if (id == 0) { // peak isn't in a footprint, so make one for it @@ -86,5 +87,10 @@ if (fp->peaks->n == 1) continue; - fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN); + // make sure the peaks are sorted in a way consistent with our cull process + if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE) { + fp->peaks = psArraySort(fp->peaks, pmPeaksSortBySmoothFluxDescend); + } else { + fp->peaks = psArraySort(fp->peaks, pmPeaksSortByRawFluxDescend); + } // XXX check for an assert on duplicates (I don't think they can happen, but Index: trunk/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 30621) +++ trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 31153) @@ -38,5 +38,8 @@ const float nsigma_delta, // how many sigma above local background a peak needs to be to survive const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance - const float min_threshold) { // minimum permitted coll height + const float min_threshold, // minimum permitted coll height + const float max_threshold, + const bool useSmoothedImage + ) { // maximum permitted coll height assert (img != NULL); assert (img->type.type == PS_TYPE_F32); assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); @@ -44,5 +47,5 @@ assert (fp != NULL); - if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do + if (fp->peaks == NULL || fp->peaks->n < 2) { // nothing to do return PS_ERR_NONE; } @@ -53,6 +56,5 @@ psImage *subImg = psImageSubset((psImage *)img, subRegion); - psImage *subWt = psImageSubset((psImage *)weight, subRegion); - assert (subImg != NULL && subWt != NULL); + psAssert (subImg, "trouble making local subimage"); psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32); @@ -68,4 +70,10 @@ psArrayAdd (brightPeaks, 128, fp->peaks->data[0]); + // XXX test point + // pmPeak *myPeak = fp->peaks->data[0]; + // if ((fabs(myPeak->x - 205) < 100) && (fabs(myPeak->y - 806) < 100)) { + // fprintf (stderr, "test peak\n"); + // } + // allocate the peakFootprint and peakFPSpans containers -- these are re-used by pmFootprintsFindAtPoint to minimize allocs in this function pmFootprint *peakFootprint = pmFootprintAlloc(fp->nspans, subImg); @@ -77,73 +85,236 @@ psImage *subMask = psImageCopy(NULL, subImg, PS_TYPE_IMAGE_MASK); - // The brightest peak is always safe; go through other peaks trying to cull them - for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop - const pmPeak *peak = fp->peaks->data[i]; - int x = peak->x - subImg->col0; - int y = peak->y - subImg->row0; - // - // Find the level nsigma below the peak that must separate the peak - // from any of its friends - // - assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); - - // const float stdev = sqrt(subWt->data.F32[y][x]); - // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev; - - const float stdev = sqrt(subWt->data.F32[y][x]); - const float flux = subImg->data.F32[y][x]; - const float fStdev = fabs(stdev/flux); - const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); - // if flux is negative, careful with fStdev - - float threshold = flux - nsigma_delta*stdev_pad; - - if (isnan(threshold)) { - // min_threshold is assumed to be below the detection threshold, - // so all the peaks are pmFootprint, and this isn't the brightest - continue; - } - - // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? - if (threshold > subImg->data.F32[y][x]) { - threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; - } - - if (threshold < min_threshold) { - threshold = min_threshold; - } - - // init peakFootprint here? - pmFootprintsFindAtPoint(peakFootprint, peakFPSpans, subImg, subMask, threshold, brightPeaks, peak->y, peak->x); - if (peakFPSpans->nStartSpans > 2000) { - // dumpfootprints(peakFootprint, peakFPSpans); - // fprintf (stderr, "big footprint %d : %d\n", peakFootprint->nspans, peakFPSpans->nStartSpans); - fprintf (stderr, "big test footprint: %f %f to %f %f (%d pix)\n", peakFootprint->bbox.x0, peakFootprint->bbox.y0, peakFootprint->bbox.x1, peakFootprint->bbox.y1, peakFootprint->npix); - } - - // at this point brightPeaks only has the peaks brighter than the current - - // we set the IDs to either 1 (in peak) or 0 (not in peak) - pmSetFootprintID (idImg, peakFootprint, IN_PEAK); - - // If this peak has not already been assigned to a source, then we can look for any - // brighter peaks within its footprint. Check if any of the previous (brighter) peaks - // are within the footprint of this peak If so, the current peak is bogus; drop it. - bool keep = true; - for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) { - const pmPeak *peak2 = fp->peaks->data[j]; - int x2 = peak2->x - subImg->col0; - int y2 = peak2->y - subImg->row0; - if (idImg->data.S32[y2][x2] == IN_PEAK) { - // There's a brighter peak within the footprint above threshold; so cull our initial peak - keep = false; - } - } - if (!keep) { - psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak"); // we should not drop any already assigned peaks. - continue; - } - - psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); + // if we have many peaks in a large footprint, we waste a lot of time generating nearly identical footprints + // here we attempt to cull peaks drawing a single footprint for all peaks in some threshold range + // fprintf (stderr, "footprint: %d x %d : %d pix, %d peaks\n", subImg->numCols, subImg->numRows, fp->npix, (int) fp->peaks->n); + if ((fp->npix > 30000) && (fp->peaks->n > 10)) { + + // max flux is above threshold for brightest peak + pmPeak *maxPeak = fp->peaks->data[0]; + float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux; + + // we have a relationship between the bin and the threshold of: + // threshold = 0.25 beta^2 bin^2 + minThreshold + // thus, the max bin is: sqrt(4.0*(maxThreshold - minThreshold)/ALPHA^2) + +# define ALPHA 0.1 + + float beta = nsigma_delta * ALPHA; + float beta2 = PS_SQR(beta); + int nBins = sqrt(4.0*(maxFlux - min_threshold)/beta2) + 10; // let's be extra generous + + // create a vector to store the threshold bins used for each peak + psVector *threshbins = psVectorAlloc(fp->peaks->n, PS_TYPE_S32); + threshbins->data.S32[0] = -1; // we skip the first peak + + /// create a vector to track if a peak has been tried or not: + psVector *peaktried = psVectorAlloc(fp->peaks->n, PS_TYPE_U8); + psVectorInit(peaktried, 0); + peaktried->data.U8[0] = true; // we skip the first peak + + psVector *threshbounds = psVectorAlloc(nBins, PS_TYPE_F32); + for (int i = 0; i < nBins; i++) { + threshbounds->data.F32[i] = 0.25*beta2*PS_SQR(i) + min_threshold; + } + psAssert(threshbounds->data.F32[threshbounds->n-1] > maxFlux, "upper limit does not include max flux"); + + psHistogram *threshist = psHistogramAllocGeneric(threshbounds); + + // assign the peaks to the histogram bins based on their nominal thresholsd + for (int i = 1; i < fp->peaks->n; i++) { + const pmPeak *peak = fp->peaks->data[i]; + float flux = useSmoothedImage ? peak->smoothFlux : peak->rawFlux; + float stdev = useSmoothedImage ? peak->smoothFluxStdev : peak->rawFluxStdev; + + // if flux is negative, careful with fStdev + const float fStdev = fabs(stdev/flux); + const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); + + float threshold = flux - nsigma_delta*stdev_pad; + psAssert(!isnan(threshold), "impossible"); + + if (threshold <= min_threshold) { + threshist->nums->data.F32[0] += 1.0; + threshbins->data.S32[i] = 0; + continue; + } + int bin = sqrt(4.0*(threshold - min_threshold)/beta2); + psAssert(bin >= 0, "impossible bin"); + + bin = PS_MIN(bin, threshist->nums->n - 1); + threshist->nums->data.F32[bin] += 1.0; + + // record the bin selected for this peak + threshbins->data.S32[i] = bin; + } + + // XXX TEST: did we assign correctly +# if (0) + int nPeaks = 1; // we don't put the brightest in the histogram + for (int i = 0; i < threshist->nums->n; i++) { + if (threshist->nums->data.F32[i] > 0) { + fprintf (stderr, "%f : %f : %d\n", threshist->bounds->data.F32[i], threshist->bounds->data.F32[i+1], (int)threshist->nums->data.F32[i]); + nPeaks += threshist->nums->data.F32[i]; + } + } + fprintf (stderr, "%d peaks vs %d in histogram\n", (int) fp->peaks->n, nPeaks); +# endif + + // XXX for the moment, we will use the alternate cull for all peaks -- it might be + // faster to use the standard process for the peaks for which the threshold bin + // contains < N sources (N ~ 5-10?) + + // loop over the threshold bins from brightest to faintest + for (int i = threshist->nums->n - 1; i >= 0; i--) { + if (threshist->nums->data.F32[i] == 0) continue; + + // we are going to examine the footprints at this threshold + float threshold = 0.5*(threshist->bounds->data.F32[i] + threshist->bounds->data.F32[i+1]); + + // generate all footprints corresponding to this threshold + psArray *myFP = pmFootprintsFind(subImg, threshold, 5); + if (!myFP) { + psWarning ("missing footprint?"); + continue; + } + if (!myFP->n) { + psWarning ("empty footprint?"); + continue; + } + + // an array to track if the footprint has a brighter peak or not + psVector *found = psVectorAlloc(myFP->n + 1, PS_TYPE_U8); + psVectorInit(found, 0); + int nFound = 0; + + // set IDs to distinguish the footprints + psImageInit(idImg, 0); + pmSetFootprintArrayIDsForImage(idImg, myFP, true); + + // check which footprints contain already-accepted (brighter) peaks + // (we can give up if/when we found a peak for all footprints) + for (int j = 0; (j < brightPeaks->n) && (nFound < found->n); j++) { + const pmPeak *peak = brightPeaks->data[j]; + int x = peak->x - subImg->col0; + int y = peak->y - subImg->row0; + int myID = idImg->data.S32[y][x]; + psAssert(myID >= 0, "impossible"); + psAssert(myID < found->n, "impossible"); + + if (myID == 0) continue; // bright peak is not in a footprint + if (found->data.U8[myID]) continue; // we already know this footprint contains a peak + found->data.U8[myID] = true; + nFound ++; + } + + // check the peaks from this threshold bin: if they land in a footprint which has + // been found, we should drop that peak. otherwise, keep it + for (int j = 0; j < fp->peaks->n; j++) { + pmPeak *peak = fp->peaks->data[j]; + if (peak->assigned) continue; // peak is already claimed by a source -- don't cull + + // skip peaks if we've already tried them + if (peaktried->data.U8[j]) continue; + + // is this peak in the threshold bin of interest? + if (threshbins->data.S32[j] != i) continue; + + // do not try this peak again + peaktried->data.U8[j] = true; + + int x = peak->x - subImg->col0; + int y = peak->y - subImg->row0; + int myID = idImg->data.S32[y][x]; + psAssert(myID < found->n, "impossible"); + + // a peak in this threshold bin without a valid footprint comes from a region + // with only a handful of pixels (1 or more from the peak itself). It probably + // cannot be joined to a neighbor + if (myID == 0) { + psArrayAdd (brightPeaks, 128, peak); + continue; + } + + // keep this peak if found->data.U8[myID] == false (no brighter peak in the footprint) + if (!found->data.U8[myID]) { + // fprintf (stderr, "keeping %d: %d,%d\n", j, peak->x, peak->y); + psArrayAdd (brightPeaks, 128, peak); + continue; + } + // fprintf (stderr, "skipping %d: %d,%d\n", j, peak->x, peak->y); + } + psFree (myFP); + psFree (found); + } + psFree (threshist); + psFree (threshbounds); + psFree (threshbins); + psFree (peaktried); + + } else { + + // The brightest peak is always safe; go through other peaks trying to cull them + for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop + const pmPeak *peak = fp->peaks->data[i]; + + float flux = useSmoothedImage ? peak->smoothFlux : peak->rawFlux; + float stdev = useSmoothedImage ? peak->smoothFluxStdev : peak->rawFluxStdev; + + // if flux is negative, careful with fStdev + const float fStdev = fabs(stdev/flux); + const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); + + float threshold = flux - nsigma_delta*stdev_pad; + + if (isnan(threshold)) { + // min_threshold is assumed to be below the detection threshold, + // so all the peaks are pmFootprint, and this isn't the brightest + continue; + } + + // just in case, force the threshold below the peak source flux + if (threshold > flux) { + threshold = flux - 10*FLT_EPSILON; + } + + if (threshold < min_threshold) { + threshold = min_threshold; + } + if (threshold > max_threshold) { + threshold = max_threshold; + } + + pmFootprintsFindAtPoint(peakFootprint, peakFPSpans, subImg, subMask, threshold, brightPeaks, peak->y, peak->x); + // if (peakFPSpans->nStartSpans > 2000) { + // // dumpfootprints(peakFootprint, peakFPSpans); + // // fprintf (stderr, "big footprint %d : %d\n", peakFootprint->nspans, peakFPSpans->nStartSpans); + // // fprintf (stderr, "big test footprint: %f %f to %f %f (%d pix)\n", peakFootprint->bbox.x0, peakFootprint->bbox.y0, peakFootprint->bbox.x1, peakFootprint->bbox.y1, peakFootprint->npix); + // } + + // at this point brightPeaks only has the peaks brighter than the current + + // we set the IDs to either 1 (in peak) or 0 (not in peak) + pmSetFootprintID (idImg, peakFootprint, IN_PEAK); + + // If this peak has not already been assigned to a source, then we can look for any + // brighter peaks within its footprint. Check if any of the previous (brighter) peaks + // are within the footprint of this peak If so, the current peak is bogus; drop it. + bool keep = true; + for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) { + // const pmPeak *peak2 = fp->peaks->data[j]; XXX isn't this an error? we only care about the kept brighter peak, right? + const pmPeak *peak2 = brightPeaks->data[j]; + int x2 = peak2->x - subImg->col0; + int y2 = peak2->y - subImg->row0; + if (idImg->data.S32[y2][x2] == IN_PEAK) { + // There's a brighter peak within the footprint above threshold; so cull our initial peak + keep = false; + } + } + if (!keep) { + psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak"); // we should not drop any already assigned peaks. + continue; + } + psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); + } } @@ -153,5 +324,4 @@ psFree(idImg); psFree(subImg); - psFree(subWt); psFree(subMask); psFree(peakFootprint); Index: trunk/psModules/src/objects/pmFootprintFindAtPoint.c =================================================================== --- trunk/psModules/src/objects/pmFootprintFindAtPoint.c (revision 30621) +++ trunk/psModules/src/objects/pmFootprintFindAtPoint.c (revision 31153) @@ -31,4 +31,8 @@ * * This is the guts of pmFootprintsFindAtPoint + * + * This function is/was ill-defined if pixel values are NAN. we should either treat NAN as > + * threshold or < threshold, but the current (r29004) code is ambiguous. + * EAM : change code so NAN is always > threshold */ bool pmFootprintSpansBuild(pmFootprint *fp, // the footprint that we're building @@ -103,5 +107,6 @@ for (int j = x0 - 1; j >= -1; j--) { double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); - if ((maskRow[j] & PM_STARTSPAN_DETECTED) || pixVal < threshold) { + bool belowThreshold = (pixVal < threshold) && isfinite(pixVal); + if ((maskRow[j] & PM_STARTSPAN_DETECTED) || belowThreshold) { if (j < x0 - 1) { // we found some pixels above threshold nx0 = j + 1; @@ -119,5 +124,6 @@ for (int j = nx0 + 1; j <= numCols; j++) { double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); - if ((maskRow[j] & PM_STARTSPAN_DETECTED) || pixVal < threshold) { + bool belowThreshold = (pixVal < threshold) && isfinite(pixVal); + if ((maskRow[j] & PM_STARTSPAN_DETECTED) || belowThreshold) { nx1 = j - 1; break; @@ -146,10 +152,12 @@ for (int j = nx1 + 1; j <= x1 + 1; j++) { double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); - if (!(maskRow[j] & PM_STARTSPAN_DETECTED) && pixVal >= threshold) { + bool aboveThreshold = (pixVal >= threshold) || !isfinite(pixVal); + if (!(maskRow[j] & PM_STARTSPAN_DETECTED) && aboveThreshold) { int sx0 = j++; // span that we're working on is sx0:sx1 int sx1 = -1; // We know that if we got here, we'll also set sx1 for (; j <= numCols; j++) { double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); - if ((maskRow[j] & PM_STARTSPAN_DETECTED) || pixVal < threshold) { // end of span + bool belowThreshold = (pixVal < threshold) && isfinite(pixVal); + if ((maskRow[j] & PM_STARTSPAN_DETECTED) || belowThreshold) { // end of span sx1 = j; break; @@ -306,5 +314,6 @@ for (i = col; i >= 0; i--) { pixVal = F32 ? imgRowF32[i] : imgRowS32[i]; - if ((maskRow[i] & PM_STARTSPAN_DETECTED) || pixVal < threshold) { + bool belowThreshold = (pixVal < threshold) && isfinite(pixVal); + if ((maskRow[i] & PM_STARTSPAN_DETECTED) || belowThreshold) { break; } @@ -313,5 +322,6 @@ for (i = col; i < numCols; i++) { pixVal = F32 ? imgRowF32[i] : imgRowS32[i]; - if ((maskRow[i] & PM_STARTSPAN_DETECTED) || pixVal < threshold) { + bool belowThreshold = (pixVal < threshold) && isfinite(pixVal); + if ((maskRow[i] & PM_STARTSPAN_DETECTED) || belowThreshold) { break; } Index: trunk/psModules/src/objects/pmFootprintIDs.c =================================================================== --- trunk/psModules/src/objects/pmFootprintIDs.c (revision 30621) +++ trunk/psModules/src/objects/pmFootprintIDs.c (revision 31153) @@ -66,5 +66,6 @@ if (footprints->n == 0) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, "You didn't provide any footprints"); + // XXX this was an error, but is an allowed condition -- NULL returned image means no footprints for any peaks + // psError(PS_ERR_BAD_PARAMETER_SIZE, true, "You didn't provide any footprints"); return NULL; } Index: trunk/psModules/src/objects/pmModel.c =================================================================== --- trunk/psModules/src/objects/pmModel.c (revision 30621) +++ trunk/psModules/src/objects/pmModel.c (revision 31153) @@ -41,4 +41,5 @@ psFree(tmp->params); psFree(tmp->dparams); + psFree(tmp->covar); psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__); } @@ -66,4 +67,5 @@ tmp->chisqNorm = NAN; tmp->nDOF = 0; + tmp->nPar = 0; tmp->nPix = 0; tmp->nIter = 0; @@ -71,4 +73,5 @@ tmp->flags = PM_MODEL_STATUS_NONE; tmp->residuals = NULL; // do not free: the model does not own this memory + tmp->covar = NULL; tmp->isPCM = false; Index: trunk/psModules/src/objects/pmModel.h =================================================================== --- trunk/psModules/src/objects/pmModel.h (revision 30621) +++ trunk/psModules/src/objects/pmModel.h (revision 31153) @@ -33,4 +33,5 @@ psVector *params; ///< Paramater values. psVector *dparams; ///< Parameter errors. + psImage *covar; ///< Optional covariance matrix float chisq; ///< Fit chi-squared. float chisqNorm; ///< re-normalized fit chi-squared. @@ -38,5 +39,6 @@ float magErr; ///< integrated model magnitude error int nPix; ///< number of pixels used for fit - int nDOF; ///< number of degrees of freedom + int nPar; ///< number of parameters in fit + int nDOF; ///< number of degrees of freedom (nDOF = nPix - nPar) int nIter; ///< number of iterations to reach min pmModelStatus flags; ///< model status flags Index: trunk/psModules/src/objects/pmModelUtils.c =================================================================== --- trunk/psModules/src/objects/pmModelUtils.c (revision 30621) +++ trunk/psModules/src/objects/pmModelUtils.c (revision 31153) @@ -107,2 +107,56 @@ return true; } + +bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments) { + + psEllipseMoments emoments; + emoments.x2 = moments->Mxx; + emoments.xy = moments->Mxy; + emoments.y2 = moments->Myy; + + // force the axis ratio to be < 20.0 + psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); + + if (!isfinite(axes.major)) return false; + if (!isfinite(axes.minor)) return false; + if (!isfinite(axes.theta)) return false; + + psEllipseShape shape = psEllipseAxesToShape (axes); + + if (!isfinite(shape.sx)) return false; + if (!isfinite(shape.sy)) return false; + if (!isfinite(shape.sxy)) return false; + + // set the shape parameters + *Sxx = PS_MAX(0.5, M_SQRT2*shape.sx); + *Syy = PS_MAX(0.5, M_SQRT2*shape.sy); + *Sxy = shape.sxy; + + return true; +} + +bool pmModelSetNorm (float *Io, pmSource *source) { + + *Io = source->peak->rawFlux; + if (!isfinite(*Io) && !source->moments) return false; + + *Io = source->moments->Peak; + if (!isfinite(*Io)) return false; + + return true; +} + +bool pmModelSetPosition (float *Xo, float *Yo, pmSource *source) { + + bool useMoments = pmSourcePositionUseMoments(source); + + if (useMoments) { + *Xo = source->moments->Mx; + *Yo = source->moments->My; + } else { + *Xo = source->peak->xf; + *Yo = source->peak->yf; + } + + return true; +} Index: trunk/psModules/src/objects/pmModelUtils.h =================================================================== --- trunk/psModules/src/objects/pmModelUtils.h (revision 30621) +++ trunk/psModules/src/objects/pmModelUtils.h (revision 31153) @@ -42,5 +42,7 @@ ); - +bool pmModelSetPosition (float *Xo, float *Yo, pmSource *source); +bool pmModelSetNorm (float *Io, pmSource *source); +bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments); /// @} Index: trunk/psModules/src/objects/pmMoments.h =================================================================== --- trunk/psModules/src/objects/pmMoments.h (revision 30621) +++ trunk/psModules/src/objects/pmMoments.h (revision 31153) @@ -51,10 +51,12 @@ int nPixels; ///< Number of pixels used. - float KronFlux; ///< Kron flux (flux in 2.5*Mrf) + float KronCore; ///< flux in r < 1.0*Mrf + float KronCoreErr; ///< error on flux in r < 1.0*Mrf + + float KronFlux; ///< Kron flux (flux in r < 2.5*Mrf) float KronFluxErr; ///< Kron flux error - float KronFinner; ///< Kron flux (flux in 2.5*Mrf) - float KronFouter; ///< Kron flux (flux in 2.5*Mrf) - + float KronFinner; ///< Kron flux (flux in 1.0*Mrf < r < 2.5*Mrf) + float KronFouter; ///< Kron flux (flux in 2.5*Mrf < r < 4.0*Mrf) } pmMoments; Index: trunk/psModules/src/objects/pmPCMdata.c =================================================================== --- trunk/psModules/src/objects/pmPCMdata.c (revision 30621) +++ trunk/psModules/src/objects/pmPCMdata.c (revision 31153) @@ -254,5 +254,6 @@ pcm->nPix = nPix; - pcm->nDOF = nPix - nParams - 1; + pcm->nPar = nParams; + pcm->nDOF = nPix - nParams; return pcm; @@ -341,5 +342,6 @@ return false; } - pcm->nDOF = pcm->nPix - nParams - 1; + pcm->nPar = nParams; + pcm->nDOF = pcm->nPix - nParams; // has the source pixel window changed? Index: trunk/psModules/src/objects/pmPCMdata.h =================================================================== --- trunk/psModules/src/objects/pmPCMdata.h (revision 30621) +++ trunk/psModules/src/objects/pmPCMdata.h (revision 31153) @@ -33,4 +33,5 @@ psMinConstraint *constraint; int nPix; + int nPar; int nDOF; } pmPCMdata; Index: trunk/psModules/src/objects/pmPSF.c =================================================================== --- trunk/psModules/src/objects/pmPSF.c (revision 30621) +++ trunk/psModules/src/objects/pmPSF.c (revision 31153) @@ -430,4 +430,9 @@ return NAN; } + + // get the model full-width at half-max + float fwhmMajor = 2*model->modelRadius (model->params, 0.5); + +# if (0) psF32 *params = model->params->data.F32; // Model parameters psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO); // Ellipse axes @@ -439,3 +444,9 @@ return fwhm; -} +# else + + psFree(model); + + return fwhmMajor; +# endif +} Index: trunk/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryFitPSF.c (revision 30621) +++ trunk/psModules/src/objects/pmPSFtryFitPSF.c (revision 31153) @@ -109,5 +109,5 @@ // This function calculates the psf and aperture magnitudes - status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal); // raw PSF mag, AP mag + status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag if (!status || isnan(source->apMag)) { psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask Index: trunk/psModules/src/objects/pmPeaks.c =================================================================== --- trunk/psModules/src/objects/pmPeaks.c (revision 30621) +++ trunk/psModules/src/objects/pmPeaks.c (revision 31153) @@ -115,4 +115,5 @@ XXX: Macro this. *****************************************************************************/ +# if (0) static bool isItInThisRegion(const psRegion valid, psS32 x, @@ -130,4 +131,5 @@ return(false); } +# endif /****************************************************************************** @@ -148,7 +150,10 @@ tmp->x = x; tmp->y = y; - tmp->value = value; - tmp->flux = value; - tmp->SN = 0; + tmp->detValue = value; + tmp->rawFlux = value; // set this by default: it is up to the user to supply a better value + tmp->rawFluxStdev = NAN; + tmp->smoothFlux = value; // set this by default: it is up to the user to supply a better value + tmp->smoothFluxStdev = NAN; + // tmp->SN = 0; tmp->xf = x; tmp->yf = y; @@ -170,9 +175,8 @@ -// psSort comparison function for peaks +// psSort comparison functions for peaks // XXX: Add error-checking for NULL args -int pmPeaksCompareAscend (const void **a, const void **b) -{ - psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__); +int pmPeaksSortByDetValueAscend (const void **a, const void **b) +{ pmPeak *A = *(pmPeak **)a; pmPeak *B = *(pmPeak **)b; @@ -180,21 +184,14 @@ psF32 diff; - diff = A->value - B->value; + diff = A->detValue - B->detValue; if (diff < FLT_EPSILON) { - psTrace("psModules.objects", 10, "---- %s(-1) end ----\n", __func__); return (-1); } else if (diff > FLT_EPSILON) { - psTrace("psModules.objects", 10, "---- %s(+1) end ----\n", __func__); return (+1); } - psTrace("psModules.objects", 10, "---- %s(0) end ----\n", __func__); return (0); } - -// psSort comparison function for peaks -// XXX: Add error-checking for NULL args -int pmPeaksCompareDescend (const void **a, const void **b) -{ - psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__); +int pmPeaksSortByDetValueDescend (const void **a, const void **b) +{ pmPeak *A = *(pmPeak **)a; pmPeak *B = *(pmPeak **)b; @@ -202,47 +199,104 @@ psF32 diff; - diff = A->value - B->value; + diff = A->detValue - B->detValue; if (diff < FLT_EPSILON) { - psTrace("psModules.objects", 10, "---- %s(+1) end ----\n", __func__); return (+1); } else if (diff > FLT_EPSILON) { - psTrace("psModules.objects", 10, "---- %s(-1) end ----\n", __func__); return (-1); } - psTrace("psModules.objects", 10, "---- %s(0) end ----\n", __func__); return (0); } - -// sort by SN (descending) -int pmPeakSortBySN (const void **a, const void **b) +int pmPeaksSortByRawFluxAscend (const void **a, const void **b) { pmPeak *A = *(pmPeak **)a; pmPeak *B = *(pmPeak **)b; - psF32 fA = A->SN; - psF32 fB = B->SN; - if (isnan (fA)) fA = 0; - if (isnan (fB)) fB = 0; - - psF32 diff = fA - fB; - if (diff > FLT_EPSILON) return (-1); - if (diff < FLT_EPSILON) return (+1); + psF32 diff; + + diff = A->rawFlux - B->rawFlux; + if (diff < FLT_EPSILON) { + return (-1); + } else if (diff > FLT_EPSILON) { + return (+1); + } return (0); } - -// sort by Y (ascending) -int pmPeakSortByY (const void **a, const void **b) +int pmPeaksSortByRawFluxDescend (const void **a, const void **b) { pmPeak *A = *(pmPeak **)a; pmPeak *B = *(pmPeak **)b; - psF32 fA = A->y; - psF32 fB = B->y; - - psF32 diff = fA - fB; - if (diff > FLT_EPSILON) return (+1); - if (diff < FLT_EPSILON) return (-1); + psF32 diff; + + diff = A->rawFlux - B->rawFlux; + if (diff < FLT_EPSILON) { + return (+1); + } else if (diff > FLT_EPSILON) { + return (-1); + } return (0); } +int pmPeaksSortBySmoothFluxAscend (const void **a, const void **b) +{ + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 diff; + + diff = A->smoothFlux - B->smoothFlux; + if (diff < FLT_EPSILON) { + return (-1); + } else if (diff > FLT_EPSILON) { + return (+1); + } + return (0); +} +int pmPeaksSortBySmoothFluxDescend (const void **a, const void **b) +{ + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 diff; + + diff = A->smoothFlux - B->smoothFlux; + if (diff < FLT_EPSILON) { + return (+1); + } else if (diff > FLT_EPSILON) { + return (-1); + } + return (0); +} + +// // sort by SN (descending) +// int pmPeakSortBySN (const void **a, const void **b) +// { +// pmPeak *A = *(pmPeak **)a; +// pmPeak *B = *(pmPeak **)b; +// +// psF32 fA = A->flux; +// psF32 fB = B->flux; +// if (isnan (fA)) fA = 0; +// if (isnan (fB)) fB = 0; +// +// psF32 diff = fA - fB; +// if (diff > FLT_EPSILON) return (-1); +// if (diff < FLT_EPSILON) return (+1); +// return (0); +// } +// +// // sort by Y (ascending) +// int pmPeakSortByY (const void **a, const void **b) +// { +// pmPeak *A = *(pmPeak **)a; +// pmPeak *B = *(pmPeak **)b; +// +// psF32 fA = A->y; +// psF32 fB = B->y; +// +// psF32 diff = fA - fB; +// if (diff > FLT_EPSILON) return (+1); +// if (diff < FLT_EPSILON) return (-1); +// return (0); +// } /****************************************************************************** @@ -554,28 +608,2 @@ return(list); } - -// return a new array of peaks which are in the valid region and below threshold -// XXX this function is unused and probably could be dropped -psArray *pmPeaksSubset( - psArray *peaks, - psF32 maxValue, - const psRegion valid) -{ - psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(peaks, NULL); - - psArray *output = psArrayAllocEmpty (200); - - psTrace ("psModules.objects", 3, "list size is %ld\n", peaks->n); - - for (int i = 0; i < peaks->n; i++) { - pmPeak *tmpPeak = (pmPeak *) peaks->data[i]; - if (tmpPeak->value > maxValue) - continue; - if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)) - continue; - psArrayAdd (output, 200, tmpPeak); - } - psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__); - return(output); -} Index: trunk/psModules/src/objects/pmPeaks.h =================================================================== --- trunk/psModules/src/objects/pmPeaks.h (revision 30621) +++ trunk/psModules/src/objects/pmPeaks.h (revision 31153) @@ -35,4 +35,5 @@ PM_PEAK_EDGE, ///< Peak on edge. PM_PEAK_FLAT, ///< Peak has equal-value neighbors. + PM_PEAK_SUSPECT_SATURATION, ///< Peak is probably saturated PM_PEAK_UNDEF ///< Undefined. } pmPeakType; @@ -45,4 +46,12 @@ * associated with the source: * + * There are 3 values which define the amplitude of the peak and which may be used to sort the + * peaks: + * * detValue - the peak in the detection image (nominally, the S/N) + * * rawFlux - the peak in the unsmoothed image + * * smoothFlux - the peak in the smoothed image + * + * For a given image, peaks do necesarily not have the same sequence for these three values. + * Depending on the analysis, it may make sense to sort by one or the other of these values */ typedef struct @@ -55,7 +64,10 @@ float dx; ///< bicube fit error on peak coord (x) float dy; ///< bicube fit error on peak coord (y) - float value; ///< level in detection image - float flux; ///< level in unsmoothed sci image - float SN; ///< S/N implied by detection level + float detValue; ///< peak flux in detection image (= S/N) + float rawFlux; ///< peak flux in unsmoothed signal image + float rawFluxStdev; ///< peak stdev in unsmoothed signal image + float smoothFlux; ///< peak flux in smoothed signal image + float smoothFluxStdev; ///< peak stdev in smoothed signal image + // float SNestimated; ///< S/N estimated from the detection image bool assigned; ///< is peak assigned to a source? pmPeakType type; ///< Description of peak. @@ -136,9 +148,10 @@ ); -int pmPeaksCompareAscend (const void **a, const void **b); -int pmPeaksCompareDescend (const void **a, const void **b); - -int pmPeakSortBySN (const void **a, const void **b); -int pmPeakSortByY (const void **a, const void **b); +int pmPeaksSortByDetValueAscend (const void **a, const void **b); +int pmPeaksSortByDetValueDescend (const void **a, const void **b); +int pmPeaksSortByRawFluxAscend (const void **a, const void **b); +int pmPeaksSortByRawFluxDescend (const void **a, const void **b); +int pmPeaksSortBySmoothFluxAscend (const void **a, const void **b); +int pmPeaksSortBySmoothFluxDescend (const void **a, const void **b); /// @} Index: trunk/psModules/src/objects/pmPhotObj.c =================================================================== --- trunk/psModules/src/objects/pmPhotObj.c (revision 30621) +++ trunk/psModules/src/objects/pmPhotObj.c (revision 31153) @@ -67,9 +67,9 @@ return false; } - if (!finite(source->peak->xf)) { + if (!isfinite(source->peak->xf)) { psError(PS_ERR_UNKNOWN, true, "NAN peak coordinate"); return false; } - if (!finite(source->peak->yf)) { + if (!isfinite(source->peak->yf)) { psError(PS_ERR_UNKNOWN, true, "NAN peak coordinate"); return false; @@ -81,7 +81,7 @@ object->x = source->peak->xf; object->y = source->peak->yf; - object->SN = source->peak->SN; + object->flux = source->peak->rawFlux; } else { - object->SN = PS_MAX(object->SN, source->peak->SN); + object->flux = PS_MAX(object->flux, source->peak->rawFlux); } psArrayAdd (object->sources, 1, source); @@ -89,12 +89,12 @@ } -// sort by SN (descending) -int pmPhotObjSortBySN (const void **a, const void **b) +// sort by flux (descending) +int pmPhotObjSortByFlux (const void **a, const void **b) { pmPhotObj *objA = *(pmPhotObj **)a; pmPhotObj *objB = *(pmPhotObj **)b; - psF32 fA = objA->SN; - psF32 fB = objB->SN; + psF32 fA = objA->flux; + psF32 fB = objB->flux; psF32 diff = fA - fB; Index: trunk/psModules/src/objects/pmPhotObj.h =================================================================== --- trunk/psModules/src/objects/pmPhotObj.h (revision 30621) +++ trunk/psModules/src/objects/pmPhotObj.h (revision 31153) @@ -34,5 +34,5 @@ float x; float y; - float SN; // max of peak->SN for all matched sources + float flux; // max of peak->rawFlux for all matched sources } pmPhotObj; @@ -41,5 +41,5 @@ bool pmPhotObjAddSource(pmPhotObj *object, pmSource *source); -int pmPhotObjSortBySN (const void **a, const void **b); +int pmPhotObjSortByFlux (const void **a, const void **b); int pmPhotObjSortByX (const void **a, const void **b); Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 30621) +++ trunk/psModules/src/objects/pmSource.c (revision 31153) @@ -171,9 +171,12 @@ // peak has the same values as the original if (in->peak != NULL) { - source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->value, in->peak->type); + source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->detValue, in->peak->type); source->peak->xf = in->peak->xf; source->peak->yf = in->peak->yf; - source->peak->flux = in->peak->flux; - source->peak->SN = in->peak->SN; + source->peak->rawFlux = in->peak->rawFlux; + source->peak->rawFluxStdev = in->peak->rawFluxStdev; + source->peak->smoothFlux = in->peak->smoothFlux; + source->peak->smoothFluxStdev = in->peak->smoothFluxStdev; + // source->peak->SN = in->peak->SN; } @@ -328,8 +331,4 @@ *****************************************************************************/ -// XXX EAM include a S/N cutoff in selecting the sources? -// XXX this function should probably accept the values, not a recipe. wrap with a -// psphot-specific function which applies the recipe values -// only apply selection to sources within specified region pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 AR_MAX) { @@ -429,4 +428,5 @@ psfClump.nSigma = stats->sampleStdev; + psfClump.nTotal = nValid; if (savedImage) { @@ -465,13 +465,13 @@ psStats *stats = NULL; - // select the single highest peak - psArraySort (peaks, pmPeaksCompareDescend); + // select the single highest peak (note that we only have detValue, not rawFlux, etc + psArraySort (peaks, pmPeaksSortByDetValueDescend); clump = peaks->data[0]; - psTrace ("psModules.objects", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->value); + psTrace ("psModules.objects", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->detValue); // XXX store the mean sigma? float meanSigma = psfClump.nSigma; - psfClump.nStars = clump->value; - psfClump.nSigma = clump->value / meanSigma; + psfClump.nStars = clump->detValue; + psfClump.nSigma = clump->detValue / meanSigma; // define section window for clump @@ -643,4 +643,22 @@ } + // check for insignificant sources or excessively low-surface brightness + float coreSN = source->moments->KronCore / source->moments->KronCoreErr; + float coreKR = source->moments->KronCore / source->moments->KronFlux; + + // XXX these values need to be in the recipe... + if (false && isfinite(coreSN) && (coreSN < 5.0)) { + source->type = PM_SOURCE_TYPE_DEFECT; + source->mode |= PM_SOURCE_MODE_DEFECT; + Ncr ++; + continue; + } + if (false && isfinite(coreKR) && (coreKR < 0.1)) { + source->type = PM_SOURCE_TYPE_DEFECT; + source->mode |= PM_SOURCE_MODE_DEFECT; + Ncr ++; + continue; + } + // likely unsaturated extended source (too large to be stellar) if (sigX > clump.X + 3*clump.dX || sigY > clump.Y + 3*clump.dY) { @@ -651,7 +669,8 @@ // the rest are probable stellar objects + // the vectors below are accumulated to give user feedback on the S/N ranges starsn_moments->data.F32[starsn_moments->n] = source->moments->SN; starsn_moments->n ++; - starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN; + starsn_peaks->data.F32[starsn_peaks->n] = sqrt(source->peak->detValue); starsn_peaks->n ++; Nstar ++; @@ -1149,12 +1168,33 @@ } +// this function decides if the source position should be based on the peak or the moments. +// this is only used if we know we should not use a model fit position (eg, no model, or no +// model yet) +bool pmSourcePositionUseMoments(pmSource *source) { + + if (!source->moments) return false; // can't if there are no moments + if (!source->moments->nPixels) return false; // can't if the moments were not measured + if (source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE) return false; // can't if the moments failed... + + if (source->mode & PM_SOURCE_MODE_SATSTAR) return true; // moments are best for SATSTARs + + float dX = source->moments->Mx - source->peak->xf; + float dY = source->moments->My - source->peak->yf; + float dR = hypot(dX, dY); + + // only use the moments position if the moment-peak offset is small or the star is saturated + if (dR > 1.5) return false; + + return true; +} + // sort by SN (descending) -int pmSourceSortBySN (const void **a, const void **b) +int pmSourceSortByFlux (const void **a, const void **b) { pmSource *A = *(pmSource **)a; pmSource *B = *(pmSource **)b; - psF32 fA = (A->peak == NULL) ? 0 : A->peak->SN; - psF32 fB = (B->peak == NULL) ? 0 : B->peak->SN; + psF32 fA = (A->peak == NULL) ? 0 : A->peak->rawFlux; + psF32 fB = (B->peak == NULL) ? 0 : B->peak->rawFlux; if (isnan (fA)) fA = 0; if (isnan (fB)) fB = 0; @@ -1166,4 +1206,23 @@ } +// sort by SN (descending) +int pmSourceSortByParentFlux (const void **a, const void **b) +{ + pmSource *Ao = *(pmSource **)a; + pmSource *Bo = *(pmSource **)b; + pmSource *A = Ao->parent; + pmSource *B = Bo->parent; + + psF32 fA = (A->peak == NULL) ? 0 : A->peak->rawFlux; + psF32 fB = (B->peak == NULL) ? 0 : B->peak->rawFlux; + if (isnan (fA)) fA = 0; + if (isnan (fB)) fB = 0; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (-1); + if (diff < FLT_EPSILON) return (+1); + return (0); +} + // sort by Y (ascending) int pmSourceSortByY (const void **a, const void **b) @@ -1201,4 +1260,21 @@ pmSource *A = *(pmSource **)a; pmSource *B = *(pmSource **)b; + + int iA = A->seq; + int iB = B->seq; + + int diff = iA - iB; + if (diff > 0) return (+1); + if (diff < 0) return (-1); + return (0); +} + +// sort by Seq (ascending) +int pmSourceSortByParentSeq (const void **a, const void **b) +{ + pmSource *Ao = *(pmSource **)a; + pmSource *Bo = *(pmSource **)b; + pmSource *A = Ao->parent; + pmSource *B = Bo->parent; int iA = A->seq; Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 30621) +++ trunk/psModules/src/objects/pmSource.h (revision 31153) @@ -121,4 +121,5 @@ float dY; int nStars; + int nTotal; float nSigma; } @@ -286,8 +287,12 @@ bool pmSourceCachePSF (pmSource *source, psImageMaskType maskVal); -int pmSourceSortBySN (const void **a, const void **b); +bool pmSourcePositionUseMoments(pmSource *source); + int pmSourceSortByY (const void **a, const void **b); int pmSourceSortByX (const void **a, const void **b); int pmSourceSortBySeq (const void **a, const void **b); +int pmSourceSortByParentSeq (const void **a, const void **b); +int pmSourceSortByFlux (const void **a, const void **b); +int pmSourceSortByParentFlux (const void **a, const void **b); pmSourceMode pmSourceModeFromString (const char *name); Index: trunk/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 31153) @@ -55,4 +55,5 @@ psFree(radial->flux); psFree(radial->fluxErr); + psFree(radial->fluxStdev); psFree(radial->fill); } @@ -65,4 +66,5 @@ radial->flux = NULL; radial->fluxErr = NULL; + radial->fluxStdev = NULL; radial->fill = NULL; return radial; Index: trunk/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 30621) +++ trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 31153) @@ -23,4 +23,5 @@ typedef struct { psVector *flux; // fluxes measured at above radii + psVector *fluxStdev; // scatter (standard deviation) of flux psVector *fluxErr; // formal error on the fluxes (sqrt\sum(variance)) psVector *fill; // angles corresponding to above radial profiles @@ -63,4 +64,5 @@ float petrosianR50; float petrosianR50Err; + float petrosianFill; } pmSourceExtendedPars; Index: trunk/psModules/src/objects/pmSourceFitModel.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceFitModel.c (revision 31153) @@ -40,4 +40,5 @@ #include "pmSourceDiffStats.h" #include "pmSource.h" +#include "pmSourcePhotometry.h" #include "pmSourceFitModel.h" @@ -59,4 +60,5 @@ opt->maxChisqDOF = NAN; opt->poissonErrors = true; + opt->saveCovariance = false; return opt; @@ -231,12 +233,24 @@ psTrace ("psModules.objects", 4, "%f +/- %f", params->data.F32[i], dparams->data.F32[i]); } + if (options->saveCovariance) { + model->covar = psMemIncrRefCounter(covar); + } + model->nIter = myMin->iter; + model->nPar = nParams; + psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); // save the resulting chisq, nDOF, nIter - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nPix = y->n; - model->nDOF = y->n - nParams; - model->chisqNorm = model->chisq / model->nDOF; + // NOTE: if (!options->poissonErrors) chisq will be wrong : recalculate + if (options->poissonErrors) { + model->chisq = myMin->value; + model->nPix = y->n; + model->nDOF = y->n - model->nPar; + model->chisqNorm = model->chisq / model->nDOF; + } else { + pmSourceChisqUnsubtracted (source, model, maskVal); + } + + // set the model success or failure status model->flags |= PM_MODEL_STATUS_FITTED; if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE; Index: trunk/psModules/src/objects/pmSourceFitModel.h =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.h (revision 30621) +++ trunk/psModules/src/objects/pmSourceFitModel.h (revision 31153) @@ -31,5 +31,7 @@ float maxChisqDOF; ///< convergence criterion float weight; ///< use this weight for constant-weight fits + float covarFactor; ///< covariance factor for calculating the chisq bool poissonErrors; ///< use poisson errors for fits? + bool saveCovariance; } pmSourceFitOptions; Index: trunk/psModules/src/objects/pmSourceFitPCM.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitPCM.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceFitPCM.c (revision 31153) @@ -38,4 +38,5 @@ #include "pmSourceDiffStats.h" #include "pmSource.h" +#include "pmSourcePhotometry.h" #include "pmSourceFitModel.h" #include "pmPCMdata.h" @@ -70,4 +71,8 @@ psTrace ("psModules.objects", 4, "%f +/- %f", params->data.F32[i], dparams->data.F32[i]); } + if (fitOptions->saveCovariance) { + psFree(pcm->modelConv->covar); + pcm->modelConv->covar = psMemIncrRefCounter(covar); + } psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); @@ -79,11 +84,19 @@ } } + pcm->modelConv->nIter = myMin->iter; + pcm->modelConv->nPar = pcm->nPar; // save the resulting chisq, nDOF, nIter - pcm->modelConv->chisq = myMin->value; - pcm->modelConv->nIter = myMin->iter; - pcm->modelConv->nPix = pcm->nPix; - pcm->modelConv->nDOF = pcm->nDOF; - pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF; + if (fitOptions->poissonErrors) { + pcm->modelConv->chisq = myMin->value; + pcm->modelConv->nPix = pcm->nPix; + pcm->modelConv->nDOF = pcm->nDOF; + pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF; + } else { + // xxx this is wrong because it does not convolve with the psf + pmSourceChisqUnsubtracted (source, pcm->modelConv, maskVal); + } + + // set the model success or failure status pcm->modelConv->flags |= PM_MODEL_STATUS_FITTED; if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE; Index: trunk/psModules/src/objects/pmSourceFitSet.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitSet.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceFitSet.c (revision 31153) @@ -39,4 +39,5 @@ #include "pmSourceDiffStats.h" #include "pmSource.h" +#include "pmSourcePhotometry.h" #include "pmSourceFitModel.h" @@ -296,7 +297,8 @@ // set the model parameters for this fit set -bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, - const psVector *param, pmSource *source, - psMinimization *myMin, int nPix, bool fitStatus) +bool pmSourceFitSetValues (pmSourceFitSetData *set, + const psVector *dparam, const psVector *param, const psImage *covar, + pmSource *source, psMinimization *myMin, int nPix, + bool fitStatus, pmSourceFitOptions *options, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(set, false); @@ -311,4 +313,5 @@ int n = 0; + int nStart = 0; for (int i = 0; i < set->paramSet->n; i++) { @@ -320,11 +323,29 @@ psTrace ("psModules.objects", 4, "%f +/- %f", param->data.F32[n], dparam->data.F32[n]); } + if (options->saveCovariance) { + // we only save the covar matrix for this object with itself (ignore cross terms between objects) + model->covar = psImageAlloc(model->params->n, model->params->n, PS_TYPE_F32); + for (int ix = 0; ix < model->params->n; ix++) { + for (int iy = 0; iy < model->params->n; iy++) { + model->covar->data.F32[iy][ix] = covar->data.F32[nStart+iy][nStart+ix]; + } + } + } + nStart += model->params->n; psTrace ("psModules.objects", 4, " src %d", i); + + model->nIter = myMin->iter; + // model->nPar is set by pmSourceFitSetMasks // save the resulting chisq, nDOF, nIter // these are not unique for any one source - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nDOF = nPix - model->params->n; + if (options->poissonErrors) { + model->chisq = myMin->value; + model->nPix = nPix; + model->nDOF = nPix - model->nPar; + model->chisqNorm = model->chisq / model->nDOF; + } else { + pmSourceChisqUnsubtracted (source, model, maskVal); + } // set the model success or failure status @@ -380,4 +401,5 @@ for (int i = 0; i < set->paramSet->n; i++) { psVector *paramOne = set->paramSet->data[i]; + pmModel *modelOne = set->modelSet->data[i]; switch (mode) { @@ -387,4 +409,5 @@ if (j == PM_PAR_I0) continue; constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + j] = 1; + modelOne->nPar = 1; } break; @@ -396,4 +419,5 @@ if (j == PM_PAR_I0) continue; constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + j] = 1; + modelOne->nPar = 3; } break; @@ -401,4 +425,5 @@ // EXT model fits all params (except sky) constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + PM_PAR_SKY] = 1; + modelOne->nPar = paramOne->n - 1; break; default: @@ -550,5 +575,5 @@ } -// parameter errors from the covariance matrix + // parameter errors from the covariance matrix psVector *dparams = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32); for (int i = 0; i < dparams->n; i++) { @@ -558,5 +583,5 @@ } -// get the Gauss-Newton distance for fixed model parameters + // get the Gauss-Newton distance for fixed model parameters if (constraint->paramMask != NULL) { psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32); @@ -580,5 +605,5 @@ } - pmSourceFitSetValues (thisSet, dparams, params, source, myMin, y->n, fitStatus); + pmSourceFitSetValues (thisSet, dparams, params, covar, source, myMin, y->n, fitStatus, options, maskVal); psTrace ("psModules.objects", 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nPix: %ld\n", onPic, fitStatus, myMin->iter, myMin->value, y->n); Index: trunk/psModules/src/objects/pmSourceFitSet.h =================================================================== --- trunk/psModules/src/objects/pmSourceFitSet.h (revision 30621) +++ trunk/psModules/src/objects/pmSourceFitSet.h (revision 31153) @@ -40,5 +40,9 @@ bool pmSourceFitSetJoin (psVector *deriv, psVector *param, pmSourceFitSetData *set); bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param); -bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, const psVector *param, pmSource *source, psMinimization *myMin, int nPix, bool fitStatus); + +bool pmSourceFitSetValues (pmSourceFitSetData *set, + const psVector *dparam, const psVector *param, const psImage *covar, + pmSource *source, psMinimization *myMin, int nPix, + bool fitStatus, pmSourceFitOptions *options, psImageMaskType maskVal); psF32 pmSourceFitSetFunction(psVector *deriv, const psVector *param, const psVector *x); @@ -57,5 +61,5 @@ psArray *modelSet, ///< model to be fitted pmSourceFitOptions *options, ///< define options for fitting process - psImageMaskType maskVal ///< Vale to mask + psImageMaskType maskVal ///< Value to mask ); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 31153) @@ -49,4 +49,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -62,29 +63,9 @@ PS_ASSERT_PTR_NON_NULL(extname, false); - int i; psArray *table; psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 errMag, chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); + // if the sequence is defined, write these in seq order; otherwise @@ -94,5 +75,5 @@ if (source->seq == -1) { // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); } else { sources = psArraySort (sources, pmSourceSortBySeq); @@ -100,9 +81,16 @@ } + short nImageOverlap; + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); + table = psArrayAllocEmpty (sources->n); // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (i = 0; i < sources->n; i++) { + for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; @@ -115,55 +103,10 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - errMag = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); pmSourceDiffStats diffStats; @@ -175,48 +118,44 @@ row = psMetadataAlloc (); psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfFlux); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", chisq); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float mxx = source->moments ? source->moments->Mxx : NAN; - float mxy = source->moments ? source->moments->Mxy : NAN; - float myy = source->moments ? source->moments->Myy : NAN; - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); - - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); @@ -331,14 +270,10 @@ // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; + source->peak->rawFlux = peakFlux; + source->peak->smoothFlux = peakFlux; source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position source->peak->dx = dPAR[PM_PAR_XPOS]; source->peak->dy = dPAR[PM_PAR_YPOS]; - if (isfinite (source->errMag) && (source->errMag > 0.0)) { - source->peak->SN = 1.0 / source->errMag; - } else { - source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N - } source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); @@ -353,4 +288,7 @@ source->moments = pmMomentsAlloc (); + source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); @@ -395,5 +333,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -568,5 +506,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 31153) @@ -49,4 +49,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -62,45 +63,8 @@ PS_ASSERT_PTR_NON_NULL(extname, false); - int i; psArray *table; psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); - - // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors - float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ"); - if (!status1) { - fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1"); - if (!status1) { - fwhmMajor = 5.0; // XXX just a guess! - } - } - float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN"); - if (!status1) { - fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2"); - if (!status1) { - fwhmMinor = 5.0; // XXX just a guess! - } - } // if the sequence is defined, write these in seq order; otherwise @@ -110,5 +74,5 @@ if (source->seq == -1) { // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); } else { sources = psArraySort (sources, pmSourceSortBySeq); @@ -118,7 +82,14 @@ table = psArrayAllocEmpty (sources->n); + short nImageOverlap; + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); + // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (i = 0; i < sources->n; i++) { + for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; @@ -131,65 +102,10 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // in linear-fit mode, there is no error on the centroid - xErr = fwhmMajor * source->errMag / 2.35; - yErr = fwhmMinor * source->errMag / 2.35; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - } else { - bool useMoments = true; - useMoments = (useMoments && source->moments); // can't if there are no moments - useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured - useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... - - if (useMoments) { - xPos = source->moments->Mx; - yPos = source->moments->My; - xErr = fwhmMajor * source->errMag / 2.35; - yErr = fwhmMinor * source->errMag / 2.35; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - } - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); pmSourceDiffStats diffStats; @@ -201,10 +117,10 @@ row = psMetadataAlloc (); psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); @@ -214,57 +130,44 @@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in real aperture", source->apMagRaw); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", chisq); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float mxx = source->moments ? source->moments->Mxx : NAN; - float mxy = source->moments ? source->moments->Mxy : NAN; - float myy = source->moments ? source->moments->Myy : NAN; - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); - - float Mrf = source->moments ? source->moments->Mrf : NAN; - float Mrh = source->moments ? source->moments->Mrh : NAN; - float Krf = source->moments ? source->moments->KronFlux : NAN; - float dKrf = source->moments ? source->moments->KronFluxErr : NAN; - - float Kinner = source->moments ? source->moments->KronFinner : NAN; - float Kouter = source->moments ? source->moments->KronFouter : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", dKrf); - - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", Kouter); - - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.Kinner); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", moments.Kouter); + + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); @@ -385,14 +288,10 @@ // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; + source->peak->rawFlux = peakFlux; + source->peak->smoothFlux = peakFlux; source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position source->peak->dx = dPAR[PM_PAR_XPOS]; source->peak->dy = dPAR[PM_PAR_YPOS]; - if (isfinite (source->errMag) && (source->errMag > 0.0)) { - source->peak->SN = 1.0 / source->errMag; - } else { - source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N - } source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); @@ -408,4 +307,7 @@ source->moments = pmMomentsAlloc (); + source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); @@ -456,5 +358,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -629,5 +531,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 31153) @@ -50,4 +50,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -66,34 +67,14 @@ psArray *table; psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 errMag, chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); - - // if the sequence is defined, write these in seq order; otherwise - // write them in S/N order: + + // if the sequence is defined, write these in seq order; otherwise write them in S/N order. + // Careful: if we are working with child sources, then we need to sort by the parent info, + // not our info if (sources->n > 0) { - pmSource *source = (pmSource *) sources->data[0]; + pmSource *source = sources->data[0]; if (source->seq == -1) { - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); } else { sources = psArraySort (sources, pmSourceSortBySeq); @@ -103,189 +84,89 @@ table = psArrayAllocEmpty (sources->n); -# if (0) - // we use this just to define the output vectors (which must be present for all objects) - bool status = false; - psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - psAssert (radMax, "this must have been defined and tested earlier!"); - psAssert (radMax->n, "this must have been defined and tested earlier!"); - psAssert (radMin->n == radMax->n, "inconsistent annular bins"); - - // write the radial profile apertures to header - for (int i = 0; i < radMax->n; i++) { - sprintf (keyword1, "RMIN_%02d", i); - sprintf (keyword2, "RMAX_%02d", i); - psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); - psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); - } -# endif + short nImageOverlap; + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - - // If source->seq is -1, source was generated in this analysis. If source->seq is - // not -1, source was read from elsewhere: in the latter case, preserve the source - // ID. source.seq is used instead of source.id since the latter is a const - // generated on Alloc, and would thus be wrong for read in sources. + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + // If source->seq is -1, source is unique and generated in this analysis. If + // source->seq is not -1, source was read from elsewhere or tied to other source (eg + // from another image): in the latter case, preserve the source ID. source.seq is used + // instead of source.id since the latter is a const generated on Alloc, and would thus + // be wrong for read in sources. if (source->seq == -1) { source->seq = i; } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - errMag = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); row = psMetadataAlloc (); psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", chisq); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float Mxx = source->moments ? source->moments->Mxx : NAN; - float Mxy = source->moments ? source->moments->Mxy : NAN; - float Myy = source->moments ? source->moments->Myy : NAN; - - float Mrf = source->moments ? source->moments->Mrf : NAN; - float Mrh = source->moments ? source->moments->Mrh : NAN; - float Krf = source->moments ? source->moments->KronFlux : NAN; - float dKrf = source->moments ? source->moments->KronFluxErr : NAN; - - float Kinner = source->moments ? source->moments->KronFinner : NAN; - float Kouter = source->moments ? source->moments->KronFouter : NAN; - - float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN; - float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN; - float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN; - float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", M_c3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", M_s3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", M_c4); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", M_s4); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", dKrf); - - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Kouter); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); - -# if (0) - // XXX if we have raw radial apertures, write them out here - psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVectorInit (radFlux, NAN); - psVectorInit (radFluxErr, NAN); - psVectorInit (radFill, NAN); - if (!source->radial) goto empty_annuli; - if (!source->radial->flux) goto empty_annuli; - if (!source->radial->fill) goto empty_annuli; - psAssert (source->radial->flux->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (source->radial->fill->n <= radFlux->n, "inconsistent vector lengths"); - - // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux - for (int j = 0; j < source->radial->flux->n; j++) { - radFlux->data.F32[j] = source->radial->flux->data.F32[j]; - radFluxErr->data.F32[j] = source->radial->fluxErr->data.F32[j]; - radFill->data.F32[j] = source->radial->fill->data.F32[j]; - } - - empty_annuli: - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); - psFree (radFlux); - psFree (radFluxErr); - psFree (radFill); -# endif // XXX not sure how to get this : need to load Nimages with weight? @@ -406,14 +287,10 @@ // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; + source->peak->rawFlux = peakFlux; + source->peak->smoothFlux = peakFlux; source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position source->peak->dx = dPAR[PM_PAR_XPOS]; source->peak->dy = dPAR[PM_PAR_YPOS]; - if (isfinite (source->errMag) && (source->errMag > 0.0)) { - source->peak->SN = 1.0 / source->errMag; - } else { - source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N - } source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); @@ -429,4 +306,7 @@ source->moments = pmMomentsAlloc (); + source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); @@ -500,5 +380,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -522,5 +402,9 @@ // we write out all sources, regardless of quality. the source flags tell us the state for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // skip sources without measurements @@ -556,5 +440,5 @@ } psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO", PS_DATA_F32, "Axial Ratio of radial profile", AxialRatio); - psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA", PS_DATA_F32, "Angle of radial profile ellipse", AxialTheta); + psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA", PS_DATA_F32, "Angle of radial profile ellipse", AxialTheta); // Petrosian measurements @@ -567,21 +451,23 @@ float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", mag); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", magErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", mag); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", magErr); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill); } else { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", NAN); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", NAN); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", NAN); } } @@ -672,10 +558,15 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams int nParamMax = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + if (source->modelFits == NULL) continue; for (int j = 0; j < source->modelFits->n; j++) { @@ -691,5 +582,8 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // XXX if no model fits are saved, write out modelEXT? @@ -747,14 +641,21 @@ if (k < model->params->n) { - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]); } else { - psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", NAN); } } - // XXX other parameters which may be set. - // XXX flag / value to define the model - // XXX write out the model type, fit status flags - + // optionally, write out the covariance matrix values + // XXX do I need to pad this to match the biggest covar matrix? + if (model->covar) { + for (int iy = 0; iy < model->covar->numCols; iy++) { + for (int ix = iy; ix < model->covar->numCols; ix++) { + snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix); + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]); + + } + } + } psArrayAdd (table, 100, row); psFree (row); @@ -790,5 +691,5 @@ bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { - + bool status = false; psArray *table; psMetadata *row; @@ -796,4 +697,10 @@ char keyword1[80], keyword2[80]; + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) { + psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n"); + return true; + } + // create a header to hold the output data psMetadata *outhead = psMetadataAlloc (); @@ -803,5 +710,4 @@ // we use this just to define the output vectors (which must be present for all objects) - bool status = false; psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); @@ -818,12 +724,9 @@ } + // the FWHM values are available if we measured a psf-matched convolved set psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); - if (!fwhmValues) { - psError (PM_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout"); - return false; - } // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -832,11 +735,16 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // skip sources without radial aper measurements (or insufficient) - if (source->radialAper == NULL) continue; - psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set"); - - for (int entry = 0; entry < fwhmValues->n; entry++) { + if (source->radialAper == NULL) continue; + + // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set"); + + for (int entry = 0; entry < source->radialAper->n; entry++) { // choose the convolved EXT model, if available, otherwise the simple one @@ -844,10 +752,5 @@ assert (radialAper); - bool useMoments = true; - useMoments = (useMoments && source->moments); // can't if there are no moments - useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured - useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... - - if (useMoments) { + if (pmSourcePositionUseMoments(source)) { xPos = source->moments->Mx; yPos = source->moments->My; @@ -863,10 +766,15 @@ psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER", 0, "Center of aperture measurements", xPos); psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER", 0, "Center of aperture measurements", yPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "FWHM of matched PSF", fwhmValues->data.F32[entry]); + if (fwhmValues) { + psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "FWHM of matched PSF", fwhmValues->data.F32[entry]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "image is not FWHM-matched", NAN); + } // XXX if we have raw radial apertures, write them out here - psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32); psVectorInit (radFlux, NAN); psVectorInit (radFluxErr, NAN); @@ -879,15 +787,18 @@ // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux for (int j = 0; j < radialAper->flux->n; j++) { - radFlux->data.F32[j] = radialAper->flux->data.F32[j]; - radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j]; - radFill->data.F32[j] = radialAper->fill->data.F32[j]; + radFlux->data.F32[j] = radialAper->flux->data.F32[j]; + radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j]; + radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j]; + radFill->data.F32[j] = radialAper->fill->data.F32[j]; } write_annuli: - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation", radFluxStdev); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); psFree (radFlux); psFree (radFluxErr); + psFree (radFluxStdev); psFree (radFill); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 31153) @@ -49,4 +49,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -62,27 +63,6 @@ psArray *table; psMetadata *row; - int i; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 errMag, chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); // if the sequence is defined, write these in seq order; otherwise @@ -92,5 +72,5 @@ if (source->seq == -1) { // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); } else { sources = psArraySort (sources, pmSourceSortBySeq); @@ -100,7 +80,14 @@ table = psArrayAllocEmpty (sources->n); + short nImageOverlap; + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); + // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (i = 0; i < sources->n; i++) { + for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; @@ -113,94 +100,45 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - errMag = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); row = psMetadataAlloc (); psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", outputs.ra); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", outputs.dec); + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", chisq); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float mxx = source->moments ? source->moments->Mxx : NAN; - float mxy = source->moments ? source->moments->Mxy : NAN; - float myy = source->moments ? source->moments->Myy : NAN; - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); @@ -317,14 +255,10 @@ // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; + source->peak->rawFlux = peakFlux; + source->peak->smoothFlux = peakFlux; source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position source->peak->dx = dPAR[PM_PAR_XPOS]; source->peak->dy = dPAR[PM_PAR_YPOS]; - if (isfinite (source->errMag) && (source->errMag > 0.0)) { - source->peak->SN = 1.0 / source->errMag; - } else { - source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N - } source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); @@ -339,4 +273,7 @@ source->moments = pmMomentsAlloc (); + source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); @@ -371,5 +308,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -549,5 +486,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 31153) @@ -49,4 +49,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -62,26 +63,6 @@ psArray *table; psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 errMag, chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); // if the sequence is defined, write these in seq order; otherwise @@ -91,5 +72,5 @@ if (source->seq == -1) { // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); } else { sources = psArraySort (sources, pmSourceSortBySeq); @@ -98,4 +79,11 @@ table = psArrayAllocEmpty (sources->n); + + short nImageOverlap; + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state @@ -112,94 +100,45 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // in linear-fit mode, there is no error on the centroid - xErr = source->peak->dx; - yErr = source->peak->dy; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - errMag = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); row = psMetadataAlloc (); psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", chisq); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float mxx = source->moments ? source->moments->Mxx : NAN; - float mxy = source->moments ? source->moments->Mxy : NAN; - float myy = source->moments ? source->moments->Myy : NAN; - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", myy); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); @@ -322,14 +261,10 @@ // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; + source->peak->rawFlux = peakFlux; + source->peak->smoothFlux = peakFlux; source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position source->peak->dx = dPAR[PM_PAR_XPOS]; source->peak->dy = dPAR[PM_PAR_YPOS]; - if (isfinite (source->errMag) && (source->errMag > 0.0)) { - source->peak->SN = 1.0 / source->errMag; - } else { - source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N - } source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); @@ -344,4 +279,7 @@ source->moments = pmMomentsAlloc (); + source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); @@ -392,5 +330,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -611,5 +549,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 31153) @@ -49,4 +49,5 @@ #include "pmSourceIO.h" +#include "pmSourceOutputs.h" // panstarrs-style FITS table output (header + table in 1st extension) @@ -62,42 +63,6 @@ psArray *table; psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - psF32 chisq, apRadius; - psS32 nPix, nDOF; pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - bool status1 = false; - bool status2 = false; - float magOffset = NAN; - float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); - - // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors - float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ"); - if (!status1) { - fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1"); - if (!status1) { - fwhmMajor = 5.0; // XXX just a guess! - } - } - float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN"); - if (!status1) { - fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2"); - if (!status1) { - fwhmMinor = 5.0; // XXX just a guess! - } - } // if the sequence is defined, write these in seq order; otherwise @@ -107,5 +72,5 @@ if (source->seq == -1) { // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); } else { sources = psArraySort (sources, pmSourceSortBySeq); @@ -114,4 +79,11 @@ table = psArrayAllocEmpty (sources->n); + + short nImageOverlap; + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state @@ -128,73 +100,19 @@ } - // no difference between PSF and non-PSF model - pmModel *model = source->modelPSF; - - if (model != NULL) { - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - xErr = fwhmMajor * source->errMag / 2.35; - yErr = fwhmMinor * source->errMag / 2.35; - } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); - } else { - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - } - chisq = model->chisq; - nDOF = model->nDOF; - nPix = model->nPix; - apRadius = source->apRadius; - } else { - bool useMoments = true; - useMoments = (useMoments && source->moments); // can't if there are no moments - useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured - useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... - - if (useMoments) { - xPos = source->moments->Mx; - yPos = source->moments->My; - xErr = fwhmMajor * source->errMag / 2.35; - yErr = fwhmMinor * source->errMag / 2.35; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - xErr = source->peak->dx; - yErr = source->peak->dy; - } - axes.major = NAN; - axes.minor = NAN; - axes.theta = NAN; - chisq = NAN; - nDOF = 0; - nPix = 0; - apRadius = NAN; - } - - float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; - psS16 nImageOverlap = 1; - - psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - float posAngle = 0.0; - float pltScale = 0.0; - pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); row = psMetadataAlloc (); psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); @@ -203,66 +121,48 @@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); - - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); + + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); // NOTE: RA & DEC (both double) need to be on an 8-byte boundary... - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", ptSky.r*PS_DEG_RAD); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", ptSky.d*PS_DEG_RAD); - - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); + + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", chisq); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); - - // distinguish moments measure from window vs S/N > XX ?? - float Mxx = source->moments ? source->moments->Mxx : NAN; - float Mxy = source->moments ? source->moments->Mxy : NAN; - float Myy = source->moments ? source->moments->Myy : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", Myy); - - float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN; - float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN; - float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN; - float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", M_c3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", M_s3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", M_c4); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", M_s4); - - float Mrf = source->moments ? source->moments->Mrf : NAN; - float Mrh = source->moments ? source->moments->Mrh : NAN; - float Krf = source->moments ? source->moments->KronFlux : NAN; - float dKrf = source->moments ? source->moments->KronFluxErr : NAN; - - float Kinner = source->moments ? source->moments->KronFinner : NAN; - float Kouter = source->moments ? source->moments->KronFouter : NAN; - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", dKrf); - - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", Kouter); - - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); + + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); + // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.KronCore); + // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR", PS_DATA_F32, "Kron Error (in 1.0 R1)", moments.KronCoreErr); + psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); @@ -384,14 +284,17 @@ // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; + source->peak->rawFlux = peakFlux; + source->peak->smoothFlux = peakFlux; source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position source->peak->dx = dPAR[PM_PAR_XPOS]; source->peak->dy = dPAR[PM_PAR_YPOS]; - if (isfinite (source->errMag) && (source->errMag > 0.0)) { - source->peak->SN = 1.0 / source->errMag; - } else { - source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N - } + + // we no longer sort by S/N, only flux + // if (isfinite (source->errMag) && (source->errMag > 0.0)) { + // source->peak->SN = 1.0 / source->errMag; + // } else { + // source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N + // } source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); @@ -407,4 +310,7 @@ source->moments = pmMomentsAlloc (); + source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); @@ -477,5 +383,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -649,5 +555,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams Index: trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 31153) @@ -91,5 +91,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -136,5 +136,5 @@ float calMag = source->psfMag + magOffset; - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; + float peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN; psS16 nImageOverlap = 1; @@ -294,5 +294,6 @@ source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; + source->peak->rawFlux = peakFlux; + source->peak->smoothFlux = peakFlux; source->peak->dx = dPAR[PM_PAR_XPOS]; source->peak->dy = dPAR[PM_PAR_YPOS]; @@ -357,5 +358,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -576,5 +577,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 31153) @@ -101,5 +101,5 @@ } - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; + float peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN; psS16 nImageOverlap = 1; psS32 ID = 0; // XXX need to figure out how to generate this @@ -220,5 +220,6 @@ source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; + source->peak->rawFlux = peakFlux; + source->peak->smoothFlux = peakFlux; source->peak->dx = dPAR[PM_PAR_XPOS]; source->peak->dy = dPAR[PM_PAR_YPOS]; Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 31153) @@ -73,5 +73,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -117,5 +117,5 @@ } - float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; + float peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN; psS16 nImageOverlap = 1; @@ -263,5 +263,6 @@ source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; + source->peak->rawFlux = peakFlux; + source->peak->smoothFlux = peakFlux; source->peak->dx = dPAR[PM_PAR_XPOS]; source->peak->dy = dPAR[PM_PAR_YPOS]; @@ -307,5 +308,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); table = psArrayAllocEmpty (sources->n); @@ -480,5 +481,5 @@ // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams Index: trunk/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_RAW.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_RAW.c (revision 31153) @@ -126,5 +126,5 @@ source[0].apMag, source[0].type, source[0].mode, logChi, logChiNorm, - source[0].peak->SN, + source[0].peak->rawFlux, source[0].apRadius, source[0].pixWeightNotBad, @@ -188,5 +188,5 @@ log10(model[0].chisq/model[0].nDOF), log10(model[0].chisqNorm/model[0].nDOF), - source[0].peak->SN, + source[0].peak->rawFlux, source[0].apRadius, source[0].pixWeightNotBad, @@ -234,5 +234,5 @@ fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %8.1f %7.1f %7.1f %7.1f %4d %2d\n", - source->peak->x, source->peak->y, source->peak->value, + source->peak->x, source->peak->y, source->peak->detValue, moment->Mx, moment->My, moment->Mxx, moment->Myy, @@ -270,5 +270,5 @@ continue; fprintf (f, "%5d %5d %8.1f %7.1f %7.1f %6.3f %6.3f %10.1f %7.1f %7.1f %7.1f %4d %2d %#5x\n", - source->peak->x, source->peak->y, source->peak->value, + source->peak->x, source->peak->y, source->peak->detValue, source->moments->Mx, source->moments->My, source->moments->Mxx, source->moments->Myy, @@ -299,8 +299,8 @@ continue; fprintf (f, "%5d %5d %7.1f %7.2f %7.2f %7.2f\n", - peak->x, peak->y, peak->value, peak->SN, peak->xf, peak->yf); - } - fclose (f); - return true; -} - + peak->x, peak->y, peak->detValue, peak->rawFlux, peak->xf, peak->yf); + } + fclose (f); + return true; +} + Index: trunk/psModules/src/objects/pmSourceIO_SMPDATA.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 31153) @@ -205,5 +205,4 @@ source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; sources->data[i] = source; Index: trunk/psModules/src/objects/pmSourceMasks.h =================================================================== --- trunk/psModules/src/objects/pmSourceMasks.h (revision 30621) +++ trunk/psModules/src/objects/pmSourceMasks.h (revision 31153) @@ -46,4 +46,9 @@ PM_SOURCE_MODE2_DIFF_WITH_DOUBLE = 0x00000002, ///< diff source matched to positive detections in both images PM_SOURCE_MODE2_MATCHED = 0x00000004, ///< diff source matched to positive detections in both images + + PM_SOURCE_MODE2_ON_SPIKE = 0x00000008, ///< > 25% of (PSF-weighted) pixels land on diffraction spike + PM_SOURCE_MODE2_ON_STARCORE = 0x00000010, ///< > 25% of (PSF-weighted) pixels land on starcore + PM_SOURCE_MODE2_ON_BURNTOOL = 0x00000020, ///< > 25% of (PSF-weighted) pixels land on burntool + PM_SOURCE_MODE2_ON_CONVPOOR = 0x00000040, ///< > 25% of (PSF-weighted) pixels land on convpoor } pmSourceMode2; Index: trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- trunk/psModules/src/objects/pmSourceMoments.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceMoments.c (revision 31153) @@ -82,4 +82,6 @@ psF32 Sum = 0.0; psF32 Var = 0.0; + psF32 SumCore = 0.0; + psF32 VarCore = 0.0; psF32 R2 = PS_SQR(radius); psF32 minSN2 = PS_SQR(minSN); @@ -248,4 +250,7 @@ int nKronPix = 0; + int nCorePix = 0; + int nInner = 0; + int nOuter = 0; Sum = Var = 0.0; float SumInner = 0.0; @@ -255,5 +260,5 @@ psF32 yDiff = row - yCM; - if (fabs(yDiff) > radKron) continue; + if (fabs(yDiff) > radKouter) continue; psF32 *vPix = source->pixels->data.F32[row]; @@ -272,5 +277,5 @@ psF32 xDiff = col - xCM; - if (fabs(xDiff) > radKron) continue; + if (fabs(xDiff) > radKouter) continue; // radKron is just a function of (xDiff, yDiff) @@ -293,16 +298,28 @@ } + // use sigma (fixed by psf) not a radKron based value + if (r < sigma) { + SumCore += pDiff; + VarCore += wDiff; + nCorePix ++; + } + if ((r > radKinner) && (r < radKron)) { SumInner += pDiff; + nInner ++; } if ((r > radKron) && (r < radKouter)) { SumOuter += pDiff; + nOuter ++; } } } - source->moments->KronFlux = Sum; - source->moments->KronFinner = SumInner; - source->moments->KronFouter = SumOuter; - source->moments->KronFluxErr = sqrt(Var); + // *** should I rescale these fluxes by pi R^2 / nNpix? + source->moments->KronCore = SumCore * M_PI * PS_SQR(sigma) / nCorePix; + source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix; + source->moments->KronFlux = Sum * M_PI * PS_SQR(radKron) / nKronPix; + source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix; + source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron) - PS_SQR(radKinner)) / nInner; + source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) - PS_SQR(radKron)) / nOuter; psTrace ("psModules.objects", 4, "Mrf: %f KronFlux: %f Mxx: %f Mxy: %f Myy: %f Mxxx: %f Mxxy: %f Mxyy: %f Myyy: %f Mxxxx: %f Mxxxy: %f Mxxyy: %f Mxyyy: %f Mxyyy: %f\n", @@ -313,5 +330,5 @@ psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f My: %f Sum: %f Mxx: %f Mxy: %f Myy: %f sky: %f Npix: %d\n", - source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx, source->moments->My, Sum, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, sky, source->moments->nPixels); + source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx, source->moments->My, Sum, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, sky, source->moments->nPixels); return(true); Index: trunk/psModules/src/objects/pmSourceOutputs.c =================================================================== --- trunk/psModules/src/objects/pmSourceOutputs.c (revision 31153) +++ trunk/psModules/src/objects/pmSourceOutputs.c (revision 31153) @@ -0,0 +1,181 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" + +#include "pmSourceIO.h" +#include "pmSourceOutputs.h" + +bool pmSourceOutputsCommonValues (short *nImageOverlap, float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header) { + + pmFPA *fpa = readout->parent->parent->parent; + + bool status1 = false; + bool status2 = false; + float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); + float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); + if (!isfinite(zeropt)) { + zeropt = psMetadataLookupF32 (&status2, header, "ZPT_OBS"); + } + if (status1 && status2 && (exptime > 0.0)) { + *magOffset = zeropt + 2.5*log10(exptime); + } + *zeroptErr = psMetadataLookupF32 (&status2, header, "ZPT_ERR"); + + // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors + *fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ"); + if (!status1) { + *fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1"); + if (!status1) { + *fwhmMajor = 5.0; // XXX just a guess! + } + } + *fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN"); + if (!status1) { + *fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2"); + if (!status1) { + *fwhmMinor = 5.0; // XXX just a guess! + } + } + + *nImageOverlap = psMetadataLookupS32 (&status2, header, "NINPUTS"); + return true; + +} + +// what is the correct postion? +// * if we have a PSF model fit, use PM_PAR_XPOS,YPOS for X_PSF,Y_PSF +// * if we do not have a model: +// ** if we have moments: +// *** if the star is saturated, use the moments +// *** if the moments and peak agree to < DR, use the moments +// *** otherwise, use the peak + +bool pmSourceOutputsSetValues (pmSourceOutputs *outputs, pmSource *source, pmChip *chip, float fwhmMajor, float fwhmMinor, float magOffset) { + + psF32 *PAR, *dPAR; + psEllipseAxes axes; + + // no difference between PSF and non-PSF model + pmModel *model = source->modelPSF; + + if (model != NULL) { + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + outputs->xPos = PAR[PM_PAR_XPOS]; + outputs->yPos = PAR[PM_PAR_YPOS]; + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { + outputs->xErr = dPAR[PM_PAR_XPOS]; + outputs->yErr = dPAR[PM_PAR_YPOS]; + } else { + outputs->xErr = fwhmMajor * source->errMag / 2.35; + outputs->yErr = fwhmMinor * source->errMag / 2.35; + } + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { + axes = pmPSF_ModelToAxes (PAR, 20.0); + outputs->psfMajor = axes.major; + outputs->psfMinor = axes.minor; + outputs->psfTheta = axes.theta*PS_DEG_RAD; + } else { + outputs->psfMajor = NAN; + outputs->psfMinor = NAN; + outputs->psfTheta = NAN; + } + outputs->chisq = model->chisq; + outputs->nDOF = model->nDOF; + outputs->nPix = model->nPix; + outputs->apRadius = source->apRadius; + } else { + bool useMoments = pmSourcePositionUseMoments(source); + + if (useMoments) { + outputs->xPos = source->moments->Mx; + outputs->yPos = source->moments->My; + outputs->xErr = fwhmMajor * source->errMag / 2.35; + outputs->yErr = fwhmMinor * source->errMag / 2.35; + } else { + outputs->xPos = source->peak->xf; + outputs->yPos = source->peak->yf; + outputs->xErr = source->peak->dx; + outputs->yErr = source->peak->dy; + } + outputs->psfMajor = NAN; + outputs->psfMinor = NAN; + outputs->psfTheta = NAN; + outputs->chisq = NAN; + outputs->nDOF = 0; + outputs->nPix = 0; + outputs->apRadius = NAN; + } + + outputs->calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; + outputs->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN; + + psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + float posAngle = 0.0; + float pltScale = 0.0; + pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, outputs->xPos, outputs->yPos); + + outputs->posAngle = posAngle*PS_DEG_RAD; + outputs->pltScale = pltScale*PS_DEG_RAD*3600.0; + + outputs->ra = ptSky.r*PS_DEG_RAD; + outputs->dec = ptSky.d*PS_DEG_RAD; + + return true; +} + +bool pmSourceOutputsSetMoments (pmSourceOutputsMoments *moments, pmSource *source) { + + // distinguish moments measure from window vs S/N > XX ?? + moments->Mxx = source->moments ? source->moments->Mxx : NAN; + moments->Mxy = source->moments ? source->moments->Mxy : NAN; + moments->Myy = source->moments ? source->moments->Myy : NAN; + moments->M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN; + moments->M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN; + moments->M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN; + moments->M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN; + moments->Mrf = source->moments ? source->moments->Mrf : NAN; + moments->Mrh = source->moments ? source->moments->Mrh : NAN; + moments->Krf = source->moments ? source->moments->KronFlux : NAN; + moments->dKrf = source->moments ? source->moments->KronFluxErr : NAN; + moments->Kinner = source->moments ? source->moments->KronFinner : NAN; + moments->Kouter = source->moments ? source->moments->KronFouter : NAN; + moments->KronCore = source->moments ? source->moments->KronCore : NAN; + moments->KronCoreErr = source->moments ? source->moments->KronCoreErr : NAN; + + return true; +} Index: trunk/psModules/src/objects/pmSourceOutputs.h =================================================================== --- trunk/psModules/src/objects/pmSourceOutputs.h (revision 31153) +++ trunk/psModules/src/objects/pmSourceOutputs.h (revision 31153) @@ -0,0 +1,62 @@ +/* @file pmSourceOutputs.h + * @brief functions to perform common I/O conversions + * + * @author EAM, IfA + * + * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-16 22:30:50 $ + * Copyright 2011 Institute for Astronomy, University of Hawaii + * + */ + +# ifndef PM_SOURCE_OUTPUTS_H +# define PM_SOURCE_OUTPUTS_H + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +typedef struct { + float xPos; + float yPos; + float xErr; + float yErr; + double ra; + double dec; + float psfMajor; + float psfMinor; + float psfTheta; + float chisq; + float apRadius; + int nPix; + int nDOF; + float calMag; + float peakMag; + float posAngle; + float pltScale; +} pmSourceOutputs; + +typedef struct { + float Mxx; + float Mxy; + float Myy; + float M_c3; + float M_s3; + float M_c4; + float M_s4; + float Mrf; + float Mrh; + float Krf; + float dKrf; + float Kinner; + float Kouter; + float KronCore; + float KronCoreErr; +} pmSourceOutputsMoments; + +bool pmSourceOutputsCommonValues (short *nImageOverlap, float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header); + +bool pmSourceOutputsSetValues (pmSourceOutputs *outputs, pmSource *source, pmChip *chip, float fwhmMajor, float fwhmMinor, float magOffset); + +bool pmSourceOutputsSetMoments (pmSourceOutputsMoments *moments, pmSource *source); + +# endif Index: trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.c (revision 30621) +++ trunk/psModules/src/objects/pmSourcePhotometry.c (revision 31153) @@ -23,4 +23,5 @@ #include "pmFPAMaskWeight.h" +#include "pmConfigMask.h" #include "pmTrend2D.h" #include "pmResiduals.h" @@ -49,10 +50,26 @@ static float AP_MIN_SN = 0.0; -bool pmSourceMagnitudesInit (psMetadata *config) -{ - PS_ASSERT_PTR_NON_NULL(config, false); +// make this a bit more clever and dynamic +static psImageMaskType maskSuspect = 0; +static psImageMaskType maskSpike = 0; +static psImageMaskType maskStarCore = 0; +static psImageMaskType maskBurntool = 0; +static psImageMaskType maskConvPoor = 0; + +bool pmSourceMagnitudesInit (pmConfig *config, psMetadata *recipe) +{ + PS_ASSERT_PTR_NON_NULL(recipe, false); bool status; - float limit = psMetadataLookupF32 (&status, config, "AP_MIN_SN"); + // we are going to test specially against these poor values + if (config) { + maskSpike = pmConfigMaskGet("SPIKE", config); + maskStarCore = pmConfigMaskGet("STARCORE", config); + maskBurntool = pmConfigMaskGet("BURNTOOL", config); + maskConvPoor = pmConfigMaskGet("CONV.POOR", config); + maskSuspect = maskSpike | maskStarCore | maskBurntool | maskConvPoor; + } + + float limit = psMetadataLookupF32 (&status, recipe, "AP_MIN_SN"); if (status) { AP_MIN_SN = limit; @@ -77,5 +94,5 @@ // XXX masked region should be (optionally) elliptical // if mode is PM_SOURCE_PHOT_PSFONLY, we skip all other magnitudes -bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal) +bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal, float radius) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -166,5 +183,5 @@ // measure the contribution of included pixels if (mode & PM_SOURCE_PHOT_WEIGHT) { - pmSourcePixelWeight (&source->pixWeightNotBad, &source->pixWeightNotPoor, model, source->maskObj, maskVal, markVal); + pmSourcePixelWeight (source, model, source->maskObj, maskVal, radius); } @@ -342,8 +359,10 @@ // return source aperture magnitude -bool pmSourcePixelWeight (float *pixWeightNotBad, float *pixWeightNotPoor, pmModel *model, psImage *mask, psImageMaskType maskVal, psImageMaskType markVal) -{ - PS_ASSERT_PTR_NON_NULL(pixWeightNotBad, false); - PS_ASSERT_PTR_NON_NULL(pixWeightNotPoor, false); +bool pmSourcePixelWeight (pmSource *source, pmModel *model, psImage *mask, psImageMaskType maskVal, float radius) +{ + PS_ASSERT_PTR_NON_NULL(source, false); + source->pixWeightNotBad = NAN; + source->pixWeightNotPoor = NAN; + PS_ASSERT_PTR_NON_NULL(mask, false); PS_ASSERT_PTR_NON_NULL(model, false); @@ -355,10 +374,14 @@ float value; + float spikeSum = 0; + float starcoreSum = 0; + float burntoolSum = 0; + float convpoorSum = 0; + int Xo, Yo, dP; int dX, DX, NX; int dY, DY, NY; - *pixWeightNotBad = 0.0; - *pixWeightNotPoor = 0.0; + float radius2 = PS_SQR(radius); // we only care about the value of the object model, not the local sky @@ -387,11 +410,19 @@ NY = mask->numRows; + psImageMaskType maskBad = maskVal; + maskBad &= ~maskSuspect; + // measure modelSum and validSum. this function is applied to a sources' subimage. the // value of DX is chosen (see above) to cover the full possible size of the subimage if it // were not by an edge; ie, if the source is cut in half by an image edge, we correctly // count the virtual pixels off the edge in normalizing the value of the pixWeight + + // we skip any pixels [real or virtual] outside of the specified radius (nominally the aperture radius) for (int ix = -DX; ix < DX + 1; ix++) { + if (ix > radius) continue; int mx = ix + dX; for (int iy = -DY; iy < DY + 1; iy++) { + if (iy > radius) continue; + if (ix*ix + iy*iy > radius2) continue; int my = iy + dY; @@ -409,16 +440,53 @@ if (my >= NY) continue; - if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)) { + // count pixels which are masked only with bad pixels + if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskBad)) { notBadSum += value; } - if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & ~markVal)) { + + // count pixels which are masked with an mask bit (bad or poor) + if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)) { notPoorSum += value; } + + // count pixels which are masked with an mask bit (bad or poor) + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskSpike) { + spikeSum += value; + } + // count pixels which are masked with an mask bit (bad or poor) + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskStarCore) { + starcoreSum += value; + } + // count pixels which are masked with an mask bit (bad or poor) + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskBurntool) { + burntoolSum += value; + } + // count pixels which are masked with an mask bit (bad or poor) + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskConvPoor) { + convpoorSum += value; + } } } psFree (coord); - *pixWeightNotBad = notBadSum / modelSum; - *pixWeightNotPoor = notPoorSum / modelSum; + source->pixWeightNotBad = notBadSum / modelSum; + source->pixWeightNotPoor = notPoorSum / modelSum; + + if ((spikeSum/modelSum) > 0.25) { + source->mode2 |= PM_SOURCE_MODE2_ON_SPIKE; + } + if ((starcoreSum/modelSum) > 0.25) { + source->mode2 |= PM_SOURCE_MODE2_ON_STARCORE; + } + if ((burntoolSum/modelSum) > 0.25) { + source->mode2 |= PM_SOURCE_MODE2_ON_BURNTOOL; + } + if ((convpoorSum/modelSum) > 0.25) { + source->mode2 |= PM_SOURCE_MODE2_ON_CONVPOOR; + } + + if (isfinite(source->pixWeightNotBad) && isfinite(source->pixWeightNotPoor)) { + psAssert (source->pixWeightNotBad <= source->pixWeightNotPoor, "error: all bad pixels should also be poor"); + } return (true); @@ -427,5 +495,5 @@ # define FLUX_LIMIT 3.0 -// measure stats that may be using in difference images for distinguishing real sources from bad residuals +// measure stats that may be used in difference images for distinguishing real sources from bad residuals bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) { @@ -673,6 +741,6 @@ # endif -// determine chisq, etc for linear normalization-only fit -bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor) +// determine chisq, nPix, nDOF, chisqNorm : model->nPar must be set +bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(model, false); @@ -689,15 +757,66 @@ if (variance->data.F32[j][i] <= 0) continue; - dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]); + dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i]; Npix ++; } } model->nPix = Npix; - model->nDOF = Npix - 1; + model->nDOF = Npix - model->nPar; model->chisq = dC; + model->chisqNorm = dC / model->nDOF; return (true); } + +// return source aperture magnitude +bool pmSourceChisqUnsubtracted (pmSource *source, pmModel *model, psImageMaskType maskVal) +{ + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(model, false); + + float dC = 0.0; + int Npix = 0; + + // the model function returns the source flux at a position + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + psVector *params = model->params; + psImage *image = source->pixels; + psImage *mask = source->maskObj; + psImage *variance = source->variance; + + int dX = image->col0; + int dY = image->row0; + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + + // skip pixels which are masked + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue; + + if (variance->data.F32[iy][ix] <= 0) continue; + + coord->data.F32[0] = (psF32) ix + dX + 0.5; + coord->data.F32[1] = (psF32) iy + dY + 0.5; + + // for the full model, add all points + float value = model->modelFunc (NULL, params, coord); + + // fprintf (stderr, "%d, %d : %f, %f : %f - %f : %f\n", + // ix, iy, coord->data.F32[0], coord->data.F32[1], image->data.F32[iy][ix], value, dC); + + dC += PS_SQR (image->data.F32[iy][ix] - value) / variance->data.F32[iy][ix]; + Npix ++; + } + } + model->nPix = Npix; + model->nDOF = Npix - model->nPar; + model->chisq = dC; + model->chisqNorm = dC / model->nDOF; + + psFree (coord); + return (true); +} double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal) Index: trunk/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.h (revision 30621) +++ trunk/psModules/src/objects/pmSourcePhotometry.h (revision 31153) @@ -64,10 +64,11 @@ ); -bool pmSourceMagnitudesInit (psMetadata *config); -bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal); +bool pmSourceMagnitudesInit (pmConfig *config, psMetadata *recipe); +bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal, float radius); -bool pmSourcePixelWeight (float *pixWeightNotBad, float *pixWeightNotPoor, pmModel *model, psImage *mask, psImageMaskType maskVal, psImageMaskType markVal); +bool pmSourcePixelWeight (pmSource *source, pmModel *model, psImage *mask, psImageMaskType maskVal, float radius); -bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor); +bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal); +bool pmSourceChisqUnsubtracted (pmSource *source, pmModel *model, psImageMaskType maskVal); bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal); Index: trunk/psModules/src/objects/pmSourceVisual.c =================================================================== --- trunk/psModules/src/objects/pmSourceVisual.c (revision 30621) +++ trunk/psModules/src/objects/pmSourceVisual.c (revision 31153) @@ -56,13 +56,5 @@ if (!pmVisualTestLevel("psphot.psf.metric", 2)) return true; - - if (kapa1 == -1) { - kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots"); - if (kapa1 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + if (!pmVisualInitWindow (&kapa1, "pmSource:plots")) return false; KapaClearSections (kapa1); @@ -140,16 +132,9 @@ if (!pmVisualTestLevel("psphot.psf.subpix", 3)) return true; - - if (kapa1 == -1) { - kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots"); - if (kapa1 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + if (!pmVisualInitWindow (&kapa1, "pmSource:plots")) return false; KapaClearSections (kapa1); KapaInitGraph (&graphdata); + section.bg = KapaColorByName ("none"); // XXX probably should be 'none' int n; @@ -302,13 +287,5 @@ if (!pmVisualTestLevel("psphot.psf.fits", 2)) return true; - - if (kapa2 == -1) { - kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + if (!pmVisualInitWindow (&kapa2, "pmSource:images")) return false; // create images 1/10 scale: @@ -385,13 +362,5 @@ if (!source->pixels) return false; if (!source->modelFlux) return false; - - if (kapa2 == -1) { - kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + if (!pmVisualInitWindow (&kapa2, "pmSource:images")) return false; // KapaClearSections (kapa2); @@ -428,16 +397,9 @@ if (!plotPSF) return true; if (!pmVisualTestLevel("psphot.psf.resid", 2)) return true; - - if (kapa1 == -1) { - kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots"); - if (kapa1 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + if (!pmVisualInitWindow (&kapa1, "pmSource:plots")) return false; KapaClearPlots (kapa1); KapaInitGraph (&graphdata); + section.bg = KapaColorByName ("none"); // XXX probably should be 'none' float Xmin = +1e32;