Changeset 27840 for branches/simtest_nebulous_branches/psModules/src/objects/pmPSF_IO.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects/pmPSF_IO.c (modified) (35 diffs)

branches/simtest_nebulous_branches/psModules/src/objects/pmPSF_IO.c

@@ -28 +28 @@
 #include "pmConfig.h"
 #include "pmDetrendDB.h"
+#include "pmErrorCodes.h"
 
 #include "pmHDU.h"
@@ -49 +50 @@
 #include "pmSource.h"
 #include "pmModelClass.h"
+#include "pmModelUtils.h"
 #include "pmSourceIO.h"
 
@@ -121 +123 @@
     if (view->chip == -1) {
         if (!pmPSFmodelWriteFPA(fpa, view, file, config)) {
-            psError(PS_ERR_IO, false, "Failed to write PSF for fpa");
+            psError(psErrorCodeLast(), false, "Failed to write PSF for fpa");
             return false;
         }
@@ -134 +136 @@
     if (view->cell == -1) {
         if (!pmPSFmodelWriteChip (chip, view, file, config)) {
-            psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+            psError(psErrorCodeLast(), false, "Failed to write PSF for chip");
             return false;
         }
@@ -140 +142 @@
     }
 
-    psError(PS_ERR_IO, false, "PSF must be written at the chip level");
+    psError(PM_ERR_CONFIG, true, "PSF must be written at the chip level");
     return false;
 }
@@ -157 +159 @@
         thisView->chip = i;
         if (!pmPSFmodelWriteChip (chip, thisView, file, config)) {
-            psError(PS_ERR_IO, false, "Failed to write PSF for %dth chip", i);
+            psError(psErrorCodeLast(), false, "Failed to write PSF for %dth chip", i);
             psFree(thisView);
             return false;
@@ -172 +174 @@
     PS_ASSERT_PTR_NON_NULL(chip, false);
 
-    if (!pmPSFmodelWrite (chip->analysis, view, file, config)) {
-        psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+    if (!pmPSFmodelWrite(chip->analysis, view, file, config)) {
+        psError(psErrorCodeLast(), false, "Failed to write PSF for chip");
         return false;
     }
@@ -196 +198 @@
     char *headName, *tableName, *residName;
 
-    if (!analysis) return false;
+    if (!analysis) {
+        psError(PM_ERR_PROG, true, "No analysis metadata for chip.");
+        return false;
+    }
 
     // select the current recipe
     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT");
     if (!recipe) {
-        psError(PS_ERR_UNKNOWN, false, "missing recipe %s", "PSPHOT");
+        psError(PM_ERR_CONFIG, false, "missing recipe %s", "PSPHOT");
         return false;
     }
@@ -212 +217 @@
     // get the current header
     pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // Suitable FPA for writing
+    if (!fpa) {
+        psError(psErrorCodeLast(), false, "Unable to get FPA for writing.");
+        return false;
+    }
     pmHDU *hdu = psMemIncrRefCounter(pmFPAviewThisHDU(view, fpa));
     psFree(fpa);
     if (!hdu) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to find HDU");
+        psError(PM_ERR_CONFIG, false, "Unable to find HDU");
         return false;
     }
@@ -233 +242 @@
         psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
         if (!menu) {
-            psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
+            psError(PM_ERR_CONFIG, true, "missing EXTNAME.RULES in camera.config");
             psFree(hdu);
             return false;
@@ -241 +250 @@
         rule = psMetadataLookupStr(&status, menu, "PSF.HEAD");
         if (!rule) {
-            psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.HEAD in EXTNAME.RULES in camera.config");
+            psError(PM_ERR_CONFIG, false, "missing entry for PSF.HEAD in EXTNAME.RULES in camera.config");
             psFree(hdu);
             return false;
@@ -250 +259 @@
         rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
         if (!rule) {
-            psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
+            psError(PM_ERR_CONFIG, false, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
             psFree (headName);
             psFree(hdu);
@@ -260 +269 @@
         rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
         if (!rule) {
-            psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
+            psError(PM_ERR_CONFIG, false, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
             psFree (headName);
             psFree (tableName);
@@ -267 +276 @@
         }
         residName = pmFPAfileNameFromRule (rule, file, view);
+
+        // EXTNAME for psf image
+        // rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
+        // if (!rule) {
+        //     psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
+        //     psFree (headName);
+        //     psFree (tableName);
+        //     psFree(hdu);
+        //     return false;
+        // }
+        // residName = pmFPAfileNameFromRule (rule, file, view);
     }
 
@@ -282 +302 @@
         }
 
-        psFitsWriteBlank (file->fits, hdu->header, headName);
+        if (!psFitsWriteBlank(file->fits, hdu->header, headName)) {
+            psError(psErrorCodeLast(), false, "Unable to write PSF PHU.");
+            psFree(hdu);
+            return false;
+        }
         psTrace ("pmFPAfile", 5, "wrote ext head %s (type: %d)\n", file->filename, file->type);
         file->header = hdu->header;
@@ -292 +316 @@
     pmPSF *psf = psMetadataLookupPtr (&status, analysis, "PSPHOT.PSF");
     if (!psf) {
-        psError(PS_ERR_UNKNOWN, true, "missing PSF for this analysis metadata");
+        psError(PM_ERR_PROG, true, "missing PSF for this analysis metadata");
         psFree (tableName);
         psFree (residName);
@@ -310 +334 @@
         psMetadataAddBool (header, PS_LIST_TAIL, "ERR_PAR",  0, "Use Poisson errors in fits?", psf->poissonErrorsParams);
 
-        int nPar = pmModelClassParameterCount (psf->type)    ;
+        int nPar = pmModelClassParameterCount (psf->type);
         psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
 
@@ -321 +345 @@
         pmPSFClump psfClump;
 
-        // we now save clump parameters for each region : need to save all of those
-        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
-        psMetadataAddS32 (header, PS_LIST_TAIL, "PSF_CLN", PS_META_REPLACE, "number of psf clump regions", nRegions);
-        for (int i = 0; i < nRegions; i++) {
-            char regionName[64];
-            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
-            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
-
-            psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   assert (status);
-            psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   assert (status);
-            psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  assert (status);
-            psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  assert (status);
-
-            char key[16];
-            snprintf (key, 9, "CLX_%03d", i);
-            psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.X);
-            snprintf (key, 9, "CLY_%03d", i);
-            psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.Y);
-            snprintf (key, 9, "CLDX_%03d", i);
-            psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dX);
-            snprintf (key, 9, "CLDY_%03d", i);
-            psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dY);
-        }
+        // we now save clump parameters for each region : need to save all of those
+        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
+        psMetadataAddS32 (header, PS_LIST_TAIL, "PSF_CLN", PS_META_REPLACE, "number of psf clump regions", nRegions);
+        for (int i = 0; i < nRegions; i++) {
+            char regionName[64];
+            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
+            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
+
+            psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   assert (status);
+            psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   assert (status);
+            psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  assert (status);
+            psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  assert (status);
+
+            char key[16];
+            snprintf (key, 9, "CLX_%03d", i);
+            psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.X);
+            snprintf (key, 9, "CLY_%03d", i);
+            psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump center", psfClump.Y);
+            snprintf (key, 9, "CLDX_%03d", i);
+            psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dX);
+            snprintf (key, 9, "CLDY_%03d", i);
+            psMetadataAddF32 (header, PS_LIST_TAIL, key, PS_META_REPLACE, "psf clump size", psfClump.dY);
+        }
 
         // save the dimensions of each parameter
@@ -424 +448 @@
         // 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);
+            psError(PM_ERR_PROG, true, "No PSF data to write.");
+            psFree(tableName);
+            psFree(residName);
+            psFree(psfTable);
+            psFree(header);
+            return false;
         } 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);
+            if (!psFitsWriteTable(file->fits, header, psfTable, tableName)) {
+                psError(psErrorCodeLast(), false, "Error writing psf table data %s\n", tableName);
                 psFree (tableName);
                 psFree (residName);
@@ -447 +475 @@
         if (psf->residuals == NULL) {
             // set some header keywords to make it clear there are no residuals?
-            psFitsWriteBlank (file->fits, header, residName);
+            if (!psFitsWriteBlank(file->fits, header, residName)) {
+                psError(psErrorCodeLast(), false, "Unable to write blank PSF residual image.");
+                psFree(residName);
+                psFree(header);
+                return false;
+            }
             psFree (residName);
             psFree (header);
@@ -458 +491 @@
         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
+        // write the residuals as planes of the image
+        psArray *images = psArrayAllocEmpty (1);
+        psArrayAdd (images, 1, psf->residuals->Ro);  // z = 0 is Ro
+
         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, residName);
-            psFree (images);
-        } else {
-            // this call creates an extension with NAXIS3 = 1
-            psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, residName);
-        }
+        }
+
+        // note that all N plane are implicitly of the same type, so we convert the mask
+        if (psf->residuals->mask) {
+            psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
+            psArrayAdd (images, 1, mask);
+            psFree (mask);
+        }
+
+        // psFitsWriteImageCube (file->fits, header, images, residName);
+        // psFree (images);
+
+        if (!psFitsWriteImageCube (file->fits, header, images, residName)) {
+            psError(psErrorCodeLast(), false, "Unable to write PSF residuals.");
+            psFree(images);
+            psFree(residName);
+            psFree(header);
+            return false;
+        }
+        psFree (images);
         psFree (residName);
+        psFree (header);
+    }
+
+    // write a representation of the psf model
+    {
+        psMetadata *header = psMetadataAlloc ();
+
+        if (0) {
+            // set some header keywords to make it clear there are no residuals?
+            if (!psFitsWriteBlank (file->fits, header, residName)) {
+                psError(psErrorCodeLast(), false, "Unable to write blank PSF residuals.");
+                psFree(residName);
+                psFree(header);
+                return false;
+            }
+            psFree (residName);
+            psFree (header);
+            return true;
+        }
+
+        int DX = 65;
+        int DY = 65;
+
+        psImage *psfMosaic = psImageAlloc (DX, DY, PS_TYPE_F32);
+        psImageInit (psfMosaic, 0.0);
+
+        pmModel *modelRef = pmModelAlloc(psf->type);
+
+        // use the center of the center pixel of the image
+        float xc = 0.5*psf->fieldNx;
+        float yc = 0.5*psf->fieldNy;
+
+        // assign the x and y coords to the image center
+        // create an object with center intensity of 1000
+        modelRef->params->data.F32[PM_PAR_SKY] = 0;
+        modelRef->params->data.F32[PM_PAR_I0] = 1.000;
+        modelRef->params->data.F32[PM_PAR_XPOS] = xc;
+        modelRef->params->data.F32[PM_PAR_YPOS] = yc;
+
+        // create modelPSF from this model
+        pmModel *model = pmModelFromPSF (modelRef, psf);
+        if (model) {
+            // place the reference object in the image center
+            pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, 0.0, 0.0);
+            psFree (model);
+
+            if (false) {
+                // this call creates an extension with NAXIS3 = 3
+                psArray *images = psArrayAllocEmpty (3);
+                psArrayAdd (images, 1, psfMosaic);
+                // psArrayAdd (images, 1, psfModel);
+                // psArrayAdd (images, 1, psfModel);
+
+                if (!psFitsWriteImageCube (file->fits, header, images, "PSF_MODEL")) {
+                    psError(psErrorCodeLast(), false, "Unable to write PSF representation.");
+                    psFree(images);
+                    psFree(psfMosaic);
+                    psFree(modelRef);
+                    psFree(header);
+                    return false;
+                }
+                psFree (images);
+            } else {
+                // this call creates an extension with NAXIS3 = 1
+                // XXX need to replace PSF_MODEL with rule-based name like residName
+                if (!psFitsWriteImage(file->fits, header, psfMosaic, 0, "PSF_MODEL")) {
+                    psError(psErrorCodeLast(), false, "Unable to write PSF representation.");
+                    psFree(psfMosaic);
+                    psFree(modelRef);
+                    psFree(header);
+                    return false;
+                }
+            }
+        }
+
+        psFree (psfMosaic);
+        psFree (modelRef);
         psFree (header);
     }
@@ -523 +645 @@
     // find the FPA phu
     pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // Suitable FPA for writing
+    if (!fpa) {
+        psError(psErrorCodeLast(), false, "Unable to build FPA to write.");
+        return false;
+    }
     pmHDU *phu = psMemIncrRefCounter(pmFPAviewThisPHU(view, fpa));
     psFree(fpa);
@@ -533 +659 @@
         psMetadataCopy (outhead, phu->header);
     } else {
-        pmConfigConformHeader (outhead, file->format);
+        if (!pmConfigConformHeader (outhead, file->format)) {
+            psError(psErrorCodeLast(), false, "Unable to conform header of PSF PHU.");
+            psFree(phu);
+            return false;
+        }
     }
     psFree(phu);
 
     psMetadataAddBool (outhead, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true);
-    psFitsWriteBlank (file->fits, outhead, "");
+    if (!psFitsWriteBlank (file->fits, outhead, "")) {
+        psError(psErrorCodeLast(), false, "Unable to write PHU for PSF.");
+        psFree(outhead);
+        return false;
+    }
     file->wrote_phu = true;
 
@@ -568 +702 @@
     }
 
-    psError(PS_ERR_IO, false, "PSF must be read at the chip level");
+    psError(PM_ERR_CONFIG, true, "PSF must be read at the chip level");
     return false;
 }
@@ -595 +729 @@
 
     if (!pmPSFmodelRead (chip->analysis, view, file, config)) {
-        psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+        psError(psErrorCodeLast(), false, "Failed to write PSF for chip");
         return false;
     }
@@ -618 +752 @@
     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT");
     if (!recipe) {
-        psError(PS_ERR_UNKNOWN, false, "missing recipe %s", "PSPHOT");
+        psError(PM_ERR_CONFIG, false, "missing recipe %s", "PSPHOT");
         return false;
     }
@@ -625 +759 @@
     psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
     if (!menu) {
-        psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
+        psError(PM_ERR_CONFIG, true, "missing EXTNAME.RULES in camera.config");
         return false;
     }
@@ -631 +765 @@
     rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
     if (!rule) {
-        psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
+        psError(PM_ERR_CONFIG, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
         return false;
     }
@@ -638 +772 @@
     rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
     if (!rule) {
-        psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
+        psError(PM_ERR_CONFIG, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
         return false;
     }
@@ -646 +780 @@
     // 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);
+    if (!psFitsMoveExtName(file->fits, tableName)) {
+        psError(psErrorCodeLast(), false, "cannot find data extension %s in %s", tableName, file->filename);
+        return false;
     }
 
     // load the PSF model table header
     header = psFitsReadHeader (NULL, file->fits);
-    if (!header) psAbort("cannot read table header");
+    if (!header) {
+        psError(psErrorCodeLast(), false, "Cannot read PSF table header.");
+        return false;
+    }
 
     pmPSFOptions *options = pmPSFOptionsAlloc();
@@ -667 +805 @@
     int nRegions = psMetadataLookupS32 (&status, header, "PSF_CLN");
     if (!status) {
-        // read old-style psf clump data
-
-        char regionName[64];
-        snprintf (regionName, 64, "PSF.CLUMP.REGION.000");
-        psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
-
-        psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
-        if (!regionMD) {
-            regionMD = psMetadataAlloc();
-            psMetadataAddMetadata (recipe, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
-            psFree (regionMD);
-        }
-
-        // 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 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
-        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
-        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
-        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
+        // read old-style psf clump data
+
+        char regionName[64];
+        snprintf (regionName, 64, "PSF.CLUMP.REGION.000");
+        psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
+
+        psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
+        if (!regionMD) {
+            regionMD = psMetadataAlloc();
+            psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
+            psFree (regionMD);
+        }
+
+        // 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 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
+        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
     } else {
-        psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", nRegions);
-
-        for (int i = 0; i < nRegions; i++) {
-            char key[10];
-            char regionName[64];
-            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
-
-            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
-            if (!regionMD) {
-                regionMD = psMetadataAlloc();
-                psMetadataAddMetadata (recipe, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
-                psFree (regionMD);
-            }
-
-            // psf clump data
-            pmPSFClump psfClump;
-
-            snprintf (key, 9, "CLX_%03d", i);
-            psfClump.X  = psMetadataLookupF32 (&status, header, key);  assert(status);
-            snprintf (key, 9, "CLY_%03d", i);
-            psfClump.Y  = psMetadataLookupF32 (&status, header, key);  assert(status);
-            snprintf (key, 9, "CLDX_%03d", i);
-            psfClump.dX = psMetadataLookupF32 (&status, header, key);  assert(status);
-            snprintf (key, 9, "CLDY_%03d", i);
-            psfClump.dY = psMetadataLookupF32 (&status, header, key);  assert(status);
-
-            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
-            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
-            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
-            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
-        }
+        psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", nRegions);
+
+        for (int i = 0; i < nRegions; i++) {
+            char key[10];
+            char regionName[64];
+            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
+
+            psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName);
+            if (!regionMD) {
+                regionMD = psMetadataAlloc();
+                psMetadataAddMetadata (analysis, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD);
+                psFree (regionMD);
+            }
+
+            // psf clump data
+            pmPSFClump psfClump;
+
+            snprintf (key, 9, "CLX_%03d", i);
+            psfClump.X  = psMetadataLookupF32 (&status, header, key);  assert(status);
+            snprintf (key, 9, "CLY_%03d", i);
+            psfClump.Y  = psMetadataLookupF32 (&status, header, key);  assert(status);
+            snprintf (key, 9, "CLDX_%03d", i);
+            psfClump.dX = psMetadataLookupF32 (&status, header, key);  assert(status);
+            snprintf (key, 9, "CLDY_%03d", i);
+            psfClump.dY = psMetadataLookupF32 (&status, header, key);  assert(status);
+
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
+        }
     }
 
@@ -764 +902 @@
         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);
+            psError(PM_ERR_PROG, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i);
             return false;
         }
@@ -798 +936 @@
     // read the raw table data
     psArray *table = psFitsReadTable (file->fits);
+    if (!table) {
+        psError(psErrorCodeLast(), false, "Unable to read PSF table.");
+        psFree(header);
+        return false;
+    }
 
     // fill in the matching psf->params entries
@@ -842 +985 @@
     // 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);
+        psError(psErrorCodeLast(), false, "Cannot find PSF data extension %s in %s",
+                imageName, file->filename);
+        return false;
     }
 
     header = psFitsReadHeader (NULL, file->fits);
+    if (!header) {
+        psError(psErrorCodeLast(), false, "Unable to read PSF header.");
+        return false;
+    }
     int Naxis = psMetadataLookupS32 (&status, header, "NAXIS");
     if (Naxis != 0) {
@@ -865 +1014 @@
 
         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
-        }
-        // XXX notice that we are not saving the resid->mask
+        if (!psFitsReadImageBuffer(psf->residuals->Ro, file->fits, fullImage, 0)) {
+            psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+            return false;
+        }
+
+        // note that all N plane are implicitly of the same type, so we convert the mask
+        psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
+        psImageInit (psf->residuals->mask, 0);
+        psImageInit (psf->residuals->Rx, 0.0);
+        psImageInit (psf->residuals->Ry, 0.0);
+        switch (Nz) {
+          case 1: // Ro only
+            break;
+          case 2: // Ro and mask
+            if (!psFitsReadImageBuffer(mask, file->fits, fullImage, 1)) {
+                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+                return false;
+            }
+            psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
+            break;
+          case 3: // Ro, Rx and Ry, no mask
+            if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) {
+                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+                return false;
+            }
+            if (!psFitsReadImageBuffer(psf->residuals->Ry, file->fits, fullImage, 2)) {
+                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+                return false;
+            }
+            break;
+          case 4: // Ro, Rx, Ry, and mask:
+            if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) {
+                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+                return false;
+            }
+            if (!psFitsReadImageBuffer(psf->residuals->Ry, file->fits, fullImage, 2)) {
+                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+                return false;
+            }
+            if (!psFitsReadImageBuffer(mask, file->fits, fullImage, 3)) {
+                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+                return false;
+            }
+            psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
+            break;
+        }
+        psFree (mask);
     }