Changeset 13034

Timestamp:

Apr 25, 2007, 3:20:29 PM (19 years ago)

Author:

magnier

Message:

incorporating updates from eam_02_branch (cached models, pmPSF FITS IO)

Location:

trunk/psModules

Files:

• src/camera/pmFPACopy.c (modified) (2 diffs)
• src/camera/pmFPAfile.c (modified) (4 diffs)
• src/camera/pmFPAfileIO.c (modified) (6 diffs)
• src/concepts/pmConcepts.c (modified) (1 diff)
• src/concepts/pmConceptsUpdate.c (modified) (2 diffs)
• src/objects/pmModel.h (modified) (4 diffs)
• src/objects/pmPSF.c (modified) (3 diffs)
• src/objects/pmPSF_IO.c (modified) (4 diffs)
• src/objects/pmPSFtry.c (modified) (4 diffs)
• src/objects/pmPeaks.c (modified) (2 diffs)
• src/objects/pmSource.c (modified) (14 diffs)
• src/objects/pmSource.h (modified) (4 diffs)
• src/objects/pmSourceFitModel.c (modified) (5 diffs)
• src/objects/pmSourceIO.c (modified) (2 diffs)
• src/objects/pmSourceIO_CMF.c (modified) (2 diffs)
• src/objects/pmSourceIO_CMP.c (modified) (2 diffs)
• src/objects/pmSourceIO_OBJ.c (modified) (2 diffs)
• src/objects/pmSourceIO_PS1_DEV_0.c (modified) (2 diffs)
• src/objects/pmSourceIO_SMPDATA.c (modified) (2 diffs)
• src/objects/pmSourceIO_SX.c (modified) (2 diffs)
• src/objects/pmSourcePhotometry.c (modified) (10 diffs)
• src/objects/pmSourcePhotometry.h (modified) (2 diffs)
• src/objects/pmSourceSky.c (modified) (5 diffs)
• test/objects/tap_pmSourceFitModel.c (modified) (1 diff)
• test/objects/tap_pmSourceFitModel_Delta.c (modified) (1 diff)

trunk/psModules/src/camera/pmFPACopy.c

@@ -211 +211 @@
         if (mdok && trimsec && !psRegionIsNaN(*trimsec)) {
             *trimsec = psImageBinningSetRuffRegion(binning, *trimsec);
+            // force integer pixels : truncate x0, roundup x1:
+            trimsec->x0 = (int)trimsec->x0;
+            if (trimsec->x1 > (int)trimsec->x1) {
+                trimsec->x1 = (int)trimsec->x1 + 1;
+            } else {
+                trimsec->x1 = (int)trimsec->x1;
+            }           
+            trimsec->y0 = (int)trimsec->y0;
+            if (trimsec->y1 > (int)trimsec->y1) {
+                trimsec->y1 = (int)trimsec->y1 + 1;
+            } else {
+                trimsec->y1 = (int)trimsec->y1;
+            }           
         }
         psList *biassecs = psMetadataLookupPtr(&mdok, target->concepts, "CELL.BIASSEC"); // The bias sections
@@ -219 +232 @@
                 if (!psRegionIsNaN(*biassec)) {
                     *biassec = psImageBinningSetRuffRegion(binning, *biassec);
+                    // force integer pixels : truncate x0, roundup x1:
+                    biassec->x0 = (int)biassec->x0;
+                    if (biassec->x1 > (int)biassec->x1) {
+                        biassec->x1 = (int)biassec->x1 + 1;
+                    } else {
+                        biassec->x1 = (int)biassec->x1;
+                    }           
+                    biassec->y0 = (int)biassec->y0;
+                    if (biassec->y1 > (int)biassec->y1) {
+                        biassec->y1 = (int)biassec->y1 + 1;
+                    } else {
+                        biassec->y1 = (int)biassec->y1;
+                    }           
                 }
             }

trunk/psModules/src/camera/pmFPAfile.c

@@ -17 +17 @@
 #include "pmFPAfile.h"
 #include "pmFPACopy.h"
+#include "pmConcepts.h"
 
 static void pmFPAfileFree(pmFPAfile *file)
@@ -326 +327 @@
     PS_ASSERT_PTR_NON_NULL(view, false);
 
+    // XXX this should be smarter (ie, only copy concepts from the current chips)
+    // but such a call is needed, so re-copy stuff rather than no copy
+    pmFPACopyConcepts (out, in);
+
     // pmFPAWrite takes care of all PHUs as needed
-    if ((view->chip == -1) || in->hdu) {
+    if (view->chip == -1) {
         status = pmFPACopyStructure (out, in, xBin, yBin);
         return status;
@@ -338 +343 @@
     pmChip *outChip = out->chips->data[view->chip];
 
-    if ((view->cell == -1) || inChip->hdu) {
+    if (view->cell == -1) {
         status = pmChipCopyStructure (outChip, inChip, xBin, yBin);
         return status;
-    }
+    } 
     if (view->cell >= inChip->cells->n) {
         psError(PS_ERR_IO, true, "Requested cell == %d>= inChip->cells->n == %ld",
@@ -350 +355 @@
     pmCell *outCell = outChip->cells->data[view->cell];
 
-    if ((view->readout == -1) || inCell->hdu) {
-        status = pmCellCopyStructure (outCell, inCell, xBin, yBin);
-        return status;
-    }
-    psError(PS_ERR_UNKNOWN, true, "Returning false");
-    return false;
+    status = pmCellCopyStructure (outCell, inCell, xBin, yBin);
+    return status;
 }
 

trunk/psModules/src/camera/pmFPAfileIO.c

@@ -434 +434 @@
       case PM_FPA_FILE_FRINGE:
       case PM_FPA_FILE_CMF:
+      case PM_FPA_FILE_PSF:
         psTrace ("psModules.camera", 5, "closing %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout);
         status = psFitsClose (file->fits);
@@ -447 +448 @@
       case PM_FPA_FILE_OBJ:
       case PM_FPA_FILE_CMP:
-      case PM_FPA_FILE_PSF:
       case PM_FPA_FILE_JPEG:
       case PM_FPA_FILE_KAPA:
@@ -500 +500 @@
       case PM_FPA_FILE_CMP:
       case PM_FPA_FILE_CMF:
+      case PM_FPA_FILE_PSF:
         psTrace ("psModules.camera", 5, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name);
         return true;
-      case PM_FPA_FILE_PSF:
       case PM_FPA_FILE_JPEG:
       case PM_FPA_FILE_KAPA:
@@ -638 +638 @@
       case PM_FPA_FILE_FRINGE:
       case PM_FPA_FILE_CMF:
+      case PM_FPA_FILE_PSF:
         psTrace ("psModules.camera", 5, "opening %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout);
         file->fits = psFitsOpen (file->filename, mode);
@@ -672 +673 @@
       case PM_FPA_FILE_CMP:
       case PM_FPA_FILE_RAW:
-      case PM_FPA_FILE_PSF:
       case PM_FPA_FILE_JPEG:
       case PM_FPA_FILE_KAPA:
@@ -751 +751 @@
 
 // XXX this function is only called from pmFPAfileWrite
+// XXX for each data type, there should be a function which writes the PHU, if needed
 bool pmFPAfileWritePHU(pmFPAfile *file, const pmFPAview *view, const pmConfig *config)
 {

trunk/psModules/src/concepts/pmConcepts.c

@@ -888 +888 @@
 
 
+// XXX EAM : Paul, please handle the biassec independently so the 
+// psList is copied without a verbose warning.  when this is fixed, you can
+// convert the trace back to a log (psMetadata.c:429)
+
 bool pmFPACopyConcepts(pmFPA *target, const pmFPA *source)
 {

trunk/psModules/src/concepts/pmConceptsUpdate.c

@@ -44 +44 @@
             psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // Trim section
             *trimsec = psImageBinningSetRuffRegion (binning, *trimsec);
+            // force integer pixels : truncate x0, roundup x1:
+            trimsec->x0 = (int)trimsec->x0;
+            if (trimsec->x1 > (int)trimsec->x1) {
+                trimsec->x1 = (int)trimsec->x1 + 1;
+            } else {
+                trimsec->x1 = (int)trimsec->x1;
+            }           
+            trimsec->y0 = (int)trimsec->y0;
+            if (trimsec->y1 > (int)trimsec->y1) {
+                trimsec->y1 = (int)trimsec->y1 + 1;
+            } else {
+                trimsec->y1 = (int)trimsec->y1;
+            }           
             psMetadataRemoveKey(cell->concepts, "CELL.TRIMSEC.UPDATE");
         }
@@ -54 +67 @@
             while ((biassec = psListGetAndIncrement(biassecsIter))) {
                 *biassec = psImageBinningSetRuffRegion (binning, *biassec);
+                // force integer pixels : truncate x0, roundup x1:
+                biassec->x0 = (int)biassec->x0;
+                if (biassec->x1 > (int)biassec->x1) {
+                    biassec->x1 = (int)biassec->x1 + 1;
+                } else {
+                    biassec->x1 = (int)biassec->x1;
+                }               
+                biassec->y0 = (int)biassec->y0;
+                if (biassec->y1 > (int)biassec->y1) {
+                    biassec->y1 = (int)biassec->y1 + 1;
+                } else {
+                    biassec->y1 = (int)biassec->y1;
+                }               
             }
             psFree(biassecsIter);

trunk/psModules/src/objects/pmModel.h

@@ -5 +5 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-04-21 19:47:14 $
+ * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-04-26 01:20:29 $
+ *
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -26 +27 @@
     PM_MODEL_BADARGS   ///< model fit called with invalid args
 } pmModelStatus;
+
+typedef enum {
+    PM_MODEL_OP_NONE    = 0x00,
+    PM_MODEL_OP_FUNC    = 0x01,
+    PM_MODEL_OP_RES0    = 0x02,
+    PM_MODEL_OP_RES1    = 0x04,
+    PM_MODEL_OP_FULL   = 0x07,
+    PM_MODEL_OP_SKY     = 0x08,
+    PM_MODEL_OP_CENTER = 0x10,
+    PM_MODEL_OP_NORM    = 0x20,
+} pmModelOpMode;
 
 /** pmModel data structure
@@ -91 +103 @@
  */
 bool pmModelAdd(
-    psImage *image,   ///< The output image (float)
-    psImage *mask,   ///< The image pixel mask (valid == 0)
-    pmModel *model,   ///< The input pmModel
-    bool center,   ///< A boolean flag that determines whether pixels are centered
-    bool sky    ///< A boolean flag that determines if the sky is subtracted
+    psImage *image,                     ///< The output image (float)
+    psImage *mask,                      ///< The image pixel mask (valid == 0)
+    pmModel *model,                     ///< The input pmModel
+    pmModelOpMode mode               ///< mode to control how the model is added into the image
 );
-
 
 /** pmModelSub()
@@ -110 +120 @@
  */
 bool pmModelSub(
-    psImage *image,   ///< The output image (float)
-    psImage *mask,   ///< The image pixel mask (valid == 0)
-    pmModel *model,   ///< The input pmModel
-    bool center,   ///< A boolean flag that determines whether pixels are centered
-    bool sky    ///< A boolean flag that determines if the sky is subtracted
+    psImage *image,                     ///< The output image (float)
+    psImage *mask,                      ///< The image pixel mask (valid == 0)
+    pmModel *model,                     ///< The input pmModel
+    pmModelOpMode mode               ///< mode to control how the model is added into the image
 );
 

trunk/psModules/src/objects/pmPSF.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -424 +424 @@
     // place the reference object in the image center
     // no need to mask the source here
-    pmModelAdd (image, NULL, model, false, false);
+    // XXX should we measure this for the analytical model only or the full model?
+    pmModelAdd (image, NULL, model, PM_MODEL_OP_FULL);
 
     // loop over a range of source fluxes
@@ -436 +437 @@
         psImageKeepCircle (mask, xc, yc, radius, "OR", PM_MASK_MARK);
         pmSourcePhotometryAper (&apMag, model, image, mask);
+
+        // XXX since we re-mask on each pass, this could be dropped.
         psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_U8(PM_MASK_MARK));
 

trunk/psModules/src/objects/pmPSF_IO.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -44 +44 @@
 #include "pmSourceIO.h"
 
+bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    pmFPA *fpa = file->fpa;
+
+    if (view->chip == -1) {
+        if (!pmFPAWritePSFmodel (fpa, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for fpa");
+            return false;
+        }
+        return true;
+    }
+
+    if (view->chip >= fpa->chips->n) {
+        return false;
+    }
+    pmChip *chip = fpa->chips->data[view->chip];
+
+    if (view->cell == -1) {
+        if (!pmChipWritePSFmodel (chip, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+            return false;
+        }
+        return true;
+    }
+
+    if (view->cell >= chip->cells->n) {
+        return false;
+    }
+    pmCell *cell = chip->cells->data[view->cell];
+
+    if (view->readout == -1) {
+        if (!pmCellWritePSFmodel (cell, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for cell");
+            return false;
+        }
+        return true;
+    }
+
+    if (view->readout >= cell->readouts->n) {
+        return false;
+    }
+    pmReadout *readout = cell->readouts->data[view->readout];
+
+    if (!pmReadoutWritePSFmodel (readout, view, file, config)) {
+        psError(PS_ERR_IO, false, "Failed to write PSF for readout");
+        return false;
+    }
+    return true;
+}
+
+// read in all chip-level PSFmodel files for this FPA
+bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < fpa->chips->n; i++) {
+
+        pmChip *chip = fpa->chips->data[i];
+        if (!pmChipWritePSFmodel (chip, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for %dth chip", i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// read in all cell-level PSFmodel files for this chip
+bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < chip->cells->n; i++) {
+
+        pmCell *cell = chip->cells->data[i];
+        if (!pmCellWritePSFmodel (cell, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for %dth cell", i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// read in all readout-level PSFmodel files for this cell
+bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < cell->readouts->n; i++) {
+
+        pmReadout *readout = cell->readouts->data[i];
+        if (!pmReadoutWritePSFmodel (readout, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for %dth readout", i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, 
+// and header + image for the PSF residual images
+bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+    bool status;
+
+    // select the psf of interest
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    if (!psf) {
+        psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout");
+        return false;
+    }
+
+    // lookup the EXTNAME values used for table data and residual image segments
+    char *rule = NULL;
+
+    // Menu of EXTNAME rules
+    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
+    if (!menu) {
+        psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
+        return false;
+    }
+    // EXTNAME for table data
+    rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
+    if (!rule) {
+        psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
+        return false;
+    }
+    char *tableName = pmFPAfileNameFromRule (rule, file, view);
+
+    // write the PSF model parameters in a FITS table
+
+    // we need to write a header for the table, 
+    psMetadata *header = psMetadataAlloc();
+
+    char *modelName = pmModelGetType (psf->type);
+    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);
+
+    int nPar = pmModelParameterCount (psf->type)    ;
+    psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
+
+    // save the dimensions of each parameter
+    for (int i = 0; i < nPar; i++) {
+        char name[9];
+        psPolynomial2D *poly = psf->params_NEW->data[i];
+        if (poly == NULL) continue;
+        snprintf (name, 9, "PAR%02d_NX", i);
+        psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX);
+        snprintf (name, 9, "PAR%02d_NY", i);
+        psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY);
+    }
+
+    // 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);
+
+    // XXX can we drop this now?
+    psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
+
+    // build a FITS table of the PSF parameters
+    psArray *psfTable = psArrayAllocEmpty (100);
+    for (int i = 0; i < nPar; i++) {
+        psPolynomial2D *poly = psf->params_NEW->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);
+            }
+        }
+    }
+    
+    // write an empty FITS segment if we have no PSF information
+    if (psfTable->n == 0) {
+        // XXX this is probably an error (if we have a PSF, how do we have no data?)
+        psFitsWriteBlank (file->fits, header, tableName);
+    } else {
+        psTrace ("pmFPAfile", 5, "writing psf data %s\n", tableName);
+        if (!psFitsWriteTable (file->fits, header, psfTable, tableName)) {
+            psError(PS_ERR_IO, false, "writing psf table data %s\n", tableName);
+            psFree (tableName);
+            psFree (psfTable);
+            psFree (header);
+            return false;
+        }
+    }
+    psFree (tableName);
+    psFree (psfTable);
+    psFree (header);
+
+    // EXTNAME for residual images
+    rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
+    if (!rule) {
+        psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
+        return false;
+    }
+    char *imageName = pmFPAfileNameFromRule (rule, file, view);
+
+    // write the residual images (3D) 
+    header = psMetadataAlloc ();
+    if (psf->residuals == NULL) {
+        // set some header keywords to make it clear there are no residuals?
+        psFitsWriteBlank (file->fits, header, imageName);
+        psFree (imageName);
+        psFree (header);
+        return true;
+    }
+
+    psMetadataAddS32 (header, PS_LIST_TAIL, "XBIN",    0, "", psf->residuals->xBin);
+    psMetadataAddS32 (header, PS_LIST_TAIL, "YBIN",    0, "", psf->residuals->yBin);
+    psMetadataAddS32 (header, PS_LIST_TAIL, "XCENTER", 0, "", psf->residuals->xCenter);
+    psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter);
+
+    // write the residuals as three planes of the image
+    // this call creates an extension with NAXIS3 = 3
+    if (psf->residuals->Rx) {
+        // this call creates an extension with NAXIS3 = 3
+        psArray *images = psArrayAllocEmpty (3);
+        psArrayAdd (images, 1, psf->residuals->Ro);
+        psArrayAdd (images, 1, psf->residuals->Rx);
+        psArrayAdd (images, 1, psf->residuals->Ry);
+
+        psFitsWriteImageCube (file->fits, header, images, imageName);
+
+        // psFitsWriteImage(file->fits, header, psf->residuals->Ro, 3, imageName);
+        // psFitsUpdateImage(file->fits, psf->residuals->Rx, 0, 0, 1);
+        // psFitsUpdateImage(file->fits, psf->residuals->Ry, 0, 0, 2);
+    } else {
+        // this call creates an extension with NAXIS3 = 1
+        psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, imageName);
+    }
+  
+    psFree (imageName);
+    psFree (header);
+    return true;
+
+    // XXX save the growth curve
+    // XXX save the ApResid fit
+
+# if (0)
+    // build a FITS table of the fit to the Aperture Residuals
+    psArray *apresTable = psArrayAllocEmpty (100);
+    psPolynomial4D *poly = psf->ApTrend;
+    for (int ix = 0; ix < poly->nX; ix++) {
+        for (int iy = 0; iy < poly->nY; iy++) {
+
+            row = psMetadataAlloc ();
+            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);
+        }
+    }
+# endif
+}
+
+// XXX add in error handling
+bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    pmFPA *fpa = file->fpa;
+
+    if (view->chip == -1) {
+        pmFPAReadPSFmodel (fpa, view, file, config);
+        return true;
+    }
+
+    if (view->chip >= fpa->chips->n) {
+        return false;
+    }
+    pmChip *chip = fpa->chips->data[view->chip];
+
+    if (view->cell == -1) {
+        pmChipReadPSFmodel (chip, view, file, config);
+        return true;
+    }
+
+    if (view->cell >= chip->cells->n) {
+        return false;
+    }
+    pmCell *cell = chip->cells->data[view->cell];
+
+    if (view->readout == -1) {
+        pmCellReadPSFmodel (cell, view, file, config);
+        return true;
+    }
+
+    if (view->readout >= cell->readouts->n) {
+        return false;
+    }
+    pmReadout *readout = cell->readouts->data[view->readout];
+
+    pmReadoutReadPSFmodel (readout, view, file, config);
+    return true;
+}
+
+// read in all chip-level PSFmodel files for this FPA
+bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < fpa->chips->n; i++) {
+
+        pmChip *chip = fpa->chips->data[i];
+        pmChipReadPSFmodel (chip, view, file, config);
+    }
+    return true;
+}
+
+// read in all cell-level PSFmodel files for this chip
+bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < chip->cells->n; i++) {
+
+        pmCell *cell = chip->cells->data[i];
+        pmCellReadPSFmodel (cell, view, file, config);
+    }
+    return true;
+}
+
+// read in all readout-level PSFmodel files for this cell
+bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < cell->readouts->n; i++) {
+
+        pmReadout *readout = cell->readouts->data[i];
+        pmReadoutReadPSFmodel (readout, view, file, config);
+    }
+    return true;
+}
+
+// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, 
+// and header + image for the PSF residual images
+bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+    bool status; 
+    char *rule = NULL;
+    psMetadata *header = NULL;
+
+    // Menu of EXTNAME rules
+    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
+    if (!menu) {
+        psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
+        return false;
+    }
+    // EXTNAME for table data
+    rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
+    if (!rule) {
+        psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
+        return false;
+    }
+    char *tableName = pmFPAfileNameFromRule (rule, file, view);
+    // EXTNAME for residual images
+    rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
+    if (!rule) {
+        psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
+        return false;
+    }
+    char *imageName = pmFPAfileNameFromRule (rule, file, view);
+
+    // move fits pointer to table and read header
+    // advance to the table data extension
+    // since we have read the IMAGE header, the TABLE header should exist
+    if (!psFitsMoveExtName (file->fits, tableName)) {
+        psAbort("cannot find data extension %s in %s", tableName, file->filename);
+    }
+
+    // load the PSF model table header
+    header = psFitsReadHeader (NULL, file->fits);
+    if (!header) psAbort("cannot read table header");
+
+    // load the PSF model parameters from the FITS table
+    char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME");
+    pmModelType type = pmModelSetType (modelName);
+
+    bool poissonErrors = psMetadataLookupBool (&status, header, "POISSON");
+
+    // we determine the PSF parameter polynomials from the MD-defined polynomials
+    pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL);
+
+    // check the number of expected parameters
+    int nPar = psMetadataLookupS32 (&status, header, "PSF_NPAR");
+    if (nPar != pmModelParameterCount (psf->type))
+        psAbort("mismatch model par count");
+
+    // load the 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);
+        snprintf (name, 9, "PAR%02d_NY", i);
+        int nYorder = psMetadataLookupS32 (&status, header, name);
+        psf->params_NEW->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder);
+    }
+
+    // other required information describing the PSF
+    psf->ApResid   = psMetadataLookupF32 (&status, header, "AP_RESID");
+    psf->dApResid  = psMetadataLookupF32 (&status, header, "AP_ERROR");
+    psf->chisq     = psMetadataLookupF32 (&status, header, "CHISQ");
+    psf->nPSFstars = psMetadataLookupS32 (&status, header, "NSTARS");
+
+    // XXX can we drop this now?
+    psf->skyBias   = psMetadataLookupF32 (&status, header, "SKY_BIAS");
+
+    // read the raw table data
+    psArray *table = psFitsReadTable (file->fits);
+
+    // fill in the matching psf->params entries
+    for (int i = 0; i > table->n; i++) {
+        psMetadata *row = table->data[i];
+        int iPar = psMetadataLookupS32 (&status, row, "MODEL_TERM");
+        int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
+        int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
+        // XXX sanity check here
+
+        psPolynomial2D *poly = psf->params_NEW->data[iPar];
+
+        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);
+    psFree (table);
+
+    // move fits pointer to residual image and read header
+    // advance to the table data extension
+    // since we have read the IMAGE header, the TABLE header should exist
+    if (!psFitsMoveExtName (file->fits, imageName)) {
+        psAbort("cannot find data extension %s in %s", imageName, file->filename);
+    }
+
+    header = psFitsReadHeader (NULL, file->fits);
+    int Naxis = psMetadataLookupS32 (&status, header, "NAXIS");
+    if (Naxis != 0) {
+
+        int Nx = psMetadataLookupS32 (&status, header, "NAXIS1");
+        int Ny = psMetadataLookupS32 (&status, header, "NAXIS2");
+        int Nz = psMetadataLookupS32 (&status, header, "NAXIS3");
+
+        int xBin  = psMetadataLookupS32 (&status, header, "XBIN");
+        int yBin  = psMetadataLookupS32 (&status, header, "YBIN");
+
+        int xSize = Nx / xBin;
+        int ySize = Ny / yBin;
+
+        psf->residuals = pmResidualsAlloc (xSize, ySize, xBin, yBin);
+
+        psf->residuals->xCenter = psMetadataLookupS32 (&status, header, "XCENTER");
+        psf->residuals->yCenter = psMetadataLookupS32 (&status, header, "YCENTER");
+
+        psRegion fullImage = {0, 0, 0, 0};
+        psFitsReadImageBuffer(psf->residuals->Ro, file->fits, fullImage, 0); // Desired pixels
+        if (Nz > 1) {
+            assert (Nz == 3);
+            psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1); // Desired pixels
+            psFitsReadImageBuffer(psf->residuals->Ry, file->fits, fullImage, 2); // Desired pixels
+        }
+    }
+
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
+    psFree (psf);
+
+    psFree (tableName);
+    psFree (imageName);
+
+    return true;
+}
+
+/************ old support functions, deprecate? **************/
+
 psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf)
 {
@@ -124 +608 @@
 }
 
-bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    pmFPA *fpa = file->fpa;
-
-    if (view->chip == -1) {
-        pmFPAWritePSFmodel (fpa, view, file, config);
-        return true;
-    }
-
-    if (view->chip >= fpa->chips->n) {
-        return false;
-    }
-    pmChip *chip = fpa->chips->data[view->chip];
-
-    if (view->cell == -1) {
-        pmChipWritePSFmodel (chip, view, file, config);
-        return true;
-    }
-
-    if (view->cell >= chip->cells->n) {
-        return false;
-    }
-    pmCell *cell = chip->cells->data[view->cell];
-
-    if (view->readout == -1) {
-        pmCellWritePSFmodel (cell, view, file, config);
-        return true;
-    }
-
-    if (view->readout >= cell->readouts->n) {
-        return false;
-    }
-    pmReadout *readout = cell->readouts->data[view->readout];
-
-    pmReadoutWritePSFmodel (readout, view, file, config);
-    return true;
-}
-
-// read in all chip-level PSFmodel files for this FPA
-bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < fpa->chips->n; i++) {
-
-        pmChip *chip = fpa->chips->data[i];
-        pmChipWritePSFmodel (chip, view, file, config);
-    }
-    return true;
-}
-
-// read in all cell-level PSFmodel files for this chip
-bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < chip->cells->n; i++) {
-
-        pmCell *cell = chip->cells->data[i];
-        pmCellWritePSFmodel (cell, view, file, config);
-    }
-    return true;
-}
-
-// read in all readout-level PSFmodel files for this cell
-bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < cell->readouts->n; i++) {
-
-        pmReadout *readout = cell->readouts->data[i];
-        pmReadoutWritePSFmodel (readout, view, file, config);
-    }
-    return true;
-}
-
 // read in all readout-level Objects files for this cell
-bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmReadoutWritePSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
     bool status;
@@ -228 +637 @@
 }
 
-
-
-bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    pmFPA *fpa = file->fpa;
-
-    if (view->chip == -1) {
-        pmFPAReadPSFmodel (fpa, view, file, config);
-        return true;
-    }
-
-    if (view->chip >= fpa->chips->n) {
-        return false;
-    }
-    pmChip *chip = fpa->chips->data[view->chip];
-
-    if (view->cell == -1) {
-        pmChipReadPSFmodel (chip, view, file, config);
-        return true;
-    }
-
-    if (view->cell >= chip->cells->n) {
-        return false;
-    }
-    pmCell *cell = chip->cells->data[view->cell];
-
-    if (view->readout == -1) {
-        pmCellReadPSFmodel (cell, view, file, config);
-        return true;
-    }
-
-    if (view->readout >= cell->readouts->n) {
-        return false;
-    }
-    pmReadout *readout = cell->readouts->data[view->readout];
-
-    pmReadoutReadPSFmodel (readout, view, file, config);
-    return true;
-}
-
-// read in all chip-level PSFmodel files for this FPA
-bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < fpa->chips->n; i++) {
-
-        pmChip *chip = fpa->chips->data[i];
-        pmChipReadPSFmodel (chip, view, file, config);
-    }
-    return true;
-}
-
-// read in all cell-level PSFmodel files for this chip
-bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < chip->cells->n; i++) {
-
-        pmCell *cell = chip->cells->data[i];
-        pmCellReadPSFmodel (cell, view, file, config);
-    }
-    return true;
-}
-
-// read in all readout-level PSFmodel files for this cell
-bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < cell->readouts->n; i++) {
-
-        pmReadout *readout = cell->readouts->data[i];
-        pmReadoutReadPSFmodel (readout, view, file, config);
-    }
-    return true;
-}
-
 // read in all readout-level Objects files for this cell
-bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmReadoutReadPSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
 

trunk/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -99 +99 @@
 {
     bool status;
-    float x;
-    float y;
     int Next = 0;
     int Npsf = 0;
@@ -120 +118 @@
             return NULL;
         }
-        x = source->peak->x;
-        y = source->peak->y;
-
-        // set temporary object mask and fit object
+
+        // set object mask to define valid pixels
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", PM_MASK_MARK);
+
         // fit model as EXT, not PSF
-
-        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
         status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT);
-        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
 
         // exclude the poor fits
@@ -158 +153 @@
         source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
         source->modelPSF->radiusFit = RADIUS;
-        x = source->peak->x;
-        y = source->peak->y;
-
-        // set the mask and fit the PSF model
-        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
+
+        // set object mask to define valid pixels
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", PM_MASK_MARK);
+
+        // fit the PSF model to the source
         status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF);
-        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
 
         // skip poor fits

trunk/psModules/src/objects/pmPeaks.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-12 18:56:52 $
+ *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -134 +134 @@
     tmp->flux = 0;
     tmp->SN = 0;
+    tmp->xf = x;
+    tmp->yf = y;
 
     tmp->type = type;

trunk/psModules/src/objects/pmSource.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -39 +39 @@
     psFree(tmp->pixels);
     psFree(tmp->weight);
-    psFree(tmp->mask);
+    psFree(tmp->maskObj);
+    psFree(tmp->maskView);
+    psFree(tmp->modelFlux);
     psFree(tmp->moments);
     psFree(tmp->modelPSF);
@@ -47 +49 @@
 }
 
+// free only the pixel data associated with this source
 void pmSourceFreePixels(pmSource *source)
 {
@@ -55 +58 @@
     psFree (source->pixels);
     psFree (source->weight);
-    psFree (source->mask);
+    psFree (source->maskObj);
+    psFree (source->maskView);
+    psFree (source->modelFlux);
 
     source->pixels = NULL;
     source->weight = NULL;
-    source->mask = NULL;
+    source->maskObj = NULL;
+    source->maskView = NULL;
+    source->modelFlux = NULL;
     return;
 }
@@ -76 +83 @@
     source->pixels = NULL;
     source->weight = NULL;
-    source->mask = NULL;
+    source->maskObj = NULL;
+    source->maskView = NULL;
+    source->modelFlux = NULL;
     source->moments = NULL;
     source->blends = NULL;
@@ -114 +123 @@
     // the original values.
     pmSource *source = pmSourceAlloc ();
-    source->peak = psMemIncrRefCounter (in->peak);
-    source->pixels = psMemIncrRefCounter (in->pixels);
-    source->weight = psMemIncrRefCounter (in->weight);
-    source->mask = psMemIncrRefCounter (in->mask);
-    source->moments = psMemIncrRefCounter (in->moments);
-    source->blends = NULL;
-    source->modelPSF = NULL;
-    source->modelEXT = NULL;
+
+    // this is actually the same peak as the original, is this the correct way to handle this?
+    source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->value, in->peak->type);
+    source->peak->xf = in->peak->xf;
+    source->peak->yf = in->peak->yf;
+    source->peak->flux = in->peak->flux;
+    source->peak->SN = in->peak->SN;
+
+    // copy the values in the moments structure
+    source->moments  =  pmMomentsAlloc();
+    *source->moments = *in->moments;
+
+    // These images are all views to the parent.  
+    // We want a new view, but pointing at the same pixels.
+    source->pixels   = psImageCopyView (NULL, in->pixels);
+    source->weight   = psImageCopyView (NULL, in->weight);
+    source->maskView = psImageCopyView (NULL, in->maskView);
+
+    // the maskObj is a unique mask array; create a new mask image
+    if (in->maskObj) {
+        source->maskObj = psImageCopy (NULL, in->maskObj, PS_TYPE_MASK);
+    }
+
     source->type = in->type;
     source->mode = in->mode;
-    psMemSetDeallocator(source, (psFreeFunc) sourceFree);
-
-    // default values are NAN
-    source->psfMag = NAN;
-    source->extMag = NAN;
-    source->errMag = NAN;
-    source->apMag  = NAN;
-    source->sky    = NAN;
-    source->skyErr = NAN;
-    source->pixWeight = NAN;
 
     return(source);
@@ -151 +165 @@
     srcRegion = psRegionForImage (readout->image, srcRegion);
 
-    mySource->pixels = psMemIncrRefCounter(psImageSubset(readout->image, srcRegion));
-    mySource->weight = psMemIncrRefCounter(psImageSubset(readout->weight, srcRegion));
-    mySource->mask   = psMemIncrRefCounter(psImageSubset(readout->mask,  srcRegion));
-    mySource->region = srcRegion;
+    // these images are subset images of the equivalent parents
+    mySource->pixels   = psMemIncrRefCounter(psImageSubset(readout->image, srcRegion));
+    mySource->weight   = psMemIncrRefCounter(psImageSubset(readout->weight, srcRegion));
+    mySource->maskView = psMemIncrRefCounter(psImageSubset(readout->mask,  srcRegion));
+    mySource->region   = srcRegion;
+
+    // the object mask is a copy, and used to define the source pixels
+    mySource->maskObj = psImageCopy (NULL, mySource->maskView, PS_TYPE_MASK);
 
     return true;
@@ -183 +201 @@
     extend |= (mySource->pixels == NULL);
     extend |= (mySource->weight == NULL);
-    extend |= (mySource->mask == NULL);
+    extend |= (mySource->maskObj == NULL);
+    extend |= (mySource->maskView == NULL);
 
     // extend = true;
@@ -190 +209 @@
         psFree (mySource->pixels);
         psFree (mySource->weight);
-        psFree (mySource->mask);
-
-        mySource->pixels = psMemIncrRefCounter(psImageSubset(readout->image,  newRegion));
-        mySource->weight = psMemIncrRefCounter(psImageSubset(readout->weight, newRegion));
-        mySource->mask   = psMemIncrRefCounter(psImageSubset(readout->mask,   newRegion));
-        mySource->region = newRegion;
-    }
-
+        psFree (mySource->maskView);
+
+        mySource->pixels   = psMemIncrRefCounter(psImageSubset(readout->image,  newRegion));
+        mySource->weight   = psMemIncrRefCounter(psImageSubset(readout->weight, newRegion));
+        mySource->maskView = psMemIncrRefCounter(psImageSubset(readout->mask,   newRegion));
+        mySource->region   = newRegion;
+
+        // re-copy the main mask pixels.  NOTE: the user will need to reset the object mask
+        // pixels (eg, with psImageKeepCircle)
+        mySource->maskObj = psImageCopy (mySource->maskObj, mySource->maskView, PS_TYPE_MASK);
+        
+        // drop the old modelFlux pixels and force the user to re-create
+        psFree (mySource->modelFlux);
+        mySource->modelFlux = NULL;
+    }
     return extend;
 }
@@ -207 +233 @@
 
 // XXX EAM include a S/N cutoff in selecting the sources?
+// XXX this function should probably accept the values, not a recipe. wrap with a 
+// psphot-specific function which applies the recipe values
 pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *recipe)
 {
@@ -458 +486 @@
         // XXX EAM : can we use the value of SATURATE if mask is NULL?
         inner = psRegionForSquare (source->peak->x, source->peak->y, 2);
-        inner = psRegionForImage (source->mask, inner);
-        int Nsatpix = psImageCountPixelMask (source->mask, inner, PM_MASK_SAT);
+        inner = psRegionForImage (source->maskView, inner);
+        int Nsatpix = psImageCountPixelMask (source->maskView, inner, PM_MASK_SAT);
 
         // saturated star (size consistent with PSF or larger)
@@ -572 +600 @@
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
@@ -622 +650 @@
         psF32 *vPix = source->pixels->data.F32[row];
         psF32 *vWgt = source->weight->data.F32[row];
-        psU8  *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row];
+        psU8  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.U8[row];
 
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -719 +747 @@
 }
 
-pmModel *pmSourceSelectModel (pmSource *source)
-{
+// construct a realization of the source model 
+bool pmSourceCacheModel (pmSource *source) {
+
+    // select appropriate model
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+        
+    // if we already have a cached image, re-use that memory
+    source->modelFlux = psImageCopy (source->modelFlux, source->pixels, PS_TYPE_F32);
+    psImageInit (source->modelFlux, 0.0);
+
+    // in some places (psphotEnsemble), we need a normalized version
+    // in others, we just want the model.  which is more commonly used?
+    // modelFlux always has unity normalization (I0 = 1.0)
+    pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM);
+    return true;
+}
+
+// should we call pmSourceCacheModel if it does not exist?
+bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add) {
+
+    bool status;
+
+    if (add) {
+        psTrace ("psphot", 3, "replacing object at %f,%f\n", source->peak->xf, source->peak->yf);
+    } else {
+        psTrace ("psphot", 3, "removing object at %f,%f\n", source->peak->xf, source->peak->yf);
+    }
+
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+
+    if (source->modelFlux) {
+        // add in the pixels from the modelFlux image
+        int dX = source->modelFlux->col0 - source->pixels->col0;
+        int dY = source->modelFlux->row0 - source->pixels->row0;
+        assert (dX >= 0);
+        assert (dY >= 0);
+        assert (dX + source->modelFlux->numCols <= source->pixels->numCols);
+        assert (dY + source->modelFlux->numRows <= source->pixels->numRows);
+        
+        // modelFlux has unity normalization
+        float Io = model->params->data.F32[PM_PAR_I0];
+        if (mode & PM_MODEL_OP_NORM) {
+            Io = 1.0;
+        }
+
+        psU8 **mask = NULL;
+        if (source->maskObj) {
+            mask = source->maskObj->data.U8;
+        }
+
+        // XXX need to respect the source and model masks?
+        for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
+            int oy = iy + dY;
+            for (int ix = 0; ix < source->modelFlux->numCols; ix++) {
+                int ox = ix + dX;
+                if (mask && mask[iy][ix]) continue;
+                float value = Io*source->modelFlux->data.F32[iy][ix];
+                if (add) {
+                    source->pixels->data.F32[oy][ox] += value;
+                } else {
+                    source->pixels->data.F32[oy][ox] -= value;
+                }
+            }
+        }
+        return true;
+    } 
+
+    if (add) {
+        status = pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
+    } else {
+        status = pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
+    }
+
+    return true;
+}
+
+bool pmSourceAdd (pmSource *source, pmModelOpMode mode) {
+    bool status = pmSourceOp (source, mode, true);
+    return status;
+}
+
+bool pmSourceSub (pmSource *source, pmModelOpMode mode) {
+    bool status = pmSourceOp (source, mode, false);
+    return status;
+}
+
+// given a source, which model is currently appropriate?
+// choose PSF or EXT based on source->type, but fall back on PSF
+// if the EXT model is NULL
+pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source)
+{
+
+    pmModel *model;
+
+    if (isPSF) {
+      *isPSF = false;
+    }
+
     switch (source->type) {
     case PM_SOURCE_TYPE_STAR:
-        return source->modelPSF;
+        model = source->modelPSF;
+        if (model == NULL)
+            return NULL;
+        if (isPSF) {
+            *isPSF = true;
+        }
+        return model;
 
     case PM_SOURCE_TYPE_EXTENDED:
-        return source->modelEXT;
+        model = source->modelEXT;
+        if (!model && source->modelPSF) {
+            if (isPSF) {
+                *isPSF = true;
+            }
+            return source->modelPSF;
+        }
+        return (model);
+        break;
 
     default:

trunk/psModules/src/objects/pmSource.h

@@ -3 +3 @@
  * @author EAM, IfA; GLG, MHPCC
  *
- * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-01-24 02:54:15 $
+ * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-04-26 01:20:29 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -13 +13 @@
 /// @addtogroup Objects Object Detection / Analysis Functions
 /// @{
-
-/**
- * In the object analysis process, we will use specific mask values to mark the
- * image pixels. The following structure defines the relevant mask values.
- *
- * XXX: This is probably a bad solution: we will want to set mask values
- * outside of the PSPHOT code.  Perhaps we can set up a registered set of mask
- * values with specific meanings that other functions can add to or define?
- 
- * XXX We will only use the PM_MASK_xxx mask values
-typedef enum {
-    PM_SOURCE_MASK_CLEAR     = 0x00,
-    PM_SOURCE_MASK_INVALID   = 0x01,
-    PM_SOURCE_MASK_SATURATED = 0x02,
-    PM_SOURCE_MASK_MARKED    = 0x08,
-} psphotMaskValues;
-*/
 
 /** pmSourceType enumeration
@@ -75 +58 @@
 typedef struct
 {
-    const int id;   ///< Unique ID for object
-    pmPeak *peak;   ///< Description of peak pixel.
-    psImage *pixels;   ///< Rectangular region including object pixels.
-    psImage *weight;   ///< Image variance.
-    psImage *mask;   ///< Mask which marks pixels associated with objects.
-    pmMoments *moments;   ///< Basic moments measure for the object.
-    pmModel *modelPSF;   ///< PSF Model fit (parameters and type)
-    pmModel *modelEXT;   ///< EXT (floating) Model fit (parameters and type).
-    pmSourceType type;   ///< Best identification of object.
-    pmSourceMode mode;   ///< Best identification of object.
+    const int id;                       ///< Unique ID for object
+    pmPeak *peak;                       ///< Description of peak pixel.
+    psImage *pixels;                    ///< Rectangular region including object pixels.
+    psImage *weight;                    ///< Image variance.
+    psImage *modelFlux;                 ///< cached copy of the model for this source
+    psImage *maskObj;                   ///< unique mask for this object which marks included pixels associated with objects.
+    psImage *maskView;                  ///< view into global image mask for this object region
+    pmMoments *moments;                 ///< Basic moments measure for the object.
+    pmModel *modelPSF;                  ///< PSF Model fit (parameters and type)
+    pmModel *modelEXT;                  ///< EXT (floating) Model fit (parameters and type).
+    pmSourceType type;                  ///< Best identification of object.
+    pmSourceMode mode;                  ///< Best identification of object.
     psArray *blends;
-    float psfMag;   ///< calculated from flux in modelPsf
-    float extMag;   ///< calculated from flux in modelEXT
-    float errMag;   ///< error in psfMag OR extMag (depending on type)
-    float apMag;   ///< apMag corresponding to psfMag or extMag (depending on type)
-    float pixWeight;   // model-weighted coverage of valid pixels
-    psRegion region;   // area on image covered by selected pixels
-    float sky, skyErr;   //?< The sky and its error at the center of the object
+    float psfMag;                       ///< calculated from flux in modelPsf
+    float extMag;                       ///< calculated from flux in modelEXT
+    float errMag;                       ///< error in psfMag OR extMag (depending on type)
+    float apMag;               ///< apMag corresponding to psfMag or extMag (depending on type)
+    float pixWeight;                    // model-weighted coverage of valid pixels
+    psRegion region;                    // area on image covered by selected pixels
+    float sky, skyErr;                  //?< The sky and its error at the center of the object
 }
 pmSource;
@@ -225 +210 @@
 );
 
-
-// select the model used for this source
-pmModel *pmSourceSelectModel (pmSource *source);
+pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source);
+
+bool pmSourceAdd (pmSource *source, pmModelOpMode mode);
+bool pmSourceSub (pmSource *source, pmModelOpMode mode);
+
+bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add);
+bool pmSourceCacheModel (pmSource *source);
 
 /// @}

trunk/psModules/src/objects/pmSourceFitModel.c

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -57 +57 @@
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
     PS_ASSERT_PTR_NON_NULL(source->weight, false);
 
@@ -77 +77 @@
         for (psS32 j = 0; j < source->pixels->numCols; j++) {
             // skip masked points
-            if (source->mask->data.U8[i][j]) {
+            if (source->maskObj->data.U8[i][j]) {
                 continue;
             }
@@ -334 +334 @@
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
     PS_ASSERT_PTR_NON_NULL(source->weight, false);
 
@@ -354 +354 @@
         for (psS32 j = 0; j < source->pixels->numCols; j++) {
             // skip masked points
-            if (source->mask->data.U8[i][j]) {
+            if (source->maskObj->data.U8[i][j]) {
                 continue;
             }

trunk/psModules/src/objects/pmSourceIO.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -402 +402 @@
 
 // Given a FITS file pointer, read the table of object data
+// XXX add in error handling
 bool pmFPAviewReadObjects (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {

trunk/psModules/src/objects/pmSourceIO_CMF.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -56 +56 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+        // no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
         if (model == NULL)
             continue;

trunk/psModules/src/objects/pmSourceIO_CMP.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -107 +107 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+        // no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
         if (model == NULL)
             continue;

trunk/psModules/src/objects/pmSourceIO_OBJ.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -61 +61 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+        // no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
         if (model == NULL)
             continue;

trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -63 +63 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+        // no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
 
         if (model != NULL) {

trunk/psModules/src/objects/pmSourceIO_SMPDATA.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -62 +62 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+        // no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
         if (model != NULL) {
             PAR = model->params->data.F32;

trunk/psModules/src/objects/pmSourceIO_SX.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -57 +57 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+        // no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
         if (model == NULL)
             continue;

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -3 +3 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -106 +106 @@
     // replace source flux
     // XXX need to be certain which model and size of mask for prior subtraction
+    // XXX full model or just analytical?
+    // XXX use pmSourceAdd instead?
     if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
-        pmModelAdd (source->pixels, source->mask, model, false, false);
+        pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
     }
 
@@ -180 +182 @@
     // set aperture mask circle to model radius
     // XXX use a different radius for apertures and fits...
-    psImageKeepCircle (source->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
+    // XXX can I remove this?  the source should have the mask defined when it is constructed or
+    // when the fit / aperture radius is changed...
+    psImageKeepCircle (source->maskObj, x, y, model->radiusFit, "OR", PM_MASK_MARK);
 
     // measure the weight of included pixels
     // XXX is this supposed to use the weight or the flux?
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->mask);
+        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj);
     }
 
     // measure object aperture photometry
-    status = pmSourcePhotometryAper  (&source->apMag, model, flux, source->mask);
+    status = pmSourcePhotometryAper  (&source->apMag, model, flux, source->maskObj);
 
     // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
@@ -206 +210 @@
 
     // unmask aperture
-    psImageKeepCircle (source->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
+    // XXX can I remove this?  this will probably break things downstream...
+    psImageKeepCircle (source->maskObj, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
 
     // if source was originally subtracted, re-subtract object, leave local sky
+    // XXX replace with pmSourceSub...
     if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
-        pmModelSub (source->pixels, source->mask, model, false, false);
+        pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
     }
 
@@ -363 +369 @@
 }
 
+# if (0)
 double pmSourceCrossProduct (const pmSource *Mi,
                              const pmSource *Mj,
@@ -379 +386 @@
     const psImage *Wi = Mi->weight;
 
-    const psImage *Ti = Mi->mask;
-    const psImage *Tj = Mj->mask;
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
 
     Xs = PS_MAX (Pi->col0, Pj->col0);
@@ -435 +442 @@
     const psImage *Wi = Mi->weight;
 
-    const psImage *Ti = Mi->mask;
-    const psImage *Tj = Mj->mask;
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
 
     Xs = PS_MAX (Pi->col0, Pj->col0);
@@ -483 +490 @@
     const psImage *Pi = Mi->pixels;
     const psImage *Wi = Mi->weight;
-    const psImage *Ti = Mi->mask;
+    const psImage *Ti = Mi->maskObj;
 
     // note that this is addressing the same image pixels,
@@ -521 +528 @@
     return flux;
 }
+# endif
 
 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight)
@@ -543 +551 @@
 }
 
-// given a source, which model is currently appropriate?
-// choose PSF or EXT based on source->type, but fall back on PSF
-// if the EXT model is NULL
-pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source)
-{
-
-    pmModel *model;
-
-    switch (source->type) {
-    case PM_SOURCE_TYPE_STAR:
-        model = source->modelPSF;
-        if (model == NULL)
-            return NULL;
-        if (isPSF)
-            *isPSF = true;
-        return model;
-
-    case PM_SOURCE_TYPE_EXTENDED:
-        model = source->modelEXT;
-        if (isPSF)
-            *isPSF = false;
-        if (model == NULL) {
-            model = source->modelPSF;
-            if (isPSF)
-                *isPSF = true;
-        }
-        if (model == NULL) {
-            if (isPSF)
-                *isPSF = FALSE;
-            return NULL;
-        }
-        return (model);
-        break;
-
-    default:
-        if (isPSF)
-            *isPSF = false;
-        return NULL;
-    }
-    return NULL;
-}
+
+double pmSourceModelWeight(const pmSource *Mi,
+                      int term,
+                      const bool unweighted_sum) // should the cross product be weighted?
+{
+    double flux = 0, wt = 0, factor = 0;
+
+    const psImage *Pi = Mi->modelFlux;
+    const psImage *Wi = Mi->weight;
+    const psImage *Ti = Mi->maskObj;
+
+    for (int yi = 0; yi < Pi->numRows; yi++) {
+        for (int xi = 0; xi < Pi->numCols; xi++) {
+            if (Ti->data.U8[yi][xi])
+                continue;
+            if (!unweighted_sum) {
+                wt = Wi->data.F32[yi][xi];
+                if (wt == 0)
+                    continue;
+            }
+
+            switch (term) {
+            case 0:
+                factor = 1;
+                break;
+            case 1:
+                factor = xi + Pi->col0;
+                break;
+            case 2:
+                factor = yi + Pi->row0;
+                break;
+            default:
+                psAbort("invalid term for pmSourceWeight");
+            }
+
+            if (unweighted_sum) {
+                flux += (factor * Pi->data.F32[yi][xi]);
+            } else {
+                flux += (factor * Pi->data.F32[yi][xi]) / wt;
+            }
+        }
+    }
+    return flux;
+}
+
+double pmSourceModelDotModel (const pmSource *Mi,
+                              const pmSource *Mj,
+                              const bool unweighted_sum) // should the cross product be weighted?
+{
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    double flux, wt;
+
+    const psImage *Pi = Mi->modelFlux;
+    const psImage *Pj = Mj->modelFlux;
+
+    const psImage *Wi = Mi->weight;
+
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
+
+    Xs = PS_MAX (Pi->col0, Pj->col0);
+    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
+
+    Ys = PS_MAX (Pi->row0, Pj->row0);
+    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
+
+    xIs = Xs - Pi->col0;
+    xJs = Xs - Pj->col0;
+    yIs = Ys - Pi->row0;
+    yJs = Ys - Pj->row0;
+
+    xIe = Xe - Pi->col0;
+    yIe = Ye - Pi->row0;
+
+    // note that weight is addressing the same image pixels
+    flux = 0;
+    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
+        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
+            if (Ti->data.U8[yi][xi])
+                continue;
+            if (Tj->data.U8[yj][xj])
+                continue;
+
+            // XXX skip the nonsense weight pixels?
+            if (unweighted_sum) {
+                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
+            } else {
+                wt = Wi->data.F32[yi][xi];
+                if (wt > 0) {
+                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
+                }
+            }
+        }
+    }
+    return flux;
+}
+
+double pmSourceDataDotModel (const pmSource *Mi,
+                             const pmSource *Mj,
+                             const bool unweighted_sum) // should the cross product be weighted?
+{
+
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    double flux, wt;
+
+    const psImage *Pi = Mi->pixels;
+    const psImage *Pj = Mj->modelFlux;
+
+    const psImage *Wi = Mi->weight;
+
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
+
+    Xs = PS_MAX (Pi->col0, Pj->col0);
+    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
+
+    Ys = PS_MAX (Pi->row0, Pj->row0);
+    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
+
+    xIs = Xs - Pi->col0;
+    xJs = Xs - Pj->col0;
+    yIs = Ys - Pi->row0;
+    yJs = Ys - Pj->row0;
+
+    xIe = Xe - Pi->col0;
+    yIe = Ye - Pi->row0;
+
+    // note that weight is addressing the same image pixels,
+    flux = 0;
+    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
+        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
+            if (Ti->data.U8[yi][xi])
+                continue;
+            if (Tj->data.U8[yj][xj])
+                continue;
+
+            // XXX skip the nonsense weight pixels?
+            if (unweighted_sum) {
+                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
+            } else {
+                wt = Wi->data.F32[yi][xi];
+                if (wt > 0) {
+                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
+                }
+            }
+        }
+    }
+    return flux;
+}
+

trunk/psModules/src/objects/pmSourcePhotometry.h

@@ -4 +4 @@
  * @author EAM, IfA; GLG, MHPCC
  *
- * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-01-24 02:54:15 $
+ * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-04-26 01:20:29 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -50 +50 @@
 bool pmSourceMagnitudesInit (psMetadata *config);
 bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode);
-double pmSourceCrossProduct(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
-double pmSourceCrossWeight(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
 bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *image, psImage *mask);
 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight);
-pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source);
 
-double pmSourceWeight(const pmSource *Mi, int term, const bool unweighted_sum);
+
+double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
+double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
+double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum);
+
+// retire these:
+// double pmSourceCrossProduct(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
+// double pmSourceCrossWeight(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
+// double pmSourceWeight(const pmSource *Mi, int term, const bool unweighted_sum);
 
 /// @}

trunk/psModules/src/objects/pmSourceSky.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -49 +49 @@
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_IMAGE_NON_NULL(source->pixels, false);
-    PS_ASSERT_IMAGE_NON_NULL(source->mask, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->maskObj, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_INT_POSITIVE(Radius, false);
@@ -61 +61 @@
 
     psImage *image = source->pixels;
-    psImage *mask  = source->mask;
+    psImage *mask  = source->maskObj;
     pmPeak *peak  = source->peak;
     psRegion srcRegion;
@@ -100 +100 @@
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_IMAGE_NON_NULL(source->weight, false);
-    PS_ASSERT_IMAGE_NON_NULL(source->mask, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->maskObj, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_INT_POSITIVE(Radius, false);
@@ -112 +112 @@
 
     psImage *image = source->weight;
-    psImage *mask  = source->mask;
+    psImage *mask  = source->maskObj;
     pmPeak *peak  = source->peak;
     psRegion srcRegion;

trunk/psModules/test/objects/tap_pmSourceFitModel.c

@@ -136 +136 @@
 
     // create an image with the model, and add noise: gain is 1, subtracted sky is 100, readnoise is 5
-    pmModelAdd (source->pixels, source->mask, source->modelEXT, false, false);
+    pmModelAdd (source->pixels, source->mask, source->modelEXT, PM_MODEL_OP_FULL);
     int npix = 0;
     for (int j = 0; j < source->pixels->numRows; j++) {

trunk/psModules/test/objects/tap_pmSourceFitModel_Delta.c

@@ -54 +54 @@
 
     // create an image with the model, and add noise: gain is 1, subtracted sky is 100, readnoise is 5
-    pmModelAdd (source->pixels, source->mask, source->modelEXT, false, false);
+    pmModelAdd (source->pixels, source->mask, source->modelEXT, PM_MODEL_OP_FULL);
     int npix = 0;
     for (int j = 0; j < source->pixels->numRows; j++) {