Changeset 15516

Timestamp:

Nov 8, 2007, 12:59:04 PM (19 years ago)

Author:

eugene

Message:

accepting updates from HEAD

Location:

branches/eam_branch_20071023/psModules/src

Files:

• camera/pmFPAConstruct.c (modified) (3 diffs)
• camera/pmFPAFlags.c (modified) (3 diffs)
• camera/pmFPAMosaic.c (modified) (8 diffs)
• camera/pmHDU.c (modified) (2 diffs)
• camera/pmHDUGenerate.c (modified) (10 diffs)
• camera/pmReadoutFake.c (modified) (2 diffs)
• concepts/pmConceptsAverage.h (modified) (1 diff)
• concepts/pmConceptsStandard.c (modified) (1 diff)
• config/pmConfig.c (modified) (1 diff)
• config/pmConfig.h (modified) (1 diff)
• config/pmConfigCamera.c (modified) (1 diff)
• config/pmConfigRecipes.c (modified) (23 diffs)
• imcombine/pmReadoutCombine.c (modified) (1 diff)
• imcombine/pmStack.c (modified) (1 diff)
• imcombine/pmStackReject.c (modified) (3 diffs)
• imcombine/pmSubtraction.h (modified) (7 diffs)
• imcombine/pmSubtractionKernels.c (modified) (4 diffs)
• imcombine/pmSubtractionKernels.h (modified) (1 diff)
• imcombine/pmSubtractionMatch.c (modified) (17 diffs)
• imcombine/pmSubtractionMatch.h (modified) (3 diffs)
• imcombine/pmSubtractionParams.c (modified) (19 diffs)
• imcombine/pmSubtractionParams.h (modified) (1 diff)
• imcombine/pmSubtractionStamps.c (modified) (21 diffs)
• imcombine/pmSubtractionStamps.h (modified) (5 diffs)
• objects/pmTrend2D.c (modified) (9 diffs)

branches/eam_branch_20071023/psModules/src/camera/pmFPAConstruct.c

@@ -45 +45 @@
     psMetadata *cellData = psMetadataLookupMetadata(&status, cells, cellName); // The data for the particular cell
     if (!status || !cellData) {
-        psLogMsg(__func__, PS_LOG_WARN, "Unable to find specs for cell %s: ignored\n", cellName);
+        psWarning("Unable to find specs for cell %s: ignored\n", cellName);
     }
 
@@ -304 +304 @@
         // Put in the cell data
         if (newCell->config) {
-            psLogMsg(__func__, PS_LOG_WARN, "Overwriting cell data in chip\n");
+            psWarning("Overwriting cell data in chip\n");
             psFree(newCell->config); // Make way!
         }
@@ -1304 +1304 @@
         const char *chipName = componentsItem->name; // Name of the chip
         if (componentsItem->type != PS_DATA_STRING) {
-            psLogMsg(__func__, PS_LOG_WARN, "Element %s in FPA within the camera configuration is not of "
+            psWarning("Element %s in FPA within the camera configuration is not of "
                      "type STR (type=%x) --- ignored.\n", chipName, componentsItem->type);
             continue;

branches/eam_branch_20071023/psModules/src/camera/pmFPAFlags.c

@@ -136 +136 @@
 
     for (int i = 0; i < fpa->chips->n; i++) {
-        pmChip *chip = fpa->chips->data[i];
-        if (chip == NULL) continue;
-        if (chip->data_exists) return true;
+        pmChip *chip = fpa->chips->data[i];
+        if (chip == NULL) continue;
+        if (chip->data_exists) return true;
     }
     return false;
@@ -170 +170 @@
 
     if (view->chip == -1) {
-        bool exists = pmFPACheckDataStatus (fpa);
-        return exists;
+        bool exists = pmFPACheckDataStatus (fpa);
+        return exists;
     }
 
     if (view->chip >= fpa->chips->n) {
-        psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n);
-        return false;
+        psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n);
+        return false;
     }
     pmChip *chip = fpa->chips->data[view->chip];
 
     if (view->cell == -1) {
-        bool exists = pmChipCheckDataStatus (chip);
-        return exists;
+        bool exists = pmChipCheckDataStatus (chip);
+        return exists;
     }
 
     if (view->cell >= chip->cells->n) {
-        psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n);
-        return false;
+        psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n);
+        return false;
     }
     pmCell *cell = chip->cells->data[view->cell];
 
     if (view->readout == -1) {
-        bool exists = pmCellCheckDataStatus (cell);
-        return exists;
+        bool exists = pmCellCheckDataStatus (cell);
+        return exists;
     }
 
     if (view->readout >= cell->readouts->n) {
-        psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouds->n == %ld", view->readout, cell->readouts->n);
-        return false;
+        psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouds->n == %ld", view->readout, cell->readouts->n);
+        return false;
     }
     pmReadout *readout = cell->readouts->data[view->readout];
@@ -290 +290 @@
     psArray *chips = fpa->chips;        // Component chips
     if (chips == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: fpa->chips == NULL\n");
+        psWarning("WARNING: fpa->chips == NULL\n");
         return(0);
     }
     if ((chipNum >= chips->n) || (NULL == (pmChip *) chips->data[chipNum])) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: the specified chip (%d) does not exist.\n", chipNum);
+        psWarning("WARNING: the specified chip (%d) does not exist.\n", chipNum);
         return(0);
     }

branches/eam_branch_20071023/psModules/src/camera/pmFPAMosaic.c

@@ -704 +704 @@
     psArray *readouts = cell->readouts; // The array of readouts
     if (readouts->n > 1) {
-        psLogMsg(__func__, PS_LOG_WARN, "Cell contains more than one readout (%ld) --- only the first will "
-                 "be mosaicked.\n", readouts->n);
+        psWarning("Cell contains more than one readout (%ld) --- only the first will be mosaicked.\n",
+                  readouts->n);
     }
     pmReadout *readout = readouts->data[0]; // The only readout we'll bother with
@@ -754 +754 @@
     int x0Target = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.X0");
     if (!mdok) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.X0 is not set for the target cell; assuming 0.\n");
+        psWarning("CELL.X0 is not set for the target cell; assuming 0.\n");
         FIX_CONCEPT(targetCell->concepts, "CELL.X0", S32, 0);
     }
     int y0Target = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.Y0");
     if (!mdok) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.Y0 is not set for the target cell; assuming 0.\n");
+        psWarning("CELL.Y0 is not set for the target cell; assuming 0.\n");
         FIX_CONCEPT(targetCell->concepts, "CELL.Y0", S32, 0);
     }
     int xParityTarget = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.XPARITY");
     if (!mdok || (xParityTarget != -1 && xParityTarget != 1)) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.XPARITY is not set for the target cell; assuming 1.\n");
+        psWarning("CELL.XPARITY is not set for the target cell; assuming 1.\n");
         FIX_CONCEPT(targetCell->concepts, "CELL.XPARITY", S32, 1);
         xParityTarget = 1;
@@ -770 +770 @@
     int yParityTarget = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.YPARITY");
     if (!mdok || (yParityTarget != -1 && yParityTarget != 1)) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.YPARITY is not set for the target cell; assuming 1.\n");
+        psWarning("CELL.YPARITY is not set for the target cell; assuming 1.\n");
         FIX_CONCEPT(targetCell->concepts, "CELL.YPARITY", S32, 1);
         yParityTarget = 1;
@@ -863 +863 @@
     int x0Target = psMetadataLookupS32(&mdok, targetChip->concepts, "CHIP.X0");
     if (!mdok) {
-        psLogMsg(__func__, PS_LOG_WARN, "CHIP.X0 is not set for the target chip; assuming 0.\n");
+        psWarning("CHIP.X0 is not set for the target chip; assuming 0.\n");
         FIX_CONCEPT(targetChip->concepts, "CHIP.X0", S32, 0);
     }
     int y0Target = psMetadataLookupS32(&mdok, targetChip->concepts, "CHIP.Y0");
     if (!mdok) {
-        psLogMsg(__func__, PS_LOG_WARN, "CHIP.Y0 is not set for the target chip; assuming 0.\n");
+        psWarning("CHIP.Y0 is not set for the target chip; assuming 0.\n");
         FIX_CONCEPT(targetChip->concepts, "CHIP.Y0", S32, 0);
     }
     x0Target += psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.X0");
     if (!mdok) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.X0 is not set for the target cell; assuming 0.\n");
+        psWarning("CELL.X0 is not set for the target cell; assuming 0.\n");
         FIX_CONCEPT(targetCell->concepts, "CELL.X0", S32, 0);
     }
     y0Target += psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.Y0");
     if (!mdok) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.Y0 is not set for the target cell; assuming 0.\n");
+        psWarning("CELL.Y0 is not set for the target cell; assuming 0.\n");
         FIX_CONCEPT(targetCell->concepts, "CELL.Y0", S32, 0);
     }
     int xParityChipTarget = psMetadataLookupS32(&mdok, targetChip->concepts, "CHIP.XPARITY");
     if (!mdok || (xParityChipTarget != -1 && xParityChipTarget != 1)) {
-        psLogMsg(__func__, PS_LOG_WARN, "CHIP.XPARITY is not set for the target chip; assuming 1.\n");
+        psWarning("CHIP.XPARITY is not set for the target chip; assuming 1.\n");
         FIX_CONCEPT(targetChip->concepts, "CHIP.XPARITY", S32, 1);
         xParityChipTarget = 1;
@@ -889 +889 @@
     int yParityChipTarget = psMetadataLookupS32(&mdok, targetChip->concepts, "CHIP.YPARITY");
     if (!mdok || (yParityChipTarget != -1 && yParityChipTarget != 1)) {
-        psLogMsg(__func__, PS_LOG_WARN, "CHIP.YPARITY is not set for the target chip; assuming 1.\n");
+        psWarning("CHIP.YPARITY is not set for the target chip; assuming 1.\n");
         FIX_CONCEPT(targetChip->concepts, "CHIP.YPARITY", S32, 1);
         yParityChipTarget = 1;
@@ -895 +895 @@
     int xParityCellTarget = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.XPARITY");
     if (!mdok || (xParityCellTarget != -1 && xParityCellTarget != 1)) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.XPARITY is not set for the target cell; assuming 1.\n");
+        psWarning("CELL.XPARITY is not set for the target cell; assuming 1.\n");
         FIX_CONCEPT(targetCell->concepts, "CELL.XPARITY", S32, 1);
         xParityCellTarget = 1;
@@ -901 +901 @@
     int yParityCellTarget = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.YPARITY");
     if (!mdok || (yParityCellTarget != -1 && yParityCellTarget != 1)) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.YPARITY is not set for the target cell; assuming 1.\n");
+        psWarning("CELL.YPARITY is not set for the target cell; assuming 1.\n");
         FIX_CONCEPT(targetCell->concepts, "CELL.YPARITY", S32, 1);
         yParityCellTarget = 1;
@@ -1334 +1334 @@
     pmHDU *sourceHDU = pmHDUGetHighest(source, firstSourceChip, firstSourceCell); // The HDU for the source
     if (!sourceHDU) {
-        psLogMsg(__func__, PS_LOG_WARN, "Unable to find HDU in source FPA; unable to copy headers.\n");
+        psWarning("Unable to find HDU in source FPA; unable to copy headers.\n");
         return false;
     }
     pmHDU *targetHDU = pmHDUGetHighest(target, targetChip, targetCell); // The HDU for the target
     if (!targetHDU) {
-        psLogMsg(__func__, PS_LOG_WARN, "Unable to find HDU in target FPA; unable to copy headers.\n");
+        psWarning("Unable to find HDU in target FPA; unable to copy headers.\n");
         return false;
     }

branches/eam_branch_20071023/psModules/src/camera/pmHDU.c

@@ -119 +119 @@
 
     if (*images) {
-        psLogMsg(__func__, PS_LOG_WARN, "HDU %s has already been read --- overwriting.\n", hdu->extname);
+        psWarning("HDU %s has already been read --- overwriting.\n", hdu->extname);
         psFree(*images);                // Blow away anything existing
     }
@@ -167 +167 @@
 
     if (!images && !hdu->header) {
-        psLogMsg(__func__, PS_LOG_WARN, "Nothing to write for HDU %s\n", hdu->extname);
+        psWarning("Nothing to write for HDU %s\n", hdu->extname);
         return false;
     }

branches/eam_branch_20071023/psModules/src/camera/pmHDUGenerate.c

@@ -134 +134 @@
     psMetadataItem *biassecItem = psMetadataLookup(cell->concepts, "CELL.BIASSEC"); // Bias sections
     if (!biassecItem) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.BIASSEC has not been initialised in cell --- "
-                 "ignored.\n");
+        psWarning("CELL.BIASSEC has not been initialised in cell --- ignored.\n");
         return false;
     }
@@ -146 +145 @@
     if (!mdok || (cellreaddir != 1 && cellreaddir != 2)) {
         // Probably unnecessary, but just in case....
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR is not set in cell --- ignored.\n");
+        psWarning("CELL.READDIR is not set in cell --- ignored.\n");
         return false;
     }
@@ -152 +151 @@
         *readdir = cellreaddir;
     } else if (*readdir != cellreaddir) {
-        psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR does not match read direction for HDU --- ignored.\n");
+        psWarning("CELL.READDIR does not match read direction for HDU --- ignored.\n");
         return false;
     }
@@ -255 +254 @@
         psMetadataItem *trimsecItem = psMetadataLookup(cell->concepts, "CELL.TRIMSEC"); // Item with trimsec
         if (!trimsecItem || trimsecItem->type != PS_DATA_REGION) {
-            psLogMsg(__func__, PS_LOG_WARN, "CELL.TRIMSEC has not been initialised in cell --- "
-                     "ignored.\n");
+            psWarning("CELL.TRIMSEC has not been initialised in cell --- ignored.\n");
             continue;
         }
@@ -264 +262 @@
         if (!mdok || (cellreaddir != 1 && cellreaddir != 2)) {
             // Probably unnecessary, but just in case....
-            psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR is not set in cell --- ignored.\n");
+            psWarning("CELL.READDIR is not set in cell --- ignored.\n");
             continue;
         }
@@ -270 +268 @@
             readdir = cellreaddir;
         } else if (readdir != cellreaddir) {
-            psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR for cells within the HDU do not match!\n");
+            psWarning("CELL.READDIR for cells within the HDU do not match!\n");
             cellreaddir = readdir;
         }
@@ -282 +280 @@
         if (readout->mask &&
                 (readout->mask->numCols != image->numCols || readout->mask->numRows != image->numRows)) {
-            psLogMsg(__func__, PS_LOG_WARN, "Image and mask have different sizes (%dx%d vs %dx%d)!\n",
+            psWarning("Image and mask have different sizes (%dx%d vs %dx%d)!\n",
                      image->numCols, image->numRows, readout->mask->numCols, readout->mask->numRows);
         }
         if (readout->weight &&
                 (readout->weight->numCols != image->numCols || readout->weight->numRows != image->numRows)) {
-            psLogMsg(__func__, PS_LOG_WARN, "Image and weight have different sizes (%dx%d vs %dx%d)!\n",
+            psWarning("Image and weight have different sizes (%dx%d vs %dx%d)!\n",
                      image->numCols, image->numRows, readout->weight->numCols, readout->weight->numRows);
         }
@@ -318 +316 @@
 
     if (previous != current) {
-        psLogMsg(__func__, PS_LOG_WARN, "Images within the HDU are of different types "
-                 "(%x vs %x) --- promoting\n", previous, current);
+        psWarning("Images within the HDU are of different types (%x vs %x) --- promoting\n",
+                  previous, current);
         return PS_MAX(previous, current);
     }
@@ -336 +334 @@
     if (source->numCols != region->x1 - region->x0 || source->numRows != region->y1 - region->y0) {
         psString regionString = psRegionToString(*region);
-        psLogMsg(__func__, PS_LOG_WARN, "Image size (%dx%d) does not match region (%s).\n",
+        psWarning("Image size (%dx%d) does not match region (%s).\n",
                  source->numCols, source->numRows, regionString);
         psFree(regionString);
@@ -475 +473 @@
 
                 if (biassecs->n != readout->bias->n) {
-                    psLogMsg(__func__, PS_LOG_WARN, "Number of bias sections (%ld) and number of biases (%ld)"
-                             " do not match.\n", biassecs->n, readout->bias->n);
+                    psWarning("Number of bias sections (%ld) and number of biases (%ld) do not match.\n",
+                              biassecs->n, readout->bias->n);
                 }
                 psListIterator *biasIter = psListIteratorAlloc(readout->bias, PS_LIST_HEAD, false); // Iteratr

branches/eam_branch_20071023/psModules/src/camera/pmReadoutFake.c

@@ -35 +35 @@
     pmReadout *readout = pmReadoutAlloc(NULL); // Output readout
     readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+    psImageInit(readout->image, 0);
 
     int numSources = sources->n;          // Number of stars
@@ -86 +87 @@
         fakeModel->params->data.F32[PM_PAR_I0] = powf(10.0, -0.4 * source->psfMag) / flux0;
 
+        psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n",
+                fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS],
+                fakeModel->params->data.F32[PM_PAR_I0]);
+
         pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
         fakeSource->peak = pmPeakAlloc(x, y, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);

branches/eam_branch_20071023/psModules/src/concepts/pmConceptsAverage.h

@@ -4 +4 @@
  * @author Paul Price, IfA
  *
- * @version $Revision: 1.9.10.1 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-10-29 01:40:54 $
+ * @version $Revision: 1.9.10.2 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-11-08 22:59:04 $
  * Copyright 2005-2006 Institute for Astronomy, University of Hawaii
  */

branches/eam_branch_20071023/psModules/src/concepts/pmConceptsStandard.c

@@ -375 +375 @@
     assert(pattern);
 
-    psList *biassecs; // List of bias sections
+    psList *biassecs = NULL; // List of bias sections
 
     switch (concept->type) {

branches/eam_branch_20071023/psModules/src/config/pmConfig.c

@@ -661 +661 @@
             if (!cameraReadFormats(item->data.md, item->name)) {
                 psWarning("Unable to read formats for camera %s: removed.\n", item->name);
+                psErrorStackPrint(stderr, "errors from read failure\n");
+                psErrorClear();
                 psMetadataRemoveKey(cameras, item->name);
+                continue;
             }
             if (!cameraReadCalibrations(item->data.md, item->name)) {
                 psWarning("Unable to read calibrations for camera %s: removed.\n", item->name);
+                psErrorStackPrint(stderr, "errors from read failure\n");
+                psErrorClear();
                 psMetadataRemoveKey(cameras, item->name);
+                continue;
             }
         }

branches/eam_branch_20071023/psModules/src/config/pmConfig.h

@@ -5 +5 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.31.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-10-29 01:40:54 $
+ *  @version $Revision: 1.31.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-11-08 22:59:04 $
  *
  *  Copyright 2005-2006 Institute for Astronomy, University of Hawaii

branches/eam_branch_20071023/psModules/src/config/pmConfigCamera.c

@@ -464 +464 @@
     while ((formatsItem = psMetadataGetAndIncrement(formatsIter))) {
         assert(formatsItem->type == PS_DATA_METADATA); // We should have read it by now!
-        psMetadata *format = formatsItem->data.V; // The camera format
+        psMetadata *format = formatsItem->data.md; // The camera format
 
         // Add a new RULE which uniquely describes the mosaicked format.  this is needed so

branches/eam_branch_20071023/psModules/src/config/pmConfigRecipes.c

@@ -5 +5 @@
 #include <stdio.h>
 #include <string.h>
-#include <strings.h>            /* for strn?casecmp */
+#include <strings.h>            /* for strn?casecmp */
 #include <unistd.h>
 #include <libgen.h>
@@ -33 +33 @@
         options = psMetadataAlloc ();
         success = psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "OPTIONS",  PS_DATA_METADATA, "", options);
-        assert (success); // type mismatch : OPTIONS already defined but wrong type
-        psFree (options); // drop extra reference
+        assert (success); // type mismatch : OPTIONS already defined but wrong type
+        psFree (options); // drop extra reference
     }
 
@@ -44 +44 @@
         recipe = psMetadataAlloc();
         success = psMetadataAddPtr(options, PS_LIST_TAIL, recipeName,  PS_DATA_METADATA, "", recipe);
-        assert (success); // type mismatch : OPTIONS already defined but wrong type
-        psFree (recipe);  // drop extra reference
+        assert (success); // type mismatch : OPTIONS already defined but wrong type
+        psFree (recipe);  // drop extra reference
     }
     return recipe;
@@ -67 +67 @@
     // master recipe files in the site:recipe location when they are built.
     if (config->site && (source & PM_RECIPE_SOURCE_SITE)) {
-        if (!loadRecipeSite(&status, config, config->site)) {
+        if (!loadRecipeSite(&status, config, config->site)) {
             psError(PS_ERR_IO, false, "Failed to read recipes from site config");
-            return false;
-        } 
-        psTrace ("psModules.config", 3, "read recipes from site config");
+            return false;
+        }
+        psTrace ("psModules.config", 3, "read recipes from site config");
     }
 
@@ -81 +81 @@
         if (!loadRecipeCamera(&status, config, config->camera)) {
             psError(PS_ERR_IO, false, "Failed to read recipes from camera config");
-            return false;
-        } 
-        if (status) {
-            psTrace ("psModules.config", 3, "read recipes from camera config");
-        } else {
-            psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by camera config");
-        }
+            return false;
+        }
+        if (status) {
+            psTrace ("psModules.config", 3, "read recipes from camera config");
+        } else {
+            psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by camera config");
+        }
     }
 
@@ -94 +94 @@
         if (!loadRecipeFromArguments(&status, config)) {
             psError(PS_ERR_IO, false, "Failed to read recipes from command-line arguments");
-            return false;
-        } 
-        if (status) {
-            psTrace ("psModules.config", 3, "read recipes from command-line arguments");
-        } else {
-            psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied on command-line arguments");
-        }
+            return false;
+        }
+        if (status) {
+            psTrace ("psModules.config", 3, "read recipes from command-line arguments");
+        } else {
+            psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied on command-line arguments");
+        }
         if (!loadRecipeSymbols(&status, config)) {
             psError(PS_ERR_IO, false, "Failed to read recipes from symbolic references");
-            return false;
-        } 
-        if (status) {
-            psTrace ("psModules.config", 3, "read recipes from symbolic references");
-        } else {
-            psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by symbolic reference");
-        }
+            return false;
+        }
+        if (status) {
+            psTrace ("psModules.config", 3, "read recipes from symbolic references");
+        } else {
+            psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by symbolic reference");
+        }
         if (!loadRecipeOptions(&status, config)) {
             psError(PS_ERR_IO, false, "Failed to read recipes from symbolic references");
-            return false;
-        } 
-        if (status) {
+            return false;
+        }
+        if (status) {
             psTrace ("psModules.config", 3, "read recipes from command-line arguments");
         } else {
@@ -135 +135 @@
         options = psMetadataAlloc();
         success = psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "OPTIONS",  PS_DATA_METADATA, "Command-line options specified with -D", options);
-        assert (success); // type mismatch : OPTIONS already defined but wrong type
-        psFree (options); // drop extra reference
+        assert (success); // type mismatch : OPTIONS already defined but wrong type
+        psFree (options); // drop extra reference
     }
 
@@ -176 +176 @@
             recipe = psMetadataAlloc();
             success = psMetadataAddPtr(options, PS_LIST_TAIL, recipeName,  PS_DATA_METADATA, "", recipe);
-            assert (success); // type mismatch : recipe already defined but wrong type
-            psFree (recipe); // drop extra reference
+            assert (success); // type mismatch : recipe already defined but wrong type
+            psFree (recipe); // drop extra reference
         }
 
@@ -216 +216 @@
 bool pmConfigLoadRecipeArguments (int *argc, char **argv, pmConfig *config)
 {
-    bool success; 
+    bool success;
 
     psMetadata *recipes = psMetadataLookupMetadata(&success, config->arguments, "RECIPES");
@@ -222 +222 @@
         recipes = psMetadataAlloc();
         success = psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "RECIPES",  PS_DATA_METADATA, "", recipes);
-        assert (success);
-        psFree (recipes);
+        assert (success);
+        psFree (recipes);
     }
 
@@ -271 +271 @@
 // Load the recipe files for SITE : REQUIRED
 static bool loadRecipeSite(bool *status,
-                           pmConfig *config, // The configuration into which to read the recipes
-                           psMetadata *source // The source configuration, from which to read the filenames
+                           pmConfig *config, // The configuration into which to read the recipes
+                           psMetadata *source // The source configuration, from which to read the filenames
     )
 {
@@ -300 +300 @@
         if (fileItem->type != PS_DATA_STRING) {
             psError(PS_ERR_IO, true, "%s in site configuration RECIPES is not of type STR", fileItem->name);
-            return false;
+            return false;
         }
 
         // Read the recipe file. if we fail on a file, give a warning, but continue
         psMetadata *recipe = NULL;
-        if (!pmConfigFileRead(&recipe, fileItem->data.V, "recipe")) {
-            psError(PS_ERR_IO, false, "Failed to read recipe file %s listed in site configuration\n", (char *)fileItem->data.V);
-            return false;
+        if (!pmConfigFileRead(&recipe, fileItem->data.str, "recipe")) {
+            psError(PS_ERR_IO, false, "Failed to read recipe file %s listed in site configuration\n",
+                    fileItem->data.str);
+            return false;
         }
 
         // if this named recipe exists, supplement it
-        psMetadataAdd(config->recipes, PS_LIST_TAIL, fileItem->name, PS_DATA_METADATA | PS_META_REPLACE,
-                      fileItem->comment, recipe);
+        psMetadataAdd(config->recipes, PS_LIST_TAIL, fileItem->name, PS_DATA_METADATA | PS_META_REPLACE,
+                      fileItem->comment, recipe);
         psFree(recipe);  // Drop reference
     }
@@ -326 +327 @@
 // for sourceType == SITE | CAMERA, RECIPES contains a list of files to be read (pmConfigFileRead)
 static bool loadRecipeCamera(bool *status, // status variable
-                             pmConfig *config, // The configuration into which to read the recipes
-                             psMetadata *source // The source configuration, from which to read the filenames
+                             pmConfig *config, // The configuration into which to read the recipes
+                             psMetadata *source // The source configuration, from which to read the filenames
     )
 {
@@ -345 +346 @@
     psMetadata *recipes = psMetadataLookupMetadata(&success, source, "RECIPES"); // The list of recipes
     if (!recipes) {
-        psTrace ("psModules.config", 3, "RECIPES not found in the camera configuration\n");
+        psTrace ("psModules.config", 3, "RECIPES not found in the camera configuration\n");
         return true;
     }
@@ -355 +356 @@
     psMetadataItem *fileItem = NULL;    // MD item containing the filename, from recipe iteration
     while ((fileItem = psMetadataGetAndIncrement(recipesIter))) {
-        char *recipeName = fileItem->name;
-        char *recipeFile = fileItem->data.V;
-
-        psTrace("psModules.config", 3, "Supplementing %s from %s within camera configuration.\n", recipeName, recipeFile);
+        char *recipeName = fileItem->name;
+        char *recipeFile = fileItem->data.str;
+
+        psTrace("psModules.config", 3, "Supplementing %s from %s within camera configuration.\n", recipeName, recipeFile);
 
         // type mismatch is a serious error
@@ -368 +369 @@
         psMetadata *recipe = NULL;
         if (!pmConfigFileRead(&recipe, recipeFile, "recipe")) {
-            psError(PS_ERR_IO, false, "Failed to read recipe file %s listed in camera configuration\n", recipeFile);
-            return false;
+            psError(PS_ERR_IO, false, "Failed to read recipe file %s listed in camera configuration\n",
+                    recipeFile);
+            return false;
         }
 
         // the named recipe must exist; supplement it
-        psMetadata *current = psMetadataLookupMetadata(NULL, config->recipes, recipeName);
-        if (!current) {
-            psError(PS_ERR_IO, false, "Failed to find recipe for %s in master recipe list", recipeName);
-            psFree(recipe);  // Drop reference
-            return false;
-        }
-
-        if (!psMetadataUpdate(current, recipe)) {
-            psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", recipeName);
-            psFree(recipe);  // Drop reference
-            return false;
-        }
+        psMetadata *current = psMetadataLookupMetadata(NULL, config->recipes, recipeName);
+        if (!current) {
+            psError(PS_ERR_IO, false, "Failed to find recipe for %s in master recipe list", recipeName);
+            psFree(recipe);  // Drop reference
+            return false;
+        }
+
+        if (!psMetadataUpdate(current, recipe)) {
+            psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", recipeName);
+            psFree(recipe);  // Drop reference
+            return false;
+        }
         psFree(recipe);  // Drop reference
     }
@@ -397 +399 @@
 // entries for an existing recipe metadata
 static bool loadRecipeFromArguments(bool *status,
-                                    pmConfig *config // The configuration into which to read the recipes
+                                    pmConfig *config // The configuration into which to read the recipes
     )
 {
@@ -425 +427 @@
         // increment the ref counter to protect the data
 
-        psMetadata *recipe = item->data.V; // Recipe of interest
+        psMetadata *recipe = item->data.md; // Recipe of interest
 
         // if this named recipe exists, supplement it
         psMetadata *current = psMetadataLookupMetadata(NULL, config->recipes, item->name);
         if (!current) {
-            psError(PS_ERR_IO, false, "Failed to find recipe for %s in master recipe list", item->name);
-            psFree(recipe);  // Drop reference
-            return false;
-        }
-        psTrace("psModules.config", 3, "Supplementing %s from arguments.\n", item->name);
+            psError(PS_ERR_IO, false, "Failed to find recipe for %s in master recipe list", item->name);
+            psFree(recipe);  // Drop reference
+            return false;
+        }
+        psTrace("psModules.config", 3, "Supplementing %s from arguments.\n", item->name);
 
         if (!psMetadataUpdate (current, recipe)) {
-            psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", item->name);
-            return false;
-        }
+            psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", item->name);
+            return false;
+        }
     }
     psFree(recipesIter);
@@ -449 +451 @@
 // entries for an existing recipe metadata
 static bool loadRecipeSymbols(bool *status,
-                              pmConfig *config // The configuration into which to read the recipes
+                              pmConfig *config // The configuration into which to read the recipes
     )
 {
@@ -465 +467 @@
     while ((item = psMetadataGetAndIncrement(iter))) {
         assert(item->type == PS_DATA_STRING); // It should be this type: we put it in ourselves
-        const char *sourceName = item->data.V; // The name of the symbolic reference
+        const char *sourceName = item->data.str; // The name of the symbolic reference
         const char *targetName = item->name;
-        psTrace("psModules.config", 3, "Supplementing %s from %s.\n", targetName, sourceName);
-
-        // the target recipe must exist; select it
-        psMetadata *targetMD = psMetadataLookupMetadata(&found, config->recipes, targetName);
-        if (!targetMD) {
-            psError(PS_ERR_IO, true, "Failed to find recipe for %s in master recipe list", targetName);
-            return false;
-        }
+        psTrace("psModules.config", 3, "Supplementing %s from %s.\n", targetName, sourceName);
+
+        // the target recipe must exist; select it
+        psMetadata *targetMD = psMetadataLookupMetadata(&found, config->recipes, targetName);
+        if (!targetMD) {
+            psError(PS_ERR_IO, true, "Failed to find recipe for %s in master recipe list", targetName);
+            return false;
+        }
 
         // search for sourceName : it may be in config->recipes or target MD
         psMetadata *sourceMD = NULL;
-        sourceMD = psMetadataLookupMetadata(&found, config->recipes, sourceName);
+        sourceMD = psMetadataLookupMetadata(&found, config->recipes, sourceName);
         if (!sourceMD) {
-            sourceMD = psMetadataLookupMetadata(&found, targetMD, sourceName);
-            if (!sourceMD) {
-                psError(PS_ERR_IO, false, "Selected symbolic name %s does not exist in recipes", sourceName);
-                return false;
-            }
-        }
-
-        if (!psMetadataUpdate(targetMD, sourceMD)) {
-            psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", targetName);
-            return false;
-        }
+            sourceMD = psMetadataLookupMetadata(&found, targetMD, sourceName);
+            if (!sourceMD) {
+                psError(PS_ERR_IO, false, "Selected symbolic name %s does not exist in recipes", sourceName);
+                return false;
+            }
+        }
+
+        if (!psMetadataUpdate(targetMD, sourceMD)) {
+            psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", targetName);
+            return false;
+        }
     }
     psFree(iter);
@@ -501 +503 @@
 // entries for an existing recipe metadata
 static bool loadRecipeOptions(bool *status,
-                              pmConfig *config // The configuration into which to read the recipes
+                              pmConfig *config // The configuration into which to read the recipes
     )
 {
@@ -524 +526 @@
     psMetadataItem *item = NULL;    // MD item containing the filename, from recipe iteration
     while ((item = psMetadataGetAndIncrement(recipesIter))) {
-        char *recipeName = item->name;
-        psMetadata *recipe = item->data.V;
+        char *recipeName = item->name;
+        psMetadata *recipe = item->data.md;
 
         // type mismatch is a serious error
@@ -535 +537 @@
         psMetadata *current = psMetadataLookupMetadata(NULL, config->recipes, recipeName);
         if (!current) {
-            psError(PS_ERR_IO, false, "Selected recipe %s is not found in camera recipe", recipeName);
-            return false;
-        }
+            psError(PS_ERR_IO, false, "Selected recipe %s is not found in camera recipe", recipeName);
+            return false;
+        }
         if (!psMetadataUpdate (current, recipe)) {
-            psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", recipeName);
-            return false;
-        }
+            psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", recipeName);
+            return false;
+        }
     }
     psFree(recipesIter);

branches/eam_branch_20071023/psModules/src/imcombine/pmReadoutCombine.c

@@ -315 +315 @@
             // Combination
             if (!psVectorStats(stats, pixels, errors, mask, 1)) {
-                psError(PS_ERR_UNKNOWN, false, "Error in statistics.");
-                psFree(index);
-                psFree(pixels);
-                psFree(mask);
-                psFree(weights);
-                psFree(errors);
-                psFree(stats);
-                psFree(invScale);
-                return false;
-            }
-
-            outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
-            if (!isfinite(outputImage[yOut][xOut])) {
+                // Can't do much about it, but it's not worth worrying about
+                psErrorClear();
+                outputImage[yOut][xOut] = NAN;
                 outputMask[yOut][xOut] = params->blank;
-            }
-            if (params->weights) {
-                float stdev = psStatsGetValue(stats, combineStdev);
-                outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance
-                // XXXX this is not the correct formal error.
-                // also, the weighted mean is not obviously the correct thing here
+                if (params->weights) {
+                    outputWeight[yOut][xOut] = NAN;
+                }
+            } else {
+                outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
+                if (!isfinite(outputImage[yOut][xOut])) {
+                    outputMask[yOut][xOut] = params->blank;
+                }
+                if (params->weights) {
+                    float stdev = psStatsGetValue(stats, combineStdev);
+                    outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance
+                    // XXXX this is not the correct formal error.
+                    // also, the weighted mean is not obviously the correct thing here
+                }
             }
         }

branches/eam_branch_20071023/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.13.8.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-10-29 01:40:54 $
+ *  @version $Revision: 1.13.8.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-11-08 22:59:04 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii

branches/eam_branch_20071023/psModules/src/imcombine/pmStackReject.c

@@ -52 +52 @@
 
     // Convolve the image with the kernel --- we're basically applying a matched filter and then thresholding
-    psImage *convolved = NULL;          // Convolved image
+    pmReadout *convRO = pmReadoutAlloc(NULL); // Readout with convolved image
+    pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
+    inRO->image = image;
     for (int i = 0; i < numRegions; i++) {
         psRegion *region = regions->data[i]; // Region of interest
         psVector *solution = solutions->data[i]; // Solution of interest
-        if (!pmSubtractionConvolve(&convolved, NULL, NULL, image, NULL, NULL, 0,
-                                   region, solution, kernels, true)) {
+        if (!pmSubtractionConvolve(convRO, inRO, NULL, NULL, 0, region, solution, kernels,
+                                   PM_SUBTRACTION_MODE_1, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
-            psFree(convolved);
-            psFree(image);
+            psFree(convRO);
+            psFree(inRO);
             return NULL;
         }
@@ -73 +75 @@
         if (!kernel) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
-            psFree(convolved);
-            psFree(image);
+            psFree(convRO);
+            psFree(inRO);
             return NULL;
         }
@@ -85 +87 @@
         psFree(kernel);
 
-        psImage *subConv = psImageSubset(convolved, *region); // Sub-image of convolved image
+        psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image
         psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32));
     }
-    psFree(image);
+    psFree(inRO);
+    psImage *convolved = psMemIncrRefCounter(convRO->image);
+    psFree(convRO);
 
     // Threshold the convolved image

branches/eam_branch_20071023/psModules/src/imcombine/pmSubtraction.h

@@ -6 +6 @@
  * @author GLG, MHPCC
  *
- * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-10-17 02:45:40 $
+ * @version $Revision: 1.19.2.1 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-11-08 22:59:04 $
  *
  * Copyright 2004-207 Institute for Astronomy, University of Hawaii
@@ -16 +16 @@
 
 #include <pslib.h>
-#include "pmSubtractionKernels.h"
-#include "pmSubtractionStamps.h"
+
+#include <pmHDU.h>
+#include <pmFPA.h>
+#include <pmSubtractionKernels.h>
+#include <pmSubtractionStamps.h>
 
 /// @addtogroup imcombine Image Combinations
 /// @{
 
+/// Mask values for the subtraction mask
 typedef enum {
-    PM_SUBTRACTION_MASK_CLEAR     = 0x00, // No masking
-    PM_SUBTRACTION_MASK_REF       = 0x01, // Reference image is bad
-    PM_SUBTRACTION_MASK_INPUT     = 0x02, // Input image is bad
-    PM_SUBTRACTION_MASK_CONVOLVE  = 0x04, // If convolved, would be bad
-    PM_SUBTRACTION_MASK_FOOTPRINT = 0x08, // Bad pixel within the stamp footprint
-    PM_SUBTRACTION_MASK_BORDER    = 0x10, // Image border
-    PM_SUBTRACTION_MASK_REJ       = 0x20, // Previously tried as a stamp, and rejected
+    PM_SUBTRACTION_MASK_CLEAR       = 0x00, // No masking
+    PM_SUBTRACTION_MASK_BAD_1       = 0x01, // Image 1 is bad
+    PM_SUBTRACTION_MASK_BAD_2       = 0x02, // Image 2 is bad
+    PM_SUBTRACTION_MASK_CONVOLVE_1  = 0x04, // If image 1 is convolved, would be bad
+    PM_SUBTRACTION_MASK_CONVOLVE_2  = 0x08, // If image 2 is convolved, would be bad
+    PM_SUBTRACTION_MASK_FOOTPRINT_1 = 0x10, // Bad pixel within the stamp footprint of image 1
+    PM_SUBTRACTION_MASK_FOOTPRINT_2 = 0x20, // Bad pixel within the stamp footprint of image 2
+    PM_SUBTRACTION_MASK_BORDER      = 0x40, // Image border
+    PM_SUBTRACTION_MASK_REJ         = 0x80, // Previously tried as a stamp, and rejected
 } pmSubtractionMasks;
 
@@ -45 +51 @@
 bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, ///< Stamp to convolve
                                 const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                int footprint ///< Half-size of region over which to calculate equation
+                                int footprint, ///< Half-size of region over which to calculate equation
+                                pmSubtractionMode mode ///< Mode for subtraction (which to convolve)
     );
 
@@ -51 +58 @@
 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
                                     const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                    int footprint ///< Half-size of region over which to calculate equation
+                                    int footprint, ///< Half-size of region over which to calculate equation
+                                    pmSubtractionMode mode ///< Mode for subtraction (which to convolve)
                                     );
 
@@ -59 +67 @@
                                      );
 
+/// Calculate deviations
+psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, ///< Stamps
+                                           const psVector *solution, ///< Solution vector
+                                           int footprint, ///< Half-size of stamp
+                                           const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                           pmSubtractionMode mode ///< Mode for subtraction
+    );
+
 /// Reject stamps
 int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, ///< Stamps
+                              const psVector *deviations, ///< Deviations for each stamp
                               psImage *subMask, ///< Subtraction mask
-                              const psVector *solution, ///< Solution vector
-                              int footprint, ///< Region to mask if stamp is bad
                               float sigmaRej, ///< Number of RMS deviations above zero at which to reject
-                              const pmSubtractionKernels *kernels ///< Kernel parameters
+                              int footprint ///< Half-size of stamp
     );
 
@@ -81 +96 @@
 
 /// Convolve image in preparation for subtraction
-bool pmSubtractionConvolve(psImage **outImage, ///< Output image
-                           psImage **outWeight, ///< Output weight map (or NULL)
-                           psImage **outMask, ///< Output mask (or NULL)
-                           const psImage *inImage, ///< Input image
-                           const psImage *inWeight, ///< Input weight map (or NULL)
+bool pmSubtractionConvolve(pmReadout *out, ///< Output image
+                           const pmReadout *ro1, // Input image 1
+                           const pmReadout *ro2, // Input image 2
                            const psImage *subMask, ///< Subtraction mask (or NULL)
                            psMaskType blank, ///< Mask value for blank regions
@@ -91 +104 @@
                            const psVector *solution, ///< The solution vector
                            const pmSubtractionKernels *kernels, ///< Kernel parameters
+                           pmSubtractionMode mode, ///< Mode for subtraction
                            bool useFFT  ///< Use Fast Fourier Transform for the convolution?
     );

branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.c

@@ -4 +4 @@
 
 #include <stdio.h>
-#include <string.h>
+#include <strings.h>
 #include <pslib.h>
 
@@ -488 +488 @@
     PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
     PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
-    PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
+    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(inner, size, NULL);
     PS_ASSERT_INT_NONNEGATIVE(ringsOrder, NULL);
 
@@ -572 +572 @@
                     for (int v = -size; v <= size; v++) {
                         int v2 = PS_SQR(v);   // Square of v
-                        float vPoly = power(v, vOrder); // Value of v^vOrder
+//                        float vPoly = power(v, vOrder); // Value of v^vOrder
+                        float vPoly = powf(v/(float)size, vOrder); // Value of v^vOrder
 
                         for (int u = -size; u <= size; u++) {
@@ -578 +579 @@
                             int distance2 = u2 + v2; // Distance from the centre
                             if (distance2 > lower2 && distance2 < upper2) {
-                                float uPoly = power(u, uOrder); // Value of u^uOrder
+//                                float uPoly = power(u, uOrder); // Value of u^uOrder
+                                float uPoly = powf(u/(float)size, uOrder); // Value of u^uOrder
 
                                 float polyVal = uPoly * vPoly; // Value of polynomial

branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.h

@@ -14 +14 @@
     PM_SUBTRACTION_KERNEL_RINGS,        ///< Rings Instead of the Normal Gaussian Subtraction
 } pmSubtractionKernelsType;
+
+/// Modes --- specifies which image to convolve
+typedef enum {
+    PM_SUBTRACTION_MODE_ERR,            // Error in the mode
+    PM_SUBTRACTION_MODE_TARGET,         // Convolve image 1 to match target PSF
+    PM_SUBTRACTION_MODE_1,              // Convolve image 1
+    PM_SUBTRACTION_MODE_2,              // Convolve image 2
+    PM_SUBTRACTION_MODE_UNSURE,         // Not sure yet which image to convolve so try to satisfy both
+    PM_SUBTRACTION_MODE_DUAL,           // Dual convolution
+} pmSubtractionMode;
 
 /// Kernels specification

branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.c

@@ -50 +50 @@
 
 static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read
-                      const pmReadout *reference, // Reference readout
-                      const pmReadout *input, // Input readout, or NULL to generate fake stamps
+                      const pmReadout *ro1, // Readout 1
+                      const pmReadout *ro2, // Readout 2
                       const psImage *subMask, // Mask for subtraction, or NULL
                       psImage *weight,  // Weight map
@@ -57 +57 @@
                       float threshold,  // Threshold for stamp finding
                       float stampSpacing, // Spacing between stamps
-                      float targetWidth,// Target width for fake stamps
                       int size,         // Kernel half-size
-                      int footprint     // Convolution footprint for stamps
+                      int footprint,     // Convolution footprint for stamps
+                      pmSubtractionMode mode // Mode for subtraction
     )
 {
     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
-    *stamps = pmSubtractionStampsFind(*stamps, reference->image, subMask, region, threshold, stampSpacing);
+    *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, stampSpacing, mode);
     if (!*stamps) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
@@ -71 +71 @@
     memCheck("  find stamps");
 
-    if (!input && !pmSubtractionStampsGenerate(*stamps, targetWidth, footprint, size)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps.");
-        return false;
-    }
-
-    memCheck("   generate stamps");
-
     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
-    if (!pmSubtractionStampsExtract(*stamps, reference->image, input ? input->image : NULL,
-                                    weight, footprint, size)) {
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, footprint, size)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
         return false;
@@ -94 +86 @@
 
 
-bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *reference, const pmReadout *input,
+bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *ro1, const pmReadout *ro2,
                         int footprint, float regionSize, float stampSpacing, float threshold,
-                        const psArray *sources, const char *stampsName, float targetWidth,
+                        const psArray *sources, const char *stampsName,
                         pmSubtractionKernelsType type, int size, int spatialOrder,
                         const psVector *isisWidths, const psVector *isisOrders,
                         int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs,
                         int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad,
-                        psMaskType maskBlank, float badFrac)
+                        psMaskType maskBlank, float badFrac, pmSubtractionMode mode)
 {
     PS_ASSERT_PTR_NON_NULL(convolved, false);
-    PS_ASSERT_PTR_NON_NULL(reference, false);
-    PS_ASSERT_IMAGE_NON_NULL(reference->image, false);
-    PS_ASSERT_IMAGE_TYPE(reference->image, PS_TYPE_F32, false);
-    if (reference->mask) {
-        PS_ASSERT_IMAGE_NON_NULL(reference->mask, false);
-        PS_ASSERT_IMAGE_TYPE(reference->mask, PS_TYPE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(reference->mask, reference->image, false);
-    }
-    if (reference->weight) {
-        PS_ASSERT_IMAGE_NON_NULL(reference->weight, false);
-        PS_ASSERT_IMAGE_TYPE(reference->weight, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(reference->weight, reference->image, false);
-    }
-    if (input) {
-        PS_ASSERT_IMAGE_NON_NULL(input->image, false);
-        PS_ASSERT_IMAGE_TYPE(input->image, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(input->image, reference->image, false);
-        if (input->mask) {
-            PS_ASSERT_IMAGE_NON_NULL(input->mask, false);
-            PS_ASSERT_IMAGE_TYPE(input->mask, PS_TYPE_MASK, false);
-            PS_ASSERT_IMAGES_SIZE_EQUAL(input->mask, reference->image, false);
-        }
-        if (input->weight) {
-            PS_ASSERT_IMAGE_NON_NULL(input->weight, false);
-            PS_ASSERT_IMAGE_TYPE(input->weight, PS_TYPE_F32, false);
-            PS_ASSERT_IMAGES_SIZE_EQUAL(input->weight, reference->image, false);
+    PS_ASSERT_PTR_NON_NULL(ro1, false);
+    PS_ASSERT_IMAGE_NON_NULL(ro1->image, false);
+    PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false);
+    if (ro1->mask) {
+        PS_ASSERT_IMAGE_NON_NULL(ro1->mask, false);
+        PS_ASSERT_IMAGE_TYPE(ro1->mask, PS_TYPE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->mask, ro1->image, false);
+    }
+    if (ro1->weight) {
+        PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false);
+        PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false);
+    }
+    if (ro2) {
+        PS_ASSERT_IMAGE_NON_NULL(ro2->image, false);
+        PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->image, ro1->image, false);
+        if (ro2->mask) {
+            PS_ASSERT_IMAGE_NON_NULL(ro2->mask, false);
+            PS_ASSERT_IMAGE_TYPE(ro2->mask, PS_TYPE_MASK, false);
+            PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->mask, ro1->image, false);
+        }
+        if (ro2->weight) {
+            PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false);
+            PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false);
+            PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false);
         }
     } else if (!stampsName && !sources) {
@@ -144 +136 @@
     }
     // stampsName may be anything
-    // targetWidth can be just about anything (except maybe negative, but it can be NAN)
     // We'll check kernel type when we allocate the kernels
     PS_ASSERT_INT_POSITIVE(size, false);
@@ -196 +187 @@
     }
 
-    psImage *inImage = NULL, *inMask = NULL; // Input image, mask, weight
-    if (input) {
-        inImage = input->image;
-        inMask = input->mask;
-    }
-
     // Where does our weight map come from?
     psImage *weight = NULL;             // Weight image to use
-    if (input && input->weight) {
-        weight = input->weight;
-    } else if (reference->weight) {
-        weight = reference->weight;
-    } else if (input) {
-        weight = input->image;
+    if (ro1->weight && ro2 && ro2->weight) {
+        weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight);
+    } else if (ro1->weight) {
+        weight = psMemIncrRefCounter(ro1->weight);
+    } else if (ro2) {
+        if (ro2->weight) {
+            weight = psMemIncrRefCounter(ro2->weight);
+        } else {
+            weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image);
+        }
     } else {
-        weight = reference->image;
+        weight = psMemIncrRefCounter(ro1->image);
     }
 
@@ -221 +210 @@
     psVector *solution = NULL;          // Solution to match PSF
     pmSubtractionKernels *kernels = NULL; // Kernel basis functions
-
-    int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions
+    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
 
     memCheck("start");
 
-    subMask = pmSubtractionMask(reference->mask, inMask, maskBad, size, footprint, badFrac, useFFT);
+    subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskBad, size, footprint,
+                                badFrac, useFFT);
     if (!subMask) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
+        psFree(weight);
         return false;
     }
@@ -260 +250 @@
 
             if (sources) {
-                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing,
-                                                           input ? 0 : 2 * footprint);
+                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, mode);
             } else if (stampsName && strlen(stampsName) > 0) {
-                // Read stamps from file
-                psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
-                const char *stampFormat = input ? "%f %f" : "%f %f %f"; // Format for reading stamp file
-                psArray *stampsData = stampsData = psVectorsReadFromFile(stampsName, stampFormat);
-                psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes
-                if (!stampsData) {
-                    psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName);
-                    goto ERROR;
-                }
-                xStamp = stampsData->data[0];
-                yStamp = stampsData->data[1];
-                if (!input) {
-                    fluxStamp = stampsData->data[2];
-                }
-
-                // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
-                psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
-                psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
-
-                stamps = pmSubtractionStampsSet(xStamp, yStamp, fluxStamp, reference->image, subMask,
-                                                region, stampSpacing, input ? 0 : footprint + size);
-                psFree(stampsData);
+                stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, stampSpacing, mode);
+            }
+
+            // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
+            // doesn't matter.
+            if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing,
+                           size, footprint, mode)) {
+                goto MATCH_ERROR;
+            }
+
+            if (mode == PM_SUBTRACTION_MODE_UNSURE || mode == PM_SUBTRACTION_MODE_TARGET) {
+                pmSubtractionMode newMode = pmSubtractionOrder(stamps, footprint); // Subtraction mode
+                switch (newMode) {
+                  case PM_SUBTRACTION_MODE_1:
+                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2.");
+                    break;
+                  case PM_SUBTRACTION_MODE_2:
+                    if (mode == PM_SUBTRACTION_MODE_TARGET) {
+                        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                                "Input PSF is larger than target PSF --- can't match image.");
+                        goto MATCH_ERROR;
+                    }
+                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1.");
+                    break;
+                  default:
+                    psError(PS_ERR_UNKNOWN, false, "Unable to determine subtraction order.");
+                    goto MATCH_ERROR;
+                }
+                mode = newMode;
             }
 
             // Define kernel basis functions
             if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
-                if (!getStamps(&stamps, reference, input, subMask, weight, NULL,
-                               threshold, stampSpacing, targetWidth, size, footprint)) {
-                    goto ERROR;
-                }
                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
-                                                          stamps, footprint, optThreshold);
+                                                          stamps, footprint, optThreshold, mode);
                 if (!kernels) {
                     psErrorClear();
@@ -314 +306 @@
                 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
 
-                if (!getStamps(&stamps, reference, input, subMask, weight, region,
-                               threshold, stampSpacing, targetWidth, size, footprint)) {
-                    goto ERROR;
+                if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing,
+                               size, footprint, mode)) {
+                    goto MATCH_ERROR;
                 }
 
                 psTrace("psModules.imcombine", 3, "Calculating equation...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels, footprint)) {
+                if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto ERROR;
+                    goto MATCH_ERROR;
                 }
 
@@ -331 +323 @@
                 if (!solution) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto ERROR;
+                    goto MATCH_ERROR;
                 }
 
                 memCheck("  solve equation");
 
+                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
+                                                                        mode); // Deviations for each stamp
+                if (!deviations) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+                    goto MATCH_ERROR;
+                }
+
+                memCheck("   calculate deviations");
+
                 psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
-                numRejected = pmSubtractionRejectStamps(stamps, subMask, solution, footprint, rej, kernels);
+                numRejected = pmSubtractionRejectStamps(stamps, deviations, subMask, rej, footprint);
                 if (numRejected < 0) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
-                    goto ERROR;
-                }
+                    psFree(deviations);
+                    goto MATCH_ERROR;
+                }
+                psFree(deviations);
+
                 memCheck("  reject stamps");
             }
@@ -350 +354 @@
                 if (!solution) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto ERROR;
-                }
-                (void)pmSubtractionRejectStamps(stamps, subMask, solution, footprint, NAN, kernels);
+                    goto MATCH_ERROR;
+                }
+                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
+                                                                        mode); // Deviations for each stamp
+                if (!deviations) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+                    goto MATCH_ERROR;
+                }
+                pmSubtractionRejectStamps(stamps, deviations, subMask, NAN, footprint);
+                psFree(deviations);
             }
             psFree(stamps);
@@ -372 +383 @@
                             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
                             psFree(convKernels);
-                            goto ERROR;
+                            goto MATCH_ERROR;
                         }
 
@@ -380 +391 @@
                             psFree(kernel);
                             psFree(convKernels);
-                            goto ERROR;
+                            goto MATCH_ERROR;
                         }
                         psFree(kernel);
@@ -416 +427 @@
 
             psTrace("psModules.imcombine", 2, "Convolving...\n");
-            if (!pmSubtractionConvolve(&convolved->image, &convolved->weight, &convolved->mask,
-                                       reference->image, reference->weight, subMask, maskBlank, region,
-                                       solution, kernels, useFFT)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image.");
-                goto ERROR;
+            if (!pmSubtractionConvolve(convolved, ro1, ro2, subMask, maskBlank, region,
+                                       solution, kernels, mode, useFFT)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
+                goto MATCH_ERROR;
             }
             psFree(kernels);
@@ -461 +471 @@
     psFree(subMask);
     subMask = NULL;
+    psFree(weight);
+    weight = NULL;
 
     if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
-        goto ERROR;
+        goto MATCH_ERROR;
     }
 
@@ -472 +484 @@
     return true;
 
-ERROR:
+MATCH_ERROR:
     psFree(region);
     psFree(regionString);
@@ -479 +491 @@
     psFree(stamps);
     psFree(solution);
+    psFree(weight);
     return false;
 }
+
+// Calculate the second order moments for an image
+static float subtractionOrderMoment(const psKernel *kernel, // Image for which to measure moments
+                                    int radius   // Maximum radius
+                                    )
+{
+    assert(kernel && kernel->kernel);
+
+    int xMin = PS_MAX(kernel->xMin, -radius), xMax = PS_MIN(kernel->xMax, radius); // Bounds in x
+    int yMin = PS_MAX(kernel->yMin, -radius), yMax = PS_MIN(kernel->yMax, radius); // Bounds in y
+
+    float xCentroid = 0.0, yCentroid = 0.0; // Centroid (first moment)
+    float sum = 0.0;       // Sum (zero-th moment)
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            xCentroid += kernel->kernel[y][x] * x;
+            yCentroid += kernel->kernel[y][x] * y;
+            sum += kernel->kernel[y][x];
+        }
+    }
+    xCentroid /= sum;
+    yCentroid /= sum;
+
+    float eta20 = 0.0, eta02 = 0.0;     // Second moments
+    for (int y = yMin; y <= yMax; y++) {
+        float yDiff = y - yCentroid;
+        for (int x = xMin; x <= xMax; x++) {
+            float xDiff = x - xCentroid;
+            eta20 += PS_SQR(xDiff) * kernel->kernel[y][x];
+            eta02 += PS_SQR(yDiff) * kernel->kernel[y][x];
+        }
+    }
+
+    // Normalise to calculate the scale-invariant
+    float sum2 = PS_SQR(sum);
+    eta20 /= sum2;
+    eta02 /= sum2;
+    // eta11 /= sum2;
+
+    return eta20 + eta02;
+}
+
+#if 0
+// Calculate the deviations for a particular subtraction order
+static psVector *subtractionOrderDeviation(float *sumKernel, // Sum of the kernel
+                                           pmSubtractionStampList *stamps, // Stamps to convolve
+                                           const pmSubtractionKernels *kernels, // Kernel basis functions
+                                           int footprint, // Stamp footprint
+                                           pmSubtractionMode mode // Mode of subtraction
+                                           )
+{
+    assert(stamps);
+    assert(footprint >= 0);
+    assert(mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_2);
+
+    if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return NULL;
+    }
+
+    psVector *solution = pmSubtractionSolveEquation(NULL, stamps);
+    if (!solution) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return NULL;
+    }
+
+    if (sumKernel) {
+        float sum = 0.0;                // Sum of the kernel
+        psImage *image = pmSubtractionKernelImage(solution, kernels, 0.0, 0.0); // Image of kernel
+        for (int y = 0; y < image->numRows; y++) {
+            for (int x = 0; x < image->numCols; x++) {
+                sum += image->data.F32[y][x];
+            }
+        }
+        psFree(image);
+        *sumKernel = sum;
+    }
+
+    psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, mode);
+    psFree(solution);
+    if (!deviations) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+        return NULL;
+    }
+
+    return deviations;
+}
+#endif
+
+pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, int radius)
+{
+    PS_ASSERT_INT_POSITIVE(radius, PM_SUBTRACTION_MODE_ERR);
+
+    psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask for stamps
+    psVector *moments = psVectorAlloc(stamps->num, PS_TYPE_F32); // Moments
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            mask->data.PS_TYPE_MASK_DATA[i] = 0xff;
+            continue;
+        }
+        mask->data.PS_TYPE_MASK_DATA[i] = 0;
+        moments->data.F32[i] = subtractionOrderMoment(stamp->image1, radius) /
+            subtractionOrderMoment(stamp->image2, radius);
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+    if (!psVectorStats(stats, moments, NULL, mask, 0xff)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio.");
+        psFree(mask);
+        psFree(moments);
+        psFree(stats);
+        return PM_SUBTRACTION_MODE_ERR;
+    }
+    psFree(moments);
+    psFree(mask);
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median ratio of second moments: %lf", stats->robustMedian);
+    pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2);
+    psFree(stats);
+
+    return mode;
+}
+
+
+

branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.h

@@ -4 +4 @@
 #include <pslib.h>
 
-#include "pmHDU.h"
-#include "pmFPA.h"
-#include "pmSubtractionKernels.h"
+#include <pmHDU.h>
+#include <pmFPA.h>
+#include <pmSubtractionKernels.h>
+#include <pmSubtractionStamps.h>
 
-
+/// Match two images
 bool pmSubtractionMatch(pmReadout *convolved, ///< Output convolved data
-                        const pmReadout *reference, ///< Reference data
-                        const pmReadout *input, ///< Input data
+                        const pmReadout *ro1, ///< Image 1
+                        const pmReadout *ro2, ///< Image 2
                         // Stamp parameters
                         int footprint,  ///< Stamp half-size
@@ -19 +20 @@
                         const psArray *sources, ///< Sources for stamps
                         const char *stampsName, ///< Filename for stamps
-                        float targetWidth, ///< Width of PSF for simulated target
                         // Kernel parameters
                         pmSubtractionKernelsType type, ///< Kernel type
@@ -38 +38 @@
                         psMaskType maskBad, ///< Value to mask
                         psMaskType maskBlank, ///< Mask for blank region
-                        float badFrac   ///< Maximum fraction of bad input pixels to accept
+                        float badFrac,   ///< Maximum fraction of bad input pixels to accept
+                        pmSubtractionMode mode ///< Mode of subtraction
     );
 
+/// Determine which image to convolve
+pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, ///< Stamps that have been extracted
+                                     int footprint ///< Stamp half-size
+    );
+
+
 #endif

branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionParams.c

@@ -50 +50 @@
 #endif
 
+/// Select the appropriate convolution, given the kernel basis function and subtraction mode
+static inline psKernel *selectConvolution(const pmSubtractionStamp *stamp, // Stamp
+                                          int kernelIndex, // Index for kernel component
+                                          pmSubtractionMode mode // Mode of subtraction
+    )
+{
+    switch (mode) {
+      case PM_SUBTRACTION_MODE_1:
+        return stamp->convolutions1->data[kernelIndex];
+      case PM_SUBTRACTION_MODE_2:
+        return stamp->convolutions2->data[kernelIndex];
+      default:
+        psAbort("Unsupported subtraction mode: %x", mode);
+    }
+    return NULL;                        // Unreached
+}
+
 // Accumulate cross-term sums for a stamp
 static void accumulateCross(double *sumI, // Sum of I(x)/sigma(x)^2
@@ -55 +72 @@
                             double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2
                             const pmSubtractionStamp *stamp, // Stamp with weight
-                            const psKernel *input, // Input image, I(x)
+                            const psKernel *target, // Target stamp
                             int kernelIndex, // Index for kernel component
-                            int footprint // Size of region of interest
+                            int footprint, // Size of region of interest
+                            pmSubtractionMode mode // Mode of subtraction
     )
 {
     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
-    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
+    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
 
     for (int y = -footprint; y <= footprint; y++) {
-        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
         psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
         psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
@@ -82 +100 @@
                                    const pmSubtractionStamp *stamp, // Stamp with input and weight
                                    int kernelIndex, // Index for kernel component
-                                   int footprint // Size of region of interest
+                                   int footprint, // Size of region of interest
+                                   pmSubtractionMode mode // Mode of subtraction
     )
 {
     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
-    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
+    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
 
     for (int y = -footprint; y <= footprint; y++) {
@@ -100 +119 @@
 }
 
-static double accumulateChi2(psKernel *input, // Input stamp
+static double accumulateChi2(const psKernel *target, // Target stamp
                              pmSubtractionStamp *stamp, // Stamp with weight
                              int kernelIndex, // Index for kernel component
                              double coeff, // Coefficient of convolution
                              double bg,  // Background term
-                             int footprint // Size of region of interest
+                             int footprint, // Size of region of interest
+                             pmSubtractionMode mode // Mode of subtraction
     )
 {
     double chi2 = 0.0;
     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
-    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
+    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
 
     for (int y = -footprint; y <= footprint; y++) {
-        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
         psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
         psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
@@ -125 +145 @@
 
 // Return the initial value of chi^2
-static double initialChi2(psKernel *input, // Input stamp
+static double initialChi2(const psKernel *target, // Target stamp
                           const pmSubtractionStamp *stamp, // Stamp with weight
-                          int footprint // Size of convolution
+                          int footprint, // Size of convolution
+                          pmSubtractionMode mode // Mode of subtraction
     )
 {
     psKernel *weight = stamp->weight;   // Weight map
-    psKernel *reference = stamp->reference; // Reference stamp
+    psKernel *source;                   // Source stamp
+    switch (mode) {
+      case PM_SUBTRACTION_MODE_1:
+        source = stamp->image1;
+        break;
+      case PM_SUBTRACTION_MODE_2:
+        source = stamp->image2;
+        break;
+      default:
+        psAbort("Unsupported subtraction mode: %x", mode);
+    }
 
     double chi2 = 0.0;                  // Chi^2
     for (int y = -footprint; y <= footprint; y++) {
-        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
         psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
-        psF32 *ref = &reference->kernel[y][-footprint]; // Derference reference
+        psF32 *ref = &source->kernel[y][-footprint]; // Derference reference
         for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) {
             float diff = *in - *ref;    // Temporary value
@@ -148 +179 @@
 
 // Subtract a convolution from the input
-static void subtractConvolution(psKernel *input, // Input stamp
+static void subtractConvolution(psKernel *target, // Target stamp
                                 const pmSubtractionStamp *stamp, // Stamp with weight
                                 int kernelIndex, // Index for kernel component
                                 float coeff, // Coefficient of subtraction
                                 float bg, // Background term
-                                int footprint // Size of region of interest
-    )
-{
-    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
-
+                                int footprint, // Size of region of interest
+                                pmSubtractionMode mode // Mode of subtraction
+    )
+{
+    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
     for (int y = -footprint; y <= footprint; y++) {
-        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
         psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
         for (int x = -footprint; x <= footprint; x++, in++, conv++) {
@@ -173 +204 @@
                                                       int spatialOrder, const psVector *fwhms, int maxOrder,
                                                       const pmSubtractionStampList *stamps, int footprint,
-                                                      float tolerance)
+                                                      float tolerance, pmSubtractionMode mode)
 {
     if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) {
@@ -210 +241 @@
     // Need to save the stamp inputs --- we're changing the values!
     int numStamps = stamps->num;        // Number of stamps
-    psArray *inputs = psArrayAlloc(numStamps); // Deep copies of the inputs
+    psArray *targets = psArrayAlloc(numStamps); // Deep copies of the targets
     psVector *badStamps = psVectorAlloc(numStamps, PS_TYPE_U8); // Mark the bad stamps
     psVectorInit(badStamps, 0);
@@ -219 +250 @@
             continue;
         }
-        psKernel *input = stamp->input; // Input image of interest
-        psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image
-        inputs->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint);
+        psKernel *target;               // Target image of interest
+        switch (mode) {
+          case PM_SUBTRACTION_MODE_1:
+            target = stamp->image2;
+            break;
+          case PM_SUBTRACTION_MODE_2:
+            target = stamp->image1;
+            break;
+          default:
+            psAbort("Unsupported subtraction mode: %x", mode);
+        }
+        psImage *copy = psImageCopy(NULL, target->image, PS_TYPE_F32); // Copy of the image
+        targets->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint);
         psFree(copy);                   // Drop reference
     }
@@ -238 +279 @@
             continue;
         }
-        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
+        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
-            psFree(inputs);
+            psFree(targets);
             psFree(kernels);
             psFree(badStamps);
@@ -258 +299 @@
                     "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n",
                     i, (int)stamp->x, (int)stamp->y);
-            psFree(inputs);
+            psFree(targets);
             psFree(kernels);
             psFree(badStamps);
@@ -265 +306 @@
 
         for (int j = 0; j < numKernels; j++) {
-            accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint);
-        }
-
-        lastChi2 += initialChi2(inputs->data[i], stamp, footprint);
+            accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint, mode);
+        }
+
+        lastChi2 += initialChi2(targets->data[i], stamp, footprint, mode);
         numPixels += PS_SQR(2 * footprint + 1);
     }
@@ -297 +338 @@
                 }
                 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
-                accumulateCross(&sumI, &sumII, &sumIC, stamp, inputs->data[j], i, footprint);
+                accumulateCross(&sumI, &sumII, &sumIC, stamp, targets->data[j], i, footprint, mode);
             }
 
@@ -310 +351 @@
                 }
                 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
-                chi2 += accumulateChi2(inputs->data[j], stamp, i, coeff, bg, footprint);
+                chi2 += accumulateChi2(targets->data[j], stamp, i, coeff, bg, footprint, mode);
             }
 
@@ -328 +369 @@
         if (bestIndex == -1) {
             psError(PS_ERR_UNKNOWN, false, "Unable to find best kernel component in round %d.", iter);
-            psFree(inputs);
+            psFree(targets);
             psFree(sumC);
             psFree(sumCC);
@@ -345 +386 @@
             }
             pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
-            subtractConvolution(inputs->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint);
+            subtractConvolution(targets->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint, mode);
         }
 
@@ -361 +402 @@
         lastChi2 = bestChi2;
     }
-    psFree(inputs);
+    psFree(targets);
     psFree(sumC);
     psFree(sumCC);
@@ -401 +442 @@
                 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
                 psArray *convolutions = convNew->data[j]; // Convolutions for this stamp
-                convolutions->data[rank] = psMemIncrRefCounter(stamp->convolutions->data[i]);
+                convolutions->data[rank] = psMemIncrRefCounter(selectConvolution(stamp, i, mode));
             }
         }
@@ -420 +461 @@
         }
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        psFree(stamp->convolutions);
-        stamp->convolutions = convNew->data[i];
+        psFree(stamp->convolutions1);
+        psFree(stamp->convolutions2);
+        switch (mode) {
+          case PM_SUBTRACTION_MODE_1:
+            stamp->convolutions1 = convNew->data[i];
+            stamp->convolutions2 = NULL;
+            break;
+          case PM_SUBTRACTION_MODE_2:
+            stamp->convolutions1 = NULL;
+            stamp->convolutions2 = convNew->data[i];
+            break;
+          default:
+            psAbort("Unsupported subtraction mode: %x", mode);
+        }
     }
 

branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionParams.h

@@ -15 +15 @@
                                                       const pmSubtractionStampList *stamps, ///< Stamps
                                                       int footprint, ///< Convolution footprint for stamps
-                                                      float tolerance ///< Maximum difference in chi^2
+                                                      float tolerance, ///< Maximum difference in chi^2
+                                                      pmSubtractionMode mode // Mode for subtraction
     );
 

branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionStamps.c

@@ -4 +4 @@
 
 #include <stdio.h>
+#include <string.h>
 #include <pslib.h>
 
@@ -45 +46 @@
     psFree(stamp->matrix);
     psFree(stamp->vector);
-    psFree(stamp->reference);
-    psFree(stamp->input);
+    psFree(stamp->image1);
+    psFree(stamp->image2);
     psFree(stamp->weight);
-    psFree(stamp->convolutions);
+    psFree(stamp->convolutions1);
+    psFree(stamp->convolutions2);
 }
 
@@ -65 +67 @@
 // Is this position unmasked?
 static bool checkStampMask(int x, int y, // Coordinates of stamp
-                           const psImage *mask
+                           const psImage *mask, // Mask
+                           pmSubtractionMode mode // Mode for subtraction
                            )
 {
@@ -74 +77 @@
         return false;
     }
-    return (mask->data.PS_TYPE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER |
-                                                  PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_REJ)) ?
-        false : true;
+
+    psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ; // Mask value
+    switch (mode) {
+      case PM_SUBTRACTION_MODE_1:
+      case PM_SUBTRACTION_MODE_TARGET:
+        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1;
+        break;
+      case PM_SUBTRACTION_MODE_2:
+        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_2;
+        break;
+      case PM_SUBTRACTION_MODE_UNSURE:
+      case PM_SUBTRACTION_MODE_DUAL:
+        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1 | PM_SUBTRACTION_MASK_FOOTPRINT_2;
+        break;
+      default:
+        psAbort("Unsupported subtraction mode: %x", mode);
+    }
+
+    return (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? false : true;
 }
 
@@ -140 +159 @@
     stamp->status = PM_SUBTRACTION_STAMP_INIT;
 
-    stamp->reference = NULL;
-    stamp->input = NULL;
+    stamp->image1 = NULL;
+    stamp->image2 = NULL;
     stamp->weight = NULL;
-    stamp->convolutions = NULL;
+    stamp->convolutions1 = NULL;
+    stamp->convolutions2 = NULL;
 
     return stamp;
@@ -151 +171 @@
 pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image,
                                                 const psImage *subMask, const psRegion *region,
-                                                float threshold, float spacing)
+                                                float threshold, float spacing, pmSubtractionMode mode)
 {
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
@@ -226 +246 @@
             for (int y = subRegion->y0; y <= subRegion->y1 ; y++) {
                 for (int x = subRegion->x0; x <= subRegion->y1 ; x++) {
-                    if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxStamp) {
+                    if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) {
                         fluxStamp = image->data.F32[y][x];
                         xStamp = x;
@@ -242 +262 @@
 
             // Reset the postage stamps since we're making a new stamp
-            psFree(stamp->reference);
-            psFree(stamp->input);
+            psFree(stamp->image1);
+            psFree(stamp->image2);
             psFree(stamp->weight);
-            stamp->reference = stamp->input = stamp->weight = NULL;
-
-            stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
+            psFree(stamp->convolutions1);
+            psFree(stamp->convolutions2);
+            stamp->image1 = stamp->image2 = stamp->weight = NULL;
+            stamp->convolutions1 = stamp->convolutions2 = NULL;
+
+            stamp->status = PM_SUBTRACTION_STAMP_FOUND;
             numFound++;
             psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n",
@@ -266 +289 @@
 pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux,
                                                const psImage *image, const psImage *subMask,
-                                               const psRegion *region, float spacing, int exclusionZone)
+                                               const psRegion *region, float spacing, pmSubtractionMode mode)
+
 {
     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
@@ -289 +313 @@
     }
     PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(exclusionZone, NULL);
 
     int numStars = x->n;                // Number of stars
-    psVector *exclude = psVectorAlloc(numStars, PS_TYPE_U8); // Exclude a star?
-    psVectorInit(exclude, 0);
-
-    // Apply exclusion zone around each star --- important when we're convolving to a specified PSF; less so
-    // when we're trying to get a good subtraction.
-    if (exclusionZone > 0) {
-        psTrace("psModules.imcombine", 2, "Applying exclusion zone of %d pixels for stamps\n", exclusionZone);
-        // We use something resembling a binary search --- coordinates are sorted in the x dimension, and then
-        // to exclude stars within a nominated distance, we need only examine (i.e., calculate the
-        // 2-dimensional distance to) other stars in the list that are within that distance of the x
-        // coordinate of the star of interest.  By marking both stars when we find a violation of the
-        // exclusion zone, we need only search upwards from the x coordinate of the star of interest.
-
-        int minDistance2 = PS_SQR(exclusionZone); // Minimum square distance for other stars
-        psVector *indexes = psVectorSortIndex(NULL, x); // Indices for sorting in x
-        for (int i = 0; i < numStars - 1; i++) {
-            int iIndex = indexes->data.S32[i]; // Index for star i
-            if (exclude->data.U8[iIndex]) {
-                continue;
-            }
-            float ix = x->data.F32[iIndex], iy = y->data.F32[iIndex]; // Coordinates for star i
-            float jx, jy;                   // Coordinates for star j
-            for (int j = i + 1, jIndex = indexes->data.S32[j];
-                 j < numStars && (jx = x->data.F32[jIndex]) < ix + exclusionZone;
-                 j++, jIndex = indexes->data.S32[j]) {
-                jy = y->data.F32[jIndex];
-                if (PS_SQR(ix - jx) + PS_SQR(iy - jy) < minDistance2) {
-                    exclude->data.U8[iIndex] = 0xff;
-                    exclude->data.U8[jIndex] = 0xff;
-                    psTrace("psModules.imcombine", 5, "Excluding stamps %d,%d and %d,%d\n",
-                            (int)x->data.F32[iIndex], (int)y->data.F32[iIndex],
-                            (int)x->data.F32[jIndex], (int)y->data.F32[jIndex]);
-                    // Would 'break' here, but there might be more stamps within the exclusion zone.
-                }
-            }
-        }
-        psFree(indexes);
-    }
-
     pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
                                                                  region, spacing); // Stamp list
@@ -347 +331 @@
     // Put the stars into their appropriate subregions
     for (int i = 0; i < numStars; i++) {
-        if (exclude->data.U8[i]) {
-            // Star has been excluded
-            continue;
-        }
         float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp
         int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp
@@ -357 +337 @@
             continue;
         }
-        if (!checkStampMask(xPix, yPix, subMask)) {
+        if (!checkStampMask(xPix, yPix, subMask, mode)) {
             // Not a good stamp
             continue;
@@ -386 +366 @@
         }
     }
-    psFree(exclude);
 
     // Sort the list by flux, with the brightest last
@@ -421 +400 @@
 
 
-bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *reference, psImage *input,
+bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
                                 psImage *weight, int footprint, int kernelSize)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
-    PS_ASSERT_IMAGE_NON_NULL(reference, false);
-    PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false);
-    if (input) {
-        PS_ASSERT_IMAGE_NON_NULL(input, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);
-        PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false);
-    }
-    if (weight) {
-        PS_ASSERT_IMAGE_NON_NULL(weight, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(weight, reference, false);
-        PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
-    }
+    PS_ASSERT_IMAGE_NON_NULL(image1, false);
+    PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, false);
+    if (image2) {
+        PS_ASSERT_IMAGE_NON_NULL(image2, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(image2, image1, false);
+        PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false);
+    }
+    PS_ASSERT_IMAGE_NON_NULL(weight, false);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false);
+    PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
     PS_ASSERT_INT_NONNEGATIVE(footprint, false);
     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
 
-    int numCols = reference->numCols, numRows = reference->numRows; // Size of images
+    int numCols = image1->numCols, numRows = image1->numRows; // Size of images
     int size = kernelSize + footprint; // Size of postage stamps
-
-    if (!weight) {
-        // Use the input (or if I must, the reference) as a rough approximation to the variance map, and HOPE
-        // that it's positive.
-        weight = input ? input : reference;
-    }
 
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
+        if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_FOUND) {
             continue;
         }
@@ -468 +439 @@
         }
 
+        // Catch memory leaks --- these should have been freed and NULLed before
+        assert(stamp->image1 == NULL);
+        assert(stamp->image2 == NULL);
+        assert(stamp->weight == NULL);
+
         psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest
 
-        psImage *refSub = psImageSubset(reference, region); // Subimage with stamp
-        stamp->reference = psKernelAllocFromImage(refSub, size, size);
-        psFree(refSub);                 // Drop reference
-
-        if (input) {
-            psImage *inSub = psImageSubset(input, region); // Subimage with stamp
-            stamp->input = psKernelAllocFromImage(inSub, size, size);
-            psFree(inSub);              // Drop reference
+        psImage *sub1 = psImageSubset(image1, region); // Subimage with stamp
+        stamp->image1 = psKernelAllocFromImage(sub1, size, size);
+        psFree(sub1);                   // Drop reference
+
+        if (image2) {
+            psImage *sub2 = psImageSubset(image2, region); // Subimage with stamp
+            stamp->image2 = psKernelAllocFromImage(sub2, size, size);
+            psFree(sub2);               // Drop reference
         }
 
         psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp
         stamp->weight = psKernelAllocFromImage(wtSub, size, size);
-        psFree(wtSub);                 // Drop reference
+        psFree(wtSub);                  // Drop reference
+
+        stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
     }
 
@@ -488 +466 @@
 }
 
-
+#if 0
 bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize)
 {
@@ -517 +495 @@
         float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame
 
-        psFree(stamp->input);
-        stamp->input = psKernelAlloc(-size, size, -size, size);
-        psKernel *input = stamp->input; // Target stamp
+        psFree(stamp->image2);
+        stamp->image2 = psKernelAlloc(-size, size, -size, size);
+        psKernel *target = stamp->image2; // Target stamp
 
         // Put in a Waussian, just for fun!
@@ -525 +503 @@
             for (int u = -size; u <= size; u++) {
                 float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma));
-                input->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));
+                target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));
             }
         }
@@ -533 +511 @@
     return true;
 }
-
+#endif
 
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask,
                                                           const psRegion *region, float spacing,
-                                                          int exclusionZone)
+                                                          pmSubtractionMode mode)
 {
     PS_ASSERT_ARRAY_NON_NULL(sources, NULL);
@@ -561 +539 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region,
-                                                            spacing, exclusionZone); // Stamps to return
+                                                            spacing, mode); // Stamps to return
     psFree(x);
     psFree(y);
@@ -572 +550 @@
     return stamps;
 }
+
+
+pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *subMask,
+                                                       const psRegion *region, float spacing,
+                                                       pmSubtractionMode mode)
+{
+    PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
+    // Let pmSubtractionStampsSet take care of the rest of the assertions
+
+    const char *format = (mode == PM_SUBTRACTION_MODE_TARGET ? "%f %f %f" : "%f %f"); // Format of file
+    psArray *data = psVectorsReadFromFile(filename, format);
+    if (!data) {
+        psError(PS_ERR_IO, false, "Unable to read stamps file %s", filename);
+        return NULL;
+    }
+    psVector *x = data->data[0], *y = data->data[1]; // Stamp positions
+    psVector *flux = (mode == PM_SUBTRACTION_MODE_TARGET ? data->data[2] : NULL); // Stamp fluxes
+
+    // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
+    psBinaryOp(x, x, "-", psScalarAlloc(1.0, PS_TYPE_F32));
+    psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32));
+
+    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, spacing,
+                                                            mode);
+    psFree(data);
+
+    return stamps;
+
+}

branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionStamps.h

@@ -9 +9 @@
 typedef enum {
     PM_SUBTRACTION_STAMP_INIT,          ///< Initial state
+    PM_SUBTRACTION_STAMP_FOUND,         ///< Found a suitable source for this stamp
+    PM_SUBTRACTION_STAMP_CALCULATE,     ///< Calculate matrix and vector values for this stamp
     PM_SUBTRACTION_STAMP_USED,          ///< Use this stamp
     PM_SUBTRACTION_STAMP_REJECTED,      ///< This stamp has been rejected
-    PM_SUBTRACTION_STAMP_CALCULATE,     ///< Calculate matrix and vector values for this stamp
     PM_SUBTRACTION_STAMP_NONE           ///< No stamp in this region
 } pmSubtractionStampStatus;
@@ -54 +55 @@
     float flux;                         ///< Flux
     float xNorm, yNorm;                 ///< Normalised position
-    psKernel *reference;                ///< Reference image postage stamp
-    psKernel *input;                    ///< Input image postage stamp
+    psKernel *image1;                   ///< Reference image postage stamp
+    psKernel *image2;                   ///< Input image postage stamp
     psKernel *weight;                   ///< Weight image postage stamp
-    psArray *convolutions;              ///< Convolutions of the reference for each kernel component
+    psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component
+    psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component
     psImage *matrix;                    ///< Associated matrix
     psVector *vector;                   ///< Assoicated vector
@@ -72 +74 @@
                                                 const psRegion *region, ///< Region to search, or NULL
                                                 float threshold, ///< Threshold for stamps in the image
-                                                float spacing ///< Rough spacing for stamps
+                                                float spacing, ///< Rough spacing for stamps
+                                                pmSubtractionMode mode ///< Mode for subtraction
     );
 
 /// Set stamps based on a list of x,y
-///
-/// We may optionally apply an exclusion zone around each star --- this is important when we're convolving to
-/// a specified PSF; less so when we're only trying to get a good subtraction.
 pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp
                                                const psVector *y, ///< y coordinates for each stamp
@@ -86 +86 @@
                                                const psRegion *region, ///< Region to search, or NULL
                                                float spacing, ///< Rough spacing for stamps
-                                               int exclusionZone ///< Exclusion zone around each star
+                                               pmSubtractionMode mode ///< Mode for subtraction
     );
 
-///
+/// Set stamps based on a list of sources
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(
-    const psArray *sources, ///< Sources for each stamp
-    const psImage *subMask, ///< Mask, or NULL
-    const psRegion *region, ///< Region to search, or NULL
-    float spacing, ///< Rough spacing for stamps
-    int exclusionZone ///< Exclusion zone around each star
+    const psArray *sources,             ///< Sources for each stamp
+    const psImage *subMask,             ///< Mask, or NULL
+    const psRegion *region,             ///< Region to search, or NULL
+    float spacing,                      ///< Rough spacing for stamps
+    pmSubtractionMode mode              ///< Mode for subtraction
+    );
+
+/// Set stamps based on values in a file
+pmSubtractionStampList *pmSubtractionStampsSetFromFile(
+    const char *filename,               ///< Filename of file containing x,y (or x,y,flux) on each line
+    const psImage *subMask,             ///< Mask, or NULL
+    const psRegion *region,             ///< Region to search, or NULL
+    float spacing,                      ///< Rough spacing for stamps
+    pmSubtractionMode mode              ///< Mode for subtraction
     );
 
@@ -107 +116 @@
     );
 
-/// Generate target for stamps based on list
-bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, ///< Stamps
-                                 float fwhm, ///< FWHM for each stamp
-                                 int footprint, ///< Stamp footprint size
-                                 int size ///< Kernel half-size
-    );
-
 #endif

branches/eam_branch_20071023/psModules/src/objects/pmTrend2D.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.4.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-10-29 22:42:34 $
+ *  @version $Revision: 1.4.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-11-08 22:59:04 $
  *
  *  Copyright 2004 Institute for Astronomy, University of Hawaii
@@ -11 +11 @@
 // XXXX: Ignore (2)
 
+// XXXX: ignore
+
+// XXXX: ignore
+
 #ifdef HAVE_CONFIG_H
 #include <config.h>
@@ -32 +36 @@
 pmTrend2D *pmTrend2DAlloc (pmTrend2DMode mode, psImage *image, int nXtrend, int nYtrend, psStats *stats)
 {
-    assert (image);
+    PS_ASSERT_PTR_NON_NULL(stats, NULL);
+    if (mode == PM_TREND_MAP) {
+        PS_ASSERT_PTR_NON_NULL(image, NULL);
+    }
 
     pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
@@ -73 +80 @@
           break;
       }
-
+      // XXX: Put a more graceful error here.
       default:
         psAbort ("error");
     }
     return (trend);
+}
+
+bool psMemCheckTrend2D(psPtr ptr)
+{
+    PS_ASSERT_PTR(ptr, false);
+    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmTrend2DFree);
 }
 
@@ -123 +136 @@
 pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats)
 {
+    PS_ASSERT_PTR_NON_NULL(stats, NULL);
     pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
     psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree);
@@ -164 +178 @@
 
       default:
+        // XXX: Put a more graceful error here.
         psAbort ("error");
     }
@@ -169 +184 @@
 }
 
-bool pmTrend2DFit (pmTrend2D *trend, psVector *mask, psMaskType maskVal, psVector *x, psVector *y, psVector *f, psVector *df) {
+bool pmTrend2DFit (pmTrend2D *trend, psVector *mask, psMaskType maskVal, psVector *x,
+                   psVector *y, psVector *f, psVector *df)
+{
+    PS_ASSERT_PTR_NON_NULL(trend, false);
 
     bool status;
@@ -199 +217 @@
 }
 
-double pmTrend2DEval (pmTrend2D *trend, float x, float y) {
+double pmTrend2DEval (pmTrend2D *trend, float x, float y)
+{
+    if (!trend) return 0.0;
 
     double result;
-
-    assert (trend);
-
     switch (trend->mode) {
       case PM_TREND_POLY_ORD:
@@ -221 +238 @@
 }
 
-psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y) {
-
+psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y)
+{
+    PS_ASSERT_PTR_NON_NULL(trend, NULL);
     psVector *result;