Index: /branches/eam_branches/ipp-20100621/psModules/src/camera/pmFPAfileIO.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/camera/pmFPAfileIO.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/camera/pmFPAfileIO.c (revision 28643) @@ -24,17 +24,26 @@ #include "pmFPAfileFitsIO.h" #include "pmFPAfileFringeIO.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" #include "pmPSF.h" -#include "pmModel.h" -#include "pmSource.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" -#include "pmResiduals.h" #include "pmPSF_IO.h" + #include "pmAstrometryModel.h" #include "pmAstrometryRefstars.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/camera/pmReadoutFake.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/camera/pmReadoutFake.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/camera/pmReadoutFake.c (revision 28643) @@ -10,19 +10,26 @@ #include "pmFPA.h" -#include "pmMoments.h" + +#include "pmTrend2D.h" #include "pmResiduals.h" #include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.h" -#include "pmModel.h" -#include "pmModelClass.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 "pmSourceUtils.h" -#include "pmModelUtils.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" + #include "pmSourceGroups.h" - #include "pmReadoutFake.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/camera/pmReadoutFake.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/camera/pmReadoutFake.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/camera/pmReadoutFake.h (revision 28643) @@ -2,14 +2,14 @@ #define PM_READOUT_FAKE_H -#include -#include -#include - -#include -#include -#include -#include -#include -#include +// #include +// #include +// #include +// +// #include +// #include +// #include +// #include +// #include +// #include /// Set threading Index: /branches/eam_branches/ipp-20100621/psModules/src/config/pmConfig.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/config/pmConfig.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/config/pmConfig.c (revision 28643) @@ -655,7 +655,9 @@ if ((argNum = psArgumentGet(*argc, argv, "-visual-all"))) { pmVisualSetLevel(".", 10); + psArgumentRemove(argNum, argc, argv); } if ((argNum = psArgumentGet(*argc, argv, "-visual-levels"))) { pmVisualPrintLevels(stdout); + psArgumentRemove(argNum, argc, argv); } Index: /branches/eam_branches/ipp-20100621/psModules/src/detrend/pmMaskStats.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/detrend/pmMaskStats.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/detrend/pmMaskStats.c (revision 28643) @@ -7,5 +7,4 @@ #include -#include #include "pmHDU.h" @@ -16,4 +15,5 @@ #include "pmFPAAstrometry.h" +#include "pmMaskStats.h" #define ESCAPE { \ Index: /branches/eam_branches/ipp-20100621/psModules/src/extras/pmVisual.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/extras/pmVisual.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/extras/pmVisual.c (revision 28643) @@ -21,16 +21,24 @@ #include "pmAstrometryObjects.h" #include "pmSubtractionStamps.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 "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" #include "pmPSF.h" #include "pmPSFtry.h" -#include "pmSource.h" + #include "pmFPAExtent.h" - -# if (HAVE_KAPA) -# include -# include "pmVisual.h" -# include "pmKapaPlots.h" -# define KAPAX 700 -# define KAPAY 700 #include "pmAstrometryVisual.h" @@ -38,4 +46,12 @@ #include "pmStackVisual.h" #include "pmSourceVisual.h" + +# if (HAVE_KAPA) +# include +#include "pmVisual.h" +#include "pmKapaPlots.h" + +# define KAPAX 700 +# define KAPAY 700 //#define TESTING @@ -102,5 +118,5 @@ /* Wait up to 1.0 second for a response, then continue */ - timeout.tv_sec = 2; + timeout.tv_sec = 10; timeout.tv_usec = 0; Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmPSFEnvelope.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmPSFEnvelope.c (revision 28643) @@ -10,23 +10,26 @@ #include "pmHDU.h" #include "pmFPA.h" -#include "pmReadoutFake.h" - -#include "pmMoments.h" + +#include "pmTrend2D.h" #include "pmResiduals.h" #include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.h" -#include "pmModel.h" -#include "pmModelClass.h" -#include "pmModelUtils.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 "pmSourceUtils.h" #include "pmSourceFitModel.h" +#include "pmPSF.h" #include "pmPSFtry.h" - +#include "pmReadoutFake.h" #include "pmPSFEnvelope.h" #include "pmStackVisual.h" @@ -410,5 +413,12 @@ options->chiFluxTrend = false; // All sources have similar flux, so fitting a trend often fails - pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true); + // options which modify the behavior of the model fitting + options->fitOptions = pmSourceFitOptionsAlloc(); + options->fitOptions->nIter = SOURCE_FIT_ITERATIONS; + options->fitOptions->tol = 0.01; + options->fitOptions->poissonErrors = true; + options->fitOptions->weight = VARIANCE_VAL; + options->fitOptions->mode = PM_SOURCE_FIT_PSF; + pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT); // Important for getting a good stack target PSF Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c (revision 28643) @@ -7,18 +7,23 @@ #include +#include "pmErrorCodes.h" +#include "pmHDU.h" +#include "pmFPA.h" + // All these includes required to get stamps out of an array of pmSources -#include "pmMoments.h" +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" -#include "pmResiduals.h" -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" #include "pmModel.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" #include "pmSource.h" -#include "pmErrorCodes.h" #include "pmSubtraction.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/Makefile.am (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/Makefile.am (revision 28643) @@ -43,4 +43,5 @@ pmSourceIO_CMF_PS1_V1.c \ pmSourceIO_CMF_PS1_V2.c \ + pmSourceIO_CMF_PS1_V3.c \ pmSourceIO_CMF_PS1_SV1.c \ pmSourceIO_CMF_PS1_DV1.c \ @@ -71,5 +72,7 @@ models/pmModel_QGAUSS.c \ models/pmModel_RGAUSS.c \ - models/pmModel_SERSIC.c + models/pmModel_SERSIC.c \ + models/pmModel_EXP.c \ + models/pmModel_DEV.c pkginclude_HEADERS = \ @@ -80,4 +83,5 @@ pmPeaks.h \ pmMoments.h \ + pmModelFuncs.h \ pmModel.h \ pmModelClass.h \ @@ -111,5 +115,7 @@ models/pmModel_QGAUSS.h \ models/pmModel_RGAUSS.h \ - models/pmModel_SERSIC.h + models/pmModel_SERSIC.h \ + models/pmModel_EXP.h \ + models/pmModel_DEV.h CLEANFILES = *~ Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_DEV.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_DEV.c (revision 28643) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_DEV.c (revision 28643) @@ -0,0 +1,472 @@ +/****************************************************************************** + * this file defines the DEV source shape model. Note that these model functions are loaded + * by pmModelClass.c using 'include', and thus need no 'include' statements of their own. The + * models use a psVector to represent the set of parameters, with the sequence used to specify + * the meaning of the parameter. The meaning of the parameters may thus vary depending on the + * specifics of the model. All models which are used as a PSF representations share a few + * parameters, for which # define names are listed in pmModel.h: + + f = exp(-z^n) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + * PM_PAR_7 7 - normalized dev parameter + + note that a standard dev model uses exp(-K*(z^(1/2n) - 1). the additional elements (K, + the -1 offset) are absorbed in this model by the normalization, the exponent, and the + radial scale. We fit the elements in this form, then re-normalize them on output. + *****************************************************************************/ + +#include +#include +#include "pmHDU.h" +#include "pmFPA.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 "pmDetections.h" + +#include "pmModel_DEV.h" + +# define PM_MODEL_NPARAM 7 +# define PM_MODEL_FUNC pmModelFunc_DEV +# define PM_MODEL_FLUX pmModelFlux_DEV +# define PM_MODEL_GUESS pmModelGuess_DEV +# define PM_MODEL_LIMITS pmModelLimits_DEV +# define PM_MODEL_RADIUS pmModelRadius_DEV +# define PM_MODEL_FROM_PSF pmModelFromPSF_DEV +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_DEV +# define PM_MODEL_FIT_STATUS pmModelFitStatus_DEV +# define PM_MODEL_SET_LIMITS pmModelSetLimits_DEV + +// f = exp(-z^0.125) +# define ALPHA 0.125 +// # define ALPHA 0.25 + +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords + +// Lax parameter limits +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; + +// Moderate parameter limits +static float *paramsMinModerate = paramsMinLax; +static float *paramsMaxModerate = paramsMaxLax; + +// Strict parameter limits +static float *paramsMinStrict = paramsMinLax; +static float *paramsMaxStrict = paramsMaxLax; + +// Parameter limits to use +static float *paramsMinUse = paramsMinLax; +static float *paramsMaxUse = paramsMaxLax; +static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 }; + +static bool limitsApply = true; // Apply limits? + +psF32 PM_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *pixcoord) +{ + psF32 *PAR = params->data.F32; + + float index = 0.5 / ALPHA; + float bn = 1.9992*index - 0.3271; + float Io = exp(bn); + + psF32 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 px = X / PAR[PM_PAR_SXX]; + psF32 py = Y / PAR[PM_PAR_SYY]; + psF32 z = (PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y); + + assert (z >= 0); + + psF32 f2 = bn*pow(z,ALPHA); + psF32 f1 = Io*exp(-f2); + psF32 z0 = PAR[PM_PAR_I0]*f1; + psF32 f0 = PAR[PM_PAR_SKY] + z0; + + assert (isfinite(f2)); + assert (isfinite(f1)); + assert (isfinite(z0)); + assert (isfinite(f0)); + + if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + + // gradient is infinite for z = 0; saturate at z = 0.01 + psF32 z1 = (z < 0.01) ? z0*bn*ALPHA*pow(0.01,ALPHA - 1.0) : z0*bn*ALPHA*pow(z,ALPHA - 1.0); + + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +2.0*f1; + + assert (isfinite(z1)); + + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; + } + return (f0); +} + +// define the parameter limits +// AR_MAX is the maximum allowed axis ratio +// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2) +# define AR_MAX 20.0 +# define AR_RATIO 0.99 +bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) +{ + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); + + // we need to calculate the limits for SXY specially + float q2 = NAN; + if (nParam == PM_PAR_SXY) { + float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); + float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); + float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); + q1 = (q1 < 0.0) ? 0.0 : q1; + // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg.. Saturate at that + // angle and let f2,f1 fight it out + q2 = 0.5*sqrtf(q1); + } + + switch (mode) { + case PS_MINIMIZE_BETA_LIMIT: { + psAssert(beta, "Require beta to limit beta"); + float limit = betaUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (fabs(beta[nParam]) > fabs(limit)) { + beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit); + psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g", + nParam, beta[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MIN: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMinUse, "Require parameter limits to limit parameters"); + float limit = paramsMinUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (params[nParam] < limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MAX: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMaxUse, "Require parameter limits to limit parameters"); + float limit = paramsMaxUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (params[nParam] > limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + default: + psAbort("invalid choice for limits"); + } + psAbort("should not reach here"); + return false; +} + +// make an initial guess for parameters +// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters +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 + 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; + + return(true); +} + +psF64 PM_MODEL_FLUX (const psVector *params) +{ + float z, norm; + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; + shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; + shape.sxy = PAR[PM_PAR_SXY]; + + // Area is equivalent to 2 pi sigma^2 + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; + + // the area needs to be multiplied by the integral of f(z) + norm = 0.0; + + # define DZ 0.25 + + float f0 = 1.0; + float f1, f2; + for (z = DZ; z < 50; z += DZ) { + f1 = exp(-pow(z,ALPHA)); + z += DZ; + f2 = exp(-pow(z,ALPHA)); + norm += f0 + 4*f1 + f2; + f0 = f2; + } + norm *= DZ / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + if (flux <= 0) + return (1.0); + if (PAR[PM_PAR_I0] <= 0) + return (1.0); + if (flux >= PAR[PM_PAR_I0]) + return (1.0); + + shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; + shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; + shape.sxy = PAR[PM_PAR_SXY]; + + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + + // f = Io exp(-z^n) -> z^n = ln(Io/f) + psF64 zn = log(PAR[PM_PAR_I0] / flux); + psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA); + + psAssert (isfinite(radius), "fix this code: z should not be nan"); + return (radius); +} + +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) +{ + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + // we require these two parameters to exist + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); + + for (int i = 0; i < psf->params->n; i++) { + if (psf->params->data[i] == NULL) { + out[i] = in[i]; + } else { + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } + } + + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + if (!pmPSF_FitToModel (out, 0.1)) { + psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); + return false; + } + + // apply the model limits here: this truncates excessive extrapolation + // XXX do we need to do this still? should we put in asserts to test? + for (int i = 0; i < psf->params->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params->data[i] == NULL) + continue; + + bool status = true; + status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MIN, i, out, NULL); + status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MAX, i, out, NULL); + if (!status) { + psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", + in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + modelPSF->flags |= PM_MODEL_STATUS_LIMITS; + } + } + + return true; +} + +// construct the PSF model from the FLT model and the psf +// XXX is this sufficiently general do be a global function, not a pmModelClass function? +bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io) +{ + psF32 *PAR = model->params->data.F32; + + // we require these two parameters to exist + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); + + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_I0] = Io; + PAR[PM_PAR_XPOS] = Xo; + PAR[PM_PAR_YPOS] = Yo; + + // supply the model-fitted parameters, or copy from the input + for (int i = 0; i < psf->params->n; i++) { + if (i == PM_PAR_SKY) continue; + if (i == PM_PAR_I0) continue; + if (i == PM_PAR_XPOS) continue; + if (i == PM_PAR_YPOS) continue; + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); + } + + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + // XXX user-defined value for limit? + if (!pmPSF_FitToModel (PAR, 0.1)) { + psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); + return false; + } + + // apply the model limits here: this truncates excessive extrapolation + // XXX do we need to do this still? should we put in asserts to test? + for (int i = 0; i < psf->params->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params->data[i] == NULL) + continue; + + bool status = true; + status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL); + status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL); + if (!status) { + psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo); + model->flags |= PM_MODEL_STATUS_LIMITS; + } + } + return(true); +} + +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + status = true; + status &= (PAR[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + return status; +} + + +void PM_MODEL_SET_LIMITS(pmModelLimitsType type) +{ + switch (type) { + case PM_MODEL_LIMITS_NONE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = true; + break; + case PM_MODEL_LIMITS_IGNORE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = false; + break; + case PM_MODEL_LIMITS_LAX: + paramsMinUse = paramsMinLax; + paramsMaxUse = paramsMaxLax; + limitsApply = true; + break; + case PM_MODEL_LIMITS_MODERATE: + paramsMinUse = paramsMinModerate; + paramsMaxUse = paramsMaxModerate; + limitsApply = true; + break; + case PM_MODEL_LIMITS_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_DEV.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_DEV.h (revision 28643) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_DEV.h (revision 28643) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_DEV_H + +#include "pmModel.h" + +psF32 pmModelFunc_DEV(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_DEV(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_DEV(pmModel *model, pmSource *source); +psF64 pmModelFlux_DEV(const psVector *params); +psF64 pmModelRadius_DEV(const psVector *params, psF64 flux); +bool pmModelFromPSF_DEV(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_DEV(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_DEV(pmModel *model); +void pmModelSetLimits_DEV(pmModelLimitsType type); + +#endif Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_EXP.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_EXP.c (revision 28643) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_EXP.c (revision 28643) @@ -0,0 +1,453 @@ +/****************************************************************************** + * this file defines the Exponential (EXP) source shape model. Note that these model functions + * are loaded by pmModelClass.c using 'include', and thus need no 'include' statements of their + * own. The models use a psVector to represent the set of parameters, with the sequence used + * to specify the meaning of the parameter. The meaning of the parameters may thus vary + * depending on the specifics of the model. All models which are used as a PSF representations + * share a few parameters, for which # define names are listed in pmModel.h: + + f = exp(-sqrt(z)) (since z is r^2) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + + *****************************************************************************/ + +#include +#include +#include "pmHDU.h" +#include "pmFPA.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 "pmDetections.h" + +#include "pmModel_EXP.h" + +# define PM_MODEL_NPARAM 7 +# define PM_MODEL_FUNC pmModelFunc_EXP +# define PM_MODEL_FLUX pmModelFlux_EXP +# define PM_MODEL_GUESS pmModelGuess_EXP +# define PM_MODEL_LIMITS pmModelLimits_EXP +# define PM_MODEL_RADIUS pmModelRadius_EXP +# define PM_MODEL_FROM_PSF pmModelFromPSF_EXP +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_EXP +# define PM_MODEL_FIT_STATUS pmModelFitStatus_EXP +# define PM_MODEL_SET_LIMITS pmModelSetLimits_EXP + +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords + +// Lax parameter limits +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; + +// Moderate parameter limits +static float *paramsMinModerate = paramsMinLax; +static float *paramsMaxModerate = paramsMaxLax; + +// Strict parameter limits +static float *paramsMinStrict = paramsMinLax; +static float *paramsMaxStrict = paramsMaxLax; + +// Parameter limits to use +static float *paramsMinUse = paramsMinLax; +static float *paramsMaxUse = paramsMaxLax; +static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5}; + +static bool limitsApply = true; // Apply limits? + +psF32 PM_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *pixcoord) +{ + psF32 *PAR = params->data.F32; + + psF32 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 px = X / PAR[PM_PAR_SXX]; + psF32 py = Y / PAR[PM_PAR_SYY]; + psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; + + // XXX if the elliptical contour is defined in valid way, this step should not be required. + // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values + // for negative values of z + // XXX use an assert here to force the elliptical parameters to be correctly determined + // if (z < 0) z = 0; + assert (z >= 0); + + psF32 f2 = sqrt(z); + psF32 f1 = exp(-f2); + psF32 z0 = PAR[PM_PAR_I0]*f1; + psF32 f0 = PAR[PM_PAR_SKY] + z0; + + assert (isfinite(f2)); + assert (isfinite(f1)); + assert (isfinite(z0)); + assert (isfinite(f0)); + + if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +f1; + + // gradient is infinite for z = 0; saturate at z = 0.01 + // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i] + psF32 z1 = (z < 0.01) ? 0.5*z0/sqrt(0.01) : 0.5*z0/sqrt(z); + + // XXX dampen SXX and SYY as in GAUSS? + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; + } + return (f0); +} + +// define the parameter limits +// AR_MAX is the maximum allowed axis ratio +// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2) +# define AR_MAX 20.0 +# define AR_RATIO 0.99 +bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) +{ + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); + + // we need to calculate the limits for SXY specially + float q2 = NAN; + if (nParam == PM_PAR_SXY) { + float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); + float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); + float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); + q1 = (q1 < 0.0) ? 0.0 : q1; + // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg.. Saturate at that + // angle and let f2,f1 fight it out + q2 = 0.5*sqrtf(q1); + } + + switch (mode) { + case PS_MINIMIZE_BETA_LIMIT: { + psAssert(beta, "Require beta to limit beta"); + float limit = betaUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (fabs(beta[nParam]) > fabs(limit)) { + beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit); + psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g", + nParam, beta[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MIN: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMinUse, "Require parameter limits to limit parameters"); + float limit = paramsMinUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (params[nParam] < limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MAX: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMaxUse, "Require parameter limits to limit parameters"); + float limit = paramsMaxUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (params[nParam] > limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + default: + psAbort("invalid choice for limits"); + } + psAbort("should not reach here"); + return false; +} + +// make an initial guess for parameters +// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters +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; + + 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; + return(true); +} + +psF64 PM_MODEL_FLUX (const psVector *params) +{ + float z, norm; + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; + shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; + shape.sxy = PAR[PM_PAR_SXY]; + + // Area is equivalent to 2 pi sigma^2 + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; + + // the area needs to be multiplied by the integral of f(z) = exp(-sqrt(z)) [0 to infinity] + norm = 0.0; + + # define DZ 0.25 + + float f0 = 1.0; + float f1, f2; + for (z = DZ; z < 50; z += DZ) { + f1 = exp(-sqrt(z)); + z += DZ; + f2 = exp(-sqrt(z)); + norm += f0 + 4*f1 + f2; + f0 = f2; + } + norm *= DZ / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psEllipseShape shape; + + psF32 *PAR = params->data.F32; + + if (flux <= 0) + return (1.0); + if (PAR[PM_PAR_I0] <= 0) + return (1.0); + if (flux >= PAR[PM_PAR_I0]) + return (1.0); + + shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; + shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; + shape.sxy = PAR[PM_PAR_SXY]; + + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + + // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f) + psF64 zn = log(PAR[PM_PAR_I0] / flux); + 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]); + return (radius); +} + +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) +{ + + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + // we require these two parameters to exist + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); + + for (int i = 0; i < psf->params->n; i++) { + if (psf->params->data[i] == NULL) { + out[i] = in[i]; + } else { + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } + } + + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + if (!pmPSF_FitToModel (out, 0.1)) { + psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); + return false; + } + + // apply the model limits here: this truncates excessive extrapolation + // XXX do we need to do this still? should we put in asserts to test? + for (int i = 0; i < psf->params->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params->data[i] == NULL) + continue; + + bool status = true; + status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MIN, i, out, NULL); + status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MAX, i, out, NULL); + if (!status) { + psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", + in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + modelPSF->flags |= PM_MODEL_STATUS_LIMITS; + } + } + + return true; +} + +// construct the PSF model from the FLT model and the psf +// XXX is this sufficiently general do be a global function, not a pmModelClass function? +bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io) +{ + psF32 *PAR = model->params->data.F32; + + // we require these two parameters to exist + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); + + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_I0] = Io; + PAR[PM_PAR_XPOS] = Xo; + PAR[PM_PAR_YPOS] = Yo; + + // supply the model-fitted parameters, or copy from the input + for (int i = 0; i < psf->params->n; i++) { + if (i == PM_PAR_SKY) continue; + if (i == PM_PAR_I0) continue; + if (i == PM_PAR_XPOS) continue; + if (i == PM_PAR_YPOS) continue; + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); + } + + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + // XXX user-defined value for limit? + if (!pmPSF_FitToModel (PAR, 0.1)) { + psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); + return false; + } + + // apply the model limits here: this truncates excessive extrapolation + // XXX do we need to do this still? should we put in asserts to test? + for (int i = 0; i < psf->params->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params->data[i] == NULL) + continue; + + bool status = true; + status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL); + status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL); + if (!status) { + psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo); + model->flags |= PM_MODEL_STATUS_LIMITS; + } + } + return(true); +} + +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + status = true; + status &= (PAR[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + return status; +} + + +void PM_MODEL_SET_LIMITS(pmModelLimitsType type) +{ + switch (type) { + case PM_MODEL_LIMITS_NONE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = true; + break; + case PM_MODEL_LIMITS_IGNORE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = false; + break; + case PM_MODEL_LIMITS_LAX: + paramsMinUse = paramsMinLax; + paramsMaxUse = paramsMaxLax; + limitsApply = true; + break; + case PM_MODEL_LIMITS_MODERATE: + paramsMinUse = paramsMinModerate; + paramsMaxUse = paramsMaxModerate; + limitsApply = true; + break; + case PM_MODEL_LIMITS_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_EXP.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_EXP.h (revision 28643) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_EXP.h (revision 28643) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_EXP_H + +#include "pmModel.h" + +psF32 pmModelFunc_EXP(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_EXP(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_EXP(pmModel *model, pmSource *source); +psF64 pmModelFlux_EXP(const psVector *params); +psF64 pmModelRadius_EXP(const psVector *params, psF64 flux); +bool pmModelFromPSF_EXP(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_EXP(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_EXP(pmModel *model); +void pmModelSetLimits_EXP(pmModelLimitsType type); + +#endif Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_GAUSS.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_GAUSS.c (revision 28643) @@ -21,11 +21,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.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 "pmPeaks.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 "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_GAUSS.h" +# define PM_MODEL_NPARAM 7 # define PM_MODEL_FUNC pmModelFunc_GAUSS # define PM_MODEL_FLUX pmModelFlux_GAUSS @@ -83,4 +103,5 @@ dPAR[PM_PAR_XPOS] = q*(2*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); dPAR[PM_PAR_YPOS] = q*(2*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + // the extra factor of 2 below is needed to avoid excessive swings dPAR[PM_PAR_SXX] = +4.0*q*px*px/PAR[PM_PAR_SXX]; @@ -102,5 +123,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -347,4 +368,5 @@ // this test is invalid if the parameters are derived // from the PSF model +// XXX how is this used? it prevents forced photometry from ever being 'successful' bool PM_MODEL_FIT_STATUS (pmModel *model) { @@ -394,12 +416,2 @@ return; } - -# undef PM_MODEL_FUNC -# undef PM_MODEL_FLUX -# undef PM_MODEL_GUESS -# undef PM_MODEL_LIMITS -# undef PM_MODEL_RADIUS -# undef PM_MODEL_FROM_PSF -# undef PM_MODEL_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_PGAUSS.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_PGAUSS.c (revision 28643) @@ -21,11 +21,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.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 "pmPeaks.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 "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_PGAUSS.h" +# define PM_MODEL_NPARAM 7 # define PM_MODEL_FUNC pmModelFunc_PGAUSS # define PM_MODEL_FLUX pmModelFlux_PGAUSS @@ -103,5 +123,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -448,12 +468,2 @@ return; } - -# undef PM_MODEL_FUNC -# undef PM_MODEL_FLUX -# undef PM_MODEL_GUESS -# undef PM_MODEL_LIMITS -# undef PM_MODEL_RADIUS -# undef PM_MODEL_FROM_PSF -# undef PM_MODEL_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_PS1_V1.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_PS1_V1.c (revision 28643) @@ -22,11 +22,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.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 "pmPeaks.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 "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_PS1_V1.h" +# define PM_MODEL_NPARAM 8 # define PM_MODEL_FUNC pmModelFunc_PS1_V1 # define PM_MODEL_FLUX pmModelFlux_PS1_V1 @@ -122,5 +142,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -222,5 +242,5 @@ PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); PAR[PM_PAR_SXY] = shape.sxy; - PAR[PM_PAR_7] = 1.0; + PAR[PM_PAR_7] = 0.5; return(true); @@ -468,14 +488,2 @@ return; } - -# undef PM_MODEL_FUNC -# undef PM_MODEL_FLUX -# undef PM_MODEL_GUESS -# undef PM_MODEL_LIMITS -# undef PM_MODEL_RADIUS -# undef PM_MODEL_FROM_PSF -# undef PM_MODEL_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS -# undef ALPHA -# undef ALPHA_M Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_QGAUSS.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_QGAUSS.c (revision 28643) @@ -22,11 +22,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.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 "pmPeaks.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 "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_QGAUSS.h" +# define PM_MODEL_NPARAM 8 # define PM_MODEL_FUNC pmModelFunc_QGAUSS # define PM_MODEL_FLUX pmModelFlux_QGAUSS @@ -123,5 +143,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -469,14 +489,2 @@ return; } - -# undef PM_MODEL_FUNC -# undef PM_MODEL_FLUX -# undef PM_MODEL_GUESS -# undef PM_MODEL_LIMITS -# undef PM_MODEL_RADIUS -# undef PM_MODEL_FROM_PSF -# undef PM_MODEL_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS -# undef ALPHA -# undef ALPHA_M Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_RGAUSS.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_RGAUSS.c (revision 28643) @@ -22,11 +22,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.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 "pmPeaks.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 "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_RGAUSS.h" +# define PM_MODEL_NPARAM 8 # define PM_MODEL_FUNC pmModelFunc_RGAUSS # define PM_MODEL_FLUX pmModelFlux_RGAUSS @@ -112,5 +132,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -209,8 +229,8 @@ 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_SXX] = PS_MAX(0.5, shape.sx); + PAR[PM_PAR_SYY] = PS_MAX(0.5, shape.sy); PAR[PM_PAR_SXY] = shape.sxy; - PAR[PM_PAR_7] = 2.25; + PAR[PM_PAR_7] = 1.5; return(true); @@ -463,12 +483,2 @@ return; } - -# undef PM_MODEL_FUNC -# undef PM_MODEL_FLUX -# undef PM_MODEL_GUESS -# undef PM_MODEL_LIMITS -# undef PM_MODEL_RADIUS -# undef PM_MODEL_FROM_PSF -# undef PM_MODEL_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_SERSIC.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_SERSIC.c (revision 28643) @@ -18,5 +18,5 @@ * PM_PAR_7 7 - normalized sersic parameter - note that a standard sersic model uses exp(-K*(z^(1/n) - 1). the additional elements (K, + note that a standard sersic model uses exp(-K*(z^(1/2n) - 1). the additional elements (K, the -1 offset) are absorbed in this model by the normalization, the exponent, and the radial scale. We fit the elements in this form, then re-normalize them on output. @@ -25,11 +25,31 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.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 "pmPeaks.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 "pmModel.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmModel_SERSIC.h" +# define PM_MODEL_NPARAM 8 # define PM_MODEL_FUNC pmModelFunc_SERSIC # define PM_MODEL_FLUX pmModelFlux_SERSIC @@ -47,5 +67,5 @@ // Lax parameter limits -static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0, 0.05 }; +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 }; static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 }; @@ -84,6 +104,10 @@ assert (z >= 0); - psF32 f2 = pow(z,PAR[PM_PAR_7]); - psF32 f1 = exp(-f2); + float index = 0.5 / PAR[PM_PAR_7]; + float bn = 1.9992*index - 0.3271; + float Io = exp(bn); + + psF32 f2 = bn*pow(z,PAR[PM_PAR_7]); + psF32 f1 = Io*exp(-f2); psF32 z0 = PAR[PM_PAR_I0]*f1; psF32 f0 = PAR[PM_PAR_SKY] + z0; @@ -98,9 +122,9 @@ // gradient is infinite for z = 0; saturate at z = 0.01 - psF32 z1 = (z < 0.01) ? z0*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0); + psF32 z1 = (z < 0.01) ? z0*bn*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*bn*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0); dPAR[PM_PAR_SKY] = +1.0; dPAR[PM_PAR_I0] = +f1; - dPAR[PM_PAR_7] = (z == 0.0) ? 0.0 : -z0*f2*log(z); + dPAR[PM_PAR_7] = (z < 0.01) ? -1.5*z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -1.5*z0*f2*log(z); assert (isfinite(z1)); @@ -111,5 +135,4 @@ dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; - dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; } @@ -127,5 +150,5 @@ return true; } - psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); // we need to calculate the limits for SXY specially @@ -219,12 +242,27 @@ if (!isfinite(shape.sxy)) return false; + // the other parameters depend on the guess for PAR_7 + if (0) { + } else { + PAR[PM_PAR_7] = 0.25; + // PAR[PM_PAR_7] = 0.125; + // PAR[PM_PAR_7] = 0.5; + } + + float index = 0.5 / PAR[PM_PAR_7]; + 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; + PAR[PM_PAR_I0] = peak->flux / Io; 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_SXX] = Sxx; + PAR[PM_PAR_SYY] = Syy; PAR[PM_PAR_SXY] = shape.sxy; - PAR[PM_PAR_7] = 0.5; return(true); @@ -254,7 +292,8 @@ float f1, f2; for (z = DZ; z < 50; z += DZ) { - f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); + // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); + f1 = exp(-pow(z,PAR[PM_PAR_7])); z += DZ; - f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); + f2 = exp(-pow(z,PAR[PM_PAR_7])); norm += f0 + 4*f1 + f2; f0 = f2; @@ -287,17 +326,10 @@ psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); - psF64 sigma = axes.major; - - psF64 limit = flux / PAR[PM_PAR_I0]; - - psF64 z = pow (-log(limit), (1.0 / PAR[PM_PAR_7])); - psAssert (isfinite(z), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]); - - psF64 radius = sigma * sqrt (2.0 * z); - psAssert (isfinite(radius), "fix this code: radius should not be nan for %f, %f", PAR[PM_PAR_7], sigma); - - if (isnan(radius)) - psAbort("error in code: radius is NaN"); - + + // f = Io exp(-z^n) -> z^n = ln(Io/f) + psF64 zn = log(PAR[PM_PAR_I0] / flux); + psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]); + + psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]); return (radius); } @@ -448,12 +480,2 @@ return; } - -# undef PM_MODEL_FUNC -# undef PM_MODEL_FLUX -# undef PM_MODEL_GUESS -# undef PM_MODEL_LIMITS -# undef PM_MODEL_RADIUS -# undef PM_MODEL_FROM_PSF -# undef PM_MODEL_PARAMS_FROM_PSF -# undef PM_MODEL_FIT_STATUS -# undef PM_MODEL_SET_LIMITS Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmDetEff.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmDetEff.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmDetEff.c (revision 28643) @@ -3,8 +3,30 @@ #endif +#include #include +#include "pmHDU.h" +#include "pmFPA.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 "pmDetections.h" #include "pmDetEff.h" - static void detEffFree(pmDetEff *de) Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmDetEff.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmDetEff.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmDetEff.h (revision 28643) @@ -1,9 +1,4 @@ #ifndef PM_DET_EFF_H #define PM_DET_EFF_H - -#include -#include - -#include "pmFPA.h" #define PM_DETEFF_ANALYSIS "DETEFF" // Location of detection efficiency on pmReadout.analysis Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprint.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprint.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprint.c (revision 28643) @@ -19,4 +19,5 @@ #include #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprint.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprint.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprint.h (revision 28643) @@ -10,8 +10,4 @@ #ifndef PM_FOOTPRINT_H #define PM_FOOTPRINT_H - -#include -#include "pmSpan.h" -#include "pmFootprintSpans.h" typedef struct { Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintArrayGrow.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintArrayGrow.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintArrayGrow.c (revision 28643) @@ -18,4 +18,5 @@ #include #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintArraysMerge.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintArraysMerge.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintArraysMerge.c (revision 28643) @@ -18,4 +18,5 @@ #include #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintAssignPeaks.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintAssignPeaks.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintAssignPeaks.c (revision 28643) @@ -18,4 +18,5 @@ #include #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintCullPeaks.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintCullPeaks.c (revision 28643) @@ -19,6 +19,6 @@ #include #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" -#include "pmFootprintSpans.h" #include "pmPeaks.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintFind.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintFind.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintFind.c (revision 28643) @@ -19,4 +19,5 @@ #include #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintFindAtPoint.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintFindAtPoint.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintFindAtPoint.c (revision 28643) @@ -19,7 +19,7 @@ #include #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" -#include "pmFootprintSpans.h" /************************************************************************************************************/ Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintIDs.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintIDs.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintIDs.c (revision 28643) @@ -19,4 +19,5 @@ #include #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintSpans.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintSpans.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintSpans.c (revision 28643) @@ -19,7 +19,7 @@ #include #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" -#include "pmFootprintSpans.h" static void pmStartSpanFree(pmStartSpan *sspan) { Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintSpans.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintSpans.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmFootprintSpans.h (revision 28643) @@ -10,7 +10,4 @@ #ifndef PM_FOOTPRINT_SPANS_H #define PM_FOOTPRINT_SPANS_H - -#include -#include "pmSpan.h" /* We define two helper structures used in building the pmFootprints: Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmGrowthCurve.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmGrowthCurve.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmGrowthCurve.c (revision 28643) @@ -20,14 +20,26 @@ #include "pmHDU.h" #include "pmFPA.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "psVectorBracket.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmGrowthCurveGenerate.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmGrowthCurveGenerate.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmGrowthCurveGenerate.c (revision 28643) @@ -24,20 +24,30 @@ #include "pmHDU.h" #include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "psVectorBracket.h" +#include "pmErrorCodes.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" -#include "pmModelUtils.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmSourcePhotometry.h" -#include "pmFPAMaskWeight.h" -#include "psVectorBracket.h" -#include "pmErrorCodes.h" pmGrowthCurve *pmGrowthCurveForPosition (psImage *image, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType markVal, float xc, float yc); Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModel.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModel.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModel.c (revision 28643) @@ -23,8 +23,14 @@ #include "pmHDU.h" #include "pmFPA.h" + +#include "pmTrend2D.h" #include "pmResiduals.h" #include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" #include "pmModel.h" #include "pmModelClass.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModel.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModel.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModel.h (revision 28643) @@ -14,7 +14,4 @@ # define PM_MODEL_H -#include -#include "pmPSF.h" - /// @addtogroup Objects Object Detection / Analysis Functions /// @{ @@ -22,68 +19,4 @@ /* pointers for the functions types below are supplied to each pmModel, and can be used by the programmer without needing to know the model class */ - -typedef enum { - PM_MODEL_STATUS_NONE = 0x00, ///< model fit not yet attempted, no other info - PM_MODEL_STATUS_FITTED = 0x01, ///< model fit completed - PM_MODEL_STATUS_NONCONVERGE = 0x02, ///< model fit did not converge - PM_MODEL_STATUS_OFFIMAGE = 0x04, ///< model fit drove out of range - PM_MODEL_STATUS_BADARGS = 0x08, ///< model fit called with invalid args - PM_MODEL_STATUS_LIMITS = 0x10 ///< model parameters hit limits -} pmModelStatus; - -typedef enum { - PM_MODEL_OP_NONE = 0x00, - PM_MODEL_OP_FUNC = 0x01, - PM_MODEL_OP_RES0 = 0x02, - PM_MODEL_OP_RES1 = 0x04, - PM_MODEL_OP_FULL = 0x07, - PM_MODEL_OP_SKY = 0x08, - PM_MODEL_OP_CENTER = 0x10, - PM_MODEL_OP_NORM = 0x20, - PM_MODEL_OP_NOISE = 0x40, -} pmModelOpMode; - -/// Parameter limit types -typedef enum { - PM_MODEL_LIMITS_NONE, ///< Apply no limits: suitable for debugging - PM_MODEL_LIMITS_IGNORE, ///< Ignore all limits: fit can go to town - PM_MODEL_LIMITS_LAX, ///< Lax limits: attempting to reproduce even bad data - PM_MODEL_LIMITS_MODERATE, ///< Moderate limits: cope with mildly bad data - PM_MODEL_LIMITS_STRICT, ///< Strict limits: insist on good quality data -} pmModelLimitsType; - -typedef struct pmModel pmModel; -typedef struct pmSource pmSource; - -// This function is the model chi-square minimization function for this model. -typedef psMinimizeLMChi2Func pmModelFunc; - -// This function sets the parameter limits for this model. -typedef psMinimizeLMLimitFunc pmModelLimits; - -// This function returns the integrated flux for the given model parameters. -typedef psF64 (*pmModelFlux)(const psVector *params); - -// This function returns the radius at which the given model and parameters -// achieves the given flux. -typedef psF64 (*pmModelRadius)(const psVector *params, double flux); - -// This function provides the model guess parameters based on the details of -// the given source. -typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); - -// This function constructs the PSF model for the given source based on the -// supplied psf and the EXT model for the object. -typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, const pmPSF *psf); - -// This function sets the model parameters based on the PSF for a given coordinate and central -// intensity -typedef bool (*pmModelParamsFromPSF)(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); - -// This function returns the success / failure status of the given model fit -typedef bool (*pmModelFitStatusFunc)(pmModel *model); - -// This function sets the parameter limits for the given model -typedef bool (*pmModelSetLimitsFunc)(pmModelLimitsType type); /** pmModel data structure @@ -122,18 +55,4 @@ pmModelSetLimitsFunc modelSetLimits; }; - -/** Symbolic names for the elements of [d]params - * Note: these are #defines not enums as a given element of [d]params - * may/will correspond to different parameters in different contexts - */ -#define PM_PAR_SKY 0 ///< Sky -#define PM_PAR_I0 1 ///< Central intensity -#define PM_PAR_XPOS 2 ///< X centre of object -#define PM_PAR_YPOS 3 ///< Y centre of object -#define PM_PAR_SXX 4 ///< shape X^2 moment -#define PM_PAR_SYY 5 ///< shape Y^2 moment -#define PM_PAR_SXY 6 ///< shape XY moment -#define PM_PAR_7 7 ///< ??? Unknown parameter -#define PM_PAR_8 8 ///< ??? Unknown parameter /** pmModelAlloc() Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelClass.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelClass.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelClass.c (revision 28643) @@ -23,15 +23,18 @@ #include "pmHDU.h" #include "pmFPA.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.h" +#include "pmModelFuncs.h" #include "pmModel.h" -#include "pmSource.h" +#include "pmModelUtils.h" #include "pmModelClass.h" + #include "pmErrorCodes.h" @@ -46,4 +49,6 @@ # include "models/pmModel_RGAUSS.h" # include "models/pmModel_SERSIC.h" +# include "models/pmModel_EXP.h" +# include "models/pmModel_DEV.h" static pmModelClass defaultModels[] = { @@ -53,5 +58,7 @@ {"PS_MODEL_PS1_V1", 8, (pmModelFunc)pmModelFunc_PS1_V1, (pmModelFlux)pmModelFlux_PS1_V1, (pmModelRadius)pmModelRadius_PS1_V1, (pmModelLimits)pmModelLimits_PS1_V1, (pmModelGuessFunc)pmModelGuess_PS1_V1, (pmModelFromPSFFunc)pmModelFromPSF_PS1_V1, (pmModelParamsFromPSF)pmModelParamsFromPSF_PS1_V1, (pmModelFitStatusFunc)pmModelFitStatus_PS1_V1, (pmModelSetLimitsFunc)pmModelSetLimits_PS1_V1 }, {"PS_MODEL_RGAUSS", 8, (pmModelFunc)pmModelFunc_RGAUSS, (pmModelFlux)pmModelFlux_RGAUSS, (pmModelRadius)pmModelRadius_RGAUSS, (pmModelLimits)pmModelLimits_RGAUSS, (pmModelGuessFunc)pmModelGuess_RGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_RGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_RGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_RGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_RGAUSS }, - {"PS_MODEL_SERSIC", 8, (pmModelFunc)pmModelFunc_SERSIC, (pmModelFlux)pmModelFlux_SERSIC, (pmModelRadius)pmModelRadius_SERSIC, (pmModelLimits)pmModelLimits_SERSIC, (pmModelGuessFunc)pmModelGuess_SERSIC, (pmModelFromPSFFunc)pmModelFromPSF_SERSIC, (pmModelParamsFromPSF)pmModelParamsFromPSF_SERSIC, (pmModelFitStatusFunc)pmModelFitStatus_SERSIC, (pmModelSetLimitsFunc)pmModelSetLimits_SERSIC } + {"PS_MODEL_SERSIC", 8, (pmModelFunc)pmModelFunc_SERSIC, (pmModelFlux)pmModelFlux_SERSIC, (pmModelRadius)pmModelRadius_SERSIC, (pmModelLimits)pmModelLimits_SERSIC, (pmModelGuessFunc)pmModelGuess_SERSIC, (pmModelFromPSFFunc)pmModelFromPSF_SERSIC, (pmModelParamsFromPSF)pmModelParamsFromPSF_SERSIC, (pmModelFitStatusFunc)pmModelFitStatus_SERSIC, (pmModelSetLimitsFunc)pmModelSetLimits_SERSIC }, + {"PS_MODEL_EXP", 7, (pmModelFunc)pmModelFunc_EXP, (pmModelFlux)pmModelFlux_EXP, (pmModelRadius)pmModelRadius_EXP, (pmModelLimits)pmModelLimits_EXP, (pmModelGuessFunc)pmModelGuess_EXP, (pmModelFromPSFFunc)pmModelFromPSF_EXP, (pmModelParamsFromPSF)pmModelParamsFromPSF_EXP, (pmModelFitStatusFunc)pmModelFitStatus_EXP, (pmModelSetLimitsFunc)pmModelSetLimits_EXP }, + {"PS_MODEL_DEV", 7, (pmModelFunc)pmModelFunc_DEV, (pmModelFlux)pmModelFlux_DEV, (pmModelRadius)pmModelRadius_DEV, (pmModelLimits)pmModelLimits_DEV, (pmModelGuessFunc)pmModelGuess_DEV, (pmModelFromPSFFunc)pmModelFromPSF_DEV, (pmModelParamsFromPSF)pmModelParamsFromPSF_DEV, (pmModelFitStatusFunc)pmModelFitStatus_DEV, (pmModelSetLimitsFunc)pmModelSetLimits_DEV }, }; Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelClass.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelClass.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelClass.h (revision 28643) @@ -28,6 +28,4 @@ # ifndef PM_MODEL_CLASS_H # define PM_MODEL_CLASS_H - -#include /// @addtogroup Objects Object Detection / Analysis Functions Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelFuncs.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelFuncs.h (revision 28643) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelFuncs.h (revision 28643) @@ -0,0 +1,119 @@ +/* @file pmModelClass.h + * + * The object model function types are desined to allow for the flexible addition of new + * object models. Every object model, with parameters represented by pmModel, has an + * associated set of functions which provide necessary support operations. A These + * functions allow the programmer to select the approriate function or property for a + * specific object model class. + * + * Every model instance belongs to a class of models, defined by the value of the + * pmModelType type entry. Various functions need access to information about each of the + * models. Some of this information varies from model to model, and may depend on the + * current parameter values or other data quantities. In order to keep the code from + * requiring the information about each model to be coded into the low-level fitting + * routines, we define a collection of functions which allow us to abstract this type of + * model-dependent information. These generic functions take the model type and return the + * corresponding function pointer for the specified model. Each model is defined by + * creating this collection of specific functions, and placing them in a single file for + * each model. We define the following structure to carry the collection of information + * about the models. + * + * @author EAM, IfA + * + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-27 03:14:57 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +# ifndef PM_MODEL_FUNCS_H +# define PM_MODEL_FUNCS_H + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +// type of model carried by the pmModel structure +typedef int pmModelType; + +typedef enum { + PM_MODEL_STATUS_NONE = 0x00, ///< model fit not yet attempted, no other info + PM_MODEL_STATUS_FITTED = 0x01, ///< model fit completed + PM_MODEL_STATUS_NONCONVERGE = 0x02, ///< model fit did not converge + PM_MODEL_STATUS_OFFIMAGE = 0x04, ///< model fit drove out of range + PM_MODEL_STATUS_BADARGS = 0x08, ///< model fit called with invalid args + PM_MODEL_STATUS_LIMITS = 0x10 ///< model parameters hit limits +} pmModelStatus; + +typedef enum { + PM_MODEL_OP_NONE = 0x00, + PM_MODEL_OP_FUNC = 0x01, + PM_MODEL_OP_RES0 = 0x02, + PM_MODEL_OP_RES1 = 0x04, + PM_MODEL_OP_FULL = 0x07, + PM_MODEL_OP_SKY = 0x08, + PM_MODEL_OP_CENTER = 0x10, + PM_MODEL_OP_NORM = 0x20, + PM_MODEL_OP_NOISE = 0x40, +} pmModelOpMode; + +/// Parameter limit types +typedef enum { + PM_MODEL_LIMITS_NONE, ///< Apply no limits: suitable for debugging + PM_MODEL_LIMITS_IGNORE, ///< Ignore all limits: fit can go to town + PM_MODEL_LIMITS_LAX, ///< Lax limits: attempting to reproduce even bad data + PM_MODEL_LIMITS_MODERATE, ///< Moderate limits: cope with mildly bad data + PM_MODEL_LIMITS_STRICT, ///< Strict limits: insist on good quality data +} pmModelLimitsType; + +/** Symbolic names for the elements of [d]params + * Note: these are #defines not enums as a given element of [d]params + * may/will correspond to different parameters in different contexts + */ +#define PM_PAR_SKY 0 ///< Sky +#define PM_PAR_I0 1 ///< Central intensity +#define PM_PAR_XPOS 2 ///< X center of object +#define PM_PAR_YPOS 3 ///< Y center of object +#define PM_PAR_SXX 4 ///< shape X^2 moment +#define PM_PAR_SYY 5 ///< shape Y^2 moment +#define PM_PAR_SXY 6 ///< shape XY moment +#define PM_PAR_7 7 ///< Model-dependent parameter +#define PM_PAR_8 8 ///< Model-dependent parameter + +/*** these prototype classes are used to define elements of the pmModelClass structure below ***/ + +typedef struct pmModel pmModel; +typedef struct pmSource pmSource; +typedef struct pmPSF pmPSF; + +// This function is the model chi-square minimization function for this model. +typedef psMinimizeLMChi2Func pmModelFunc; + +// This function sets the parameter limits for this model. +typedef psMinimizeLMLimitFunc pmModelLimits; + +// This function returns the integrated flux for the given model parameters. +typedef psF64 (*pmModelFlux)(const psVector *params); + +// This function returns the radius at which the given model and parameters +// achieves the given flux. +typedef psF64 (*pmModelRadius)(const psVector *params, double flux); + +// This function provides the model guess parameters based on the details of +// the given source. +typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); + +// This function constructs the PSF model for the given source based on the +// supplied psf and the EXT model for the object. +typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, const pmPSF *psf); + +// This function sets the model parameters based on the PSF for a given coordinate and central +// intensity +typedef bool (*pmModelParamsFromPSF)(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); + +// This function returns the success / failure status of the given model fit +typedef bool (*pmModelFitStatusFunc)(pmModel *model); + +// This function sets the parameter limits for the given model +typedef bool (*pmModelSetLimitsFunc)(pmModelLimitsType type); + +/// @} +# endif /* PM_MODEL_FUNCS_H */ Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelUtils.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelUtils.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmModelUtils.c (revision 28643) @@ -22,11 +22,25 @@ #include "pmHDU.h" #include "pmFPA.h" + +#include "pmTrend2D.h" #include "pmResiduals.h" #include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" + #include "pmErrorCodes.h" -#include "pmModelUtils.h" /***************************************************************************** Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmMoments.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmMoments.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmMoments.c (revision 28643) @@ -29,4 +29,13 @@ psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__); pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments)); + + tmp->Mrf = 0.0; + tmp->Mrh = 0.0; + tmp->KronFlux = 0.0; + tmp->KronFluxErr = 0.0; + + tmp->KronFinner = 0.0; + tmp->KronFouter = 0.0; + tmp->Mx = 0.0; tmp->My = 0.0; Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmMoments.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmMoments.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmMoments.h (revision 28643) @@ -24,4 +24,7 @@ typedef struct { + float Mrf; ///< radial first moment + float Mrh; ///< radial half moment + float Mx; ///< X-coord of centroid. float My; ///< Y-coord of centroid. @@ -47,4 +50,11 @@ float SN; ///< approx signal-to-noise int nPixels; ///< Number of pixels used. + + float KronFlux; ///< Kron flux (flux in 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) + } pmMoments; Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF.c (revision 28643) @@ -25,19 +25,28 @@ #include "pmHDU.h" #include "pmFPA.h" +#include "pmFPAMaskWeight.h" +#include "psVectorBracket.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" -#include "pmModelUtils.h" -#include "pmSourcePhotometry.h" -#include "pmFPAMaskWeight.h" -#include "psVectorBracket.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmErrorCodes.h" @@ -54,4 +63,5 @@ psFree (options->stats); + psFree (options->fitOptions); return; } @@ -65,4 +75,5 @@ options->stats = NULL; + options->fitOptions = NULL; // XXX this has to be set before calling pmPSF fit functions options->psfTrendMode = PM_TREND_NONE; Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF.h (revision 28643) @@ -14,15 +14,6 @@ # define PM_PSF_H -#include -#include "pmTrend2D.h" -#include "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmFPA.h" - /// @addtogroup Objects Object Detection / Analysis Functions /// @{ - -// type of model carried by the pmModel structure -typedef int pmModelType; /** pmPSF data structure @@ -38,5 +29,5 @@ * */ -typedef struct { +struct pmPSF { pmModelType type; ///< PSF Model in use psArray *params; ///< Model parameters (psPolynomial2D) @@ -64,5 +55,5 @@ pmGrowthCurve *growth; ///< apMag vs Radius pmResiduals *residuals; ///< normalized residual image (no spatial variation) -} pmPSF; +}; typedef struct { @@ -82,4 +73,5 @@ float apRadius; bool chiFluxTrend; // Fit a trend in Chi2 as a function of flux? + pmSourceFitOptions *fitOptions; } pmPSFOptions; Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSF_IO.c (revision 28643) @@ -38,17 +38,26 @@ #include "pmFPAfileFitsIO.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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + #include "pmPSF_IO.h" -#include "pmSource.h" -#include "pmModelClass.h" -#include "pmModelUtils.h" #include "pmSourceIO.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtry.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtry.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtry.c (revision 28643) @@ -19,19 +19,26 @@ #include "pmFPA.h" #include "pmFPAMaskWeight.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmSourceUtils.h" #include "pmPSFtry.h" -#include "pmModelClass.h" -#include "pmModelUtils.h" -#include "pmSourceFitModel.h" +#include "pmDetections.h" + #include "pmSourcePhotometry.h" #include "pmSourceVisual.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryFitEXT.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryFitEXT.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryFitEXT.c (revision 28643) @@ -19,19 +19,27 @@ #include "pmFPA.h" #include "pmFPAMaskWeight.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.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 "pmSourceUtils.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" #include "pmPSFtry.h" -#include "pmModelClass.h" -#include "pmModelUtils.h" -#include "pmSourceFitModel.h" +#include "pmDetections.h" + #include "pmSourcePhotometry.h" #include "pmSourceVisual.h" @@ -44,4 +52,7 @@ psTimerStart ("psf.fit"); + + // in this segment, we are fitting the full PSF model class (shape unconstrained) + options->fitOptions->mode = PM_SOURCE_FIT_EXT; // maskVal is used to test for rejected pixels, and must include markVal @@ -73,5 +84,5 @@ // fit model as EXT, not PSF - status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal); + status = pmSourceFitModel (source, source->modelEXT, options->fitOptions, maskVal); // clear object mask to define valid pixels Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryFitPSF.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryFitPSF.c (revision 28643) @@ -17,19 +17,26 @@ #include "pmFPA.h" #include "pmFPAMaskWeight.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmSourceUtils.h" #include "pmPSFtry.h" -#include "pmModelClass.h" -#include "pmModelUtils.h" -#include "pmSourceFitModel.h" +#include "pmDetections.h" + #include "pmSourcePhotometry.h" #include "pmSourceVisual.h" @@ -42,4 +49,7 @@ psTimerStart ("psf.fit"); + + // in this segment, we are fitting the fitted PSF model class (shape constrained) + options->fitOptions->mode = PM_SOURCE_FIT_PSF; // maskVal is used to test for rejected pixels, and must include markVal @@ -81,5 +91,5 @@ // fit the PSF model to the source - status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_NORM, maskVal); + status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal); // skip poor fits Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryMakePSF.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryMakePSF.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryMakePSF.c (revision 28643) @@ -18,19 +18,26 @@ #include "pmFPA.h" #include "pmFPAMaskWeight.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmSourceUtils.h" #include "pmPSFtry.h" -#include "pmModelClass.h" -#include "pmModelUtils.h" -#include "pmSourceFitModel.h" +#include "pmDetections.h" + #include "pmSourcePhotometry.h" #include "pmSourceVisual.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryMetric.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryMetric.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryMetric.c (revision 28643) @@ -18,19 +18,26 @@ #include "pmFPA.h" #include "pmFPAMaskWeight.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmSourceUtils.h" #include "pmPSFtry.h" -#include "pmModelClass.h" -#include "pmModelUtils.h" -#include "pmSourceFitModel.h" +#include "pmDetections.h" + #include "pmSourcePhotometry.h" #include "pmSourceVisual.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryModel.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryModel.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPSFtryModel.c (revision 28643) @@ -19,20 +19,26 @@ #include "pmFPA.h" #include "pmFPAMaskWeight.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmSourceUtils.h" #include "pmPSFtry.h" -#include "pmModelClass.h" -#include "pmModelUtils.h" -#include "pmSourceFitModel.h" -#include "pmSourcePhotometry.h" +#include "pmDetections.h" + #include "pmSourceVisual.h" @@ -49,4 +55,6 @@ pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) { + assert (options->fitOptions); + // validate the requested model name options->type = pmModelClassGetType (modelName); Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPeaks.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPeaks.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPeaks.c (revision 28643) @@ -22,4 +22,5 @@ #include #include "pmSpan.h" +#include "pmFootprintSpans.h" #include "pmFootprint.h" #include "pmPeaks.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPeaks.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPeaks.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPeaks.h (revision 28643) @@ -17,7 +17,4 @@ # ifndef PM_PEAKS_H # define PM_PEAKS_H - -#include -#include "pmFootprint.h" /// @addtogroup Objects Object Detection / Analysis Functions Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPhotObj.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPhotObj.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPhotObj.c (revision 28643) @@ -16,6 +16,27 @@ #include #include + +#include "pmHDU.h" +#include "pmFPA.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 "pmPhotObj.h" -#include "pmSource.h" + static void pmPhotObjFree (pmPhotObj *tmp) Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPhotObj.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPhotObj.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPhotObj.h (revision 28643) @@ -10,9 +10,4 @@ # ifndef PM_PHOT_OBJ_H # define PM_PHOT_OBJ_H - -#include -#include "pmPeaks.h" -#include "pmModel.h" -#include "pmMoments.h" /// @addtogroup Objects Object Detection / Analysis Functions Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.c (revision 28643) @@ -23,13 +23,20 @@ #include "pmFPA.h" #include "pmFPAMaskWeight.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.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" @@ -482,5 +489,5 @@ } psfClump.X = stats->clippedMean; - psfClump.dX = stats->clippedStdev; + psfClump.dX = hypot(stats->clippedStdev, PSF_CLUMP_GRID_SCALE); if (!psVectorStats (stats, tmpSy, NULL, NULL, 0)) { @@ -489,5 +496,5 @@ } psfClump.Y = stats->clippedMean; - psfClump.dY = stats->clippedStdev; + psfClump.dY = hypot(stats->clippedStdev, PSF_CLUMP_GRID_SCALE); psTrace ("psModules.objects", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.h (revision 28643) @@ -11,15 +11,6 @@ # define PM_SOURCE_H -#include -#include "pmPeaks.h" -#include "pmModel.h" -#include "pmMoments.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" - /// @addtogroup Objects Object Detection / Analysis Functions /// @{ - -#include /** pmSourceType enumeration Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceContour.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceContour.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceContour.c (revision 28643) @@ -23,14 +23,22 @@ #include "pmHDU.h" #include "pmFPA.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.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 "pmSourceContour.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitModel.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitModel.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitModel.c (revision 28643) @@ -23,37 +23,45 @@ #include "pmHDU.h" #include "pmFPA.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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.h" -#include "pmPSF.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 "pmModelClass.h" #include "pmSourceFitModel.h" -// save as static values so they may be set externally -static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; -static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; -static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0; -static bool PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true; - -bool pmSourceFitModelInit (float nIter, float tol, float weight, bool poissonErrors) +void pmSourceFitOptionsFree(pmSourceFitOptions *opt) { - - PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter; - PM_SOURCE_FIT_MODEL_TOLERANCE = tol; - PM_SOURCE_FIT_MODEL_WEIGHT = weight; - PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = poissonErrors; - - return true; + return; +} + +pmSourceFitOptions *pmSourceFitOptionsAlloc(void) { + + pmSourceFitOptions *opt = (pmSourceFitOptions *) psAlloc(sizeof(pmSourceFitOptions)); + psMemSetDeallocator(opt, (psFreeFunc) pmSourceFitOptionsFree); + + opt->mode = PM_SOURCE_FIT_PSF; + opt->nIter = 15; + opt->tol = 0.1; + opt->weight = 1.0; + opt->poissonErrors = true; + + return opt; } bool pmSourceFitModel (pmSource *source, pmModel *model, - pmSourceFitMode mode, + pmSourceFitOptions *options, psImageMaskType maskVal) { @@ -111,8 +119,8 @@ // as variance to avoid the bias from systematic errors here we would just use the // source sky variance - if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) { + if (options->poissonErrors) { yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j]; } else { - yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT; + yErr->data.F32[nPix] = 1.0 / options->weight; } nPix++; @@ -133,5 +141,5 @@ // set parameter mask based on fitting mode int nParams = 0; - switch (mode) { + switch (options->mode) { case PM_SOURCE_FIT_NORM: // NORM-only model fits only source normalization (Io) @@ -173,5 +181,5 @@ } - psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE); + psMinimization *myMin = psMinimizationAlloc (options->nIter, options->tol); psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitModel.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitModel.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitModel.h (revision 28643) @@ -22,10 +22,22 @@ } pmSourceFitMode; -bool pmSourceFitModelInit( - float nIter, ///< max number of allowed iterations - float tol, ///< convergence criterion - float weight, ///< use this weight for constant-weight fits - bool poissonErrors // use poisson errors for fits? -); +typedef struct { + pmSourceFitMode mode; ///< optionally fit all or a subset of parameters + float nIter; ///< max number of allowed iterations + float tol; ///< convergence criterion + float weight; ///< use this weight for constant-weight fits + bool poissonErrors; ///< use poisson errors for fits? +} pmSourceFitOptions; + +// the pmSourceFitOptions structure is used to control details of the fitting process +pmSourceFitOptions *pmSourceFitOptionsAlloc(void); + +// bool pmSourceFitModelInit( +// pmSourceFitMode mode, ///< what parameter set should be fitted? +// float nIter, ///< max number of allowed iterations +// float tol, ///< convergence criterion +// float weight, ///< use this weight for constant-weight fits +// bool poissonErrors ///< use poisson errors for fits? +// ); /** pmSourceFitModel() @@ -38,40 +50,39 @@ */ bool pmSourceFitModel( - pmSource *source, ///< The input pmSource - pmModel *model, ///< model to be fitted - pmSourceFitMode mode, ///< define parameters to be fitted + pmSource *source, ///< The input pmSource + pmModel *model, ///< model to be fitted + pmSourceFitOptions *options, ///< define parameters to be fitted psImageMaskType maskVal ///< Value to mask ); - -// initialize data for a group of object models -bool pmSourceFitSetInit (pmModelType type); - -// clear data for a group of object models -void pmSourceFitSetClear (void); - -// function used to set limits for a group of models -bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas); - -// function used to fit a group of object models -psF32 pmSourceFitSet_Function(psVector *deriv, - const psVector *params, - const psVector *x); - -/** pmSourceFitSet() - * - * Fit the requested model to the specified source. The starting guess for the model is given - * by the input source.model parameter values. The pixels of interest are specified by the - * source.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image - * data. The function returns TRUE on success or FALSE on failure. - * - */ -bool pmSourceFitSet( - pmSource *source, ///< The input pmSource - psArray *modelSet, ///< model to be fitted - pmSourceFitMode mode, ///< define parameters to be fitted - psImageMaskType maskVal ///< Vale to mask - -); +// // initialize data for a group of object models +// bool pmSourceFitSetInit (pmModelType type); +// +// // clear data for a group of object models +// void pmSourceFitSetClear (void); +// +// // function used to set limits for a group of models +// bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas); +// +// // function used to fit a group of object models +// psF32 pmSourceFitSet_Function(psVector *deriv, +// const psVector *params, +// const psVector *x); +// +// /** pmSourceFitSet() +// * +// * Fit the requested model to the specified source. The starting guess for the model is given +// * by the input source.model parameter values. The pixels of interest are specified by the +// * source.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image +// * data. The function returns TRUE on success or FALSE on failure. +// * +// */ +// bool pmSourceFitSet( +// pmSource *source, ///< The input pmSource +// psArray *modelSet, ///< model to be fitted +// pmSourceFitMode mode, ///< define parameters to be fitted +// psImageMaskType maskVal ///< Vale to mask +// +// ); /// @} Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitSet.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitSet.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitSet.c (revision 28643) @@ -22,23 +22,30 @@ #include "pmHDU.h" #include "pmFPA.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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.h" -#include "pmPSF.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 "pmModelClass.h" + #include "pmSourceFitModel.h" #include "pmSourceFitSet.h" // save as static values so they may be set externally -static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; -static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; -static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0; -static bool PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true; +// static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; +// static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; +// static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0; +// static bool PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true; /********************* Source Model Set Functions ***************************/ @@ -429,5 +436,5 @@ bool pmSourceFitSet (pmSource *source, psArray *modelSet, - pmSourceFitMode mode, + pmSourceFitOptions *options, psImageMaskType maskVal) { @@ -478,11 +485,11 @@ // as variance to avoid the bias from systematic errors here we would just use the // source sky variance - if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) { - yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j]; - } else { - yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT; - } - nPix++; - } + if (options->poissonErrors) { + yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j]; + } else { + yErr->data.F32[nPix] = 1.0 / options->weight; + } + nPix++; + } } x->n = nPix; @@ -490,49 +497,49 @@ yErr->n = nPix; - // create the FitSet for this thread and set the initial parameter guesses +// create the FitSet for this thread and set the initial parameter guesses pmSourceFitSetData *thisSet = pmSourceFitSetDataSet(modelSet); - // define param and deriv vectors for complete set of parameters +// define param and deriv vectors for complete set of parameters psVector *params = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32); - // set the param and deriv vectors based on the curent values +// set the param and deriv vectors based on the curent values pmSourceFitSetJoin (NULL, params, thisSet); - // create the minimization constraints +// create the minimization constraints psMinConstraint *constraint = psMinConstraintAlloc(); constraint->paramMask = psVectorAlloc (thisSet->nParamSet, PS_TYPE_VECTOR_MASK); constraint->checkLimits = pmSourceFitSetCheckLimits; - pmSourceFitSetMasks (constraint, thisSet, mode); - - // force the floating parameters to fall within the contraint ranges + pmSourceFitSetMasks (constraint, thisSet, options->mode); + +// force the floating parameters to fall within the contraint ranges for (int i = 0; i < params->n; i++) { - pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL); - pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL); + pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL); + pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL); } if (psTraceGetLevel("psModules.objects") >= 5) { - for (int i = 0; i < params->n; i++) { - fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]); - } + for (int i = 0; i < params->n; i++) { + fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]); + } } if (nPix < thisSet->nParamSet + 1) { - psTrace (__func__, 4, "insufficient valid pixels\n"); - psTrace("psModules.objects", 10, "---- %s() end : fail pixels ----\n", __func__); - for (int i = 0; i < modelSet->n; i++) { - pmModel *model = modelSet->data[i]; - model->flags |= PM_MODEL_STATUS_BADARGS; - } - psFree (x); - psFree (y); - psFree (yErr); - psFree (params); - psFree(constraint); - pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets - return(false); - } - - psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE); + psTrace (__func__, 4, "insufficient valid pixels\n"); + psTrace("psModules.objects", 10, "---- %s() end : fail pixels ----\n", __func__); + for (int i = 0; i < modelSet->n; i++) { + pmModel *model = modelSet->data[i]; + model->flags |= PM_MODEL_STATUS_BADARGS; + } + psFree (x); + psFree (y); + psFree (yErr); + psFree (params); + psFree(constraint); + pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets + return(false); + } + + psMinimization *myMin = psMinimizationAlloc (options->nIter, options->tol); psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); @@ -540,35 +547,35 @@ fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSetFunction); if (!fitStatus) { - psTrace("psModules.objects", 4, "Failed to fit model (%ld components)\n", modelSet->n); - } - - // parameter errors from the covariance matrix + psTrace("psModules.objects", 4, "Failed to fit model (%ld components)\n", modelSet->n); + } + +// parameter errors from the covariance matrix psVector *dparams = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32); for (int i = 0; i < dparams->n; i++) { - if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) - continue; - dparams->data.F32[i] = sqrt(covar->data.F32[i][i]); - } - - // get the Gauss-Newton distance for fixed model parameters + if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F32[i][i]); + } + +// get the Gauss-Newton distance for fixed model parameters if (constraint->paramMask != NULL) { - psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32); - psVector *altmask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK); - altmask->data.PS_TYPE_VECTOR_MASK_DATA[0] = 1; - for (int i = 1; i < dparams->n; i++) { - altmask->data.PS_TYPE_VECTOR_MASK_DATA[i] = (constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) ? 0 : 1; - } - psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSetFunction); - - for (int i = 0; i < dparams->n; i++) { - if (!constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) - continue; - // note that delta is the value *subtracted* from the parameter - // to get the new guess. for dparams to represent the direction - // of motion, we need to take -delta - dparams->data.F32[i] = -delta->data.F32[i]; - } - psFree (delta); - psFree (altmask); + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32); + psVector *altmask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK); + altmask->data.PS_TYPE_VECTOR_MASK_DATA[0] = 1; + for (int i = 1; i < dparams->n; i++) { + altmask->data.PS_TYPE_VECTOR_MASK_DATA[i] = (constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) ? 0 : 1; + } + psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSetFunction); + + for (int i = 0; i < dparams->n; i++) { + if (!constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) + continue; + // note that delta is the value *subtracted* from the parameter + // to get the new guess. for dparams to represent the direction + // of motion, we need to take -delta + dparams->data.F32[i] = -delta->data.F32[i]; + } + psFree (delta); + psFree (altmask); } Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitSet.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitSet.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitSet.h (revision 28643) @@ -56,5 +56,5 @@ pmSource *source, ///< The input pmSource psArray *modelSet, ///< model to be fitted - pmSourceFitMode mode, ///< define parameters to be fitted + pmSourceFitOptions *options, ///< define options for fitting process psImageMaskType maskVal ///< Vale to mask Index: anches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceGroup.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceGroup.h (revision 28642) +++ (revision ) @@ -1,44 +1,0 @@ -/* @file pmSourceGroup.h - * - * @author EAM, IfA - * - * @version $Revision: $ - * @date $Date: 2009-02-16 22:30:50 $ - * Copyright 2009 Institute for Astronomy, University of Hawaii - */ - -# ifndef PM_SOURCE_GROUP_H -# define PM_SOURCE_GROUP_H - -#include -#include "pmPeaks.h" -#include "pmModel.h" -#include "pmMoments.h" -#include "pmSourceExtendedPars.h" - -/// @addtogroup Objects Object Detection / Analysis Functions -/// @{ - -#include - -/** pmSourceGroup data structure - * - * A source group is a connected set of source measurements with a - * common connection. Each source in the group represents the - * detection of some astronomical object on a single each. The group - * represents the related collection of measurements. The fits are - * coupled in some way. For example, they may all have the same - * position, but independent fluxes. Or, they may have a common set - * of positions and shape parameters. Or the position in each image - * may be related by a function. - * - * XXX is thre any info that is neaded for each source group beyond that carried by the sources (besides ID)? - */ -struct pmSource { - int seq; ///< ID for output (generated on write OR set on read) - psArray *sources; -} pmSourceGroup; - - -/// @} -# endif /* PM_SOURCE_GROUP_H */ Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceGroups.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceGroups.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceGroups.c (revision 28643) @@ -6,6 +6,28 @@ #include +#include "pmHDU.h" #include "pmFPA.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 "pmDetections.h" + #include "pmSourceGroups.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceGroups.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceGroups.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceGroups.h (revision 28643) @@ -1,8 +1,4 @@ #ifndef PM_SOURCE_GROUPS_H #define PM_SOURCE_GROUPS_H - -#include - -#include "pmFPA.h" /// Groups of sources Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO.c (revision 28643) @@ -31,16 +31,25 @@ #include "pmConceptsRead.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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" +#include "pmPSFtry.h" + #include "pmDetections.h" -#include "pmSource.h" -#include "pmModelClass.h" #include "pmDetEff.h" #include "pmSourceIO.h" @@ -326,4 +335,15 @@ } +# define PM_SOURCES_WRITE(NAME,TYPE) \ + if (!strcmp (exttype, NAME)) { \ + status &= pmSourcesWrite_##TYPE(file->fits, readout, sources, file->header, outhead, dataname, recipe); \ + if (xsrcname) { \ + status &= pmSourcesWrite_##TYPE##_XSRC(file->fits, readout, sources, file->header, xsrcname, recipe); \ + } \ + if (xfitname) { \ + status &= pmSourcesWrite_##TYPE##_XFIT (file->fits, readout, sources, file->header, xfitname); \ + } \ + } + // write out all readout-level Objects files for this cell bool pmReadoutWriteObjects (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) @@ -360,4 +380,5 @@ } + // the older types (RAW, OBJ, SX, CMP) are for backwards compatibility -- deprecate eventually? switch (file->type) { case PM_FPA_FILE_RAW: @@ -518,72 +539,18 @@ psMetadataAddStr (outhead, PS_LIST_TAIL, "XFITNAME", PS_META_REPLACE, "name of XFIT table extension", xfitname); } - - // XXX these are case-sensitive since the EXTYPE is case-sensitive + + + // these are case-sensitive since the EXTYPE is case-sensitive status = true; - if (!strcmp (exttype, "SMPDATA")) { - status &= pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname); - } - if (!strcmp (exttype, "PS1_DEV_0")) { - status &= pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname); - } - if (!strcmp (exttype, "PS1_DEV_1")) { - status &= pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname); - } - if (!strcmp (exttype, "PS1_CAL_0")) { - status &= pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname); - } - if (!strcmp (exttype, "PS1_V1")) { - status &= pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname); - } - if (!strcmp (exttype, "PS1_V2")) { - status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname); - } - if (!strcmp (exttype, "PS1_SV1")) { - status &= pmSourcesWrite_CMF_PS1_SV1 (file->fits, readout, sources, file->header, outhead, dataname, recipe); - } - if (!strcmp (exttype, "PS1_DV1")) { - status &= pmSourcesWrite_CMF_PS1_DV1 (file->fits, readout, sources, file->header, outhead, dataname); - } - - if (xsrcname) { - if (!strcmp (exttype, "PS1_DEV_1")) { - status &= pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe); - } - if (!strcmp (exttype, "PS1_CAL_0")) { - status &= pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); - } - if (!strcmp (exttype, "PS1_V1")) { - status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); - } - if (!strcmp (exttype, "PS1_V2")) { - status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); - } - if (!strcmp (exttype, "PS1_SV1")) { - status &= pmSourcesWrite_CMF_PS1_SV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); - } - if (!strcmp (exttype, "PS1_DV1")) { - status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe); - } - } - if (xfitname) { - if (!strcmp (exttype, "PS1_DEV_1")) { - status &= pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname); - } - if (!strcmp (exttype, "PS1_CAL_0")) { - status &= pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname); - } - if (!strcmp (exttype, "PS1_V1")) { - status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, readout, sources, xfitname); - } - if (!strcmp (exttype, "PS1_V2")) { - status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, readout, sources, xfitname); - } - if (!strcmp (exttype, "PS1_SV1")) { - status &= pmSourcesWrite_CMF_PS1_SV1_XFIT (file->fits, readout, sources, xfitname); - } - if (!strcmp (exttype, "PS1_DV1")) { - status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, readout, sources, xfitname); - } - } + PM_SOURCES_WRITE("SMPDATA", SMPDATA); + PM_SOURCES_WRITE("PS1_DEV_0", PS1_DEV_0); + PM_SOURCES_WRITE("PS1_DEV_1", PS1_DEV_1); + PM_SOURCES_WRITE("PS1_CAL_0", PS1_CAL_0); + PM_SOURCES_WRITE("PS1_V1", CMF_PS1_V1); + PM_SOURCES_WRITE("PS1_V2", CMF_PS1_V2); + PM_SOURCES_WRITE("PS1_V3", CMF_PS1_V3); + PM_SOURCES_WRITE("PS1_SV1", CMF_PS1_SV1); + PM_SOURCES_WRITE("PS1_DV1", CMF_PS1_DV1); + psFree (outhead); psFree (exttype); Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO.h (revision 28643) @@ -24,32 +24,42 @@ bool pmSourcesWriteCMP (psArray *sources, char *filename, psMetadata *header); -bool pmSourcesWrite_SMPDATA (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); +bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config); -bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); -bool pmSourcesWrite_PS1_DEV_1_XFIT (psFits *fits, psArray *sources, char *extname); +// All of these functions need to use the same API, even if not all elements are used in a specific case +bool pmSourcesWrite_SMPDATA(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_SMPDATA_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_SMPDATA_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); -bool pmSourcesWrite_PS1_CAL_0 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_PS1_CAL_0_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); -bool pmSourcesWrite_PS1_CAL_0_XFIT (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); +bool pmSourcesWrite_PS1_DEV_0(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_PS1_DEV_0_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_PS1_DEV_0_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); +bool pmSourcesWrite_PS1_DEV_1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_PS1_DEV_1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_PS1_DEV_1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); +bool pmSourcesWrite_PS1_CAL_0(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_PS1_CAL_0_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_PS1_CAL_0_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_SV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_V1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); -bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); -bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_V2(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V2_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); -bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config); +bool pmSourcesWrite_CMF_PS1_V3(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V3_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_V3_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); + +bool pmSourcesWrite_CMF_PS1_SV1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_SV1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); + +bool pmSourcesWrite_CMF_PS1_DV1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_DV1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_DV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); psArray *pmSourcesReadCMP (char *filename, psMetadata *header); Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" @@ -47,6 +56,5 @@ // This version has elements intended for difference images & forced photometry: // diffStats entries (good for dipoles); flux + flux error (for insignificant detections) -bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, - psMetadata *imageHeader, psMetadata *tableHeader, char *extname) +bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -536,5 +544,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) +bool pmSourcesWrite_CMF_PS1_DV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) { Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" @@ -46,5 +55,5 @@ // NOTE: this output function is intended for psphotStack analysis: it includes per-psf radial fluxes -// XXX currently in the 'read' function is NOT consistent with the 'write' function (does not read radial fluxes) +// XXX currently, the 'read' function is NOT consistent with the 'write' function (does not read radial fluxes) bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) @@ -616,5 +625,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) +bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) { Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" @@ -45,6 +54,5 @@ // followed by a zero-size matrix, followed by the table data -bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, - psMetadata *imageHeader, psMetadata *tableHeader, char *extname) +bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -516,5 +524,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) +bool pmSourcesWrite_CMF_PS1_V1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) { Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" @@ -45,5 +54,5 @@ // followed by a zero-size matrix, followed by the table data -bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname) +bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -578,5 +587,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname) +bool pmSourcesWrite_CMF_PS1_V2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) { Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 28643) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 28643) @@ -0,0 +1,766 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 02:44:19 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#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" + +// panstarrs-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data + +bool pmSourcesWrite_CMF_PS1_V3 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(extname, false); + + 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 (sources->n > 0) { + pmSource *source = (pmSource *) sources->data[0]; + if (source->seq == -1) { + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + } else { + sources = psArraySort (sources, pmSourceSortBySeq); + } + } + + 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 (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. + 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); + + 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, "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_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, "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, "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, "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_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeight); + 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, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + + // XXX not sure how to get this : need to load Nimages with weight? + psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap); + psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); + + psArrayAdd (table, 100, row); + psFree (row); + + // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not + // subtracted + + // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image + // edge; 3) any pixels in the 3x3 peak region are masked; + } + + // XXX why do we make a copy here to be supplemented with the masks? why not do this in the calling function? + psMetadata *header = psMetadataCopy(NULL, tableHeader); + pmSourceMasksHeader(header); + + if (table->n == 0) { + if (!psFitsWriteBlank(fits, header, extname)) { + psError(psErrorCodeLast(), false, "Unable to write blank sources file."); + psFree(table); + psFree(header); + return false; + } + psFree(table); + psFree(header); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable(fits, header, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree(table); + psFree(header); + return false; + } + psFree(table); + psFree(header); + + return true; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_CMF_PS1_V3 (psFits *fits, psMetadata *header) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(header, false); + + bool status; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + + // define PSF model type + // XXX need to carry the extra model parameters + int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelClassGetType (PSF_NAME); + } + assert (modelType > -1); + + // We get the size of the table, and allocate the array of sources first because the table + // is large and ephemeral --- when the table gets blown away, whatever is allocated after + // the table is read blocks the free. In fact, it's better to read the table row by row. + long numSources = psFitsTableSize(fits); // Number of sources in table + psArray *sources = psArrayAlloc(numSources); // Array of sources, to return + + // convert the table to the pmSource entriesa + for (int i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %d of sources", i); + psFree(sources); + return NULL; + } + + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags + + // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); + dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); + dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); + axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); + axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); + axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); + + PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); + dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); + source->sky = PAR[PM_PAR_SKY]; + source->skyErr = dPAR[PM_PAR_SKY]; + + // XXX use these to determine PAR[PM_PAR_I0]? + source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); + source->errMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); + source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); + + // XXX this scaling is incorrect: does not include the 2 \pi AREA factor + PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; + dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN; + + pmPSF_AxesToModel (PAR, axes); + + float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); + float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; + + // 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->dx = dPAR[PM_PAR_XPOS]; + source->peak->dy = dPAR[PM_PAR_YPOS]; + source->peak->SN = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N + + source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF"); + source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); + source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); + + // note that some older versions used PSF_PROBABILITY: this was not well defined. + model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); + model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); + model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); + + source->moments = pmMomentsAlloc (); + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); + source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); + source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); + + source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); + source->moments->Mrh = psMetadataLookupF32 (&status, row, "MOMENTS_RH"); + source->moments->KronFlux = psMetadataLookupF32 (&status, row, "KRON_FLUX"); + source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR"); + + source->moments->KronFinner = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER"); + source->moments->KronFouter = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER"); + + // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data, + // we are storing enough information so the output will be consistent with the input + source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); + source->moments->Mxxy = 0.0; + source->moments->Mxyy = 0.0; + source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S"); + + source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C"); + source->moments->Mxxxy = 0.0; + source->moments->Mxxyy = 0.0; + source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S"); + source->moments->Myyyy = 0.0; + + source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); + assert (status); + + sources->data[i] = source; + psFree(row); + } + + return sources; +} + +bool pmSourcesWrite_CMF_PS1_V3_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) +{ + + bool status; + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psF32 xPos, yPos; + psF32 xErr, yErr; + int nRow = -1; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + pmChip *chip = readout->parent->parent; + pmFPA *fpa = chip->parent; + + // zero point corrections + bool status1 = false; + bool status2 = false; + float magOffset = 0.0; + 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); + } + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + table = psArrayAllocEmpty (sources->n); + + // which extended source analyses should we perform? + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + + psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + psAssert (radMin->n == radMax->n, "inconsistent annular bins"); + + // int nRadialBins = 0; + // if (doAnnuli) { + // // get the max count of radial bins + // for (int i = 0; i < sources->n; i++) { + // pmSource *source = sources->data[i]; + // if (!source->extpars) continue; + // if (!source->extpars->radProfile ) continue; + // if (!source->extpars->radProfile->binSB) continue; + // nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n); + // } + // } + + // we write out all sources, regardless of quality. the source flags tell us the state + for (int i = 0; i < sources->n; i++) { + // skip source if it is not a ext sourc + // XXX we have two places that extended source parameters are measured: + // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, + // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) + // should we require both? + + pmSource *source = sources->data[i]; + + // skip sources without measurements + if (source->extpars == NULL) continue; + + // we require a PSF model fit (ignore the real crud) + pmModel *model = source->modelPSF; + if (model == NULL) continue; + + // XXX I need to split the extended models from the extended aperture measurements + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); + + float AxialRatio = NAN; + float AxialTheta = NAN; + pmSourceExtendedPars *extpars = source->extpars; + if (extpars) { + AxialRatio = extpars->axes.minor / extpars->axes.major; + AxialTheta = extpars->axes.theta; + } + 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); + + // Petrosian measurements + // XXX insert header data: petrosian ref radius, flux ratio + // XXX check flags to see if Pet was measured + if (doPetrosian) { + pmSourceExtendedPars *extpars = source->extpars; + if (extpars) { + // XXX note that this mag is either calibrated or instrumental depending on existence of zero point + 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_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); + } 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_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); + } + } + +# if (0) + // Kron measurements + if (doKron) { + pmSourceKronValues *kron = source->extpars->kron; + if (kron) { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", kron->mag); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", kron->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); + } + } + + // Isophot measurements + // XXX insert header data: isophotal level + if (doIsophotal) { + pmSourceIsophotalValues *isophot = source->extpars->isophot; + if (isophot) { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", isophot->mag); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); + } + } +# endif + + // Flux Annuli (if we have extended source measurements, we have these. only optionally save them) + if (doAnnuli) { + psVector *radSB = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVectorInit (radSB, NAN); + psVectorInit (radFlux, NAN); + psVectorInit (radFill, NAN); + if (!source->extpars) goto empty_annuli; + if (!source->extpars->radProfile) goto empty_annuli; + if (!source->extpars->radProfile->binSB) goto empty_annuli; + psAssert (source->extpars->radProfile->binSum, "programming error"); + psAssert (source->extpars->radProfile->binFill, "programming error"); + psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->extpars->radProfile->binFill->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->extpars->radProfile->binSB->n; j++) { + radSB->data.F32[j] = source->extpars->radProfile->binSB->data.F32[j]; + radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j]; + radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j]; + } + + empty_annuli: + psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB); + psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); + psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psFree (radSB); + psFree (radFlux); + psFree (radFill); + } + if (nRow < 0) { + nRow = row->list->n; + } else { + psAssert (nRow == row->list->n, "inconsistent row lengths"); + } + psArrayAdd (table, 100, row); + psFree (row); + } + + if (table->n == 0) { + if (!psFitsWriteBlank (fits, outhead, extname)) { + psError(psErrorCodeLast(), false, "Unable to write empty sources file."); + psFree(outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + + return true; +} + +// XXX this layout is still the same as PS1_DEV_1 +bool pmSourcesWrite_CMF_PS1_V3_XFIT (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) +{ + + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + char name[64]; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + // 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]; + if (source->modelFits == NULL) continue; + for (int j = 0; j < source->modelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + assert (model); + nParamMax = PS_MAX (nParamMax, model->params->n); + } + } + + table = psArrayAllocEmpty (sources->n); + + // 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]; + + // XXX if no model fits are saved, write out modelEXT? + if (source->modelFits == NULL) continue; + + // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) + for (int j = 0; j < source->modelFits->n; j++) { + + // choose the convolved EXT model, if available, otherwise the simple one + pmModel *model = source->modelFits->data[j]; + assert (model); + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + axes = pmPSF_ModelToAxes (PAR, 20.0); + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); + + psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); + psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); + + // XXX these should be major and minor, not 'x' and 'y' + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); + + // write out the other generic parameters + for (int k = 0; k < nParamMax; k++) { + if (k == PM_PAR_I0) continue; + if (k == PM_PAR_SKY) continue; + if (k == PM_PAR_XPOS) continue; + if (k == PM_PAR_YPOS) continue; + if (k == PM_PAR_SXX) continue; + if (k == PM_PAR_SXY) continue; + if (k == PM_PAR_SYY) continue; + + snprintf (name, 64, "EXT_PAR_%02d", k); + + if (k < model->params->n) { + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); + } + } + + // XXX other parameters which may be set. + // XXX flag / value to define the model + // XXX write out the model type, fit status flags + + psArrayAdd (table, 100, row); + psFree (row); + } + } + + if (table->n == 0) { + if (!psFitsWriteBlank (fits, outhead, extname)) { + psError(psErrorCodeLast(), false, "Unable to write empty sources file."); + psFree(outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; +} Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMP.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMP.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_CMP.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_MatchedRefs.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 28643) @@ -28,16 +28,26 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" + #include "pmAstrometryObjects.h" #include "pmAstrometryWCS.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_OBJ.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_OBJ.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_OBJ.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" @@ -49,6 +58,5 @@ // XXX how do I generate the source tables which I need to send to PSPS? -bool pmSourcesWrite_PS1_CAL_0 (psFits *fits, pmReadout *readout, psArray *sources, - psMetadata *imageHeader, psMetadata *tableHeader, char *extname) +bool pmSourcesWrite_PS1_CAL_0 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_PS1_DEV_0.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" @@ -49,6 +58,5 @@ // XXX how do I generate the source tables which I need to send to PSPS? // XXX: input parameter imageHeader is never used. -bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *imageHeader, - psMetadata *tableHeader, char *extname) +bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -228,2 +236,12 @@ return (sources); } + +bool pmSourcesWrite_PS1_DEV_0_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) +{ + return true; +} + +bool pmSourcesWrite_PS1_DEV_0_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) +{ + return true; +} Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" @@ -47,8 +56,6 @@ // this output format is valid for psphot analysis of an image, and does not include calibrated // values derived in the DVO database. -// XXX how do I generate the source tables which I need to send to PSPS? - -bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, - psMetadata *imageHeader, psMetadata *tableHeader, char *extname) + +bool pmSourcesWrite_PS1_DEV_1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -281,5 +288,5 @@ } -bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { @@ -453,5 +460,5 @@ } -bool pmSourcesWrite_PS1_DEV_1_XFIT (psFits *fits, psArray *sources, char *extname) +bool pmSourcesWrite_PS1_DEV_1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) { Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_RAW.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_RAW.c (revision 28643) @@ -28,16 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" -#include "pmSourcePhotometry.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_SMPDATA.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" @@ -45,6 +54,5 @@ // followed by a zero-size matrix, followed by the table data // XXX: input parameter imageHeader is never used -bool pmSourcesWrite_SMPDATA (psFits *fits, psArray *sources, psMetadata *imageHeader, - psMetadata *tableHeader, char *extname) +bool pmSourcesWrite_SMPDATA (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -204,2 +212,12 @@ return (sources); } + +bool pmSourcesWrite_SMPDATA_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) +{ + return true; +} + +bool pmSourcesWrite_SMPDATA_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) +{ + return true; +} Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_SX.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_SX.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceIO_SX.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourceIO.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceMatch.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceMatch.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceMatch.c (revision 28643) @@ -6,6 +6,28 @@ #include #include - +#include "pmHDU.h" +#include "pmFPA.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 "pmDetections.h" + #include "pmErrorCodes.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceMatch.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceMatch.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceMatch.h (revision 28643) @@ -1,6 +1,4 @@ #ifndef PM_SOURCE_MATCH_H #define PM_SOURCE_MATCH_H - -#include /// Mask values for matched sources Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceMoments.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceMoments.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceMoments.c (revision 28643) @@ -22,16 +22,23 @@ #include #include + #include "pmHDU.h" #include "pmFPA.h" -#include "pmFPAMaskWeight.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.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" @@ -54,4 +61,7 @@ # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) +static bool beVerbose = false; +void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; } + bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal) { @@ -61,11 +71,15 @@ PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); - // use sky from moments if defined, 0.0 otherwise + // use sky from moments if defined, 0.0 otherwise + + // XXX this value comes from the sky model at the source center, and tends to over-estimate + // the sky in the vicinity of bright sources. we are better off assuming the model worked + // well: psF32 sky = 0.0; - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } else { - sky = source->moments->Sky; - } + // XXX if (source->moments == NULL) { + // XXX source->moments = pmMomentsAlloc(); + // XXX } else { + // XXX sky = source->moments->Sky; + // XXX } // First Pass: calculate the first moments (these are subtracted from the coordinates below) @@ -197,4 +211,7 @@ // Xn = SUM (x - xc)^n * (z - sky) + psF32 RF = 0.0; + psF32 RH = 0.0; + psF32 RS = 0.0; psF32 XX = 0.0; psF32 XY = 0.0; @@ -244,5 +261,6 @@ if (r > radius) continue; - psF32 pDiff = *vPix - sky; + psF32 fDiff = *vPix - sky; + psF32 pDiff = fDiff; psF32 wDiff = *vWgt; @@ -257,5 +275,7 @@ if (sigma > 0.0) { // XXX a lot of extra flops; can we do pre-calculate? - psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + // XXX we were re-calculating r2 (maybe the compiler caught this?) + // psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + psF32 z = r2 * rsigma2; assert (z >= 0.0); psF32 weight = exp(-z); @@ -266,4 +286,17 @@ Sum += pDiff; + +# if (1) +# if (0) + if (fDiff < 0) continue; +# endif + psF32 rf = r * fDiff; + psF32 rh = sqrt(r) * fDiff; + psF32 rs = fDiff; +# else + psF32 rf = r * pDiff; + psF32 rh = sqrt(r) * pDiff; + psF32 rs = pDiff; +# endif psF32 x = xDiff * pDiff; @@ -284,4 +317,8 @@ psF32 xyyy = xDiff * yyy / r2; psF32 yyyy = yDiff * yyy / r2; + + RF += rf; + RH += rh; + RS += rs; XX += xx; @@ -302,4 +339,7 @@ } + source->moments->Mrf = RF/RS; + source->moments->Mrh = RH/RS; + source->moments->Mxx = XX/Sum; source->moments->Mxy = XY/Sum; @@ -317,12 +357,71 @@ source->moments->Myyyy = YYYY/Sum; - // if (source->moments->Mxx < 0) { - // fprintf (stderr, "error: neg second moment??\n"); - // } - // if (source->moments->Myy < 0) { - // fprintf (stderr, "error: neg second moment??\n"); - // } - - psTrace ("psModules.objects", 4, "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", + // Calculate the Kron magnitude (make this block optional?) + // float radKron = 2.5*source->moments->Mrf; + float radKinner = 1.0*source->moments->Mrf; + float radKron = 2.5*source->moments->Mrf; + float radKouter = 4.0*source->moments->Mrf; + + int nKronPix = 0; + Sum = Var = 0.0; + float SumInner = 0.0; + float SumOuter = 0.0; + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 yDiff = row - yCM; + if (fabs(yDiff) > radKron) continue; + + psF32 *vPix = source->pixels->data.F32[row]; + psF32 *vWgt = source->variance->data.F32[row]; + psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; + + for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { + if (vMsk) { + if (*vMsk & maskVal) { + vMsk++; + continue; + } + vMsk++; + } + if (isnan(*vPix)) continue; + + psF32 xDiff = col - xCM; + if (fabs(xDiff) > radKron) continue; + + // radKron is just a function of (xDiff, yDiff) + psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + psF32 r = sqrt(r2); + + psF32 pDiff = *vPix - sky; + psF32 wDiff = *vWgt; + + // skip pixels below specified significance level. this is allowed, but should be + // avoided -- the over-weights the wings of bright stars compared to those of faint + // stars. + if (PS_SQR(pDiff) < minSN2*wDiff) continue; + + if (r < radKron) { + Sum += pDiff; + Var += wDiff; + nKronPix ++; + // if (beVerbose) fprintf (stderr, "mome: %d %d %f %f %f\n", col, row, sky, *vPix, Sum); + } + + if ((r > radKinner) && (r < radKron)) { + SumInner += pDiff; + } + if ((r > radKron) && (r < radKouter)) { + SumOuter += pDiff; + } + } + } + source->moments->KronFlux = Sum; + source->moments->KronFinner = SumInner; + source->moments->KronFouter = SumOuter; + source->moments->KronFluxErr = sqrt(Var); + + 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", + source->moments->Mrf, source->moments->KronFlux, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePhotometry.c (revision 28643) @@ -22,15 +22,25 @@ #include "pmFPA.h" #include "pmFPAMaskWeight.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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" -#include "pmSource.h" -#include "pmModelClass.h" +#include "pmPSFtry.h" + #include "pmSourcePhotometry.h" @@ -257,5 +267,8 @@ PS_ASSERT_PTR_NON_NULL(image, false); PS_ASSERT_PTR_NON_NULL(mask, false); - PS_ASSERT_PTR_NON_NULL(model, false); + + if (DO_SKY) { + PS_ASSERT_PTR_NON_NULL(model, false); + } float apSum = 0; @@ -271,11 +284,14 @@ psF32 **imData = image->data.F32; psImageMaskType **mkData = mask->data.PS_TYPE_IMAGE_MASK_DATA; + int nAperPix = 0; // measure apMag - for (int ix = 0; ix < image->numCols; ix++) { - for (int iy = 0; iy < image->numRows; iy++) { + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { if (mkData[iy][ix] & maskVal) continue; apSum += imData[iy][ix] - sky; + nAperPix ++; + // fprintf (stderr, "aper: %d %d %f %f %f\n", ix, iy, sky, imData[iy][ix], apSum); } } Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePlotApResid.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePlotApResid.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePlotApResid.c (revision 28643) @@ -25,15 +25,25 @@ #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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.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 "pmModel.h" #include "pmDetections.h" -#include "pmSource.h" + #include "pmSourcePlots.h" #include "pmKapaPlots.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePlotMoments.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePlotMoments.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePlotMoments.c (revision 28643) @@ -28,15 +28,24 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" #include "pmDetections.h" -#include "pmSource.h" #include "pmSourcePlots.h" #include "pmKapaPlots.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePlotPSFModel.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePlotPSFModel.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourcePlotPSFModel.c (revision 28643) @@ -28,15 +28,25 @@ #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 "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmTrend2D.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 "pmModel.h" #include "pmDetections.h" -#include "pmSource.h" #include "pmSourcePlots.h" #include "pmKapaPlots.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceSky.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceSky.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceSky.c (revision 28643) @@ -21,17 +21,25 @@ #include #include + #include "pmHDU.h" #include "pmFPA.h" -#include "pmFPAMaskWeight.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.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 "pmSourceSky.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceUtils.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceUtils.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceUtils.c (revision 28643) @@ -21,17 +21,25 @@ #include #include + #include "pmHDU.h" #include "pmFPA.h" -#include "pmFPAMaskWeight.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 "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.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 "pmSourceUtils.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceVisual.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceVisual.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceVisual.c (revision 28643) @@ -4,8 +4,28 @@ #include +#include "pmHDU.h" +#include "pmFPA.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 "pmSource.h" +#include "pmDetections.h" + #include "pmSourceVisual.h" Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSpan.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSpan.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSpan.h (revision 28643) @@ -10,6 +10,4 @@ # ifndef PM_SPAN_H # define PM_SPAN_H - -#include /// @addtogroup Objects Object Detection / Analysis Functions Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmTrend2D.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmTrend2D.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmTrend2D.h (revision 28643) @@ -12,6 +12,4 @@ # ifndef PM_TREND_2D_H # define PM_TREND_2D_H - -#include /// @addtogroup Objects Object Detection / Analysis Functions Index: /branches/eam_branches/ipp-20100621/psModules/src/psmodules.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/psmodules.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psModules/src/psmodules.h (revision 28643) @@ -114,4 +114,8 @@ // the following headers are from psModule:objects +#include +#include +#include + #include #include @@ -120,23 +124,23 @@ #include #include + +#include +#include + +#include #include #include -#include -#include -#include +#include +#include #include -#include -#include -#include +#include #include #include #include #include -#include #include #include #include #include -#include #include #include Index: /branches/eam_branches/ipp-20100621/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphot.h (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphot.h (revision 28643) @@ -211,9 +211,9 @@ bool psphotInitLimitsPSF (psMetadata *recipe, pmReadout *readout); bool psphotInitLimitsEXT (psMetadata *recipe); -bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal); -bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal); -bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal); -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal); -psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal); +bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); +bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); +bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal); +psArray *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); // functions to support simultaneous multi-source fitting @@ -252,4 +252,5 @@ bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout); bool psphotVisualShowSignificance (psImage *image, float min, float max); +bool psphotVisualShowLogSignificance (psImage *image, float min, float max); bool psphotVisualShowPeaks (pmDetections *detections); bool psphotVisualShowFootprints (pmDetections *detections); Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotBlendFit.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotBlendFit.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotBlendFit.c (revision 28643) @@ -75,5 +75,10 @@ // Define source fitting parameters for extended source fits - pmSourceFitModelInit(fitIter, fitTol, PS_SQR(skySig), poisson); + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->nIter = fitIter; + fitOptions->tol = fitTol; + fitOptions->poissonErrors = poisson; + fitOptions->weight = PS_SQR(skySig); + fitOptions->mode = PM_SOURCE_FIT_PSF; psphotInitLimitsPSF (recipe, readout); @@ -88,4 +93,5 @@ if (!sources->n) { psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + psFree (fitOptions); return true; } @@ -112,4 +118,5 @@ psArrayAdd(job->args, 1, psf); psArrayAdd(job->args, 1, newSources); // return for new sources + psArrayAdd(job->args, 1, fitOptions); // default fit options psFree (newSources); @@ -121,31 +128,7 @@ if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (fitOptions); return NULL; } - -# if (0) - { - int nfit = 0; - int npsf = 0; - int next = 0; - int nfail = 0; - psArray *newSources = psArrayAllocEmpty(16); - - if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return NULL; - } - Nfit += nfit; - Npsf += npsf; - Next += next; - Nfail += nfail; - - // add these back onto sources - for (int k = 0; k < newSources->n; k++) { - psArrayAdd (sources, 16, newSources->data[k]); - } - psFree (newSources); - } -# endif } @@ -153,4 +136,5 @@ if (!psThreadPoolWait (false)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (fitOptions); return NULL; } @@ -163,11 +147,11 @@ } else { psScalar *scalar = NULL; - scalar = job->args->data[5]; + scalar = job->args->data[6]; Nfit += scalar->data.S32; - scalar = job->args->data[6]; + scalar = job->args->data[7]; Npsf += scalar->data.S32; - scalar = job->args->data[7]; + scalar = job->args->data[8]; Next += scalar->data.S32; - scalar = job->args->data[8]; + scalar = job->args->data[9]; Nfail += scalar->data.S32; @@ -186,4 +170,5 @@ psphotSaveImage (NULL, readout->image, "image.v2.fits"); } + psFree (fitOptions); psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit); @@ -204,9 +189,10 @@ psScalar *scalar = NULL; - pmReadout *readout = job->args->data[0]; - psMetadata *recipe = job->args->data[1]; - psArray *sources = job->args->data[2]; - pmPSF *psf = job->args->data[3]; - psArray *newSources = job->args->data[4]; + pmReadout *readout = job->args->data[0]; + psMetadata *recipe = job->args->data[1]; + psArray *sources = job->args->data[2]; + pmPSF *psf = job->args->data[3]; + psArray *newSources = job->args->data[4]; + pmSourceFitOptions *fitOptions = job->args->data[5]; // bit-masks to test for good/bad pixels @@ -269,8 +255,14 @@ Nfit ++; + if (0) { + psF32 Mxx = source->moments->Mxx; + psF32 Myy = source->moments->Myy; + fprintf (stderr, "1: Mxx: %f, Myy: %f\n", Mxx, Myy); + } + // try fitting PSFs or extended sources depending on source->mode // these functions subtract the resulting fitted source if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) { + if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) { source->type = PM_SOURCE_TYPE_EXTENDED; psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); @@ -280,5 +272,5 @@ } } else { - if (psphotFitBlend (readout, source, psf, maskVal, markVal)) { + if (psphotFitBlend (readout, source, psf, fitOptions, maskVal, markVal)) { source->type = PM_SOURCE_TYPE_STAR; psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); @@ -289,4 +281,10 @@ } + if (0) { + psF32 Mxx = source->moments->Mxx; + psF32 Myy = source->moments->Myy; + fprintf (stderr, "2: Mxx: %f, Myy: %f\n", Mxx, Myy); + } + psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf); Nfail ++; @@ -298,14 +296,14 @@ // change the value of a scalar on the array (wrap this and put it in psArray.h) - scalar = job->args->data[5]; + scalar = job->args->data[6]; scalar->data.S32 = Nfit; - scalar = job->args->data[6]; + scalar = job->args->data[7]; scalar->data.S32 = Npsf; - scalar = job->args->data[7]; + scalar = job->args->data[8]; scalar->data.S32 = Next; - scalar = job->args->data[8]; + scalar = job->args->data[9]; scalar->data.S32 = Nfail; Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotChoosePSF.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotChoosePSF.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotChoosePSF.c (revision 28643) @@ -74,5 +74,5 @@ // structure to store user options defining the psf - pmPSFOptions *options = pmPSFOptionsAlloc (); + pmPSFOptions *options = pmPSFOptionsAlloc(); // load user options from the recipe. no need to check existence -- they are @@ -143,6 +143,13 @@ return false; } - pmSourceFitModelInit(fitIter, fitTol, PS_SQR(SKY_SIG), options->poissonErrorsPhotLMM); - + + // options which modify the behavior of the model fitting + options->fitOptions = pmSourceFitOptionsAlloc(); + options->fitOptions->nIter = fitIter; + options->fitOptions->tol = fitTol; + options->fitOptions->poissonErrors = options->poissonErrorsPhotLMM; + options->fitOptions->weight = PS_SQR(SKY_SIG); + options->fitOptions->mode = PM_SOURCE_FIT_PSF; + psArray *stars = psArrayAllocEmpty (sources->n); Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotExtendedSourceFits.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotExtendedSourceFits.c (revision 28643) @@ -129,4 +129,13 @@ sources = psArraySort (sources, pmSourceSortBySN); + // Define source fitting parameters for extended source fits + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->mode = PM_SOURCE_FIT_EXT; + // XXX for now, use the defaults for the rest: + // fitOptions->nIter = fitIter; + // fitOptions->tol = fitTol; + // fitOptions->poissonErrors = poisson; + // fitOptions->weight = PS_SQR(skySig); + // choose the sources of interest for (int i = 0; i < sources->n; i++) { @@ -215,5 +224,5 @@ psFree (source->modelFlux); source->modelFlux = NULL; - modelFit = psphotFitEXT (readout, source, modelType, maskVal, markVal); + modelFit = psphotFitEXT (readout, source, fitOptions, modelType, maskVal, markVal); if (!modelFit) { psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My); @@ -311,4 +320,5 @@ } } + psFree (fitOptions); psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next); Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotFindDetections.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotFindDetections.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotFindDetections.c (revision 28643) @@ -90,10 +90,6 @@ psphotVisualShowSignificance (significance, -1.0, PS_SQR(3.0*NSIGMA_PEAK)); - // XXX getting some strange results from significance image - if (0) { - psImage *lsig = (psImage *) psUnaryOp (NULL, significance, "log"); - psphotVisualShowSignificance (lsig, 0.0, 4.0); - psFree (lsig); - } + // display the log significance image + psphotVisualShowLogSignificance (significance, 0.0, 4.5); // detect the peaks in the significance image Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotFitSet.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotFitSet.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotFitSet.c (revision 28643) @@ -24,6 +24,15 @@ } + // Define source fitting parameters for extended source fits + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->mode = PM_SOURCE_FIT_EXT; + // XXX for now, use the defaults for the rest: + // fitOptions->nIter = fitIter; + // fitOptions->tol = fitTol; + // fitOptions->poissonErrors = poisson; + // fitOptions->weight = PS_SQR(skySig); + // XXX pmSourceFitSet must cache the modelFlux? - pmSourceFitSet (source, modelSet, mode, maskVal); + pmSourceFitSet (source, modelSet, fitOptions, maskVal); // write out positive object Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotReadout.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotReadout.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotReadout.c (revision 28643) @@ -53,6 +53,7 @@ } - // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are determined and saved on - // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT + // load the psf model, if suppled. FWHM_MAJ,FWHM_MIN,etc are determined and saved on + // readout->analysis. XXX Note: this function currently only works with a single + // PSPHOT.INPUT if (!psphotLoadPSF (config, view)) { // ??? need to supply 2 ? psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotSetThreads.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotSetThreads.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotSetThreads.c (revision 28643) @@ -30,5 +30,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 9); + task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 10); task->function = &psphotBlendFit_Threaded; psThreadTaskAdd(task); Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c (revision 28643) @@ -2,5 +2,4 @@ // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful -// XXX this function does not call pmSourceFitModelInit : fix this? static int NfitPSF = 0; @@ -22,7 +21,9 @@ } -bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) { +bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { float x, y, dR; + + pmSourceFitOptions options = *fitOptions; // maskVal is used to test for rejected pixels, and must include markVal @@ -31,5 +32,5 @@ // if this source is not a possible blend, just fit as PSF if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) { - bool status = psphotFitPSF (readout, source, psf, maskVal, markVal); + bool status = psphotFitPSF (readout, source, psf, fitOptions, maskVal, markVal); return status; } @@ -91,5 +92,6 @@ // fit PSF model - pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal); + options.mode = PM_SOURCE_FIT_PSF; + pmSourceFitSet (source, modelSet, &options, maskVal); // clear the circular mask @@ -154,7 +156,8 @@ } -bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) { +bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { double chiTrend; + pmSourceFitOptions options = *fitOptions; // maskVal is used to test for rejected pixels, and must include markVal @@ -171,5 +174,6 @@ // fit PSF model (set/unset the pixel mask) - pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF, maskVal); + options.mode = PM_SOURCE_FIT_PSF; + pmSourceFitModel (source, PSF, &options, maskVal); if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); @@ -213,5 +217,5 @@ } -bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) { +bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { bool okEXT, okDBL; @@ -246,5 +250,5 @@ // XXX need to handle failures better here - EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal, markVal); + EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal); if (!EXT) goto escape; if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape; @@ -254,5 +258,5 @@ // DBL will always be defined, but DBL->data[n] might not - DBL = psphotFitDBL (readout, source, maskVal, markVal); + DBL = psphotFitDBL (readout, source, fitOptions, maskVal, markVal); if (!DBL) goto escape; if (!DBL->n) goto escape; @@ -379,5 +383,5 @@ // fit a double PSF source to an extended blob -psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) { +psArray *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { float dx, dy; @@ -387,4 +391,6 @@ psEllipseMoments moments; psArray *modelSet; + + pmSourceFitOptions options = *fitOptions; NfitDBL ++; @@ -426,9 +432,12 @@ // fit PSF model (set/unset the pixel mask) - pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal); + options.mode = PM_SOURCE_FIT_PSF; + pmSourceFitSet (source, modelSet, &options, maskVal); return (modelSet); } -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) { +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) { + + pmSourceFitOptions options = *fitOptions; NfitEXT ++; @@ -436,4 +445,6 @@ // maskVal is used to test for rejected pixels, and must include markVal maskVal |= markVal; + + psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); // use the source moments, etc to guess basic model parameters @@ -448,6 +459,28 @@ } + // for sersic models, get the initial guess more carefully + // if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { + // psphotGuessSersic (); + // // test and return NULL on failure? + // } + // fit EXT (not PSF) model (set/unset the pixel mask) - pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT, maskVal); + options.mode = PM_SOURCE_FIT_EXT; + pmSourceFitModel (source, EXT, &options, maskVal); + + psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); return (EXT); } + +// A sersic model is very sensitive to the index. attempt to find the index first by grid search in just the index +// bool psphotFitSersic (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) { +// +// assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); +// +// +// +// if (!model->modelGuess(model, source)) { +// } +// +// +// } Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c (revision 28643) @@ -1,4 +1,6 @@ # include "psphotInternal.h" # include + +# define KRON 1 typedef struct { @@ -197,5 +199,10 @@ pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); + // XXX test: switch to kron flux +# if (KRON) + float apMag = -2.5*log10(source->moments->KronFlux); +# else float apMag = -2.5*log10(source->moments->Sum); +# endif float dMag = source->psfMag - apMag; @@ -343,12 +350,21 @@ pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); +# if (KRON) + float apMag = -2.5*log10(source->moments->KronFlux); +# else float apMag = -2.5*log10(source->moments->Sum); +# endif float dMag = source->psfMag - apMag; // set nSigma to include both systematic and poisson error terms // XXX the 'poisson error' contribution for size is probably wrong... - float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr); + // XXX add in a hard floor on the Ap Sys Err (to be a bit generous) + float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.025)); float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag); float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag); + + fprintf (stderr, "Mxx: %f, Myy: %f, dx: %f, dy: %f, psfMag: %f, apMag: %f, dMag: %f, errMag: %f, nSigmaMag: %f\n", + Mxx, Myy, source->peak->xf - source->moments->Mx, source->peak->yf - source->moments->My, + source->psfMag, apMag, dMag, source->errMag, nSigmaMAG); // partially-masked sources are more likely to be mis-measured PSFs Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c (revision 28643) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +void pmSourceMomentsSetVerbose(bool state); // convert detections to sources and measure their basic properties (moments, local sky, sky @@ -369,4 +370,6 @@ maskVal |= markVal; + // XXX test : pmSourceMomentsSetVerbose(true); + // threaded measurement of the sources moments int Nfail = 0; @@ -408,4 +411,17 @@ // measure basic source moments (no S/N clipping on input pixels) status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal); + // XXX moments / aperture test: + if (0) { + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, 4.0*source->moments->Mrf, "OR", markVal); + + float apMag = NAN; + pmSourcePhotometryAper (&apMag, NULL, source->pixels, source->maskObj, maskVal); + fprintf (stderr, "apMag: %f, kronMag: %f\n", apMag, -2.5*log10(source->moments->KronFlux)); + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + } if (status) { Nmoments ++; Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotVisual.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotVisual.c (revision 28642) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotVisual.c (revision 28643) @@ -228,4 +228,19 @@ psphotVisualRangeImage (kapa, image, "signif", 2, min, max); + + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualShowLogSignificance (psImage *image, float min, float max) { + + if (!pmVisualTestLevel("psphot.image.logsignif", 3)) return true; + + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + + psImage *lsig = (psImage *) psUnaryOp (NULL, image, "log"); + psphotVisualRangeImage (kapa, lsig, "log-signif", 2, min, max); + psFree (lsig); pmVisualAskUser(NULL);