Index: trunk/psModules/src/astrom/pmAstrometryUtils.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryUtils.c (revision 35726) +++ trunk/psModules/src/astrom/pmAstrometryUtils.c (revision 35768) @@ -78,7 +78,10 @@ Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo); + // Beta = (-L,-M) for the current guess Beta->data.F32[0] = -1.0 * psPolynomial2DEval (trans->x, Xo, Yo); Beta->data.F32[1] = -1.0 * psPolynomial2DEval (trans->y, Xo, Yo); - + // fprintf (stderr, " Beta: %f %f\n", Beta->data.F32[0], Beta->data.F32[1]); + + // since we want (L,M) = (0,0), Beta is also the offset if (!psMatrixGJSolve (Alpha, Beta)) { psError(PS_ERR_UNKNOWN, false, "Unable to solve for center."); Index: trunk/psModules/src/astrom/pmAstrometryVisual.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryVisual.c (revision 35726) +++ trunk/psModules/src/astrom/pmAstrometryVisual.c (revision 35768) @@ -51,6 +51,6 @@ bool pmAstromVisualClose(void) { - if (kapa1 != -1) KiiClose(kapa1); - if (kapa2 != -1) KiiClose(kapa2); + if (kapa1 != -1) KapaClose(kapa1); + if (kapa2 != -1) KapaClose(kapa2); return true; } Index: trunk/psModules/src/camera/pmReadoutFake.c =================================================================== --- trunk/psModules/src/camera/pmReadoutFake.c (revision 35726) +++ trunk/psModules/src/camera/pmReadoutFake.c (revision 35768) @@ -35,5 +35,7 @@ #include "pmReadoutFake.h" -#define MAX_AXIS_RATIO 20.0 // Maximum axis ratio for PSF model +// XXX this is now hard-wired in pmModelParamsToAxes +// #define MAX_AXIS_RATIO 20.0 // Maximum axis ratio for PSF model + #define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \ PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models @@ -52,5 +54,5 @@ psF32 *params = model->params->data.F32; // Model parameters - psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO, model->type); // Ellipse axes + psEllipseAxes axes = pmPSF_ModelToAxes(params, model->type); // Ellipse axes // Curiously, the minor axis can be larger than the major axis, so need to check. if (axes.major >= axes.minor) { Index: trunk/psModules/src/extras/Makefile.am =================================================================== --- trunk/psModules/src/extras/Makefile.am (revision 35726) +++ trunk/psModules/src/extras/Makefile.am (revision 35768) @@ -10,4 +10,5 @@ pmVisual.c \ pmVisualUtils.c \ + pmThreadTools.c \ ippStages.c \ pmCensor.c @@ -20,4 +21,5 @@ pmVisual.h \ pmVisualUtils.h \ + pmThreadTools.h \ ippDiffMode.h \ ippStages.h \ Index: trunk/psModules/src/extras/pmThreadTools.c =================================================================== --- trunk/psModules/src/extras/pmThreadTools.c (revision 35768) +++ trunk/psModules/src/extras/pmThreadTools.c (revision 35768) @@ -0,0 +1,152 @@ +/** These utility functions help threading source analysis on readouts + * @author Eugene Magnier, IfA + * @date June 05, 2013 + */ + +/* Include Files */ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#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 "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSourcePhotometry.h" +#include "pmSourceSatstar.h" +#include "pmSource.h" + +// the strategy here is to divide the image into 2x2 blocks of cells and cycle through +// the four discontiguous sets of cells, threading all within a set and blocking between +// sets + +// we divide the image region into 2*2 blocks of size Nx*Ny, the image will have +// Cx*Cy blocks so that (2Nx)Cx = numCols, (2Ny)Cy = numRows. We want to choose Cx and +// Cy so that (2Nx)Cx * (2Ny)Cy = 4 * NFILL * nThreads -- each of the four sets of cells +// has enough cells to allow NFILL cells for each thread (to better distribute heavy and +// light load cells + +// the array runs from readout->image->col0 to readout->image->col0 + readout->image->numCols + +// we save these for threaded analysis runs +static int Xo = 0; +static int Yo = 0; +static int Nx = 1; +static int Ny = 1; + +bool pmReadoutChooseCellSizes (int *Cx, int *Cy, pmReadout *readout, int nThreads) { + + int nCells = nThreads * 2*2; // number of cells in a single set + int C = sqrt(nCells) + 0.5; + + // we need to assign Cx and Cy based on the dimensionality of the image + // crude way to find most evenly balanced factors of nCells: + for (int i = C; i >= 1; i--) { + int C1 = nCells / C; + int C2 = nCells / C1; + if (C1*C2 != nCells) continue; + + if (readout->image->numRows > readout->image->numCols) { + *Cx = PS_MAX (C1, C2); + *Cy = PS_MIN (C1, C2); + } else { + *Cx = PS_MAX (C1, C2); + *Cy = PS_MIN (C1, C2); + } + + Xo = readout->image->col0; + Yo = readout->image->row0; + Nx = readout->image->numCols / (*Cx*2); + Ny = readout->image->numRows / (*Cy*2); + + return true; + } + *Cx = 1; + *Cy = 1; + + Xo = readout->image->col0; + Yo = readout->image->row0; + Nx = readout->image->numCols / (*Cx*2); + Ny = readout->image->numRows / (*Cy*2); + + return true; +} + +bool pmReadoutCoordToCell (int *group, int *cell, float x, float y, int Cx, int Cy) { + + // XXX need to handle edges + int ix = (x - Xo)/(2*Nx); + ix = PS_MAX (0, PS_MIN (ix, Cx - 1)); + + int iy = (y - Yo)/(2*Ny); + iy = PS_MAX (0, PS_MIN (iy, Cy - 1)); + + int jx = (((int)(x - Xo))%(2*Nx))/Nx; + jx = PS_MAX (0, PS_MIN (jx, Nx - 1)); + + int jy = (((int)(y - Yo))%(2*Ny))/Ny; + jy = PS_MAX (0, PS_MIN (jy, Ny - 1)); + + *group = jx + 2*jy; + *cell = ix + Cx*iy; + + return true; +} + +// we have 2x2 * Cx*Cy cells for the image. we need a function to convert an x,y +// coordinate pair into the index for these cells + +// first, how shall we number them? is there one index for all cells, or one set of +// indices for each of the 2x2 cell groups? + +psArray *pmReadoutAssignSourcesToCells (int Cx, int Cy, psArray *sources) { + + psArray *cellGroups = psArrayAlloc (4); + for (int i = 0; i < cellGroups->n; i++) { + psArray *cells = psArrayAlloc (Cx*Cy); + cellGroups->data[i] = cells; + for (int j = 0; j < cells->n; j++) { + psArray *cellSources = psArrayAllocEmpty (50); + cells->data[j] = cellSources; + } + } + + for (int i = 0; i < sources->n; i++) { + + int group = 0; + int cell = 0; + + pmSource *source = sources->data[i]; + + pmReadoutCoordToCell (&group, &cell, source->peak->xf, source->peak->yf, Cx, Cy); + + psArray *cells = cellGroups->data[group]; + psArray *cellSources = cells->data[cell]; + + psArrayAdd (cellSources, 100, source); + } + + return cellGroups; +} Index: trunk/psModules/src/extras/pmThreadTools.h =================================================================== --- trunk/psModules/src/extras/pmThreadTools.h (revision 35768) +++ trunk/psModules/src/extras/pmThreadTools.h (revision 35768) @@ -0,0 +1,14 @@ +/* @file pmVisual.h + * @brief functions to create visual diagnostics with the help of 'kapa' + * @author Chris Beaumont, IfA + * + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_THREAD_TOOLS_H +#define PM_THREAD_TOOLS_H + +psArray *pmReadoutAssignSourcesToCells (int Cx, int Cy, psArray *sources); +bool pmReadoutChooseCellSizes (int *Cx, int *Cy, pmReadout *readout, int nThreads); + +#endif //ndef PM_THREAD_TOOLS_H Index: trunk/psModules/src/imcombine/pmStackVisual.c =================================================================== --- trunk/psModules/src/imcombine/pmStackVisual.c (revision 35726) +++ trunk/psModules/src/imcombine/pmStackVisual.c (revision 35768) @@ -36,7 +36,7 @@ { if(kapa != -1) - KiiClose(kapa); + KapaClose(kapa); if(kapa2 != -1) - KiiClose(kapa2); + KapaClose(kapa2); return true; } Index: trunk/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionVisual.c (revision 35726) +++ trunk/psModules/src/imcombine/pmSubtractionVisual.c (revision 35768) @@ -56,6 +56,6 @@ bool pmSubtractionVisualClose(void) { - if(kapa1 != -1) KiiClose(kapa1); - if(kapa2 != -1) KiiClose(kapa2); + if(kapa1 != -1) KapaClose(kapa1); + if(kapa2 != -1) KapaClose(kapa2); return true; } Index: trunk/psModules/src/objects/models/pmModel_DEV.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_DEV.c (revision 35726) +++ trunk/psModules/src/objects/models/pmModel_DEV.c (revision 35768) @@ -123,12 +123,6 @@ // first, use Rmajor and index to find the central pixel flux (fraction of total flux) - psEllipseShape shape; - - shape.sx = PAR[PM_PAR_SXX]; - shape.sy = PAR[PM_PAR_SYY]; - shape.sxy = PAR[PM_PAR_SXY]; - - // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); // get the central pixel flux from the lookup table @@ -238,5 +232,6 @@ // 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); + // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy + q2 = 2.0*0.5*sqrtf(q1); } @@ -303,6 +298,5 @@ // set the shape parameters - // XXX adjust this? - if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) { return false; } @@ -330,14 +324,8 @@ psF64 PM_MODEL_FLUX (const psVector *params) { - psEllipseShape shape; - psF32 *PAR = params->data.F32; - shape.sx = PAR[PM_PAR_SXX]; - shape.sy = PAR[PM_PAR_SYY]; - shape.sxy = PAR[PM_PAR_SXY]; - - // for a non-circular DeVaucouleur, the flux of the Rmajor equivalent is scaled by the AspectRatio - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); float AspectRatio = axes.minor / axes.major; @@ -359,6 +347,4 @@ psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) { - psEllipseShape shape; - psF32 *PAR = params->data.F32; @@ -370,9 +356,6 @@ return (1.0); - shape.sx = PAR[PM_PAR_SXX]; - shape.sy = PAR[PM_PAR_SYY]; - shape.sxy = PAR[PM_PAR_SXY]; - - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); // f = Io exp(-z^n) -> z^n = ln(Io/f) @@ -382,5 +365,4 @@ psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]); - return (radius); } @@ -407,5 +389,6 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - if (!pmPSF_FitToModel (out, 0.1)) { + bool useReff = pmModelUseReff (modelPSF->type); + if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); return false; @@ -460,5 +443,6 @@ // convert to shape terms (SXX,SYY,SXY) // XXX user-defined value for limit? - if (!pmPSF_FitToModel (PAR, 0.1)) { + bool useReff = pmModelUseReff (model->type); + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); return false; Index: trunk/psModules/src/objects/models/pmModel_EXP.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_EXP.c (revision 35726) +++ trunk/psModules/src/objects/models/pmModel_EXP.c (revision 35768) @@ -115,12 +115,6 @@ // first, use Rmajor and index to find the central pixel flux (fraction of total flux) - psEllipseShape shape; - - shape.sx = PAR[PM_PAR_SXX]; - shape.sy = PAR[PM_PAR_SYY]; - shape.sxy = PAR[PM_PAR_SXY]; - - // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); // get the central pixel flux from the lookup table @@ -230,5 +224,6 @@ // 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); + // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy + q2 = 2.0*0.5*sqrtf(q1); } @@ -295,6 +290,5 @@ // set the shape parameters - // XXX adjust this? - if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) { return false; } @@ -316,14 +310,8 @@ psF64 PM_MODEL_FLUX (const psVector *params) { - psEllipseShape shape; - psF32 *PAR = params->data.F32; - shape.sx = PAR[PM_PAR_SXX]; - shape.sy = PAR[PM_PAR_SYY]; - shape.sxy = PAR[PM_PAR_SXY]; - - // for a non-circular Exponential, the flux of the Rmajor equivalent is scaled by the AspectRatio - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); float AspectRatio = axes.minor / axes.major; @@ -345,6 +333,4 @@ psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) { - psEllipseShape shape; - psF32 *PAR = params->data.F32; @@ -356,9 +342,6 @@ return (1.0); - shape.sx = PAR[PM_PAR_SXX]; - shape.sy = PAR[PM_PAR_SYY]; - shape.sxy = PAR[PM_PAR_SXY]; - - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f) @@ -392,5 +375,6 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - if (!pmPSF_FitToModel (out, 0.1)) { + bool useReff = pmModelUseReff (modelPSF->type); + if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); return false; @@ -445,5 +429,6 @@ // convert to shape terms (SXX,SYY,SXY) // XXX user-defined value for limit? - if (!pmPSF_FitToModel (PAR, 0.1)) { + bool useReff = pmModelUseReff (model->type); + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); return false; Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 35726) +++ trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 35768) @@ -129,6 +129,7 @@ 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]); + // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy + float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]); + float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]); float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); q1 = (q1 < 0.0) ? 0.0 : q1; @@ -200,5 +201,5 @@ // set the shape parameters - if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) { return false; } @@ -219,16 +220,10 @@ psF64 PM_MODEL_FLUX (const psVector *params) { - - 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]; + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); // Area is equivalent to 2 pi sigma^2 - // axes ratio < 20 - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); psF64 Area = 2.0 * M_PI * axes.major * axes.minor; @@ -242,6 +237,4 @@ psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) { - psEllipseShape shape; - psF32 *PAR = params->data.F32; @@ -253,9 +246,7 @@ 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); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); + psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux)); psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", @@ -285,12 +276,8 @@ } - // the OLD 2D model for SXY actually fitted SXY / (SXX^-2 + SYY^-2); correct here - // out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); - // 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 (out, 0.1)) { - // psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); + bool useReff = pmModelUseReff (modelPSF->type); + if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); return false; @@ -343,6 +330,6 @@ // 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)) { + bool useReff = pmModelUseReff (model->type); + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); return false; Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 35726) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 35768) @@ -129,6 +129,7 @@ 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]); + // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy + float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]); + float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]); float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); q1 = (q1 < 0.0) ? 0.0 : q1; @@ -201,5 +202,5 @@ // set the shape parameters - if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) { return false; } @@ -222,13 +223,10 @@ { 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]; - - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); + float AspectRatio = axes.minor / axes.major; @@ -262,5 +260,4 @@ { psF64 z; - psEllipseShape shape; psF32 *PAR = params->data.F32; @@ -273,10 +270,7 @@ 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]; - - // this estimates the radius assuming f(z) is roughly exp(-z) - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); + psF64 sigma = axes.major; @@ -347,10 +341,8 @@ } - // the OLD 2D model for SXY actually fitted SXY / (SXX^-2 + SYY^-2); correct here - // out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); - // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - if (!pmPSF_FitToModel (out, 0.1)) { + bool useReff = pmModelUseReff (modelPSF->type); + if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); return false; @@ -403,6 +395,6 @@ // 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)) { + bool useReff = pmModelUseReff (model->type); + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); return false; Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 35726) +++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 35768) @@ -1,9 +1,10 @@ /****************************************************************************** - * this file defines the PS1_V1 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: + * this file defines the PS1_V1 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: power-law with fitted linear term @@ -148,6 +149,6 @@ 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 f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]); + float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]); float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); q1 = (q1 < 0.0) ? 0.0 : q1; @@ -220,5 +221,5 @@ // set the shape parameters - if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) { return false; } @@ -244,13 +245,9 @@ { 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]; - - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); float AspectRatio = axes.minor / axes.major; @@ -284,5 +281,4 @@ { psF64 z; - psEllipseShape shape; psF32 *PAR = params->data.F32; @@ -292,9 +288,6 @@ 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); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); psF64 sigma = axes.major; @@ -363,5 +356,6 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - if (!pmPSF_FitToModel (out, 0.1)) { + bool useReff = pmModelUseReff (modelPSF->type); + if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); return false; @@ -416,5 +410,6 @@ // convert to shape terms (SXX,SYY,SXY) // XXX user-defined value for limit? - if (!pmPSF_FitToModel (PAR, 0.1)) { + bool useReff = pmModelUseReff (model->type); + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); return false; Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 35726) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 35768) @@ -1,9 +1,10 @@ /****************************************************************************** - * this file defines the QGAUSS source shape model (XXX need a better name!). 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 a PSF - * representations share a few parameters, for which # define names are listed in pmModel.h: + * this file defines the QGAUSS 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: power-law with fitted linear term @@ -14,6 +15,6 @@ * 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_SXX 4 - X^2 term of elliptical contour (SigmaX / sqrt(2)) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (SigmaY / sqrt(2)) * PM_PAR_SXY 6 - X*Y term of elliptical contour * PM_PAR_7 7 - amplitude of the linear component (k) @@ -138,5 +139,4 @@ # define AR_MAX 20.0 # define AR_RATIO 0.99 - bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) { @@ -149,6 +149,7 @@ 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]); + // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy + float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]); + float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]); float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); q1 = (q1 < 0.0) ? 0.0 : q1; @@ -203,5 +204,5 @@ return true; } - default: + default: psAbort("invalid choice for limits"); } @@ -221,5 +222,5 @@ // set the shape parameters - if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) { return false; } @@ -245,13 +246,9 @@ { 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]; - - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); float AspectRatio = axes.minor / axes.major; @@ -285,5 +282,4 @@ { psF64 z; - psEllipseShape shape; psF32 *PAR = params->data.F32; @@ -293,11 +289,6 @@ if (flux >= PAR[PM_PAR_I0]) return 1.0; - // if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA); - - 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); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); psF64 sigma = axes.major; @@ -307,5 +298,4 @@ return ( sigma * sqrt (2.0 * z) ); } - psF64 limit = flux / PAR[PM_PAR_I0]; @@ -367,5 +357,6 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - if (!pmPSF_FitToModel (out, 0.1)) { + bool useReff = pmModelUseReff (modelPSF->type); + if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); return false; @@ -424,6 +415,6 @@ // 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)) { + bool useReff = pmModelUseReff (model->type); + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); return false; Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 35726) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 35768) @@ -138,6 +138,7 @@ 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]); + // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy + float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]); + float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]); float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); q1 = (q1 < 0.0) ? 0.0 : q1; @@ -210,5 +211,5 @@ // set the shape parameters - if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) { return false; } @@ -234,13 +235,9 @@ { 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]; - - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); float AspectRatio = axes.minor / axes.major; @@ -274,5 +271,4 @@ { psF64 z; - psEllipseShape shape; psF32 *PAR = params->data.F32; @@ -285,9 +281,6 @@ 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); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); psF64 sigma = axes.major; @@ -357,5 +350,6 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - if (!pmPSF_FitToModel (out, 0.1)) { + bool useReff = pmModelUseReff (modelPSF->type); + if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); return false; @@ -409,6 +403,6 @@ // 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)) { + bool useReff = pmModelUseReff (model->type); + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); return false; Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 35726) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 35768) @@ -125,12 +125,6 @@ // first, use Rmajor and index to find the central pixel flux (fraction of total flux) - psEllipseShape shape; - - shape.sx = PAR[PM_PAR_SXX]; - shape.sy = PAR[PM_PAR_SYY]; - shape.sxy = PAR[PM_PAR_SXY]; - - // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); // get the central pixel flux from the lookup table @@ -250,5 +244,6 @@ // 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); + // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy + q2 = 2.0*0.5*sqrtf(q1); } @@ -347,9 +342,6 @@ axes.major = Rmajor; axes.minor = Rminor; - psEllipseShape shape = psEllipseAxesToShape (axes); - - if (!isfinite(shape.sx)) return false; - if (!isfinite(shape.sy)) return false; - if (!isfinite(shape.sxy)) return false; + + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], axes, true); float bn = 1.9992*index - 0.3271; @@ -357,11 +349,4 @@ 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_SXX] = Sxx; - PAR[PM_PAR_SYY] = Syy; - PAR[PM_PAR_SXY] = shape.sxy; - // set the model normalization (adjust for Sersic best guess) if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { @@ -381,14 +366,8 @@ psF64 PM_MODEL_FLUX (const psVector *params) { - psEllipseShape shape; - psF32 *PAR = params->data.F32; - shape.sx = PAR[PM_PAR_SXX]; - shape.sy = PAR[PM_PAR_SYY]; - shape.sxy = PAR[PM_PAR_SXY]; - - // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); float AspectRatio = axes.minor / axes.major; @@ -410,6 +389,4 @@ psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) { - psEllipseShape shape; - psF32 *PAR = params->data.F32; @@ -421,9 +398,6 @@ return (1.0); - shape.sx = PAR[PM_PAR_SXX]; - shape.sy = PAR[PM_PAR_SYY]; - shape.sxy = PAR[PM_PAR_SXY]; - - psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); // f = Io exp(-z^n) -> z^n = ln(Io/f) @@ -457,5 +431,6 @@ // the 2D PSF model fits polarization terms (E0,E1,E2) // convert to shape terms (SXX,SYY,SXY) - if (!pmPSF_FitToModel (out, 0.1)) { + bool useReff = pmModelUseReff (modelPSF->type); + if (!pmPSF_FitToModel (out, 0.1, useReff)) { psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); return false; @@ -510,5 +485,6 @@ // convert to shape terms (SXX,SYY,SXY) // XXX user-defined value for limit? - if (!pmPSF_FitToModel (PAR, 0.1)) { + bool useReff = pmModelUseReff (model->type); + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); return false; Index: trunk/psModules/src/objects/models/pmModel_TRAIL.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_TRAIL.c (revision 35726) +++ trunk/psModules/src/objects/models/pmModel_TRAIL.c (revision 35768) @@ -350,10 +350,9 @@ PAR[PM_PAR_SKY] = 0.0; - // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: - psEllipseShape psfShape; - psfShape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; - psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; - psfShape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; - psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0); + psF32 *psfPAR = source->modelPSF->params->data.F32; + bool useReff = pmModelUseReff (source->modelPSF->type); + + psEllipseAxes psfAxes; + pmModelParamsToAxes (&psfAxes, psfPAR[PM_PAR_SXX], psfPAR[PM_PAR_SXY], psfPAR[PM_PAR_SYY], useReff); psEllipseMoments emoments; @@ -369,7 +368,10 @@ if (!isfinite(axes.theta)) return false; - float size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf); - // if (size > psfAxes.major) { size -= psfAxes.major; } - //else { size = psfAxes.major; } + float size = NAN; + if (!isfinite(source->moments->Mrf)) { + size = axes.major; + } else { + size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf); + } float theta, peak; Index: trunk/psModules/src/objects/pmModelUtils.c =================================================================== --- trunk/psModules/src/objects/pmModelUtils.c (revision 35726) +++ trunk/psModules/src/objects/pmModelUtils.c (revision 35768) @@ -118,5 +118,66 @@ } -bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments) { +bool pmModelUseReff (pmModelType type) { + bool useReff = false; + useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC")); + useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV")); + useReff |= (type == pmModelClassGetType ("PS_MODEL_EXP")); + return useReff; +} + +// this function and the one below handle the two cases, where the model shape is uses R_eff or Sigma +bool pmModelAxesToParams (float *Sxx, float *Sxy, float *Syy, psEllipseAxes axes, bool useReff) { + + psEllipseShape shape = psEllipseAxesToShape (axes); + + if (!isfinite(shape.sx)) return false; + if (!isfinite(shape.sy)) return false; + if (!isfinite(shape.sxy)) return false; + + // set the shape parameters + if (useReff) { + *Sxx = PS_MAX(0.5, shape.sx); + *Syy = PS_MAX(0.5, shape.sy); + *Sxy = shape.sxy * 2.0; + } else { + *Sxx = PS_MAX(0.5, M_SQRT2*shape.sx); + *Syy = PS_MAX(0.5, M_SQRT2*shape.sy); + *Sxy = shape.sxy; + } + + return true; +} + +bool pmModelParamsToAxes (psEllipseAxes *axes, float Sxx, float Sxy, float Syy, bool useReff) { + + psEllipseShape shape; + + // set the shape parameters + if (useReff) { + shape.sx = Sxx; + shape.sy = Syy; + shape.sxy = Sxy / 2.0; + } else { + shape.sx = Sxx / M_SQRT2; + shape.sy = Syy / M_SQRT2; + shape.sxy = Sxy; + } + + if ((shape.sx == 0) || (shape.sy == 0)) { + axes->major = 0.0; + axes->minor = 0.0; + axes->theta = 0.0; + } else { + // axes ratio < 20 + // replace with maxAR argument? + *axes = psEllipseShapeToAxes (shape, 20.0); + } + + return true; +} + +// Reff says if this is a model which uses R_eff (like exp or dev) instead of Sigma +// set the parameter values SXX, SXY, SYY +bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments, bool useReff) { psEllipseMoments emoments; @@ -137,14 +198,5 @@ axes.minor *= scale; - psEllipseShape shape = psEllipseAxesToShape (axes); - - if (!isfinite(shape.sx)) return false; - if (!isfinite(shape.sy)) return false; - if (!isfinite(shape.sxy)) return false; - - // set the shape parameters - *Sxx = PS_MAX(0.5, M_SQRT2*shape.sx); - *Syy = PS_MAX(0.5, M_SQRT2*shape.sy); - *Sxy = shape.sxy; + pmModelAxesToParams (Sxx, Sxy, Syy, axes, useReff); return true; Index: trunk/psModules/src/objects/pmModelUtils.h =================================================================== --- trunk/psModules/src/objects/pmModelUtils.h (revision 35726) +++ trunk/psModules/src/objects/pmModelUtils.h (revision 35768) @@ -44,5 +44,9 @@ bool pmModelSetPosition (float *Xo, float *Yo, pmSource *source); bool pmModelSetNorm (float *Io, pmSource *source); -bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments); +bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments, bool useReff); + +bool pmModelUseReff (pmModelType type); +bool pmModelAxesToParams (float *Sxx, float *Sxy, float *Syy, psEllipseAxes axes, bool useReff); +bool pmModelParamsToAxes (psEllipseAxes *axes, float Sxx, float Sxy, float Syy, bool useReff); // XXX void pmModelSetModelVarOption (bool option); Index: trunk/psModules/src/objects/pmPCM_MinimizeChisq.c =================================================================== --- trunk/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 35726) +++ trunk/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 35768) @@ -81,8 +81,9 @@ psAbort ("programming error: no unmasked parameters to be fit\n"); } + psAssert (pcm->nPar == Beta->n, "did we set the masked parameters correctly??"); // allocate internal arrays (current vs Guess) - psImage *alpha = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32); - psVector *beta = psVectorAlloc(Beta->n, PS_TYPE_F32); + psImage *alpha = psImageAlloc(pcm->nPar, pcm->nPar, PS_TYPE_F32); + psVector *beta = psVectorAlloc(pcm->nPar, PS_TYPE_F32); psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32); @@ -90,4 +91,5 @@ psF32 lambda = 0.001; psF32 dLinear = 0.0; + psF32 nu = 2.0; # if (USE_FFT && PRE_CONVOLVE) @@ -120,6 +122,11 @@ bool done = (min->iter >= min->maxIter); while (!done) { - psTrace("psModules.objects", 5, "Iteration number %d. (max iterations is %d).\n", min->iter, min->maxIter); - psTrace("psModules.objects", 5, "Last delta is %f. stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol); + psTrace(FACILITY, 5, "Iteration number %d. (max iterations is %d).\n", min->iter, min->maxIter); + + if (min->chisqConvergence) { + psTrace(FACILITY, 5, "Last delta is %f. stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol); + } else { + psTrace(FACILITY, 5, "Last delta is %f. stop if < %f, accept if < %f\n", min->rParSigma, min->minTol*pcm->nPar, min->maxTol*pcm->nPar); + } // set a new guess for Alpha, Beta, Params @@ -140,4 +147,27 @@ p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)"); } + + // calculate the parameter change (rParDelta) and error radius (rParSigma) + // rParDelta : radius of parameter change; + // rParSigma : radius of parameter error + + // note that (before SetABX) Alpha[i][i] is the covariance matrix and + // Beta is the actual parameter change for this pass + + // note that Alpha & Beta only represent unmasked parameters, while params and Params have all + + // dParSigma = Alpha[i][i] : error (squared) on parameter i + // dParDelta = Params->data.F32[i] - params->data.F32[i] : change on parameter i + float rParSigma = 0.0; + for (int j = 0, J = 0; j < Params->n; j++) { + if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) { + continue; + } + rParSigma += PS_SQR(Params->data.F32[j] - params->data.F32[j]) / Alpha->data.F32[J][J]; + J++; + } + rParSigma = sqrt(rParSigma); + psTrace(FACILITY, 5, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter); + // fprintf (stderr, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter); // calculate Chisq for new guess, update Alpha & Beta @@ -164,27 +194,77 @@ } - /* if (Chisq < min->value) { */ + // change in chisq/nDOF since last minimum + min->lastDelta = (min->value - Chisq) / pcm->nDOF; + + // rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho) + + // XXX the old version of lambda changes: + // XXX : Madsen gives suggestion for better use of rho + // rho is positive if the new chisq is smaller if (rho >= -1e-6) { - min->lastDelta = (min->value - Chisq) / (source->pixels->numCols*source->pixels->numRows - params->n); min->value = Chisq; alpha = psImageCopy(alpha, Alpha, PS_TYPE_F32); beta = psVectorCopy(beta, Beta, PS_TYPE_F32); params = psVectorCopy(params, Params, PS_TYPE_F32); - lambda *= 0.25; // save the new convolved model image psFree (source->modelFlux); source->modelFlux = pmPCMdataSaveImage(pcm); - } else { - lambda *= 10.0; - } + } + switch (min->gainFactorMode) { + case 0: + if (rho >= -1e-6) { + lambda *= 0.25; + } else { + lambda *= 10.0; + } + break; + + case 1: + // adjust the gain ratio (lambda) based on rho + if (rho < 0.25) { + lambda *= 2.0; + } + if (rho > 0.75) { + lambda *= 0.333; + } + break; + + case 2: + if (rho > 0.0) { + lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0))); + nu = 2.0; + } else { + lambda *= nu; + nu *= 2.0; + } + break; + } min->iter++; + // ending conditions: + // 1) hard limit : too many iterations done = (min->iter >= min->maxIter); - // check for convergence: + // 2) require deltaChi > 1e-6 (ie, chisq is decreasing, but accept an insignificant change) + if (min->lastDelta < -1e-6) { + continue; + } + + // save this value in case we stop iterating + min->rParSigma = rParSigma; + + // 2) require chisqDOF < maxChisqDOF (if maxChisqDOF is not NAN) + // keep iterating regardless of rParSigma in this case float chisqDOF = Chisq / pcm->nDOF; - if (!isfinite(min->maxChisqDOF) || ((chisqDOF < min->maxChisqDOF) && isfinite(min->lastDelta))) { + if (isfinite(min->maxChisqDOF) && (chisqDOF > min->maxChisqDOF)) { + continue; + } + + // delta-chisq or rParSigma ? + if (min->chisqConvergence) { done |= (min->lastDelta < min->minTol); + } else { + done |= (rParSigma < min->minTol*pcm->nPar); } } @@ -220,7 +300,14 @@ // if the last improvement was at least as good as maxTol, accept the fit: - if (min->lastDelta <= min->maxTol) { - psTrace(FACILITY, 6, "---- end (true) ----\n"); - return(true); + if (min->chisqConvergence) { + if (min->lastDelta <= min->maxTol) { + psTrace(FACILITY, 6, "---- end (true) ----\n"); + return(true); + } + } else { + if (min->rParSigma <= min->maxTol*pcm->nPar) { + psTrace(FACILITY, 6, "---- end (true) ----\n"); + return(true); + } } psTrace(FACILITY, 6, "---- end (false) ----\n"); Index: trunk/psModules/src/objects/pmPCMdata.c =================================================================== --- trunk/psModules/src/objects/pmPCMdata.c (revision 35726) +++ trunk/psModules/src/objects/pmPCMdata.c (revision 35768) @@ -291,13 +291,10 @@ psAssert (modelPSF, "psf model must be defined"); - psEllipseShape shape; psEllipseAxes axes; - - shape.sx = modelPSF->params->data.F32[PM_PAR_SXX]; - shape.sy = modelPSF->params->data.F32[PM_PAR_SYY]; - shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; - axes = psEllipseShapeToAxes (shape, 20.0); + bool useReff = pmModelUseReff (modelPSF->type); + psF32 *PAR = modelPSF->params->data.F32; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); - float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[PM_PAR_I0]); + float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); @@ -451,13 +448,10 @@ psAssert (modelPSF, "psf model must be defined"); - psEllipseShape shape; psEllipseAxes axes; - - shape.sx = modelPSF->params->data.F32[PM_PAR_SXX]; - shape.sy = modelPSF->params->data.F32[PM_PAR_SYY]; - shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; - axes = psEllipseShapeToAxes (shape, 20.0); + bool useReff = pmModelUseReff (modelPSF->type); + psF32 *PAR = modelPSF->params->data.F32; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); - float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[PM_PAR_I0]); + float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); Index: trunk/psModules/src/objects/pmPSF.c =================================================================== --- trunk/psModules/src/objects/pmPSF.c (revision 35726) +++ trunk/psModules/src/objects/pmPSF.c (revision 35768) @@ -249,39 +249,43 @@ // Mxy = SXY * (SXX^-4 + SYY^-4 - 2 SXY ^2) +// XXX deprecated // input: model->param, output: psf->param[PM_PAR_SXY] -double pmPSF_SXYfromModel (psF32 *modelPar) -{ - PS_ASSERT_PTR_NON_NULL(modelPar, NAN); - - double SXX = modelPar[PM_PAR_SXX]; - double SYY = modelPar[PM_PAR_SYY]; - double SXY = modelPar[PM_PAR_SXY]; - - double par = SXY / PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY)); - return (par); -} - +// XXX double pmPSF_SXYfromModel (psF32 *modelPar) +// XXX { +// XXX PS_ASSERT_PTR_NON_NULL(modelPar, NAN); +// XXX +// XXX double SXX = modelPar[PM_PAR_SXX]; +// XXX double SYY = modelPar[PM_PAR_SYY]; +// XXX double SXY = modelPar[PM_PAR_SXY]; +// XXX +// XXX double par = SXY / PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY)); +// XXX return (par); +// XXX } + +// XXX deprecated // input: fitted psf->param, output: model->param[PM_PAR_SXY] -double pmPSF_SXYtoModel (psF32 *fittedPar) -{ - PS_ASSERT_PTR_NON_NULL(fittedPar, NAN); - - double SXX = fittedPar[PM_PAR_SXX]; - double SYY = fittedPar[PM_PAR_SYY]; - double fit = fittedPar[PM_PAR_SXY]; - - double SXY = fit * PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY)); - - assert (!isnan(SXY)); - - return SXY; -} - -// New Concept: the PSF modelling function fits the polarization terms e0, e1, e2: - -// convert the parameters used in the fitted source model -// to the parameters used in the 2D PSF model -// XXX this function may be invalid for SERSIC, DEV, EXP models (SQRT2 not used?) -bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis) +// XXX double pmPSF_SXYtoModel (psF32 *fittedPar) +// XXX { +// XXX PS_ASSERT_PTR_NON_NULL(fittedPar, NAN); +// XXX +// XXX double SXX = fittedPar[PM_PAR_SXX]; +// XXX double SYY = fittedPar[PM_PAR_SYY]; +// XXX double fit = fittedPar[PM_PAR_SXY]; +// XXX +// XXX double SXY = fit * PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY)); +// XXX +// XXX assert (!isnan(SXY)); +// XXX +// XXX return SXY; +// XXX } + +// The PSF modelling function fits the polarization terms e0, e1, e2: + +// the FIT is the 2D representation of the shape using polarization parameters for the elliptical contour +// the MODEL is the realized psf model for a given location + +// convert the parameters (in situ) used in the fitted source model to the parameters used in +// the 2D PSF model +bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis, bool useReff) { PS_ASSERT_PTR_NON_NULL(fittedPar, false); @@ -298,17 +302,12 @@ return false; } - psEllipseShape shape = psEllipseAxesToShape (axes); - - fittedPar[PM_PAR_SXX] = shape.sx * M_SQRT2; - fittedPar[PM_PAR_SYY] = shape.sy * M_SQRT2; - fittedPar[PM_PAR_SXY] = shape.sxy; - + + pmModelAxesToParams (&fittedPar[PM_PAR_SXX], &fittedPar[PM_PAR_SXY], &fittedPar[PM_PAR_SYY], axes, useReff); return true; } -// convert the PSF parameters used in the 2D PSF model fit into the -// parameters used in the source model -// XXX this function may be invalid for SERSIC, DEV, EXP models (SQRT2 not used?) -psEllipsePol pmPSF_ModelToFit (psF32 *modelPar) +// convert the parameters (in situ) used in the 2D PSF model fit into the parameters used in +// the source model +psEllipsePol pmPSF_ModelToFit (psF32 *modelPar, bool useReff) { // must assert non-NULL input parameter @@ -319,11 +318,8 @@ PS_ASSERT_PTR_NON_NULL(modelPar, pol); - psEllipseShape shape; - - shape.sx = modelPar[PM_PAR_SXX] / M_SQRT2; - shape.sy = modelPar[PM_PAR_SYY] / M_SQRT2; - shape.sxy = modelPar[PM_PAR_SXY]; - - pol = psEllipseShapeToPol (shape); + psEllipseAxes axes; + pmModelParamsToAxes (&axes, modelPar[PM_PAR_SXX], modelPar[PM_PAR_SXY], modelPar[PM_PAR_SYY], useReff); + + pol = psEllipseAxesToPol (axes); return pol; @@ -332,40 +328,15 @@ // convert the parameters used in the fitted source model to the psEllipseAxes representation // (major,minor,theta) -psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR, pmModelType type) -{ - psEllipseShape shape; +psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, pmModelType type) +{ psEllipseAxes axes; axes.major = NAN; axes.minor = NAN; axes.theta = NAN; - // XXX: must assert non-NULL input parameter + PS_ASSERT_PTR_NON_NULL(modelPar, axes); - bool useReff = false; - useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC")); - useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV")); - useReff |= (type == pmModelClassGetType ("PS_MODEL_EXP")); - - if (useReff) { - shape.sx = modelPar[PM_PAR_SXX]; - shape.sy = modelPar[PM_PAR_SYY]; - shape.sxy = modelPar[PM_PAR_SXY] / 2.0; - // XXX I *think* dividing by 2.0 is the right direction, but this - // needs to be checked with a real test - } else { - shape.sx = modelPar[PM_PAR_SXX] / M_SQRT2; - shape.sy = modelPar[PM_PAR_SYY] / M_SQRT2; - shape.sxy = modelPar[PM_PAR_SXY]; - } - - if ((shape.sx == 0) || (shape.sy == 0)) { - axes.major = 0.0; - axes.minor = 0.0; - axes.theta = 0.0; - } else { - // XXX this is not really consistent with the model fit range above - axes = psEllipseShapeToAxes (shape, maxAR); - } - + bool useReff = pmModelUseReff (type); + pmModelParamsToAxes (&axes, modelPar[PM_PAR_SXX], modelPar[PM_PAR_SXY], modelPar[PM_PAR_SYY], useReff); return axes; } @@ -377,27 +348,14 @@ PS_ASSERT_PTR_NON_NULL(modelPar, false); + modelPar[PM_PAR_SXX] = 0.0; + modelPar[PM_PAR_SYY] = 0.0; + modelPar[PM_PAR_SXY] = 0.0; + if ((axes.major <= 0) || (axes.minor <= 0)) { - modelPar[PM_PAR_SXX] = 0.0; - modelPar[PM_PAR_SYY] = 0.0; - modelPar[PM_PAR_SXY] = 0.0; return true; } - - psEllipseShape shape = psEllipseAxesToShape (axes); - - bool useReff = false; - useReff |= ( type == pmModelClassGetType ("PS_MODEL_SERSIC")); - useReff |= ( type == pmModelClassGetType ("PS_MODEL_DEV")); - useReff |= ( type == pmModelClassGetType ("PS_MODEL_EXP")); - - if (useReff) { - modelPar[PM_PAR_SXX] = shape.sx; - modelPar[PM_PAR_SYY] = shape.sy; - modelPar[PM_PAR_SXY] = shape.sxy * 2.0; // XXX NEED factor of 2 here for correct angle conversion - } else { - modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2; - modelPar[PM_PAR_SYY] = shape.sy * M_SQRT2; - modelPar[PM_PAR_SXY] = shape.sxy; - } + + bool useReff = pmModelUseReff (type); + pmModelAxesToParams (&modelPar[PM_PAR_SXX], &modelPar[PM_PAR_SXY], &modelPar[PM_PAR_SYY], axes, useReff); return true; } @@ -423,5 +381,6 @@ par->data.F32[PM_PAR_SXY] = sxy; - psEllipsePol pol = pmPSF_ModelToFit(par->data.F32); + bool useReff = pmModelUseReff (options->type); + psEllipsePol pol = pmPSF_ModelToFit(par->data.F32, useReff); pmTrend2D *trend = NULL; Index: trunk/psModules/src/objects/pmPSF.h =================================================================== --- trunk/psModules/src/objects/pmPSF.h (revision 35726) +++ trunk/psModules/src/objects/pmPSF.h (revision 35768) @@ -107,8 +107,8 @@ bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, pmModelType type); -bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis); +bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis, bool useReff); -psEllipsePol pmPSF_ModelToFit (psF32 *modelPar); -psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR, pmModelType type); +psEllipsePol pmPSF_ModelToFit (psF32 *modelPar, bool useReff); +psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, pmModelType type); /// Calculate FWHM value from a PSF Index: trunk/psModules/src/objects/pmPSFtry.h =================================================================== --- trunk/psModules/src/objects/pmPSFtry.h (revision 35726) +++ trunk/psModules/src/objects/pmPSFtry.h (revision 35768) @@ -99,8 +99,12 @@ /** fit EXT models to all possible psf sources */ bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal); +bool pmPSFtryFitEXT_Threaded (psThreadJob *job); bool pmPSFtryMakePSF (bool *pGoodFit, pmPSFtry *psfTry); bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal); +bool pmPSFtryFitPSF_Threaded (psThreadJob *job); + +bool pmPSFThreads (void); /** pmPSFtryMetric() Index: trunk/psModules/src/objects/pmPSFtryFitEXT.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryFitEXT.c (revision 35726) +++ trunk/psModules/src/objects/pmPSFtryFitEXT.c (revision 35768) @@ -46,9 +46,26 @@ #include "pmSourceVisual.h" +bool pmPSFThreads (void) { + + psThreadTask *task = NULL; + + task = psThreadTaskAlloc("PSF_TRY_FIT_EXT", 6); + task->function = &pmPSFtryFitEXT_Threaded; + psThreadTaskAdd(task); + psFree(task); + + task = psThreadTaskAlloc("PSF_TRY_FIT_PSF", 6); + task->function = &pmPSFtryFitPSF_Threaded; + psThreadTaskAdd(task); + psFree(task); + + return true; +} + +static int Next = 0; + // Fit an EXT model to all candidates PSF sources. // Note: this is independent of the modeled 2D variations in the PSF. bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) { - - bool status; psTimerStart ("psf.fit"); @@ -60,52 +77,74 @@ maskVal |= markVal; - int Next = 0; + Next = 0; for (int i = 0; i < psfTry->sources->n; i++) { pmSource *source = psfTry->sources->data[i]; - if (!source->moments) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no moments\n", i, source->peak->x, source->peak->y); - continue; - } - if (!source->moments->nPixels) { - psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no pixels\n", i, source->peak->x, source->peak->y); - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } - // If mask object does not exist, mark the source as bad. - // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes - // psphot to exit with a fault. - if (source->maskObj == NULL) { - psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y); - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } - source->modelEXT = pmSourceModelGuess (source, options->type, maskVal, markVal); - if (source->modelEXT == NULL) { - psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y); - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } + if (!source->moments) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no moments\n", i, source->peak->x, source->peak->y); + continue; + } + if (!source->moments->nPixels) { + psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no pixels\n", i, source->peak->x, source->peak->y); + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + continue; + } + // If mask object does not exist, mark the source as bad. + // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes + // psphot to exit with a fault. + if (source->maskObj == NULL) { + psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y); + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + continue; + } - // set object mask to define valid pixels - // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)? - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal); + source->modelEXT = pmSourceModelGuess (source, options->type, maskVal, markVal); + if (source->modelEXT == NULL) { + psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y); + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + continue; + } - // fit model as EXT, not PSF - status = pmSourceFitModel (source, source->modelEXT, options->fitOptions, maskVal); + // do some actual work on this source + psThreadJob *job = psThreadJobAlloc ("PSF_TRY_FIT_EXT"); + psArrayAdd(job->args, 1, source); + psArrayAdd(job->args, 1, psfTry); + psArrayAdd(job->args, 1, options); + + PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); - // clear object mask to define valid pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); - // exclude the poor fits - if (!status) { - psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y); - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } - Next ++; +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to create psf model."); + return false; + } +# else + if (!pmPSFtryFitEXT_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to create psf model."); + return false; + } + psFree(job); +# endif } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psError(PS_ERR_UNKNOWN, false, "failure to model psf"); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + // we have no returned data from this operation + if (job->args->n < 1) fprintf (stderr, "error with job\n"); + psFree(job); + } + psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, psfTry->sources->n); psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, psfTry->sources->n); @@ -118,2 +157,34 @@ return true; } + +bool pmPSFtryFitEXT_Threaded (psThreadJob *job) { + + pmSource *source = job->args->data[0]; + pmPSFtry *psfTry = job->args->data[1]; + pmPSFOptions *options = job->args->data[2]; + + int i = PS_SCALAR_VALUE(job->args->data[3], S32); + + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); + + // set object mask to define valid pixels + // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)? + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal); + + // fit model as EXT, not PSF + bool status = pmSourceFitModel (source, source->modelEXT, options->fitOptions, maskVal); + + // clear object mask to define valid pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + + // exclude the poor fits + if (!status) { + psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y); + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + return true; + } + Next ++; + + return true; +} Index: trunk/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryFitPSF.c (revision 35726) +++ trunk/psModules/src/objects/pmPSFtryFitPSF.c (revision 35768) @@ -43,10 +43,10 @@ #include "pmSourceVisual.h" +static int Npsf = 0; + // stage 3: Refit with fixed shape parameters. This function uses the LMM fitting, but could // be re-written to use the simultaneous linear fitting (see psphotFitSourcesLinear.c) bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) { - bool status; - psTimerStart ("psf.fit"); @@ -57,13 +57,5 @@ maskVal |= markVal; - // DEBUG code: save the PSF model fit data in detail -# ifdef DEBUG - char filename[64]; - snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy); - FILE *f = fopen (filename, "w"); - psAssert (f, "failed open"); -# endif - - int Npsf = 0; + Npsf = 0; for (int i = 0; i < psfTry->sources->n; i++) { @@ -77,56 +69,68 @@ } - // set shape for this model based on PSF + // set shape for this model based on PSF psFree (source->modelPSF); - source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf); - if (source->modelPSF == NULL) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL; - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y); - continue; - } - // PSF fit and aperture mags use different radii - source->modelPSF->fitRadius = options->fitRadius; - source->apRadius = options->apRadius; - - // set object mask to define valid pixels for PSF model fit - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal); - - // fit the PSF model to the source - status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal); - - // skip poor fits - if (!status) { - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL; - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); - continue; - } - - // set object mask to define valid pixels for APERTURE magnitude - if (options->fitRadius != options->apRadius) { - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal); + source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf); + if (source->modelPSF == NULL) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL; + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y); + return false; } - // This function calculates the psf and aperture magnitudes - status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag - if (!status || isnan(source->apMag)) { - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT; - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); - continue; - } - - // clear object mask to define valid pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask - - psfTry->fitMag->data.F32[i] = source->psfMag; - psfTry->metric->data.F32[i] = source->apMag - source->psfMag; - psfTry->metricErr->data.F32[i] = source->psfMagErr; - - // XXX this did not work: modifies shape of psf too much - // psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag; - + // do some actual work on this source + psThreadJob *job = psThreadJobAlloc ("PSF_TRY_FIT_PSF"); + psArrayAdd(job->args, 1, source); + psArrayAdd(job->args, 1, psfTry); + psArrayAdd(job->args, 1, options); + + PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); + + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to create psf model."); + return false; + } +# else + if (!pmPSFtryFitPSF_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to create psf model."); + return false; + } + psFree(job); +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psError(PS_ERR_UNKNOWN, false, "failure to model psf"); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + // we have no returned data from this operation + if (job->args->n < 1) fprintf (stderr, "error with job\n"); + psFree(job); + } + psfTry->psf->nPSFstars = Npsf; + + // DEBUG code: save the PSF model fit data in detail # ifdef DEBUG + + char filename[64]; + snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy); + FILE *f = fopen (filename, "w"); + psAssert (f, "failed open"); + + for (int i = 0; i < psfTry->sources->n; i++) { + + // skip masked sources + if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue; + + pmSource *source = psfTry->sources->data[i]; + fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n", source->peak->xf, source->peak->yf, @@ -136,12 +140,5 @@ source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_7]); -# endif - - psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n); - Npsf ++; - } - psfTry->psf->nPSFstars = Npsf; - -# ifdef DEBUG + } fclose (f); # endif @@ -159,2 +156,56 @@ return true; } + +bool pmPSFtryFitPSF_Threaded (psThreadJob *job) { + + pmSource *source = job->args->data[0]; + pmPSFtry *psfTry = job->args->data[1]; + pmPSFOptions *options = job->args->data[2]; + + int i = PS_SCALAR_VALUE(job->args->data[3], S32); + + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); + + // PSF fit and aperture mags use different radii + source->modelPSF->fitRadius = options->fitRadius; + source->apRadius = options->apRadius; + + // set object mask to define valid pixels for PSF model fit + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal); + + // fit the PSF model to the source + bool status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal); + + // skip poor fits + if (!status) { + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL; + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); + return true; + } + + // set object mask to define valid pixels for APERTURE magnitude + if (options->fitRadius != options->apRadius) { + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal); + } + + // This function calculates the psf and aperture magnitudes + status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + + if (!status || isnan(source->apMag)) { + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT; + psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); + return true; + } + + psfTry->fitMag->data.F32[i] = source->psfMag; + psfTry->metric->data.F32[i] = source->apMag - source->psfMag; + psfTry->metricErr->data.F32[i] = source->psfMagErr; + + psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n); + Npsf ++; + return true; +} Index: trunk/psModules/src/objects/pmPSFtryMakePSF.c =================================================================== --- trunk/psModules/src/objects/pmPSFtryMakePSF.c (revision 35726) +++ trunk/psModules/src/objects/pmPSFtryMakePSF.c (revision 35768) @@ -165,8 +165,5 @@ for (int i = 0; i < psf->params->n; i++) { - if (psf->params->data[i] == NULL) { - psFree(modelPSF); - continue; - } + if (psf->params->data[i] == NULL) continue; fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]); } @@ -214,5 +211,6 @@ assert (source->modelEXT); // all unmasked sources should have modelEXT - psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); + bool useReff = pmModelUseReff (source->modelEXT->type); + psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32, useReff); e0->data.F32[i] = pol.e0; Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 35726) +++ trunk/psModules/src/objects/pmSource.c (revision 35768) @@ -1145,5 +1145,4 @@ bool status; psEllipseShape oldshape; - psEllipseShape newshape; psEllipseAxes axes; @@ -1166,13 +1165,13 @@ if (!isfinite(oldI0)) return false; + bool useReff = pmModelUseReff (model->type); + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); + // increase size and height of source - axes = psEllipseShapeToAxes (oldshape, 20.0); axes.major *= SIZE; axes.minor *= SIZE; - newshape = psEllipseAxesToShape (axes); + + pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], axes, useReff); PAR[PM_PAR_I0] = FACTOR*oldI0; - PAR[PM_PAR_SXX] = newshape.sx; - PAR[PM_PAR_SYY] = newshape.sy; - PAR[PM_PAR_SXY] = newshape.sxy; psImage *target = source->variance; Index: trunk/psModules/src/objects/pmSourceFitModel.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceFitModel.c (revision 35768) @@ -62,4 +62,8 @@ opt->poissonErrors = true; opt->saveCovariance = false; + + // we default to the old algorithm + opt->gainFactorMode = 0; + opt->chisqConvergence = true; return opt; @@ -241,4 +245,6 @@ psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol); + myMin->gainFactorMode = options->gainFactorMode; + myMin->chisqConvergence = options->chisqConvergence; psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); @@ -252,4 +258,5 @@ } if (options->saveCovariance) { + psFree (model->covar); model->covar = psMemIncrRefCounter(covar); } @@ -273,5 +280,10 @@ model->flags |= PM_MODEL_STATUS_FITTED; if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE; - if (myMin->lastDelta > myMin->minTol) model->flags |= PM_MODEL_STATUS_WEAK_FIT; + + if (myMin->chisqConvergence) { + if (myMin->lastDelta > myMin->minTol) model->flags |= PM_MODEL_STATUS_WEAK_FIT; + } else { + if (myMin->rParSigma > myMin->minTol*nParams) model->flags |= PM_MODEL_STATUS_WEAK_FIT; + } // get the Gauss-Newton distance for fixed model parameters Index: trunk/psModules/src/objects/pmSourceFitModel.h =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.h (revision 35726) +++ trunk/psModules/src/objects/pmSourceFitModel.h (revision 35768) @@ -35,4 +35,6 @@ bool poissonErrors; ///< use poisson errors for fits? bool saveCovariance; + int gainFactorMode; + bool chisqConvergence; } pmSourceFitOptions; Index: trunk/psModules/src/objects/pmSourceFitPCM.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitPCM.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceFitPCM.c (revision 35768) @@ -66,4 +66,6 @@ // set up the minimization process psMinimization *myMin = psMinimizationAlloc (fitOptions->nIter, fitOptions->minTol, fitOptions->maxTol); + myMin->chisqConvergence = fitOptions->chisqConvergence; + myMin->gainFactorMode = fitOptions->gainFactorMode; psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); @@ -118,4 +120,10 @@ if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE; + if (myMin->chisqConvergence) { + if (myMin->lastDelta > myMin->minTol) pcm->modelConv->flags |= PM_MODEL_STATUS_WEAK_FIT; + } else { + if (myMin->rParSigma > myMin->minTol*pcm->nPar) pcm->modelConv->flags |= PM_MODEL_STATUS_WEAK_FIT; + } + // once we have fitted a model, we need to record that this model is a PCM model: pcm->modelConv->isPCM = true; @@ -145,9 +153,12 @@ bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) { - if (!pcm->modelConv->modelGuess(pcm->modelConv, source, maskVal, markVal)) { + if (!pcm->modelConv->modelGuess(pcm->modelConv, source, maskVal, markVal)) { return false; } return true; + // the following was an attempt to make analytical modifications to the shape terms based on the psf + // this has been replaced with a more empirical approach +# if (0) // generate copy of the model // XXX we could modify the parameter values or even the model @@ -196,3 +207,4 @@ return true; +# endif } Index: trunk/psModules/src/objects/pmSourceFitSet.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitSet.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceFitSet.c (revision 35768) @@ -568,4 +568,6 @@ psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol); + myMin->gainFactorMode = options->gainFactorMode; + myMin->chisqConvergence = options->chisqConvergence; psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); Index: trunk/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 35726) +++ trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 35768) @@ -925,5 +925,5 @@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]); } else { - psEllipseAxes axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", axes.major); psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", axes.minor); Index: trunk/psModules/src/objects/pmSourceIO_CMP.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 35768) @@ -135,5 +135,5 @@ lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky); - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + axes = pmPSF_ModelToAxes (PAR, model->type); float psfMagErr = isfinite(source->psfMagErr) ? source->psfMagErr : 999; Index: trunk/psModules/src/objects/pmSourceIO_OBJ.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_OBJ.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceIO_OBJ.c (revision 35768) @@ -92,5 +92,5 @@ } - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + axes = pmPSF_ModelToAxes (PAR, model->type); psLineInit (line); Index: trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 35768) @@ -114,5 +114,5 @@ yErr = dPAR[PM_PAR_YPOS]; if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + axes = pmPSF_ModelToAxes (PAR, model->type); } else { axes.major = NAN; @@ -623,5 +623,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + axes = pmPSF_ModelToAxes (PAR, model->type); // generate RA,DEC Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 35768) @@ -90,5 +90,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + axes = pmPSF_ModelToAxes (PAR, model->type); } else { // XXX: This code seg faults if source->peak is NULL. Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 35768) @@ -96,5 +96,5 @@ yErr = dPAR[PM_PAR_YPOS]; if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + axes = pmPSF_ModelToAxes (PAR, model->type); } else { axes.major = NAN; @@ -523,5 +523,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + axes = pmPSF_ModelToAxes (PAR, model->type); row = psMetadataAlloc (); Index: trunk/psModules/src/objects/pmSourceIO_SMPDATA.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 35768) @@ -92,5 +92,5 @@ lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky); - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + axes = pmPSF_ModelToAxes (PAR, model->type); } else { Index: trunk/psModules/src/objects/pmSourceIO_SX.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_SX.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceIO_SX.c (revision 35768) @@ -81,5 +81,5 @@ // pmSourceSextractType (source, &type, &flags); - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + axes = pmPSF_ModelToAxes (PAR, model->type); psLineInit (line); Index: trunk/psModules/src/objects/pmSourceOutputs.c =================================================================== --- trunk/psModules/src/objects/pmSourceOutputs.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceOutputs.c (revision 35768) @@ -107,5 +107,5 @@ } if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXY]) && isfinite(PAR[PM_PAR_SYY])) { - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + axes = pmPSF_ModelToAxes (PAR, model->type); outputs->psfMajor = axes.major; outputs->psfMinor = axes.minor; Index: trunk/psModules/src/objects/pmSourcePlotPSFModel.c =================================================================== --- trunk/psModules/src/objects/pmSourcePlotPSFModel.c (revision 35726) +++ trunk/psModules/src/objects/pmSourcePlotPSFModel.c (revision 35768) @@ -146,5 +146,5 @@ // force the axis ratio to be < 20.0 psEllipseAxes axes_mnt = psEllipseMomentsToAxes (moments, 20.0); - psEllipseAxes axes_psf = pmPSF_ModelToAxes (PAR, 20.0, model->type); + psEllipseAxes axes_psf = pmPSF_ModelToAxes (PAR, model->type); // moments major axis Index: trunk/psModules/src/objects/pmSourceVisual.c =================================================================== --- trunk/psModules/src/objects/pmSourceVisual.c (revision 35726) +++ trunk/psModules/src/objects/pmSourceVisual.c (revision 35768) @@ -46,5 +46,5 @@ bool pmSourceVisualClose() { if (kapa1 != -1) - KiiClose(kapa1); + KapaClose(kapa1); return true; } Index: trunk/psModules/src/psmodules.h =================================================================== --- trunk/psModules/src/psmodules.h (revision 35726) +++ trunk/psModules/src/psmodules.h (revision 35768) @@ -158,4 +158,5 @@ #include #include +#include #endif