Changeset 15000

Timestamp:

Sep 24, 2007, 11:27:58 AM (19 years ago)

Author:

eugene

Message:

update from eam_branch_20070921

Location:

trunk

Files:

• psModules/src/camera/pmFPAfileIO.c (modified) (1 diff)
• psModules/src/objects/models/pmModel_GAUSS.c (modified) (2 diffs)
• psModules/src/objects/models/pmModel_PGAUSS.c (modified) (2 diffs)
• psModules/src/objects/models/pmModel_QGAUSS.c (modified) (2 diffs)
• psModules/src/objects/models/pmModel_RGAUSS.c (modified) (2 diffs)
• psModules/src/objects/models/pmModel_SERSIC.c (modified) (2 diffs)
• psModules/src/objects/pmPSF.c (modified) (11 diffs)
• psModules/src/objects/pmPSF.h (modified) (6 diffs)
• psModules/src/objects/pmPSF_IO.c (modified) (12 diffs)
• psModules/src/objects/pmPSF_IO.h (modified) (2 diffs)
• psModules/src/objects/pmPSFtoMetadata.c (added)
• psModules/src/objects/pmPSFtry.c (modified) (20 diffs)
• psModules/src/objects/pmPSFtry.h (modified) (4 diffs)
• psModules/src/objects/pmSourcePhotometry.c (modified) (5 diffs)
• psModules/src/objects/pmTrend2D.c (modified) (10 diffs)
• psModules/src/objects/pmTrend2D.h (modified) (4 diffs)
• psphot/doc/psfmodel.txt (modified) (1 diff)
• psphot/src/models/pmModel_STRAIL.c (modified) (2 diffs)
• psphot/src/models/pmModel_TEST1.c (modified) (2 diffs)
• psphot/src/psphot.h (modified) (1 diff)
• psphot/src/psphotApResid.c (modified) (14 diffs)
• psphot/src/psphotChoosePSF.c (modified) (8 diffs)
• psphot/src/psphotRoughClass.c (modified) (2 diffs)
• psphot/src/psphotSourcePlots.c (modified) (4 diffs)

trunk/psModules/src/camera/pmFPAfileIO.c

@@ -774 +774 @@
     }
     pmFPAview *thisView = pmFPAAddSourceFromHeader (file->fpa, phu, file->format);
+    assert (thisView); // XXX we are having some trouble with input psf files not having the Cell and fpa names matching.
     psFree (thisView);
     psFree (phu);

trunk/psModules/src/objects/models/pmModel_GAUSS.c

@@ -274 +274 @@
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params->data[i];
-            out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
+            pmTrend2D *trend = psf->params->data[i];
+            out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
     }
@@ -331 +331 @@
         if (i == PM_PAR_XPOS) continue;
         if (i == PM_PAR_YPOS) continue;
-        psPolynomial2D *poly = psf->params->data[i];
-        assert (poly);
-        PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+        pmTrend2D *trend = psf->params->data[i];
+        PAR[i] = pmTrend2DEval(trend, Xo, Yo);
     }
 

trunk/psModules/src/objects/models/pmModel_PGAUSS.c

@@ -328 +328 @@
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params->data[i];
-            out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
+            pmTrend2D *trend = psf->params->data[i];
+            out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
     }
@@ -384 +384 @@
         if (i == PM_PAR_XPOS) continue;
         if (i == PM_PAR_YPOS) continue;
-        psPolynomial2D *poly = psf->params->data[i];
-        assert (poly);
-        PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+        pmTrend2D *trend = psf->params->data[i];
+        PAR[i] = pmTrend2DEval(trend, Xo, Yo);
     }
 

trunk/psModules/src/objects/models/pmModel_QGAUSS.c

@@ -358 +358 @@
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params->data[i];
-            out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
+            pmTrend2D *trend = psf->params->data[i];
+            out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
     }
@@ -412 +412 @@
         if (i == PM_PAR_XPOS) continue;
         if (i == PM_PAR_YPOS) continue;
-        psPolynomial2D *poly = psf->params->data[i];
-        assert (poly);
-        PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+        pmTrend2D *trend = psf->params->data[i];
+        PAR[i] = pmTrend2DEval(trend, Xo, Yo);
     }
 

trunk/psModules/src/objects/models/pmModel_RGAUSS.c

@@ -351 +351 @@
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params->data[i];
-            out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
+            pmTrend2D *trend = psf->params->data[i];
+            out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
     }
@@ -405 +405 @@
         if (i == PM_PAR_XPOS) continue;
         if (i == PM_PAR_YPOS) continue;
-        psPolynomial2D *poly = psf->params->data[i];
-        assert (poly);
-        PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+        pmTrend2D *trend = psf->params->data[i];
+        PAR[i] = pmTrend2DEval(trend, Xo, Yo);
     }
 

trunk/psModules/src/objects/models/pmModel_SERSIC.c

@@ -343 +343 @@
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params->data[i];
-            out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
+            pmTrend2D *trend = psf->params->data[i];
+            out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
     }
@@ -397 +397 @@
         if (i == PM_PAR_XPOS) continue;
         if (i == PM_PAR_YPOS) continue;
-        psPolynomial2D *poly = psf->params->data[i];
-        assert (poly);
-        PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+        pmTrend2D *trend = psf->params->data[i];
+        PAR[i] = pmTrend2DEval(trend, Xo, Yo);
     }
 

trunk/psModules/src/objects/pmPSF.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-21 02:46:25 $
+ *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-24 21:27:58 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -41 +41 @@
 
 /*****************************************************************************/
-/* DEFINE STATEMENTS                                                         */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* TYPE DEFINITIONS                                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* GLOBAL VARIABLES                                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* FILE STATIC VARIABLES                                                     */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* FUNCTION IMPLEMENTATION - LOCAL                                           */
-/*****************************************************************************/
-
-/*****************************************************************************/
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
 /*****************************************************************************/
+
+static void pmPSFOptionsFree (pmPSFOptions *options) {
+
+    if (!options) return;
+
+    psFree (options->stats);
+    return;
+}
+
+pmPSFOptions *pmPSFOptionsAlloc () {
+
+    pmPSFOptions *options = (pmPSFOptions *) psAlloc(sizeof(pmPSFOptions));
+    psMemSetDeallocator(options, (psFreeFunc) pmPSFOptionsFree);
+
+    options->type          = 0;
+
+    options->stats         = NULL;
+
+    options->psfTrendMode  = PM_TREND_NONE;
+    options->psfTrendNx    = 0;
+    options->psfTrendNy    = 0;
+    options->psfFieldNx    = 0;
+    options->psfFieldNy    = 0;
+    options->psfFieldXo    = 0;
+    options->psfFieldYo    = 0;
+
+    options->poissonErrorsPhotLMM = true;
+    options->poissonErrorsPhotLin = false;
+    options->poissonErrorsParams  = true;
+
+    return options;
+}
 
 /*****************************************************************************
@@ -74 +86 @@
 
     psFree (psf->ChiTrend);
+    psFree (psf->psfTrendStats);
     psFree (psf->ApTrend);
     psFree (psf->FluxScale);
@@ -94 +107 @@
  Object Normalization
  *****************************************************************************/
-pmPSF *pmPSFAlloc (pmModelType type, bool poissonErrors, psPolynomial2D *psfTrendMask)
+pmPSF *pmPSFAlloc (pmPSFOptions *options)
 {
     int Nparams;
 
     pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF));
-
-    psf->type     = type;
+    psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
+
+    psf->type     = options->type;
     psf->chisq    = 0.0;
     psf->ApResid  = 0.0;
@@ -108 +122 @@
     psf->nPSFstars  = 0;
     psf->nApResid   = 0;
-    psf->poissonErrors = poissonErrors;
+
+    psf->poissonErrorsPhotLMM = options->poissonErrorsPhotLMM;
+    psf->poissonErrorsPhotLin = options->poissonErrorsPhotLin;
+    psf->poissonErrorsParams = options->poissonErrorsParams;
 
     // the ApTrend components are (x, y).  It may be represented with a polynomial or with a
@@ -121 +138 @@
     psf->FluxScale = NULL;
 
-    if (psf->poissonErrors) {
+    if (psf->poissonErrorsPhotLMM) {
         psf->ChiTrend = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     } else {
@@ -133 +150 @@
     psf->residuals = NULL;
 
-    Nparams = pmModelClassParameterCount (type);
+    Nparams = pmModelClassParameterCount (options->type);
     if (!Nparams) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -139 +156 @@
     }
     psf->params = psArrayAlloc(Nparams);
+    
+    // save the trend stats on the psf for use in pmPSFFromPSFtry
+    psf->psfTrendStats = psMemIncrRefCounter (options->stats);
 
     // the order of the PSF parameter (X,Y) fits is determined by the psfTrendMask polynomial
@@ -147 +167 @@
     // PM_PAR_XPOS, etc)
 
-    if (psfTrendMask) {
+    psImageBinning *binning = psImageBinningAlloc();
+    binning->nXruff = options->psfTrendNx;
+    binning->nYruff = options->psfTrendNy;
+    binning->nXfine = options->psfFieldNx;
+    binning->nYfine = options->psfFieldNy;
+    psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
+    psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo);
+
+    psf->fieldNx = options->psfFieldNx;
+    psf->fieldNy = options->psfFieldNy;
+    psf->fieldXo = options->psfFieldXo;
+    psf->fieldYo = options->psfFieldYo;
+
+    // define the parameter trends
+    if (options->psfTrendMode != PM_TREND_NONE) {
         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;
-
-            psPolynomial2D *param = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, psfTrendMask->nX, psfTrendMask->nY);
-            for (int nx = 0; nx < param->nX + 1; nx++) {
-                for (int ny = 0; ny < param->nY + 1; ny++) {
-                    param->mask[nx][ny] = psfTrendMask->mask[nx][ny];
-                }
-            }
-            psf->params->data[i] = param;
+            if (i == PM_PAR_SKY) continue;
+            if (i == PM_PAR_I0) continue;
+            if (i == PM_PAR_XPOS) continue;
+            if (i == PM_PAR_YPOS) continue;
+
+            psf->params->data[i] = pmTrend2DNoImageAlloc (options->psfTrendMode, binning, options->stats);
         }
     }
-
-    psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
-    return(psf);
+    psFree (binning);
+    return psf;
 }
 
@@ -300 +323 @@
     va_start(ap, sxy);
 
-    pmModelType type = pmModelClassGetType (typeName);
-    psPolynomial2D *psfTrend = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 0, 0);
-    pmPSF *psf = pmPSFAlloc (type, true, psfTrend);
+    pmPSFOptions *options = pmPSFOptionsAlloc ();
+    options->type = pmModelClassGetType (typeName);
+    options->psfTrendMode = PM_TREND_POLY_ORD;
+    options->psfTrendNx = 0;
+    options->psfTrendNy = 0;
+
+    pmPSF *psf = pmPSFAlloc (options);
 
     psVector *par = psVectorAlloc (psf->params->n, PS_TYPE_F32);
@@ -311 +338 @@
     psEllipsePol pol = pmPSF_ModelToFit(par->data.F32);
 
+    pmTrend2D *trend = NULL;
+
     // set the psf shape parameters
-    psPolynomial2D *poly = NULL;
-    poly = psf->params->data[PM_PAR_E0];
-    poly->coeff[0][0] = pol.e0;
-
-    poly = psf->params->data[PM_PAR_E1];
-    poly->coeff[0][0] = pol.e1;
-
-    poly = psf->params->data[PM_PAR_E2];
-    poly->coeff[0][0] = pol.e2;
+    trend = psf->params->data[PM_PAR_E0];
+    trend->poly->coeff[0][0] = pol.e0;
+
+    trend = psf->params->data[PM_PAR_E1];
+    trend->poly->coeff[0][0] = pol.e1;
+
+    trend = psf->params->data[PM_PAR_E2];
+    trend->poly->coeff[0][0] = pol.e2;
 
     for (int i = PM_PAR_SXY + 1; i < psf->params->n; i++) {
-        poly = psf->params->data[i];
-        poly->coeff[0][0] = (psF32)va_arg(ap, psF64);
+        trend = psf->params->data[i];
+        trend->poly->coeff[0][0] = (psF32)va_arg(ap, psF64);
     }
     va_end(ap);
 
     psFree (par);
-    psFree (psfTrend);
+    psFree (options);
     return psf;
 }

trunk/psModules/src/objects/pmPSF.h

@@ -6 +6 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-09-21 00:05:35 $
+ * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-09-24 21:27:58 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -19 +19 @@
 // type of model carried by the pmModel structure
 typedef int pmModelType;
-
-typedef enum {
-    PM_PSF_APTREND_ERROR = 0,
-    PM_PSF_APTREND_NONE,
-    PM_PSF_APTREND_CONSTANT,
-    PM_PSF_APTREND_SKYBIAS,
-    PM_PSF_APTREND_SKYSAT,
-    PM_PSF_APTREND_XY_LIN,
-    PM_PSF_APTREND_XY_QUAD,
-    PM_PSF_APTREND_SKY_XY_LIN,
-    PM_PSF_APTREND_SKYSAT_XY_LIN,
-    PM_PSF_APTREND_ALL
-} pmPSFApTrendOptions;
 
 /** pmPSF data structure
@@ -49 +36 @@
     pmModelType type;                   ///< PSF Model in use
     psArray *params;                    ///< Model parameters (psPolynomial2D)
+    psStats *psfTrendStats;             ///< psf parameter trend clipping stats
     psPolynomial1D *ChiTrend;           ///< Chisq vs flux fit (correction for systematic errors)
     pmTrend2D *ApTrend;                 ///< ApResid vs (x,y)
@@ -59 +47 @@
     int nPSFstars;                      ///< number of stars used to measure PSF
     int nApResid;                       ///< number of stars used to measure ApResid
-    bool poissonErrors;
+    int fieldNx;
+    int fieldNy;
+    int fieldXo;
+    int fieldYo;
+    bool poissonErrorsPhotLMM;          ///< use poission errors for non-linear model fitting
+    bool poissonErrorsPhotLin;          ///< use poission errors for linear model fitting
+    bool poissonErrorsParams;           ///< use poission errors for model parameter fitting
     pmGrowthCurve *growth;              ///< apMag vs Radius
     pmResiduals *residuals;             ///< normalized residual image (no spatial variation)
 }
 pmPSF;
+
+typedef struct {
+    pmModelType   type;
+    psStats      *stats;                // psfTrend clipping stats
+    pmTrend2DMode psfTrendMode;
+    int           psfTrendNx;
+    int           psfTrendNy;
+    int           psfFieldNx;
+    int           psfFieldNy;
+    int           psfFieldXo;
+    int           psfFieldYo;
+    bool          poissonErrorsPhotLMM; ///< use poission errors for non-linear model fitting
+    bool          poissonErrorsPhotLin; ///< use poission errors for linear model fitting
+    bool          poissonErrorsParams; ///< use poission errors for model parameter fitting
+    float         radius;
+} pmPSFOptions;
 
 # define PM_PAR_E0 PM_PAR_SXX
@@ -74 +84 @@
  *
  */
-pmPSF *pmPSFAlloc(
-    pmModelType type,     // type of model for PSF
-    bool poissonErrors,
-    psPolynomial2D *psfTrendMask
-);
-
-bool pmPSFMaskApTrend (psPolynomial4D *trend, pmPSFApTrendOptions option);
-
-pmPSFApTrendOptions pmPSFApTrendOptionFromName (char *name);
+pmPSF *pmPSFAlloc (pmPSFOptions *options);
+pmPSFOptions *pmPSFOptionsAlloc ();
 
 double pmPSF_SXYfromModel (psF32 *modelPar);
@@ -91 +94 @@
 
 bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes);
+bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis);
+
+psEllipsePol pmPSF_ModelToFit (psF32 *modelPar);
 psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR);
-bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis);
-psEllipsePol pmPSF_ModelToFit (psF32 *modelPar);
 
 /// @}

trunk/psModules/src/objects/pmPSF_IO.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-21 00:05:57 $
+ *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-24 21:27:58 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -148 +148 @@
 
 // for a pmPSF supplied on the analysis metadata, we write out
+// if needed:
+//   - a PHU blank header
 // - image header        : FITS Image NAXIS = 0
-// - psf table (+header) : FITS Table
-// - psf resid (+header) : FITS Image
-// if needed, we also write out a PHU blank header
+// if (trendMode == MAP) 
+//   - psf resid (+header) : FITS Image
+// else
+//   - psf table (+header) : FITS Table
 bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
@@ -159 +162 @@
 
     if (!analysis) return false;
+
+    // select the current recipe
+    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT");
+    if (!recipe) {
+        psError(PS_ERR_UNKNOWN, false, "missing recipe %s", "PSPHOT");
+        return false;
+    }
 
     // write a PHU? (only if input image is MEF)
@@ -215 +225 @@
 
     // write out the IMAGE header segment
-    // this header block is new, write it to disk
+    // if this header block is new, write it to disk
     if (hdu->header != file->header) {
         // add EXTNAME, EXTHEAD, EXTTYPE to header
@@ -250 +260 @@
         psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName);
 
-        psMetadataAddBool (header, PS_LIST_TAIL, "POISSON",  0, "Use Poisson errors in fits?", psf->poissonErrors);
+        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LMM",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLMM);
+        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LIN",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLin);
+        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_PAR",  0, "Use Poisson errors in fits?", psf->poissonErrorsParams);
 
         int nPar = pmModelClassParameterCount (psf->type)    ;
         psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
+
+        psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS1", 0, "Image X Size", psf->fieldNx);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS2", 0, "Image Y Size", psf->fieldNy);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF1",  0, "Image X Ref",  psf->fieldXo);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF2",  0, "Image Y Ref",  psf->fieldYo);
+
+        // extract PSF Clump info
+        pmPSFClump psfClump;
+
+        psfClump.X  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X");   assert (status);
+        psfClump.Y  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y");   assert (status);
+        psfClump.dX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX");  assert (status);
+        psfClump.dY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY");  assert (status);
+
+        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLX", 0, "psf clump center", psfClump.X);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLY", 0, "psf clump center", psfClump.Y);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDX", 0, "psf clump size", psfClump.dX);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDY", 0, "psf clump size", psfClump.dY);
 
         // save the dimensions of each parameter
         for (int i = 0; i < nPar; i++) {
             char name[9];
-            psPolynomial2D *poly = psf->params->data[i];
-            if (poly == NULL) continue;
+            int nX, nY;
+
+            pmTrend2D *trend = psf->params->data[i];
+            if (trend == NULL) continue;
+
+            if (trend->mode == PM_TREND_MAP) {
+              nX = trend->map->map->numCols;
+              nY = trend->map->map->numRows;
+            } else {
+              nX = trend->poly->nX;
+              nY = trend->poly->nY;
+            }
             snprintf (name, 9, "PAR%02d_NX", i);
-            psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX);
+            psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nX);
             snprintf (name, 9, "PAR%02d_NY", i);
-            psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY);
+            psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nY);
+            snprintf (name, 9, "PAR%02d_MD", i);
+            char *modeName = pmTrend2DModeToString (trend->mode);
+            psMetadataAddStr (header, PS_LIST_TAIL, name, 0, "", modeName);
+            psFree (modeName);
         }
 
         // other required information describing the PSF
-        psMetadataAddF32  (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
-        psMetadataAddF32  (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
-        psMetadataAddF32  (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
-        psMetadataAddS32  (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);
 
         // XXX can we drop this now?
@@ -278 +322 @@
         psArray *psfTable = psArrayAllocEmpty (100);
         for (int i = 0; i < nPar; i++) {
-            psPolynomial2D *poly = psf->params->data[i];
-            if (poly == NULL) continue; // skip unset parameters (eg, XPOS)
-            for (int ix = 0; ix <= poly->nX; ix++) {
-                for (int iy = 0; iy <= poly->nY; iy++) {
-
-                    psMetadata *row = psMetadataAlloc ();
-                    psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
-                    psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
-                    psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
-                    psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
-
-                    psArrayAdd (psfTable, 100, row);
-                    psFree (row);
+            pmTrend2D *trend = psf->params->data[i];
+            if (trend == NULL) continue; // skip unset parameters (eg, XPOS)
+
+            if (trend->mode == PM_TREND_MAP) {
+                // write the image components into a table: this is needed because they may each be a different size
+                psImageMap *map = trend->map;
+                for (int ix = 0; ix < map->map->numCols; ix++) {
+                    for (int iy = 0; iy < map->map->numRows; iy++) {
+                        psMetadata *row = psMetadataAlloc ();
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
+                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
+                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
+                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
+
+                        psArrayAdd (psfTable, 100, row);
+                        psFree (row);
+                    }
+                }
+            } else {
+                // write the polynomial components into a table
+                psPolynomial2D *poly = trend->poly;
+                for (int ix = 0; ix <= poly->nX; ix++) {
+                    for (int iy = 0; iy <= poly->nY; iy++) {
+                        psMetadata *row = psMetadataAlloc ();
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
+                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
+                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
+                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
+
+                        psArrayAdd (psfTable, 100, row);
+                        psFree (row);
+                    }
                 }
             }
@@ -356 +421 @@
     // XXX save the growth curve
     // XXX save ApTrend (as image?)
+    // XXX write the ApTrend with the same API as will be used for the PSF parameters above
 
 # if (0)
@@ -474 +540 @@
     psTrace ("psModules.objects", 5, "read psf model for %s\n", file->filename);
 
+    // select the current recipe
+    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT");
+    if (!recipe) {
+        psError(PS_ERR_UNKNOWN, false, "missing recipe %s", "PSPHOT");
+        return false;
+    }
+
     // Menu of EXTNAME rules
     psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
@@ -506 +579 @@
     if (!header) psAbort("cannot read table header");
 
+    pmPSFOptions *options = pmPSFOptionsAlloc();
+
     // load the PSF model parameters from the FITS table
     char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME");
-    pmModelType type = pmModelClassGetType (modelName);
-    if (type == -1) {
+    options->type = pmModelClassGetType (modelName);
+    if (options->type == -1) {
         psError(PS_ERR_UNKNOWN, true, "invalid model name %s in psf file %s", modelName, file->filename);
         return false;
     }
 
-    bool poissonErrors = psMetadataLookupBool (&status, header, "POISSON");
+    // psf clump data
+    pmPSFClump psfClump;
+
+    psfClump.X  = psMetadataLookupF32 (&status, header, "PSF_CLX" );  assert(status);
+    psfClump.Y  = psMetadataLookupF32 (&status, header, "PSF_CLY" );  assert(status);
+    psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX");  assert(status);
+    psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY");  assert(status);
+
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.X" , 0, "psf clump center", psfClump.X);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.Y" , 0, "psf clump center", psfClump.Y);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DX", 0, "psf clump size",   psfClump.dX);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DY", 0, "psf clump size",   psfClump.dY);
+
+    options->poissonErrorsPhotLMM = psMetadataLookupBool (&status, header, "ERR_LMM");
+    options->poissonErrorsPhotLin = psMetadataLookupBool (&status, header, "ERR_LIN");
+    options->poissonErrorsParams  = psMetadataLookupBool (&status, header, "ERR_PAR");
+
+    options->psfFieldNx = psMetadataLookupS32 (&status, header, "IMAXIS1");
+    options->psfFieldNy = psMetadataLookupS32 (&status, header, "IMAXIS2");
+    options->psfFieldXo = psMetadataLookupS32 (&status, header, "IMREF1");
+    options->psfFieldYo = psMetadataLookupS32 (&status, header, "IMREF2");
+
+    psImageBinning *binning = psImageBinningAlloc();
+    binning->nXfine = options->psfFieldNx;
+    binning->nYfine = options->psfFieldNy;
 
     // we determine the PSF parameter polynomials from the MD-defined polynomials
-    pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL);
+    pmPSF *psf = pmPSFAlloc (options);
 
     // check the number of expected parameters
@@ -524 +623 @@
         psAbort("mismatch model par count");
 
-    // load the dimensions of each parameter
+    // load the trend mode and dimensions of each parameter
     for (int i = 0; i < nPar; i++) {
         char name[9];
         snprintf (name, 9, "PAR%02d_NX", i);
-        int nXorder = psMetadataLookupS32 (&status, header, name);
+        binning->nXruff = psMetadataLookupS32 (&status, header, name);
         if (!status) continue;          // not all parameters are defined
+
         snprintf (name, 9, "PAR%02d_NY", i);
-        int nYorder = psMetadataLookupS32 (&status, header, name);
+        binning->nYruff = psMetadataLookupS32 (&status, header, name);
         if (!status) {
-            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but no NY", i);
+            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but not NY", i);
             return false;
         }
-        psf->params->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder);
-    }
+
+        snprintf (name, 9, "PAR%02d_MD", i);
+        char *modeName = psMetadataLookupStr (&status, header, name);
+        if (!status) {
+            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i);
+            return false;
+        }
+        pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName);
+        if (psfTrendMode == PM_TREND_NONE) {
+            psfTrendMode = PM_TREND_POLY_ORD;
+        }
+
+        psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
+        psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo);
+        psf->params->data[i] = pmTrend2DNoImageAlloc (psfTrendMode, binning, NULL);
+    }
+    psFree (binning);
 
     // other required information describing the PSF
@@ -557 +672 @@
         int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
         int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
-        // XXX sanity check here
-
-        psPolynomial2D *poly = psf->params->data[iPar];
-        if (poly == NULL) {
-            psError(PS_ERR_UNKNOWN, true, "values for parameter %d, but missing NX", iPar);
+
+        pmTrend2D *trend = psf->params->data[iPar];
+        if (trend == NULL) {
+            psError(PS_ERR_UNKNOWN, true, "parameter %d not available", iPar);
             return false;
         }
 
-        poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
-        poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
-        poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
+        if (trend->mode == PM_TREND_MAP) {
+            psImageMap *map = trend->map;
+            assert (map);
+            assert (map->map);
+            assert (map->error);
+            assert (xPow >= 0);
+            assert (yPow >= 0);
+            assert (xPow < map->map->numCols);
+            assert (yPow < map->map->numRows);
+            map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
+            map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
+        } else {
+            psPolynomial2D *poly = trend->poly;
+            assert (poly);
+            assert (xPow >= 0);
+            assert (yPow >= 0);
+            assert (xPow <= poly->nX);
+            assert (yPow <= poly->nY);
+            poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
+            poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
+            poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
+        }
     }
     psFree (header);
@@ -618 +751 @@
 
 // create a psMetadata representation (human-readable) of a psf model
-psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf)
-{
-
-    if (metadata == NULL) {
-        metadata = psMetadataAlloc ();
-    }
-
-    char *modelName = pmModelClassGetName (psf->type);
-    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
-
-    int nPar = pmModelClassParameterCount (psf->type)    ;
-    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
-
-    for (int i = 0; i < nPar; i++) {
-        psPolynomial2D *poly = psf->params->data[i];
-        if (poly == NULL)
-            continue;
-        psPolynomial2DtoMetadata (metadata, poly, "PSF_PAR%02d", i);
-    }
-
-    // XXX fix this
-    psWarning ("APTREND is currently missing");
-    // psPolynomial4DtoMetadata (metadata, psf->ApTrend, "APTREND");
-
-    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid);
-    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid);
-    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
-
-    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq);
-    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars);
-    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_POISSON_ERRORS", PS_DATA_BOOL, "Poisson errors for fits", psf->poissonErrors);
-
-    return metadata;
-}
-
-// parse a psMetadata representation (human-readable) of a psf model
-pmPSF *pmPSFfromMetadata (psMetadata *metadata)
-{
-
-    bool status;
-    char keyword[80];
-
-    char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME");
-    pmModelType type = pmModelClassGetType (modelName);
-
-    bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS");
-    if (!status)
-        poissonErrors = true;
-
-    // we determine the PSF parameter polynomials from the MD-defined polynomials
-    pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL);
-
-    int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR");
-    if (nPar != pmModelClassParameterCount (psf->type))
-        psAbort("mismatch model par count");
-
-    // un-fitted terms, not in the Metadata, are left NULL
-    // XXX add a double-check of the expected number?
-    for (int i = 0; i < nPar; i++) {
-        sprintf (keyword, "PSF_PAR%02d", i);
-        psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword);
-        if (!status)
-            continue;
-        psPolynomial2D *poly = psPolynomial2DfromMetadata (folder);
-        psFree (psf->params->data[i]);
-        psf->params->data[i] = poly;
-    }
-
-    // load the APTREND data
-    // XXX fix this to work with pmTrend2D
-    psWarning ("APTREND is not being read");
-    # if (0)
-    sprintf (keyword, "APTREND");
-    psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword);
-    psPolynomial4D *poly = psPolynomial4DfromMetadata (folder);
-    psFree (psf->ApTrend);
-    psf->ApTrend = poly;
-    # endif
-
-    psf->ApResid = psMetadataLookupF32 (&status, metadata, "PSF_AP_RESID");
-    psf->dApResid = psMetadataLookupF32 (&status, metadata, "PSF_dAP_RESID");
-    psf->skyBias = psMetadataLookupF32 (&status, metadata, "PSF_SKY_BIAS");
-
-    psf->chisq = psMetadataLookupF32 (&status, metadata, "PSF_CHISQ");
-    psf->nPSFstars = psMetadataLookupS32 (&status, metadata, "PSF_NSTARS");
-
-    psFree (metadata);
-    return (psf);
-}
+// XXX pmPSF to/from Metadata functions were defined for 1.22 and earlier, but were dropped

trunk/psModules/src/objects/pmPSF_IO.h

@@ -6 +6 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-08-24 00:11:02 $
+ * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-09-24 21:27:58 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -33 +33 @@
 bool pmPSFmodelCheckDataStatusForChip (const pmChip *chip);
 
-psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf);
-pmPSF *pmPSFfromMetadata (psMetadata *metadata);
-
 /// @}
 # endif

trunk/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-21 00:06:57 $
+ *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-24 21:27:58 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -44 +44 @@
 static void pmPSFtryFree (pmPSFtry *test)
 {
-
-    if (test == NULL)
-        return;
+    if (test == NULL) return;
 
     psFree (test->psf);
@@ -58 +56 @@
 
 // allocate a pmPSFtry based on the desired sources and the model (identified by name)
-pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors, psPolynomial2D *psfTrendMask)
+pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options)
 {
-
-    // validate the requested model name
-    pmModelType type = pmModelClassGetType (modelName);
-    if (type == -1) {
-        psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
-        return NULL;
-    }
-
     pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
-
-    test->psf       = pmPSFAlloc (type, poissonErrors, psfTrendMask);
+    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
+
+    test->psf       = pmPSFAlloc (options);
     test->metric    = psVectorAlloc (sources->n, PS_TYPE_F32);
     test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32);
@@ -76 +67 @@
     test->mask      = psVectorAlloc (sources->n, PS_TYPE_U8);
 
+    psVectorInit (test->mask,        0);
+    psVectorInit (test->metric,    0.0);
+    psVectorInit (test->metricErr, 0.0);
+    psVectorInit (test->fitMag,    0.0);
+
     test->sources   = psArrayAlloc (sources->n);
+
     for (int i = 0; i < sources->n; i++) {
         test->sources->data[i] = pmSourceCopy (sources->data[i]);
-        test->mask->data.U8[i]  = 0;
-        test->metric->data.F32[i] = 0;
-        test->metricErr->data.F32[i] = 0;
-        test->fitMag->data.F32[i] = 0;
-    }
-
-    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
+    }
+
     return (test);
 }
@@ -99 +91 @@
 
 // generate a pmPSFtry with a copy of the test PSF sources
-pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask, bool applyWeights, psMaskType maskVal, psMaskType mark)
+pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, pmPSFOptions *options, psMaskType maskVal, psMaskType mark)
 {
     bool status;
@@ -105 +97 @@
     int Npsf = 0;
 
-    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors, psfTrendMask);
+    // validate the requested model name
+    options->type = pmModelClassGetType (modelName);
+    if (options->type == -1) {
+        psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
+        return NULL;
+    }
+
+    pmPSFtry *psfTry = pmPSFtryAlloc (sources, options);
     if (psfTry == NULL) {
         psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model");
@@ -111 +110 @@
     }
 
-    // stage 1:  fit an independent model (freeModel) to all sources
+    // stage 1:  fit an EXT model to all candidates PSF sources
     psTimerStart ("fit");
     for (int i = 0; i < psfTry->sources->n; i++) {
@@ -123 +122 @@
         }
 
-        // set object mask to define valid pixels
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", mark);
+        // set object mask to define valid pixels -- XXX not unset?
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark);
 
         // fit model as EXT, not PSF
@@ -140 +139 @@
 
     // stage 2: construct a psf (pmPSF) from this collection of model fits
-    if (!pmPSFFromPSFtry (psfTry, applyWeights)) {
+    if (!pmPSFFromPSFtry (psfTry)) {
         psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
         psFree(psfTry);
@@ -162 +161 @@
             continue;
         }
-        source->modelPSF->radiusFit = RADIUS;
-
-        // set object mask to define valid pixels
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", mark);
+        source->modelPSF->radiusFit = options->radius;
+
+        // set object mask to define valid pixels -- XXX not unset?
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark);
 
         // fit the PSF model to the source
@@ -241 +240 @@
 
     // XXX this function wants aperture radius for pmSourcePhotometry
-    pmPSFtryMetric (psfTry, RADIUS);
+    pmPSFtryMetric (psfTry, options->radius);
     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
               modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
@@ -339 +338 @@
 /*****************************************************************************
 pmPSFFromPSFtry (psfTry): build a PSF model from a collection of
-source->modelEXT entires.  The PSF ignores the first 4 (independent) model
+source->modelEXT entries.  The PSF ignores the first 4 (independent) model
 parameters and constructs a polynomial fit to the remaining as a function of
 image coordinate.
@@ -345 +344 @@
 Note: some of the array entries may be NULL (failed fits); ignore them.
  *****************************************************************************/
-bool pmPSFFromPSFtry (pmPSFtry *psfTry, bool applyWeight)
+bool pmPSFFromPSFtry (pmPSFtry *psfTry)
 {
     pmPSF *psf = psfTry->psf;
@@ -355 +354 @@
 
     psVector *dz = NULL;
-    if (applyWeight) {
+    if (psf->poissonErrorsParams) {
         dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
     }
@@ -377 +376 @@
     // more deviant than three sigma.
     // mask is currently updated for each pass. this is inconsistent?
-
-    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
 
     // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
@@ -389 +386 @@
     // threshold.
 
+    // XXX this function needs to check the fit residuals and modify Nx,Ny as needed
+
     // convert the measured source shape paramters to polarization terms
     psVector *e0   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
@@ -406 +405 @@
     // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
     // This way, the parameters masked by one of the fits will be applied to the others
-    for (int i = 0; i < stats->clipIter; i++) {
-        psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y);
-        psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", stats->clippedNvalues, e0->n);
-        psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y);
-        psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", stats->clippedNvalues, e1->n);
-        psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y);
-        psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", stats->clippedNvalues, e2->n);
-    }
-
-    // XXX temporary dump of the psf parameters
+    for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
+
+        // XXX we are using the same stats structure on each pass: do we need to re-init it?
+        pmTrend2DFit (psf->params->data[PM_PAR_E0], psfTry->mask, 0xff, x, y, e0, dz);
+        psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n);
+        pmTrend2DFit (psf->params->data[PM_PAR_E1], psfTry->mask, 0xff, x, y, e1, dz);
+        psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n);
+        pmTrend2DFit (psf->params->data[PM_PAR_E2], psfTry->mask, 0xff, x, y, e2, dz);
+        psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n);
+    }
+
+    // test dump of the psf parameters
     if (psTraceGetLevel("psModules.objects") >= 4) {
         FILE *f = fopen ("pol.dat", "w");
+        fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
         for (int i = 0; i < e0->n; i++) {
-            fprintf (f, "%f %f  :  %f %f %f  : %d\n",
+            fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
                      x->data.F32[i], y->data.F32[i],
-                     e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], psfTry->mask->data.U8[i]);
+                     e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], 
+                     pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
+                     pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]),
+                     pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]),
+                     psfTry->mask->data.U8[i]);
         }
         fclose (f);
@@ -455 +461 @@
         // fit the collection of measured parameters to the PSF 2D model
         // the mask is carried from previous steps and updated with this operation
-        // the weight is either the flux error or NULL, depending on 'applyWeights'
-        if (!psVectorClipFitPolynomial2D(psf->params->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) {
+        // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams'
+        if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, e0, NULL)) {
             psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
-            psFree(stats);
             psFree(x);
             psFree(y);
@@ -485 +490 @@
 
             for (int i = 0; i < psf->params->n; i++) {
-                if (psf->params->data[i] == NULL)
-                    continue;
+                if (psf->params->data[i] == NULL) continue;
                 fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]);
             }
@@ -496 +500 @@
     }
 
-    psFree (stats);
     psFree (x);
     psFree (y);
     psFree (z);
     psFree (dz);
-    return (true);
+    return true;
 }
 

trunk/psModules/src/objects/pmPSFtry.h

@@ -6 +6 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-08-15 20:21:18 $
+ * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-09-24 21:27:58 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -78 +78 @@
  *
  */
-pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors, psPolynomial2D *psfTrendMask);
+pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options);
 
 /** pmPSFtryModel()
@@ -87 +87 @@
  *
  */
-pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask, bool applyWeights, psMaskType maskVal, psMaskType mark);
+pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, pmPSFOptions *options, psMaskType maskVal, psMaskType mark);
 
 /** pmPSFtryMetric()
@@ -123 +123 @@
  *
  */
-bool pmPSFFromPSFtry (pmPSFtry *psfTry, bool applyWeights);
+bool pmPSFFromPSFtry (pmPSFtry *psfTry);
 
 /// @}

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -3 +3 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-21 02:46:46 $
+ *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-24 21:27:58 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -81 +81 @@
     // we must have a valid model
     model = pmSourceGetModel (&isPSF, source);
-    if (model == NULL)
-        return false;
+    if (model == NULL) {
+        psTrace ("psModules.objects", 3, "fail mag : no valid model");
+        return false;
+    }
 
     if (model->dparams->data.F32[PM_PAR_I0] > 0) {
@@ -102 +104 @@
             source->psfMag = NAN;
         }
-        fprintf (stderr, ".");
     } else {
         status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
@@ -109 +110 @@
     // measure EXT model photometry
     status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+
     // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
-    if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf) {
-        // convert to the equivalent 2D model?
+    if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) {
         source->psfMag += pmTrend2DEval (psf->ApTrend, x, y);
     }
 
-    if (!isfinite(SN) || (SN < AP_MIN_SN))
-        return false;
+    if (!isfinite(SN) || (SN < AP_MIN_SN)) {
+        psTrace ("psModules.objects", 3, "fail mag : bad SN: %f (limit: %f)", SN, AP_MIN_SN);
+        return false;
+    }
 
     // replace source flux
@@ -209 +212 @@
     status = pmSourcePhotometryAper  (&source->apMag, model, flux, source->maskObj, maskVal);
     if (!status) {
+        psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
         psErrorCode last = psErrorCodeLast();
         if (last == PM_ERR_PHOTOM) {

trunk/psModules/src/objects/pmTrend2D.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-21 00:09:18 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-24 21:27:58 $
  *
  *  Copyright 2004 Institute for Astronomy, University of Hawaii
@@ -14 +14 @@
 #endif
 
+# include <strings.h>
 # include <pslib.h>
 # include "pmTrend2D.h"
@@ -31 +32 @@
 {
     assert (image);
-    assert (stats);
 
     pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
@@ -44 +44 @@
       case PM_TREND_POLY_ORD:
         trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend);
+        // set masking somehow
+        for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
+            for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
+                if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
+                    trend->poly->mask[nx][ny] = 1;
+                } else {
+                    trend->poly->mask[nx][ny] = 0;
+                }
+            }
+        }
         break;
 
@@ -69 +79 @@
 }
 
-pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats)
+pmTrend2D *pmTrend2DNoImageAlloc (pmTrend2DMode mode, psImageBinning *binning, psStats *stats)
 {
     pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
@@ -81 +91 @@
     switch (mode) {
       case PM_TREND_POLY_ORD:
+        trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, binning->nXruff, binning->nYruff);
+        // set masking somehow
+        for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
+            for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
+                if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
+                    trend->poly->mask[nx][ny] = 1;
+                } else {
+                    trend->poly->mask[nx][ny] = 0;
+                }
+            }
+        }
+        break;
+
+      case PM_TREND_POLY_CHEB:
+        trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_CHEB, binning->nXruff, binning->nYruff);
+        break;
+
+      case PM_TREND_MAP: {
+          // binning defines the map scale relationship
+          trend->map = psImageMapNoImageAlloc (binning, stats);
+          break;
+      }
+
+      default:
+        psAbort ("error");
+    }
+    return (trend);
+}
+
+pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats)
+{
+    pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
+    psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree);
+
+    trend->map = NULL;
+    trend->poly = NULL;
+    trend->stats = psMemIncrRefCounter (stats);
+    trend->mode = mode;
+        
+    switch (mode) {
+      case PM_TREND_POLY_ORD:
         trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend);
+        // set masking somehow
+        for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
+            for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
+                if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
+                    trend->poly->mask[nx][ny] = 1;
+                } else {
+                    trend->poly->mask[nx][ny] = 0;
+                }
+            }
+        }
         break;
 
@@ -111 +172 @@
     bool status;
 
+    assert (trend);
+    assert (x);
+    assert (y);
+    assert (f);
+
     switch (trend->mode) {
       case PM_TREND_POLY_ORD:
@@ -136 +202 @@
     double result;
 
-    if (!trend) return 0.0;
+    assert (trend);
 
     switch (trend->mode) {
@@ -158 +224 @@
     psVector *result;
 
+    assert (trend);
+    assert (x);
+    assert (y);
+
     switch (trend->mode) {
       case PM_TREND_POLY_ORD:
@@ -173 +243 @@
     return result;
 }
+
+psString pmTrend2DModeToString (pmTrend2DMode mode) {
+    
+    psString name;
+
+    switch (mode) {
+      case PM_TREND_NONE:
+        name = psStringCopy ("NONE");
+        break;
+      case PM_TREND_POLY_ORD:
+        name = psStringCopy ("POLY_ORD");
+        break;
+      case PM_TREND_POLY_CHEB:
+        name = psStringCopy ("POLY_CHEB");
+        break;
+      case PM_TREND_MAP:
+        name = psStringCopy ("MAP");
+        break;
+      default:
+        psAbort ("invalid mode %d", mode);
+    }
+    return name;
+}
+
+pmTrend2DMode pmTrend2DModeFromString (psString name) {
+
+    if (!name) return PM_TREND_NONE;
+
+    if (!strcasecmp (name, "NONE")) {
+        return PM_TREND_NONE;
+    }
+    if (!strcasecmp (name, "POLY_ORD")) {
+        return PM_TREND_POLY_ORD;
+    }
+    if (!strcasecmp (name, "POLY_CHEB")) {
+        return PM_TREND_POLY_CHEB;
+    }
+    if (!strcasecmp (name, "MAP")) {
+        return PM_TREND_MAP;
+    }
+    psError (PS_ERR_UNKNOWN, true, "Unknown pmTrend2D mode %s\n", name);
+    return PM_TREND_NONE;
+}

trunk/psModules/src/objects/pmTrend2D.h

@@ -5 +5 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-09-21 00:09:18 $
+ * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-09-24 21:27:58 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -17 +17 @@
 
 typedef enum {
+    PM_TREND_NONE,
     PM_TREND_POLY_ORD,
     PM_TREND_POLY_CHEB,
@@ -33 +34 @@
 pmTrend2D *pmTrend2DAlloc (pmTrend2DMode mode, psImage *image, int nXtrend, int nYtrend, psStats *stats);
 
+pmTrend2D *pmTrend2DNoImageAlloc (pmTrend2DMode mode, psImageBinning *binning, psStats *stats);
+
 // allocate a pmTrend2D tied to an abstract field with size nXfield,nYfield 
 pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats);
@@ -41 +44 @@
 psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y);
 
+psString pmTrend2DModeToString (pmTrend2DMode mode);
+pmTrend2DMode pmTrend2DModeFromString (psString name);
+
 /// @}
 # endif

trunk/psphot/doc/psfmodel.txt

@@ +1 @@
+
+2007.09.21
+
+  there are three places where we can choose to use errors in the fits or not:
+
+  * non-linear fitting of the models to the pixel flux distribution (poissonErrorsPhotLMM)
+  * linear fitting of the models to the pixel flux distribution (poissonErrorsPhotLin)
+  * fitting of the 2D variations in the psf parameters (poissonErrorsParams)
+  * fitting of the 2D variations in the aperture residuals
+
+2007.09.20
+
+  I am upgrading the PSF model to allow the parameter variation to be
+  modeled with pmTrend2D instead of just polynomials.  I am making a
+  list of places to modify the code:
+
+pmPSFAlloc : need a method beyong psfTrendMask to carry in the psf
+options
+
+pmPSF_ModelToFit : no need to change these
+
+update pmPSFBuildSimple to set the parameters of the pmTrend, which
+ever is used.
+
+pmPSFtry.c: some significant re-work!
+
+
+
+pmPSF_IO : need new functions to save / load the trend (psImages)
 
 2006.10.27

trunk/psphot/src/models/pmModel_STRAIL.c

@@ -530 +530 @@
  
     for (int i = 4; i < 7; i++) { 
-      psPolynomial2D *poly = psf->params->data[i-4]; 
-        out[i] = psPolynomial2DEval (poly, out[2], out[3]); 
+      pmTrend2D *trend = psf->params->data[i-4]; 
+        out[i] = pmTrend2DEval (trend, out[2], out[3]); 
     } 
     return(true); 
@@ -554 +554 @@
     for (int i = 0; i < psf->params->n; i++) {
         if (i == PM_PAR_SKY) continue;
-        psPolynomial2D *poly = psf->params->data[i];
-        assert (poly);
-        PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+        pmTrend2D *trend = psf->params->data[i];
+        PAR[i] = pmTrend2DEval(trend, Xo, Yo);
     }
 

trunk/psphot/src/models/pmModel_TEST1.c

@@ -217 +217 @@
             out[i] = in[i];
         } else {            
-            psPolynomial2D *poly = psf->params->data[i];
-            out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
+            pmTrend2D *trend = psf->params->data[i];
+            out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
     }
@@ -246 +246 @@
     for (int i = 0; i < psf->params->n; i++) {
         if (i == PM_PAR_SKY) continue;
-        psPolynomial2D *poly = psf->params->data[i];
-        assert (poly);
-        PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+        pmTrend2D *trend = psf->params->data[i];
+        PAR[i] = pmTrend2DEval(trend, Xo, Yo);
     }
 

trunk/psphot/src/psphot.h

@@ -47 +47 @@
 bool            psphotReplaceAll (psArray *sources, psMaskType maskVal);
 bool            psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark);
-bool            psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, int nGroup);
+bool            psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup);
 bool            psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag);
 bool            psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background, psMaskType maskVal, psMaskType mark);

trunk/psphot/src/psphotApResid.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+# define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; }
 // measure the aperture residual statistics and 2D variations
 bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark)
@@ -63 +64 @@
         model = source->modelPSF;
 
-        if (source->type != PM_SOURCE_TYPE_STAR) continue;
-        if (source->mode &  PM_SOURCE_MODE_SATSTAR) continue;
-        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
-        if (source->mode &  PM_SOURCE_MODE_FAIL) continue;
-        if (source->mode &  PM_SOURCE_MODE_POOR) continue;
+        if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR");
+        if (source->mode &  PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR");
+        if (source->mode &  PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND");
+        if (source->mode &  PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR");
+        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
 
         // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures
@@ -73 +74 @@
         if (!pmSourceMagnitudes (source, psf, photMode, maskVal, mark)) {
             Nskip ++;
+            psTrace ("psphot", 3, "skip : bad source mag");
             continue;
         }
 
-        if (!isfinite(source->apMag) || !isfinite(source->apMag)) {
+        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
             Nfail ++;
+            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
             continue;
         }
@@ -87 +90 @@
         if ((MAX_AP_OFFSET > 0) && (fabs(dap) > MAX_AP_OFFSET)) {
             Nfail ++;
+            psTrace ("psphot", 3, "fail : bad dap %f %f", dap, MAX_AP_OFFSET);
             continue;
         }
@@ -99 +103 @@
         dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
 
+        psVectorExtend (mag,     100, 1);
         psVectorExtend (mask,    100, 1);
         psVectorExtend (xPos,    100, 1);
@@ -146 +151 @@
     }
 
+    // XXX catch error condition 
     psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag);
 
@@ -152 +158 @@
     psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
     psVector *dMagSys = (psVector *) psBinaryOp (NULL, (void *) dMag, "*", (void *) psScalarAlloc(errorScale, PS_TYPE_F32));
-
-    psphotSaveImage (NULL, psf->ApTrend->map->map, "apresid.map.fits");
 
     if (psTraceGetLevel("psphot") >= 5) {
@@ -170 +174 @@
     float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
     float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
+
     psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
     psf->dApResid = errorFloor;
@@ -206 +211 @@
 */
 
-bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, int nGroup) {
+bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) {
 
     psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
@@ -225 +230 @@
     psVector *dMSubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
     psVector *dASubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
+    psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_U8);
 
     int n = 0;
@@ -233 +239 @@
             dMSubset->data.F32[j] = dMag->data.F32[N];
             dASubset->data.F32[j] = dap->data.F32[N];
+            mkSubset->data.U8[j]  = mask->data.U8[N];
         }
         dMSubset->n = j;
         dASubset->n = j;
+        mkSubset->n = j;
 
         psStatsInit (statsS);
         psStatsInit (statsM);
 
-        psVectorStats (statsS, dASubset, NULL, NULL, 0xff);
-        psVectorStats (statsM, dMSubset, NULL, NULL, 0xff);
+        psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff);
+        psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff);
 
         dSo->data.F32[i] = statsS->robustStdev;
@@ -250 +258 @@
     psFree (dMSubset);
     psFree (dASubset);
+    psFree (mkSubset);
     
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
@@ -286 +295 @@
     }
 
+    // the mask marks the values not used to calculate the ApTrend
+    psVectorInit (mask, 0);
+
     // XXX stats structure for use by ApTrend : make parameters user setable
     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
@@ -305 +317 @@
     // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
     int nGroup = PS_MAX (3*Nx*Ny, 10);
-    psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, nGroup);
+    psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, mask, nGroup);
 
     psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);

trunk/psphot/src/psphotChoosePSF.c

@@ -11 +11 @@
     sources = psArraySort (sources, pmSourceSortBySN);
 
+    // structure to store user options defining the psf
+    pmPSFOptions *options = pmPSFOptionsAlloc ();
+
+    // load user options from the recipe. no need to check existence -- they are 
     // array to store candidate PSF stars
     int NSTARS = psMetadataLookupS32 (&status, recipe, "PSF_MAX_NSTARS");
-    PS_ASSERT (status, NULL);
+    assert (status);
 
     float PSF_MIN_DS = psMetadataLookupF32 (&status, recipe, "PSF_MIN_DS");
-    PS_ASSERT (status, NULL);
+    assert (status);
 
     // supply the measured sky variance for optional constant errors (non-poissonian)
     float SKY_SIG = psMetadataLookupF32 (&status, recipe, "SKY_SIG");
-    PS_ASSERT (status, NULL);
+    assert (status);
 
     // use poissonian errors or local-sky errors
-    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
-    PS_ASSERT (status, NULL);
+    options->poissonErrorsPhotLMM = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PHOT.LMM");
+    assert (status);
 
     // use poissonian errors or local-sky errors
-    bool PSF_PARAM_WEIGHTS = psMetadataLookupBool (&status, recipe, "PSF_PARAM_WEIGHTS");
-    PS_ASSERT (status, NULL);
+    options->poissonErrorsPhotLin = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PHOT.LIN");
+    assert (status);
+
+    // use poissonian errors or local-sky errors
+    options->poissonErrorsParams = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PARAMS");
+    assert (status);
 
     // how to model the PSF variations across the field
-    psMetadata *md = psMetadataLookupMetadata (&status, recipe, "PSF.TREND.MASK");
-    PS_ASSERT (status, NULL);
+    options->psfTrendMode = pmTrend2DModeFromString (psMetadataLookupStr (&status, recipe, "PSF.TREND.MODE"));
+    assert (status);
+    
+    options->psfTrendNx = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NX");
+    assert (status);
+
+    options->psfTrendNy = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NY");
+    assert (status);
 
     // get the fixed PSF fit radius
     // XXX check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
-    float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
-    PS_ASSERT (status, NULL);
-
-    pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), POISSON_ERRORS);
-
-    psPolynomial2D *psfTrendMask = psPolynomial2DfromMetadata (md);
-    if (!psfTrendMask) {
-        psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe");
-        return NULL;
-    }
+    options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
+    assert (status);
+
+    options->stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+
+    // dimensions of the field for which the PSF is defined
+    options->psfFieldNx = readout->image->numCols;
+    options->psfFieldNy = readout->image->numRows;
+    options->psfFieldXo = readout->image->col0;
+    options->psfFieldYo = readout->image->row0;
+
+    pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), options->poissonErrorsPhotLMM);
 
     psArray *stars = psArrayAllocEmpty (sources->n);
@@ -66 +82 @@
         psLogMsg ("psphot.choosePSF", PS_LOG_WARN, "Failed to find any PSF candidates");
         psFree (stars);
-        psFree (psfTrendMask);
         return NULL;
     }
@@ -91 +106 @@
         psMetadataItem *item = psListGetAndIncrement (iter);
         char *modelName = item->data.V;
-        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS, maskVal, mark);
+        models->data[i] = pmPSFtryModel (stars, modelName, options, maskVal, mark);
     }
 
@@ -97 +112 @@
     psFree (list);
     psFree (stars);
-    psFree (psfTrendMask);
 
     // select the best of the models
@@ -128 +142 @@
     if (psTraceGetLevel("psphot") >= 5) {
         for (int i = PM_PAR_SXX; i < try->psf->params->n; i++) {
-            psPolynomial2D *poly = try->psf->params->data[i];
-            for (int nx = 0; nx <= poly->nX; nx++) {
-                for (int ny = 0; ny <= poly->nY; ny++) {
-                    if (poly->mask[nx][ny]) continue;
-                    fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny);
-                }
-            }
+            // XXX re-write the output or some other status info
         }
     }
@@ -236 +244 @@
             // use pmModelSub because modelFlux has not been generated
             assert (source->maskObj);
-            psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK);
+            psImageKeepCircle (source->maskObj, x, y, options->radius, "OR", PM_MASK_MARK);
             pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
-            psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
+            psImageKeepCircle (source->maskObj, x, y, options->radius, "AND", PS_NOT_U8(PM_MASK_MARK));
         }
 
@@ -283 +291 @@
         pmSourcesWritePSFs (try->sources, "psfstars.dat");
         pmSourcesWriteEXTs (try->sources, "extstars.dat", false);
-        psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
-        psMetadataConfigWrite (psfData, "psfmodel.dat");
-        psFree (psfData);
+        // XXX need alternative output function
+        // psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
+        // psMetadataConfigWrite (psfData, "psfmodel.dat");
         psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n");
         exit (0);
@@ -306 +314 @@
     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
 
+    psFree (options);
     return (psf);
 }

trunk/psphot/src/psphotRoughClass.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+# define CHECK_STATUS(S,MSG) { \
+  if (!status) { \
+    psError(PSPHOT_ERR_CONFIG, false, "missing PSF Clump entry: %s\n", MSG); \
+    return false; \
+  } }
 
 // 2006.02.02 : no leaks
 bool psphotRoughClass (psArray *sources, psMetadata *recipe, const bool havePSF, psMaskType maskSat) {
 
-    static pmPSFClump   psfClump;
+    bool status;
+    static pmPSFClump psfClump;
 
     psTimerStart ("psphot");
@@ -11 +17 @@
         // determine the PSF parameters from the source moment values
         psfClump = pmSourcePSFClump (sources, recipe);
+        psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.X", 0, "psf clump center", psfClump.X);
+        psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.Y", 0, "psf clump center", psfClump.Y);
+        psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DX", 0, "psf clump center", psfClump.dX);
+        psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DY", 0, "psf clump center", psfClump.dY);
     } else {
         // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
-        bool status_x, status_y, status_a;
-        float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X");
-        float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y");
-        float ANGLE  = psMetadataLookupF32 (&status_a, recipe, "ANGLE");
-        if (!status_x | !status_y | !status_a) {
-            psError(PSPHOT_ERR_CONFIG, false, "FWHM_X, FWHM_Y, or ANGLE not defined");
-            return false;
-        }
-
-        psEllipseAxes axes;
-        axes.major = FWHM_X / (2.0*sqrt(2.0*log(2.0)));
-        axes.minor = FWHM_Y / (2.0*sqrt(2.0*log(2.0)));
-        axes.theta = ANGLE;
-        psEllipseShape shape = psEllipseAxesToShape (axes);
-
-        psfClump.X   = shape.sx;
-        psfClump.Y   = shape.sy;
-        psfClump.dX  = 0.1*shape.sx;
-        psfClump.dY  = 0.1*shape.sy;
-        // dX,dY are somewhat crudely defined, but only used to select PSF candidates.
-        // if we already have a PSF, this is not actually used...
+        psfClump.X  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X");  CHECK_STATUS (status, "PSF.CLUMP.X");
+        psfClump.Y  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y");  CHECK_STATUS (status, "PSF.CLUMP.Y");
+        psfClump.dX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX");  CHECK_STATUS (status, "PSF.CLUMP.DX");
+        psfClump.dY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY");  CHECK_STATUS (status, "PSF.CLUMP.DY");
     }
 

trunk/psphot/src/psphotSourcePlots.c

@@ -20 +20 @@
     int dY = 0;                         // height of row so far
     int NX = 20*DX;                     // full width of output image
-    int NY = 0;                         // total height of output image
+    int NY = DY;                        // total height of output image so far
 
-    // first, examine the PSF and SAT stars:
-    // - determine bounding boxes for summary image
+    // first, examine the PSF and SAT stars to set output image size:
+    // - add stamp widths until we exceed output image width,
+    // - then start a new row offset by max height
     for (int i = 0; i < sources->n; i++) {
 
@@ -31 +32 @@
         keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
         keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
-        if (!keep) continue;
+        if (!keep) {
+            psTrace ("psphot", 4, "not plotting star %d: %x", i, source->mode);
+            continue;
+        }
 
         // how does this subimage get placed into the output image?
@@ -57 +61 @@
     }
 
+    if (NY == 0) {
+        psWarning ("no PSF or SAT stars to plot? skipping.\n");
+        return false;
+    }
+
     // allocate output image
     psImage *outpos = psImageAlloc (NX, NY, PS_TYPE_F32);
     psImage *outsub = psImageAlloc (NX, NY, PS_TYPE_F32);
+
+    psImageInit (outpos, 0.0);
+    psImageInit (outsub, 0.0);
 
     int Xo = 0;                         // starting corner of next box
@@ -142 +154 @@
     psFree (outpos);
     psFree (outsub);
-    return (sources);
+    return true;
 }