Changeset 6873

Timestamp:

Apr 17, 2006, 8:10:08 AM (20 years ago)

Author:

magnier

Message:

fixed up conflicts

Location:

trunk/psModules/src

Files:

• astrom/pmAstrometry.c (modified) (1 diff)
• astrom/pmAstrometry.h (modified) (1 diff)
• astrom/pmAstrometryObjects.h (modified) (1 diff)
• config/pmConfig.c (modified) (22 diffs)
• config/pmConfig.h (modified) (1 diff)
• detrend/pmFlatField.c (modified) (1 diff)
• detrend/pmMaskBadPixels.c (modified) (1 diff)
• detrend/pmMaskBadPixelsErrors.h (modified) (1 diff)
• detrend/pmNonLinear.c (modified) (1 diff)
• imcombine/pmReadoutCombine.c (modified) (5 diffs)
• imcombine/pmReadoutCombine.h (modified) (1 diff)
• imsubtract/pmSubtractBias.c (modified) (13 diffs)
• imsubtract/pmSubtractBias.h (modified) (1 diff)
• imsubtract/pmSubtractSky.h (modified) (1 diff)
• objects/pmModelGroup.c (modified) (1 diff)
• objects/pmModelGroup.h (modified) (1 diff)
• objects/pmObjects.c (modified) (2 diffs)
• objects/pmPSF.c (modified) (8 diffs)
• objects/pmPSF.h (modified) (1 diff)
• objects/pmPSFtry.c (modified) (13 diffs)
• objects/pmPSFtry.h (modified) (1 diff)

trunk/psModules/src/astrom/pmAstrometry.c

@@ -13 +13 @@
 * XXX: Should we implement non-linear cell->chip transforms?
 *
-*  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-04-17 18:01:04 $
+*  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-04-17 18:10:08 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/astrom/pmAstrometry.h

@@ -8 +8 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-04-17 18:01:04 $
+*  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-04-17 18:10:08 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/astrom/pmAstrometryObjects.h

@@ -8 +8 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-04-17 18:01:04 $
+*  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-04-17 18:10:08 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/config/pmConfig.c

@@ -3 +3 @@
  *  @author PAP, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-10 02:43:19 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -11 +11 @@
 #include <stdio.h>
 #include <string.h>
+#include <unistd.h>
+#include <assert.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <glob.h>
 #include "pslib.h"
 #include "pmConfig.h"
 
-#define PS_SITE "PS_SITE"               // Name of the environment variable containing the site config file
-#define PS_DEFAULT_SITE "ipprc.config"  // Default site config file
+#define PS_SITE "PS_SITE"         // Name of the environment variable containing the site config file
+#define PS_DEFAULT_SITE ".ipprc"  // Default site config file
+
+static psArray *configPath = NULL;
+
+static void configFree(pmConfig *config)
+{
+    psFree(config->site);
+    psFree(config->files);
+    psFree(config->camera);
+    psFree(config->recipes);
+    psFree(config->arguments);
+    psFree(config->database);
+}
+
+pmConfig *pmConfigAlloc(void)
+{
+    pmConfig *config = psAlloc(sizeof(pmConfig));
+    (void)psMemSetDeallocator(config, (psFreeFunc)configFree);
+
+    // Initialise
+    config->site = NULL;
+    config->camera = NULL;
+    config->recipes = NULL;
+    config->arguments = NULL;
+    config->database = NULL;
+
+    // the file structure is used to carry pmFPAfiles
+    config->files = psMetadataAlloc ();
+    return config;
+}
+
+void pmConfigDone(void)
+{
+    psFree(configPath);
+}
+
 
 /** readConfig
@@ -23 +63 @@
  *
  */
-static bool readConfig(
+bool readConfig(
     psMetadata **config,                // Config to output
     const char *name,                   // Name of file
     const char *description)            // Description of file
 {
+    char *realName = NULL;
     unsigned int numBadLines = 0;
+    struct stat filestat;
 
     psLogMsg(__func__, PS_LOG_INFO, "Loading %s configuration from file %s\n",
              description, name);
-    *config = psMetadataConfigParse(NULL, &numBadLines, name, true);
+
+    uid_t uid = getuid();
+    gid_t gid = getgid();
+
+    // we try: name, path[0]/name, path[1]/name, ...
+    // find the first existing entry in the path (starting with the bare name)
+    realName = psStringCopy (name);
+    psTrace (__func__, 5, "trying %s\n", realName);
+
+    int status = stat (realName, &filestat);
+    if (status == 0) {
+        if ((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR))
+            goto found;
+        if ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP))
+            goto found;
+        if (filestat.st_mode & S_IROTH)
+            goto found;
+    }
+    psFree (realName);
+
+    if (configPath == NULL) {
+        psError(PS_ERR_IO, false, "Cannot find %s configuration file in path\n", description);
+        return false;
+    }
+
+    for (int i = 0; i < configPath->n; i++) {
+        realName = psStringCopy (configPath->data[i]);
+        psStringAppend (&realName, "/%s", name);
+        psTrace (__func__, 5, "trying %s\n", realName);
+
+        status = stat (realName, &filestat);
+        if (status == 0) {
+            if ((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR))
+                goto found;
+            if ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP))
+                goto found;
+            if (filestat.st_mode & S_IROTH)
+                goto found;
+        }
+        psFree (realName);
+    }
+
+    psError(PS_ERR_IO, false, "Cannot find %s configuration file in path\n", description);
+    return false;
+
+found:
+    *config = psMetadataConfigParse(NULL, &numBadLines, realName, true);
     if (numBadLines > 0) {
         psLogMsg(__func__, PS_LOG_WARN, "%d bad lines in %s configuration file (%s)\n",
-                 description, name);
+                 description, realName);
     }
     if (!*config) {
         psError(PS_ERR_IO, false, "Unable to read %s configuration from %s\n",
-                description, name);
+                description, realName);
+        psFree (realName);
         return false;
     }
 
+    psFree (realName);
     return true;
 }
@@ -57 +147 @@
  line.
  *****************************************************************************/
-bool pmConfigRead(
-    psMetadata **site,
-    psMetadata **camera,
-    psMetadata **recipe,
+pmConfig *pmConfigRead(
     int *argc,
-    char **argv,
-    const char *recipeName)
-{
-    /* XXX: We need clarification on what is to be done if these arguments are
-    NULL.
-        PS_ASSERT_PTR_NON_NULL(site, false);
-        PS_ASSERT_PTR_NON_NULL(*site, false);
-        PS_ASSERT_PTR_NON_NULL(camera, false);
-        PS_ASSERT_PTR_NON_NULL(*camera, false);
-        PS_ASSERT_PTR_NON_NULL(recipe, false);
-        PS_ASSERT_PTR_NON_NULL(*recipe, false);
-        PS_ASSERT_INT_POSITIVE(*argc, false);
-        PS_ASSERT_PTR_NON_NULL(argv, false);
-    */
+    char **argv)
+{
+    PS_ASSERT_INT_POSITIVE(*argc, false);
+    PS_ASSERT_PTR_NON_NULL(argv, false);
+
+    pmConfig *config = pmConfigAlloc(); // The configuration, containing site, camera and recipes
+
     //
     // The following section of code attempts to determine which file is
@@ -97 +177 @@
                      "-site command-line switch provided without the required filename --- ignored.\n");
         } else {
-            siteName = argv[argNum];
+            siteName = psStringCopy(argv[argNum]);
             psArgumentRemove(argNum, argc, argv);
         }
@@ -106 +186 @@
     if (!siteName) {
         siteName = getenv(PS_SITE);
+        if (siteName) {
+            siteName = psStringCopy (siteName);
+        }
     }
 
@@ -111 +194 @@
     // Last chance is ~/.ipprc
     //
-    bool cleanupSiteName = false; // Do I have to psFree siteName?
     if (!siteName) {
-        siteName = psStringCopy(PS_DEFAULT_SITE);
-        cleanupSiteName = true;
-    }
-
-    //
-    // We have the connfiguration filename; now we read and parse the config
+        char *home = getenv("HOME");
+        siteName = psStringCopy(home);
+        psStringAppend(&siteName, "/%s", PS_DEFAULT_SITE);
+    }
+
+
+    // We have the configuration filename; now we read and parse the config
     // file and store in psMetadata struct site.
     //
 
-    if (!readConfig(site, siteName, "site")) {
-        return false;
-    }
-    if (cleanupSiteName) {
-        psFree(siteName);
-    }
-
-
-    //
-    // Next, we do a similar thing for the recipe configuration file.  The
-    // file is read and parsed into psMetadata struct "recipe".
-    //
-    //
-    argNum = psArgumentGet(*argc, argv, "-recipe");
-    if (argNum > 0) {
-        psArgumentRemove(argNum, argc, argv);
-        if (argNum >= *argc) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "-recipe command-line switch provided without the required filename --- ignored.\n");
-        } else {
-            psArgumentRemove(argNum, argc, argv);
-            readConfig(recipe, argv[argNum], "recipe");
-        }
-    }
-    // Or, load the recipe from the camera file, if appropriate
-    if (! *recipe && *camera && recipeName) {
-        *recipe = pmConfigRecipeFromCamera(*camera, recipeName);
-    }
-
-
-    //
+    if (!readConfig(&config->site, siteName, "site")) {
+        psFree(config);
+        return NULL;
+    }
+    psFree(siteName);
+
+    // define the config-file search path (configPath)
+    if (configPath) {
+        psFree(configPath);
+        configPath = NULL;
+    }
+    char *path = psMetadataLookupStr(NULL, config->site, "PATH");
+    if (path) {
+        psList *list = psStringSplit(path, ":");
+        configPath = psListToArray(list);
+        psFree(list);
+    }
+
     // Next, we do a similar thing for the camera configuration file.  The
     // file is read and parsed into psMetadata struct "camera".
@@ -164 +234 @@
         } else {
             psArgumentRemove(argNum, argc, argv);
-            readConfig(camera, argv[argNum], "camera");
-        }
-    } else {
-        // XXX: Not sure is this is correct.
-        *camera = NULL;
+            readConfig(&config->camera, argv[argNum], "camera");
+        }
+    }
+
+    // Load the recipes from the camera file, if appropriate
+    if (! config->recipes && config->camera) {
+        pmConfigReadRecipes(config);
     }
 
@@ -183 +255 @@
     // with a call to psTimeInitialize.
     //
-    psString timeName = psMetadataLookupStr(&mdok, *site, "TIME");
+    psString timeName = psMetadataLookupStr(&mdok, config->site, "TIME");
     if (mdok && timeName) {
         psTrace(__func__, 7, "Initialising psTime with file %s\n", timeName);
@@ -195 +267 @@
     // with a call to psLogSetLevel().
     //
-    int logLevel = psMetadataLookupS32(&mdok, *site, "LOGLEVEL");
+    int logLevel = psMetadataLookupS32(&mdok, config->site, "LOGLEVEL");
     if (mdok && logLevel >= 0) {
         psTrace(__func__, 7, "Setting log level to %d\n", logLevel);
@@ -206 +278 @@
     // with a call to psLogSetFormat().
     //
-    psString logFormat = psMetadataLookupStr(&mdok, *site, "LOGFORMAT");
+    psString logFormat = psMetadataLookupStr(&mdok, config->site, "LOGFORMAT");
     if (mdok && logFormat) {
         psTrace(__func__, 7, "Setting log format to %s\n", logFormat);
@@ -218 +290 @@
     // XXX: This is not spec'ed in the SDRS.
     //
-    psString logDest = psMetadataLookupStr(&mdok, *site, "LOGDEST");
+    psString logDest = psMetadataLookupStr(&mdok, config->site, "LOGDEST");
     if (mdok && logDest) {
-        // XXX: Only stdout is provided for now; this section should be
+        // XXX: Only stdout and stderr are provided for now; this section should be
         // expanded in the future to do files, and perhaps even sockets.
-        if (strcasecmp(logDest, "STDOUT") != 0) {
-            psLogMsg(__func__, PS_LOG_WARN, "Only STDOUT is currently supported as a log destination.\n");
+        int logFD = STDIN_FILENO; // a known invalid value
+        if (!strcasecmp(logDest, "STDOUT")) {
+            logFD = STDOUT_FILENO;
+        }
+        if (!strcasecmp(logDest, "STDERR")) {
+            logFD = STDERR_FILENO;
+        }
+        if (logFD == STDIN_FILENO) {
+            psLogMsg(__func__, PS_LOG_WARN, "Only STDERR and STDOUT currently supported as a log destination.\n");
+            logFD = STDERR_FILENO;
         }
         psTrace(__func__, 7, "Setting log destination to STDOUT.\n");
-        // XXX: Use something other than "1"
-        psLogSetDestination(1);
-    }
-
+        psLogSetDestination(logFD);
+    }
 
     //
@@ -236 +314 @@
     // XXX: This is not spec'ed in the SDRS.
     //
-    psMetadata *trace = psMetadataLookupMD(&mdok, *site, "TRACE");
+    psMetadata *trace = psMetadataLookupMD(&mdok, config->site, "TRACE");
     if (mdok && trace) {
         psMetadataIterator *traceIter = psMetadataIteratorAlloc(trace, PS_LIST_HEAD, NULL); // Iterator
@@ -264 +342 @@
         psLogSetLevel(saveLogLevel);
     }
-    return(true);
-}
-
-bool pmConfigValidateCamera(
-    const psMetadata *camera,
+
+    return config;
+}
+
+
+bool pmConfigValidateCameraFormat(
+    const psMetadata *cameraFormat,
     const psMetadata *header)
 {
-    // Read the rule for that camera
+    // Read the rule for that camera format
     bool mdStatus = true;
-    psMetadata *rule = psMetadataLookupMD(&mdStatus, camera, "RULE");
+    psMetadata *rule = psMetadataLookupMD(&mdStatus, cameraFormat, "RULE");
     if (! mdStatus || ! rule) {
         psLogMsg(__func__, PS_LOG_WARN, "Unable to read rule for camera.\n");
@@ -282 +362 @@
     psMetadataIterator *ruleIter = psMetadataIteratorAlloc(rule, PS_LIST_HEAD, NULL); // Rule iterator
     psMetadataItem *ruleItem = NULL;    // Item from the metadata
-    bool match = true;                  // Does it match?
-    while ((ruleItem = psMetadataGetAndIncrement(ruleIter)) && match) {
+    while ((ruleItem = psMetadataGetAndIncrement(ruleIter))) {
         // Check for the existence of the rule
         psMetadataItem *headerItem = psMetadataLookup((psMetadata*)header, ruleItem->name);
         if (! headerItem || headerItem->type != ruleItem->type) {
-            match = false;
+            psFree(ruleIter);
+            return false;
         }
 
@@ -297 +377 @@
             if (strncmp(ruleItem->data.V, headerItem->data.V,
                         strlen(ruleItem->data.V)) != 0) {
-                match = false;
+                psFree(ruleIter);
+                return false;
             }
             break;
@@ -305 +386 @@
                     ruleItem->data.S32, headerItem->data.S32);
             if (ruleItem->data.S32 != headerItem->data.S32) {
-                match = false;
+                psFree(ruleIter);
+                return false;
             }
             break;
@@ -312 +394 @@
                     ruleItem->data.F32, headerItem->data.F32);
             if (ruleItem->data.F32 != headerItem->data.F32) {
-                match = false;
+                psFree(ruleIter);
+                return false;
             }
             break;
@@ -319 +402 @@
                     ruleItem->data.F64, headerItem->data.F64);
             if (ruleItem->data.F64 != headerItem->data.F64) {
-                match = false;
+                psFree(ruleIter);
+                return false;
             }
             break;
@@ -329 +413 @@
 
     psFree(ruleIter);
-
-    return match;
-}
-
-
-
-// Work out what camera we have, based on the FITS header and a set of
-// rules specified in the IPP configuration; return the camera configuration
-psMetadata *pmConfigCameraFromHeader(
-    const psMetadata *ipprc,            // The IPP configuration
-    const psMetadata *header)           // The FITS header
-{
-    bool mdStatus = false;  // Metadata lookup status
-    psMetadata *cameras = psMetadataLookupMD(&mdStatus, ipprc, "CAMERAS");
-    if (! mdStatus) {
-        psError(PS_ERR_IO, false, "Unable to find CAMERAS in the configuration.\n");
-        return NULL;
-    }
-    psMetadata *winner = NULL;       // The camera configuration whose rule first matches
-    //  the supplied header
-    // Iterate over the cameras
-    psMetadataIterator *iterator = psMetadataIteratorAlloc(cameras, PS_LIST_HEAD, NULL);
-    psMetadataItem *cameraItem = NULL; // Item from the metadata
-
-    while ((cameraItem = psMetadataGetAndIncrement(iterator))) {
-        // Open the camera information
-        psTrace(__func__, 3, "Inspecting camera %s (%s)\n", cameraItem->name,
-                cameraItem->comment);
-        psMetadata *camera = NULL; // The camera metadata
-        if (cameraItem->type == PS_DATA_METADATA) {
-            camera = psMemIncrRefCounter(cameraItem->data.md);
-        } else if (cameraItem->type == PS_DATA_STRING) {
-            psTrace(__func__, 5, "Reading camera configuration for %s...\n", cameraItem->name);
-            unsigned int badLines = 0;  // Number of bad lines in reading camera configuration
-            camera = psMetadataConfigParse(NULL, &badLines, cameraItem->data.V, true);
+    return true;
+}
+
+
+static bool formatFromHeader(psMetadata **format, // Format to return
+                             psMetadata *camera, // Camera configuration
+                             const psMetadata *header, // FITS header
+                             const char *cameraName // Name of camera
+                            )
+{
+    psMetadata *testFormat;
+
+    assert(format);
+    assert(camera);
+    assert(header);
+
+    bool result = false;                // Did we find the first match?
+
+    // Read the list of formats
+    bool mdok = true;                   // Status of MD lookup
+    psMetadata *formats = psMetadataLookupMD(&mdok, camera, "FORMATS"); // List of formats
+    if (!mdok || !formats) {
+        psLogMsg(__func__, PS_LOG_WARN, "Unable to find list of FORMATS in camera %s --- ignored\n",
+                 cameraName);
+        return false;
+    }
+    // Iterate over the formats
+    psMetadataIterator *formatsIter = psMetadataIteratorAlloc(formats, PS_LIST_HEAD, NULL);
+    psMetadataItem *formatsItem = NULL; // Item from formats
+    while ((formatsItem = psMetadataGetAndIncrement(formatsIter))) {
+        if (formatsItem->type != PS_DATA_STRING) {
+            psLogMsg(__func__, PS_LOG_WARN, "In camera %s, camera format %s is not of type STR --- "
+                     "ignored.\n", cameraName, formatsItem->name);
+            continue;
+        }
+        psTrace(__func__, 5, "Reading camera format for %s...\n", formatsItem->name);
+        // unsigned int badLines = 0;  // Number of bad lines in reading camera configuration
+        // Format to test against what we've got
+
+        if (!readConfig(&testFormat, formatsItem->data.V, formatsItem->name)) {
+            psLogMsg(__func__, PS_LOG_WARN, "trouble reading reading camera format %s\n", formatsItem->name);
+            psFree(testFormat);
+            continue;
+        }
+
+        # if (0)
+            psMetadata *testFormat = psMetadataConfigParse(NULL, &badLines, formatsItem->data.V, true);
+        if (badLines > 0) {
+            psLogMsg(__func__, PS_LOG_WARN, "%d bad lines encountered while reading camera"
+                     "format %s\n", badLines, formatsItem->name);
+        }
+        # endif
+
+        if (pmConfigValidateCameraFormat(testFormat, header)) {
+            if (!*format) {
+                psLogMsg(__func__, PS_LOG_INFO, "Camera %s, format %s matches header.\n", cameraName,
+                         formatsItem->name);
+                *format = psMemIncrRefCounter(testFormat);
+                result = true;
+            } else {
+                psLogMsg(__func__, PS_LOG_WARN, "Camera %s, format %s also matches header --- ignored.\n",
+                         cameraName, formatsItem->name);
+            }
+        }
+        psFree(testFormat);
+    }
+    psFree(formatsIter);
+
+    return result;
+}
+
+// Work out what camera we have, based on the FITS header and a set of rules specified in the IPP
+// configuration; return the camera configuration and format
+psMetadata *pmConfigCameraFormatFromHeader(
+    pmConfig *config,                   // The configuration
+    const psMetadata *header           // The FITS header
+)
+{
+    psMetadata *format = NULL;          // The winning format
+    bool mdok = false;                  // Metadata lookup status
+    psMetadata *testCamera = NULL;
+
+    // If we don't know what sort of camera we have, we try all that we know
+    if (! config->camera) {
+        psMetadata *cameras = psMetadataLookupMD(&mdok, config->site, "CAMERAS");
+        if (! mdok) {
+            psError(PS_ERR_IO, false, "Unable to find CAMERAS in the configuration.\n");
+            return false;
+        }
+
+        // Iterate over the cameras
+        psMetadataIterator *camerasIter = psMetadataIteratorAlloc(cameras, PS_LIST_HEAD, NULL);
+        psMetadataItem *camerasItem = NULL; // Item from the metadata
+        while ((camerasItem = psMetadataGetAndIncrement(camerasIter))) {
+            // Open the camera information
+            psTrace(__func__, 3, "Inspecting camera %s (%s)\n", camerasItem->name, camerasItem->comment);
+            if (camerasItem->type != PS_DATA_STRING) {
+                psLogMsg(__func__, PS_LOG_WARN, "Camera configuration for %s in CAMERAS is not of type STR "
+                         "--- ignored.\n", camerasItem->name);
+                continue;
+            }
+
+            psTrace(__func__, 5, "Reading camera configuration for %s...\n", camerasItem->name);
+            // unsigned int badLines = 0;  // Number of bad lines in reading camera configuration
+            // Camera to test against what we've got:
+
+            if (!readConfig(&testCamera, camerasItem->data.V, camerasItem->name)) {
+                psLogMsg(__func__, PS_LOG_WARN, "trouble reading reading camera configuration %s\n", camerasItem->name);
+                psFree(testCamera);
+                continue;
+            }
+
+            # if (0)
+                psMetadata *testCamera = psMetadataConfigParse(NULL, &badLines, camerasItem->data.V, true);
             if (badLines > 0) {
                 psLogMsg(__func__, PS_LOG_WARN, "%d bad lines encountered while reading camera"
-                         "configuration %s\n", badLines, cameraItem->name);
-            }
-        }
-
-        if (! camera) {
-            psLogMsg(__func__, PS_LOG_WARN, "Unable to interpret camera configuration for %s (%s)\n",
-                     cameraItem->name, cameraItem->comment);
-            continue;
-        }
-
-        if (pmConfigValidateCamera(camera, header)) {
-            if (! winner) {
-                // This is the first match
-                winner = psMemIncrRefCounter(camera);
-                psLogMsg(__func__, PS_LOG_INFO, "FITS header matches camera %s\n",
-                         cameraItem->name);
-            } else {
-                // We have a duplicate match
-                psLogMsg(__func__, PS_LOG_WARN, "Additional camera found that matches the rules: %s\n",
-                         cameraItem->name);
-            }
-        } // Done inspecting the camera
-
-        psFree(camera);
-
-    } // Done looking at all cameras
-    if (! winner) {
-        psError(PS_ERR_IO, true, "Unable to find an camera that matches input FITS header!\n");
-    }
-
-    psFree(iterator);
-    return winner;
-}
-
-psMetadata *pmConfigRecipeFromCamera(
-    const psMetadata *camera,
-    const char *recipeName)
-{
-    PS_ASSERT_PTR_NON_NULL(camera, false);
-    PS_ASSERT_PTR_NON_NULL(recipeName, false);
-
-    psMetadata *recipe = NULL;          // Recipe to read
+                         "configuration %s\n", badLines, camerasItem->name);
+            }
+            # endif
+
+            if (! testCamera) {
+                psLogMsg(__func__, PS_LOG_WARN, "Unable to interpret camera configuration for %s (%s) --- "
+                         "ignored\n", camerasItem->name, camerasItem->comment);
+                continue;
+            }
+
+            if (formatFromHeader(&format, testCamera, header, camerasItem->name)) {
+                config->camera = psMemIncrRefCounter(testCamera);
+            }
+            psFree(testCamera);
+        } // Done looking at all cameras
+        psFree(camerasIter);
+
+        if (! config->camera) {
+            psError(PS_ERR_IO, true, "Unable to find a camera that matches input FITS header!\n");
+            return NULL;
+        }
+
+        // Now we have the camera, we can read the recipes
+        if (!config->recipes) {
+            pmConfigReadRecipes(config);
+        }
+
+        return format;
+    }
+
+    // Otherwise, try the specific camera
+    if (! formatFromHeader(&format, config->camera, header, "specified camera")) {
+        psError(PS_ERR_IO, true, "Unable to find a format with the specified camera that matches the "
+                "given header.\n");
+        return NULL;
+    }
+    return format;
+}
+
+bool pmConfigReadRecipes(
+    pmConfig *config
+)
+{
+    PS_ASSERT_PTR_NON_NULL(config, false);
+    PS_ASSERT_PTR_NON_NULL(config->camera, false);
+
+    if (!config->recipes) {
+        config->recipes = psMetadataAlloc();
+    }
+
     bool mdok = true;                   // Status of MD lookup
-    psMetadata *recipes = psMetadataLookupMD(&mdok, camera, "RECIPES"); // The list of recipes
+    psMetadata *recipes = psMetadataLookupMD(&mdok, config->camera, "RECIPES"); // The list of recipes
     if (! mdok || ! recipes) {
         psLogMsg(__func__, PS_LOG_WARN, "RECIPES in the camera configuration file is not of type METADATA\n");
-    } else {
-        psString recipeFileName = psMetadataLookupStr(&mdok, recipes, recipeName);
-        (void)readConfig(&recipe, recipeFileName, "recipe");
-    }
-
-    return recipe;
+        return false;
+    }
+    // Go through the recipes and load each one
+    psMetadataIterator *recipesIter = psMetadataIteratorAlloc(recipes, PS_LIST_HEAD, NULL); // Iterator
+    psMetadataItem *recipesItem = NULL; // Item from iteration
+    while ((recipesItem = psMetadataGetAndIncrement(recipesIter))) {
+        if (recipesItem->type != PS_DATA_STRING) {
+            psLogMsg(__func__, PS_LOG_WARN, "Filename for recipe %s isn't of type STR --- ignored.\n",
+                     recipesItem->name);
+            continue;
+        }
+
+        psMetadata *recipe = NULL;      // Recipe from file
+        if (readConfig(&recipe, recipesItem->data.V, "recipe")) {
+            psMetadataAdd(config->recipes, PS_LIST_TAIL, recipesItem->name,
+                          PS_DATA_METADATA | PS_META_REPLACE, recipesItem->comment, recipe);
+        }
+        psFree(recipe);                 // Drop reference
+    }
+    psFree(recipesIter);
+
+    return true;
 }
 
@@ -423 +610 @@
 pmConfigDB(*site)
  
+XXX: this should allow the option of having NO database server, if chosen by config
 XXX: What should we use for the Database namespace in the call to psDBInit()?
 This is currently NULL.
  *****************************************************************************/
 
-#ifdef OMIT_PSDB
-psDB *pmConfigDB(psMetadata *site)
-{
-    psError(PS_ERR_UNKNOWN, true, "pslib was built without psDB support");
-    return NULL;
-}
-#else
 psDB *pmConfigDB(
     psMetadata *site)
@@ -442 +623 @@
     psBool mdStatus03 = false;
 
+    // XXX leaky strings
     psString dbServer = psMetadataLookupStr(&mdStatus01, site, "DBSERVER");
     psString dbUsername = psMetadataLookupStr(&mdStatus02, site, "DBUSER");
     psString dbPassword = psMetadataLookupStr(&mdStatus03, site, "DBPASSWORD");
     psString dbName = psMetadataLookupStr(&mdStatus01, site, "DBNAME");
-
     if (!(mdStatus01 & mdStatus02 & mdStatus03)) {
         psLogMsg(__func__, PS_LOG_WARN, "Could not determine database server name, userID, and password from site metadata.\n");
-        return NULL;
-    }
-
-
-
-    psDB *dbh = psDBInit(dbServer, dbUsername, dbPassword, dbName);
-    psFree(dbServer);
-    psFree(dbUsername);
-    psFree(dbPassword);
-    psFree(dbName);
-
-    if (!dbh) {
-        psError(PS_ERR_UNKNOWN, false, "database connection failed");
-    }
-
-    return dbh;
-}
-#endif
+        return(NULL);
+    }
+
+    return(psDBInit(dbServer, dbUsername, dbPassword, dbName));
+}
+
+
+bool pmConfigConformHeader(psMetadata *header, // Header to conform
+                           const psMetadata *format // Camera format
+                          )
+{
+    bool mdok = true;                   // Status of MD lookup
+    psMetadata *rules = psMetadataLookupMD(&mdok, format, "RULE"); // How to identify this format
+    if (!mdok || !rules) {
+        psError(PS_ERR_IO, false, "Unable to find RULE in camer format.\n");
+        return false;
+    }
+
+    psMetadataIterator *rulesIter = psMetadataIteratorAlloc(rules, PS_LIST_HEAD, NULL); // Iterator for rules
+    psMetadataItem *rulesItem = NULL;   // Item from iteration
+    while ((rulesItem = psMetadataGetAndIncrement(rulesIter))) {
+        psMetadataItem *newItem = psMetadataItemCopy(rulesItem); // Copy of item
+        psMetadataAddItem(header, newItem, PS_LIST_TAIL, PS_META_REPLACE);
+        psFree(newItem);                // Drop reference
+    }
+    psFree(rulesIter);
+
+    return true;
+}
+
+// given the 'file' and 'list' words, find the arguments associated with these words
+// and interpret them as lists of files.  return an array of the resulting filenames.
+psArray *pmConfigFileSets (int *argc, char **argv, char *file, char *list)
+{
+
+    int Narg;
+
+    // we load all input files onto a psArray, to be parsed later
+    psArray *input = psArrayAlloc (16);
+    input->n = 0;
+
+    // load the list of filenames the supplied file (may be a glob: "file*.fits")
+    if ((Narg = psArgumentGet (*argc, argv, file))) {
+        glob_t globList;
+        psArgumentRemove (Narg, argc, argv);
+        globList.gl_offs = 0;
+        glob (argv[Narg], 0, NULL, &globList);
+        for (int i = 0; i < globList.gl_pathc; i++) {
+            char *filename = psStringCopy (globList.gl_pathv[i]);
+            psArrayAdd (input, 16, filename);
+            psFree (filename);
+        }
+        psArgumentRemove (Narg, argc, argv);
+    }
+
+    // load the list from the supplied text file
+    if ((Narg = psArgumentGet (*argc, argv, list))) {
+        int nItems;
+        char line[1024]; // XXX limits the list lines to 1024 chars
+        char word[1024];
+        char *filename;
+
+        psArgumentRemove (Narg, argc, argv);
+        FILE *f = fopen (argv[Narg], "r");
+        if (f == NULL) {
+            psAbort ("psphot", "unable to open specified list file");
+        }
+        while (fgets (line, 1024, f) != NULL) {
+            nItems = sscanf (line, "%s", word);
+            switch (nItems) {
+            case 0:
+                break;
+            case 1:
+                filename = psStringCopy (word);
+                psArrayAdd (input, 16, filename);
+                psFree (filename);
+                break;
+            default:
+                // rigid format, no comments allowed?
+                psAbort ("pmConfig", "error parsing input list file");
+                break;
+            }
+        }
+        psArgumentRemove (Narg, argc, argv);
+    }
+
+    return input;
+}
+
+bool pmConfigFileSetsMD (psMetadata *metadata, int *argc, char **argv, char *name, char *file, char *list)
+{
+
+    psArray *files = pmConfigFileSets (argc, argv, file, list);
+    if (files->n == 0) {
+        psFree (files);
+        return false;
+    }
+
+    psMetadataAddPtr (metadata, PS_LIST_TAIL, name,  PS_DATA_ARRAY, "", files);
+    psFree (files);
+    return true;
+}

trunk/psModules/src/config/pmConfig.h

@@ -3 +3 @@
  *  @author PAP, IfA
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/detrend/pmFlatField.c

@@ -24 +24 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/detrend/pmMaskBadPixels.c

@@ -24 +24 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/detrend/pmMaskBadPixelsErrors.h

@@ -12 +12 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/detrend/pmNonLinear.c

@@ -10 +10 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/imcombine/pmReadoutCombine.c

@@ -5 +5 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-04 01:01:33 $
+ *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -36 +36 @@
 }
 
+static void pmCombineParamsFree (pmCombineParams *params)
+{
+
+    if (params == NULL)
+        return;
+
+    psFree (params->stats);
+    return;
+}
+
+pmCombineParams *pmCombineParamsAlloc (psStatsOptions statsOptions)
+{
+
+    pmCombineParams *params = psAlloc (sizeof(pmCombineParams));
+    psMemSetDeallocator(params, (psFreeFunc) pmCombineParamsFree);
+
+    params->stats = psStatsAlloc (statsOptions);
+    params->maskVal = 0;
+    params->fracHigh = 0.25;
+    params->fracHigh = 0.25;
+    params->nKeep = 3;
+
+    return (params);
+}
+
 /******************************************************************************
 XXX: Must add support for S16 and S32 types.  F32 currently supported.
  *****************************************************************************/
 psImage *pmReadoutCombine(psImage *output,
-                          const psList *inputs,
-                          psCombineParams *params,
+                          const psArray *inputs,
                           const psVector *zero,
                           const psVector *scale,
+                          pmCombineParams *params,
                           bool applyZeroScale,
                           psF32 gain,
                           psF32 readnoise)
+{
+    PS_ASSERT_PTR_NON_NULL(inputs, NULL);
+    PS_ASSERT_PTR_NON_NULL(params, NULL);
+    PS_ASSERT_PTR_NON_NULL(params->stats, NULL);
+    if (zero != NULL) {
+        PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL);
+        //        PS_ASSERT_VECTOR_TYPE_S16_S32_F32(zero, NULL);
+    }
+    if (scale != NULL) {
+        PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL);
+        //        PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL);
+    }
+    if ((zero != NULL) && (scale != NULL)) {
+        PS_ASSERT_VECTOR_TYPE_EQUAL(zero, scale, NULL);
+        // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL);
+    }
+
+    psStats *stats = params->stats;
+    psS32 maxInputCols = 0;
+    psS32 maxInputRows = 0;
+    psS32 minInputCols = PS_MAX_S32;
+    psS32 minInputRows = PS_MAX_S32;
+    pmReadout *tmpReadout = NULL;
+    psS32 tmpI;
+    psElemType outputType = PS_TYPE_F32;
+
+    if (DetermineNumBits(stats->options) != 1) {
+        psError(PS_ERR_UNKNOWN, true,
+                "Multiple statistical options have been requested.  Returning NULL.\n");
+        return(NULL);
+    }
+
+    // We step through each readout in the input image list.  If any readout is
+    // NULL, empty, or has the wrong type, we generate an error and return
+    // NULL.  We determine how big of an output image is needed to combine
+    // these input images.  We do this by taking the
+    //     max(readout->col0 + readout->numCols + image->col0 + image->numCols)
+    //     max(readout->row0 + readout->numRows + image->row0 + image->numRows)
+    //
+    for (int i = 0; i < inputs->n; i++) {
+        tmpReadout = inputs->data[i];
+        PS_ASSERT_READOUT_NON_NULL(tmpReadout, output);
+        PS_ASSERT_READOUT_NON_EMPTY(tmpReadout, output);
+        PS_ASSERT_READOUT_TYPE(tmpReadout, PS_TYPE_F32, output);
+
+        minInputRows = PS_MIN(minInputRows, (tmpReadout->row0 + tmpReadout->image->row0));
+        tmpI = tmpReadout->row0 +
+               tmpReadout->image->row0 +
+               tmpReadout->image->numRows;
+        maxInputRows = PS_MAX(maxInputRows, tmpI);
+
+        minInputCols = PS_MIN(minInputCols, (tmpReadout->col0 + tmpReadout->image->col0));
+        tmpI = tmpReadout->col0 +
+               tmpReadout->image->col0 +
+               tmpReadout->image->numCols;
+        maxInputCols = PS_MAX(maxInputCols, tmpI);
+    }
+
+    // We ensure that the zero vector is of the proper size.
+    if (zero != NULL) {
+        PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL);
+        if (zero->n < inputs->n) {
+            psError(PS_ERR_UNKNOWN, true, "zero vector has incorrect size (%d).  Returning NULL.\n", zero->n);
+            return(NULL);
+        } else if (zero->n > inputs->n) {
+            // XXX EAM : abort on this condition? is probably an error
+            psLogMsg(__func__, PS_LOG_WARN,
+                     "WARNING: the zero vector too many elements (%d)\n", zero->n);
+        }
+    }
+
+    // We ensure that the scale vector is of the proper size.
+    if (scale != NULL) {
+        PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL);
+        if (scale->n < inputs->n) {
+            psError(PS_ERR_UNKNOWN, true, "scale vector has incorrect size (%d).  Returning NULL.\n", scale->n);
+            return(NULL);
+        } else if (scale->n > inputs->n) {
+            // XXX EAM : abort on this condition? is probably an error
+            psLogMsg(__func__, PS_LOG_WARN,
+                     "WARNING: the scale vector has too many elements (%d)\n", scale->n);
+        }
+    }
+
+    // At this point, the following variables have been computed:
+    // maxInputRows: the largest input row value, in output image space.
+    // maxInputCols: the largest input column value, in output image space.
+    // minInputRows: the smallest input row value, in output image space.
+    // minInputCols: the smallest input column value, in output image space.
+    //
+    if (output == NULL) {
+        output = psImageAlloc(maxInputCols-minInputCols, maxInputRows-minInputRows, outputType);
+        *(psS32 *) &(output->col0) = minInputCols;
+        *(psS32 *) &(output->row0) = minInputRows;
+    } else {
+
+        // XXX EAM : recycle the existing output image?  why would we care about the existing pixels?
+        PS_ASSERT_IMAGE_TYPE(output, PS_TYPE_F32, NULL);
+        if (((output->col0 + output->numCols) < maxInputCols) ||
+                ((output->row0 + output->numRows) < maxInputRows)) {
+            psError(PS_ERR_UNKNOWN, true,
+                    "Output image (%d, %d) is too small to hold combined images.  Returning NULL.\n",
+                    output->row0 + output->numRows,
+                    output->col0 + output->numCols);
+            return(NULL);
+        }
+
+        // reset output origin using logic of above
+        *(psS32 *) &(output->col0) = minInputCols;
+        *(psS32 *) &(output->row0) = minInputRows;
+    }
+
+    psVector *tmpPixels     = psVectorAlloc(inputs->n, PS_TYPE_F32);
+    psVector *tmpPixelsKeep = psVectorAlloc(inputs->n, PS_TYPE_F32);
+    psVector *outRowLower   = psVectorAlloc(inputs->n, PS_TYPE_U32);
+    psVector *outRowUpper   = psVectorAlloc(inputs->n, PS_TYPE_U32);
+    psVector *outColLower   = psVectorAlloc(inputs->n, PS_TYPE_U32);
+    psVector *outColUpper   = psVectorAlloc(inputs->n, PS_TYPE_U32);
+
+    // For each input readout, we store the min/max pixel indices for that readout, in detector coordinates,
+    // in the psVectors (outRowLower, outColLower, outRowUpper, outColUpper).
+    for (int i = 0; i < inputs->n; i++) {
+        tmpReadout = (pmReadout *) inputs->data[i];
+        outRowLower->data.U32[i] = tmpReadout->row0 + tmpReadout->image->row0;
+        outColLower->data.U32[i] = tmpReadout->col0 + tmpReadout->image->col0;
+        outRowUpper->data.U32[i] = tmpReadout->row0 +
+                                   tmpReadout->image->row0 +
+                                   tmpReadout->image->numRows;
+        outColUpper->data.U32[i] = tmpReadout->col0 +
+                                   tmpReadout->image->col0 +
+                                   tmpReadout->image->numCols;
+    }
+
+    // We loop through each pixel in the output image.  We loop through each
+    // input readout.  We determine if that output pixel is contained in the
+    // image from that readout.  If so, we save it in psVector tmpPixels.
+    // If not, we set a mask for that element in tmpPixels.  Then, we mask off
+    // pixels not between fracLow and fracHigh.  Then we call the vector
+    // stats routine on those pixels/mask.  Then we set the output pixel value
+    // to the result of the stats call.
+
+    int nx, ny;
+    int nKeep, nMin;
+    float keepFrac = 1.0 - params->fracLow - params->fracHigh;
+    float value = 0;
+    psF32 *saveVector = tmpPixelsKeep->data.F32;
+
+    for (int j = output->row0; j < (output->row0 + output->numRows) ; j++) {
+        if (j % 10 == 0)
+            fprintf (stderr, ".");
+        for (int i = output->col0; i < (output->col0 + output->numCols) ; i++) {
+            int nPix = 0;
+            for (int r = 0; r < inputs->n; r++) {
+                tmpReadout = (pmReadout *) inputs->data[r];
+
+                // psTrace (__func__, 6, "[%d][%d]: [%d][%d] to [%d][%d]\n", i, j, outColLower->data.U32[r], outRowLower->data.U32[r], outColUpper->data.U32[r], outRowUpper->data.U32[r]);
+                if (i <  outColLower->data.U32[r])
+                    continue;
+                if (i >= outColUpper->data.U32[r])
+                    continue;
+                if (j <  outRowLower->data.U32[r])
+                    continue;
+                if (j >= outRowUpper->data.U32[r])
+                    continue;
+
+                nx = i - (tmpReadout->col0 + tmpReadout->image->col0);
+                ny = j - (tmpReadout->row0 + tmpReadout->image->row0);
+
+                if (tmpReadout->mask != NULL) {
+                    if (tmpReadout->mask->data.U8[ny][nx] && params->maskVal)
+                        continue;
+                }
+
+                tmpPixels->data.F32[nPix] = tmpReadout->image->data.F32[ny][nx];
+                // psTrace (__func__, 6, "readout[%d], image [%d][%d] is %f\n", r, i, j, tmpPixels->data.F32[nPix]);
+                nPix ++;
+            }
+            tmpPixels->n = nPix;
+
+            // are there enough valid pixels to apply fracLow,fracHigh?
+            nKeep = nPix * keepFrac;
+            if (nKeep >= params->nKeep) {
+                psVectorSort (tmpPixels, tmpPixels);
+                nMin = nPix * params->fracLow;
+                tmpPixelsKeep->data.F32 = &tmpPixels->data.F32[nMin];
+                tmpPixelsKeep->n = nKeep;
+            } else {
+                tmpPixelsKeep->data.F32 = tmpPixels->data.F32;
+                tmpPixelsKeep->n = nPix;
+            }
+
+            // tmpPixelsKeep is already sorted.  sample mean and median are very easy
+            if (stats->options & PS_STAT_SAMPLE_MEAN) {
+                value = 0;
+                for (int r = 0; r < tmpPixelsKeep->n; r++) {
+                    value += tmpPixelsKeep->data.F32[r];
+                }
+                if (tmpPixelsKeep->n == 0) {
+                    value = 0;
+                } else {
+                    value = value / tmpPixelsKeep->n;
+                }
+            }
+            if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
+                int r = tmpPixelsKeep->n / 2;
+                if (tmpPixelsKeep->n == 0) {
+                    value = 0;
+                    goto got_value;
+                }
+                if (tmpPixelsKeep->n % 2 == 1) {
+                    int r = 0.5*tmpPixelsKeep->n;
+                    value = tmpPixelsKeep->data.F32[r];
+                    goto got_value;
+                }
+                if (tmpPixelsKeep->n % 2 == 0) {
+                    value = 0.5*(tmpPixelsKeep->data.F32[r] +
+                                 tmpPixelsKeep->data.F32[r-1]);
+                    goto got_value;
+                }
+            }
+got_value:
+            output->data.F32[j-output->row0][i-output->col0] = value;
+        }
+    }
+    tmpPixelsKeep->data.F32 = saveVector;
+
+    psFree(tmpPixels);
+    psFree(tmpPixelsKeep);
+    psFree(outRowLower);
+    psFree(outRowUpper);
+    psFree(outColLower);
+    psFree(outColUpper);
+
+    return(output);
+}
+
+/******************************************************************************
+XXX: Must add support for S16 and S32 types.  F32 currently supported.
+ *****************************************************************************/
+psImage *pmReadoutCombine_OLD(psImage *output,
+                              const psList *inputs,
+                              pmCombineParams *params,
+                              const psVector *zero,
+                              const psVector *scale,
+                              bool applyZeroScale,
+                              psF32 gain,
+                              psF32 readnoise)
 {
     PS_ASSERT_PTR_NON_NULL(inputs, NULL);
@@ -418 +690 @@
     psRegion minRegion;
     psRegion maxRegion;
-    psStats *minStats = psStatsAlloc(PS_STAT_FITTED_MEAN);
-    psStats *maxStats = psStatsAlloc(PS_STAT_FITTED_MEAN);
+    psStats *minStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+    psStats *maxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     psStats *diffStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     psVector *diffs = psVectorAlloc(fringePoints->n, PS_TYPE_F32);
@@ -454 +726 @@
         }
 
-        fp->midValue = 0.5 * (maxStats->fittedMean + minStats->fittedMean);
-        fp->delta = maxStats->fittedMean - minStats->fittedMean;
+        fp->midValue = 0.5 * (maxStats->robustMedian + minStats->robustMedian);
+        fp->delta = maxStats->robustMedian - minStats->robustMedian;
         diffs->data.F32[i] = fp->delta;
     }
@@ -464 +736 @@
     psFree(diffs);
     if (diffStats == NULL) {
-        psError(PS_ERR_UNKNOWN, true, "Could not determine fitted median of the differences.\n");
+        psError(PS_ERR_UNKNOWN, true, "Could not determine robust median of the differences.\n");
         return(NULL);
     }
     return(diffStats);
 }
-
-
 
 /**

trunk/psModules/src/imcombine/pmReadoutCombine.h

@@ -5 +5 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/imsubtract/pmSubtractBias.c

@@ +1 @@
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
+// one that was being worked on.
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
 /** @file  pmSubtractBias.c
  *
@@ -6 +11 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-04 01:01:33 $
+ *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  *
  */
-/*****************************************************************************/
-/* INCLUDE FILES                                                             */
-/*****************************************************************************/
-#include <stdio.h>
-#include <math.h>
-#include <string.h>
-#include "pslib.h"
+
 #if HAVE_CONFIG_H
 #include <config.h>
 #endif
+
+#include <assert.h>
 #include "pmSubtractBias.h"
 
-/*****************************************************************************/
-/* DEFINE STATEMENTS                                                         */
-/*****************************************************************************/
+#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2
+#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3
+
+
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+
+
 // XXX: put these in psConstants.h
-void PS_POLY1D_PRINT(
-    psPolynomial1D *poly)
+void PS_POLY1D_PRINT(psPolynomial1D *poly)
 {
     printf("-------------- PS_POLY1D_PRINT() --------------\n");
@@ -57 +62 @@
 }\
 
-/*****************************************************************************/
-/* TYPE DEFINITIONS                                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* GLOBAL VARIABLES                                                          */
-/*****************************************************************************/
-psS32 currentId = 0;                // XXX: remove
-psS32 memLeaks = 0;                 // XXX: remove
-//PRINT_MEMLEAKS(8); XXX
-/*****************************************************************************/
-/* FILE STATIC VARIABLES                                                     */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* FUNCTION IMPLEMENTATION - LOCAL                                           */
-/*****************************************************************************/
+
+void overscanOptionsFree(pmOverscanOptions *options)
+{
+    psFree(options->stat);
+    psFree(options->poly);
+    psFree(options->spline);
+}
+
+pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat)
+{
+    pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions));
+    psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree);
+
+    // Inputs
+    opts->single = single;
+    opts->fitType = fitType;
+    opts->order = order;
+    opts->stat = psMemIncrRefCounter(stat);
+
+    // Outputs
+    opts->poly = NULL;
+    opts->spline = NULL;
+
+    return opts;
+}
+
 
 /******************************************************************************
-psSubtractFrame(): this routine will take as input the pmReadout for the input
-image and a pmReadout for the bias image.  The bias image is subtracted in
-place from the input image.  We assume that sizes and types are checked
-elsewhere.
- 
-XXX: Verify that the image and readout offsets are being used the right way.
- 
-XXX: Ensure that it does the correct thing with image size.
+psSubtractFrame(): this routine will take as input a readout for the input
+image and a readout for the bias image.  The bias image is subtracted in
+place from the input image.
 *****************************************************************************/
-static pmReadout *SubtractFrame(
-    pmReadout *in,
-    const pmReadout *bias)
-{
-    // XXX: When did the ->row0 and ->col0 offsets get coded?
-    for (psS32 i=0;i<in->image->numRows;i++) {
-        for (psS32 j=0;j<in->image->numCols;j++) {
-            in->image->data.F32[i][j]-=
-                bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0];
-
-            if ((in->mask != NULL) && (bias->mask != NULL)) {
-                (in->mask->data.U8[i][j])|=
-                    bias->mask->data.U8[i+in->row0-bias->row0][j+in->col0-bias->col0];
+static bool SubtractFrame(pmReadout *in,// Input readout
+                          const pmReadout *sub, // Readout to be subtracted from input
+                          float scale   // Scale to apply before subtracting
+                         )
+{
+    assert(in);
+    assert(sub);
+
+    psImage *inImage  = in->image;      // The input image
+    psImage *inMask   = in->mask;       // The input mask
+    psImage *subImage = sub->image;     // The image to be subtracted
+    psImage *subMask  = sub->mask;      // The mask for the subtraction image
+
+    // Offsets of the cells
+    int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0");
+    int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0");
+    int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0");
+    int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0");
+
+    if ((inImage->numCols + x0in - x0sub) > subImage->numCols) {
+        psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction.\n");
+        return false;
+    }
+    if ((inImage->numRows + y0in - y0sub) > subImage->numRows) {
+        psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction.\n");
+        return false;
+    }
+
+    if (scale == 1.0) {
+        for (int i = 0; i < inImage->numRows; i++) {
+            for (int j = 0; j < inImage->numCols; j++) {
+                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub];
+                if (inMask && subMask) {
+                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
+                }
             }
         }
-    }
-
-    return(in);
-}
-
-
-/******************************************************************************
-psSubtractDarkFrame(): this routine will take as input the pmReadout for the
-input image and a pmReadout for the dark image.  The dark image is scaled and
-subtracted in place from the input image.
- 
-XXX: Verify that the image and readout offsets are being used the right way.
- 
-XXX: Ensure that it does the correct thing with image size.
-*****************************************************************************/
-static pmReadout *SubtractDarkFrame(
-    pmReadout *in,
-    const pmReadout *dark,
-    psF32 scale)
-{
-    // XXX: When did the ->row0 and ->col0 offsets get coded?
-    if (fabs(scale) > FLT_EPSILON) {
-        for (psS32 i=0;i<in->image->numRows;i++) {
-            for (psS32 j=0;j<in->image->numCols;j++) {
-                in->image->data.F32[i][j]-=
-                    (scale * dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0]);
-
-                if ((in->mask != NULL) && (dark->mask != NULL)) {
-                    (in->mask->data.U8[i][j])|=
-                        dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0];
+    } else {
+        for (int i = 0; i < inImage->numRows; i++) {
+            for (int j = 0; j < inImage->numCols; j++) {
+                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale;
+                if (inMask && subMask) {
+                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
                 }
             }
         }
-    } else {
-        for (psS32 i=0;i<in->image->numRows;i++) {
-            for (psS32 j=0;j<in->image->numCols;j++) {
-                in->image->data.F32[i][j]-=
-                    dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0];
-
-                if ((in->mask != NULL) && (dark->mask != NULL)) {
-                    (in->mask->data.U8[i][j])|=
-                        dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0];
-                }
-            }
-        }
-    }
-
-    return(in);
-}
-
+    }
+
+    return true;
+}
+
+
+#if 0
 /******************************************************************************
 ImageSubtractScalar(): subtract a scalar from the input image.
  
-XXX: Is there a psLib function for this?
+XXX: Use a psLib function for this.
+ 
+XXX: This should
  *****************************************************************************/
-static psImage *ImageSubtractScalar(
-    psImage *image,
-    psF32 scalar)
+static psImage *ImageSubtractScalar(psImage *image,
+                                    psF32 scalar)
 {
     for (psS32 i=0;i<image->numRows;i++) {
@@ -166 +164 @@
     return(image);
 }
+#endif
 
 /******************************************************************************
@@ -179 +178 @@
     psStatsOptions opt = 0;
 
-    /*
-        if (stat->options & PS_STAT_ROBUST_MODE) {
-            if (numOptions == 0) {
-                opt = PS_STAT_ROBUST_MODE;
-            }
-            numOptions++;
-        }
-    */
     if (stat->options & PS_STAT_ROBUST_MEDIAN) {
         if (numOptions == 0) {
@@ -194 +185 @@
     }
 
-    if (stat->options & PS_STAT_FITTED_MEAN) {
-        if (numOptions == 0) {
-            opt = PS_STAT_FITTED_MEAN;
-        }
-        numOptions++;
-    }
-
     if (stat->options & PS_STAT_CLIPPED_MEAN) {
         if (numOptions == 0) {
@@ -222 +206 @@
 
     if (numOptions == 0) {
-        psError(PS_ERR_UNKNOWN,true, "No allowable statistics options have been specified.\n");
+        psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n");
     }
     if (numOptions != 1) {
@@ -231 +215 @@
 }
 
-/******************************************************************************
-Polynomial1DCopy(): This private function copies the members of the existing
-psPolynomial1D "in" into the existing psPolynomial1D "out".  The previous
-members of the existing psPolynomial1D "out" are psFree'ed.
- *****************************************************************************/
-static psBool Polynomial1DCopy(
-    psPolynomial1D *out,
-    psPolynomial1D *in)
-{
-    psFree(out->coeff);
-    psFree(out->coeffErr);
-    psFree(out->mask);
-
-    out->type = in->type;
-    out->nX = in->nX;
-
-    out->coeff = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64));
-    // XXX: use memcpy
-    for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
-        out->coeff[i] = in->coeff[i];
-    }
-
-    out->coeffErr = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64));
-    // XXX: use memcpy
-    for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
-        out->coeffErr[i] = in->coeffErr[i];
-    }
-
-    out->mask = (psMaskType *) psAlloc((in->nX + 1) * sizeof(psMaskType));
-    // XXX: use memcpy
-    for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
-        out->mask[i] = in->mask[i];
-    }
-
-    return(true);
-}
-
-/******************************************************************************
-Polynomial1DDup(): This private function duplicates and then returns the input
-psPolynomial1D "in".
- *****************************************************************************/
-static psPolynomial1D *Polynomial1DDup(
-    psPolynomial1D *in)
-{
-    psPolynomial1D *out = psPolynomial1DAlloc(in->type, in->nX);
-    Polynomial1DCopy(out, in);
-    return(out);
-}
-
-
-/******************************************************************************
-SplineCopy(): This private function copies the members of the existing
-psSpline in into the existing psSpline out.
- *****************************************************************************/
-static psBool SplineCopy(
-    psSpline1D *out,
-    psSpline1D *in)
-{
-    PS_ASSERT_PTR_NON_NULL(out, false);
-    PS_ASSERT_PTR_NON_NULL(in, false);
-
-    for (psS32 i = 0 ; i < out->n ; i++) {
-        psFree(out->spline[i]);
-    }
-    psFree(out->spline);
-    psFree(out->knots);
-    psFree(out->p_psDeriv2);
-
-    out->n = in->n;
-    out->spline = (psPolynomial1D **) psAlloc(in->n * sizeof(psPolynomial1D *));
-    for (psS32 i = 0 ; i < in->n ; i++) {
-        out->spline[i] = Polynomial1DDup(in->spline[i]);
-    }
-
-    // XXX: use psVectorCopy if they get it working.
-    out->knots = psVectorAlloc(in->knots->n, in->knots->type.type);
-    out->knots->n = out->knots->nalloc;
-    for (psS32 i = 0 ; i < in->knots->n ; i++) {
-        out->knots->data.F32[i] = in->knots->data.F32[i];
-    }
-    /*
-        out->knots = psVectorCopy(out->knots, in->knots, in->knots->type.type);
-    */
-
-    out->p_psDeriv2 = (psF32 *) psAlloc((in->n + 1) * sizeof(psF32));
-    // XXX: use memcpy
-    for (psS32 i = 0 ; i < (in->n + 1) ; i++) {
-        out->p_psDeriv2[i] = in->p_psDeriv2[i];
-    }
-
-    return(true);
-}
-
+
+
+#if 0
 /******************************************************************************
 ScaleOverscanVector(): this routine takes as input an arbitrary vector,
@@ -330 +224 @@
     PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data.
     PM_FIT_SPLINE: fit splines to the input vector data.
-The resulting spline or polynomial is set in the fitSpec argument.
+XXX: Doesn't it make more sense to do polynomial interpolation on a few
+elements of the input vector, rather than fit a polynomial to the entire
+vector?
  *****************************************************************************/
-static psVector *ScaleOverscanVector(
-    psVector *overscanVector,
-    psS32 n,
-    void *fitSpec,
-    pmFit fit)
+static psVector *ScaleOverscanVector(psVector *overscanVector,
+                                     psS32 n,
+                                     void *fitSpec,
+                                     pmFit fit)
 {
     psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4,
             "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n);
+    //    PS_VECTOR_PRINT_F32(overscanVector);
 
     if (NULL == overscanVector) {
@@ -347 +243 @@
     // Allocate the new vector.
     psVector *newVec = psVectorAlloc(n, PS_TYPE_F32);
-    newVec->n = newVec->nalloc;
+
     //
     // If the new vector is the same size as the old, simply copy the data.
     //
     if (n == overscanVector->n) {
-        return(psVectorCopy(newVec, overscanVector, PS_TYPE_F32));
-    }
+        for (psS32 i = 0 ; i < n ; i++) {
+            newVec->data.F32[i] = overscanVector->data.F32[i];
+        }
+        return(newVec);
+    }
+    psPolynomial1D *myPoly;
+    psSpline1D *mySpline;
     psF32 x;
-
+    psS32 i;
     if (fit == PM_FIT_POLYNOMIAL) {
         // Fit a polynomial to the old overscan vector.
-        psPolynomial1D *myPoly = (psPolynomial1D *) fitSpec;
+        myPoly = (psPolynomial1D *) fitSpec;
         PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
-        PS_ASSERT_POLY1D(myPoly, NULL);
         myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
         if (myPoly == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the psVector.\n");
+            psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n");
             return(NULL);
         }
@@ -370 +270 @@
         // of the old vector, use the fitted polynomial to determine the
         // interpolated value at that point, and set the new vector.
-        for (psS32 i=0;i<n;i++) {
+        for (i=0;i<n;i++) {
             x = ((psF32) i) * ((psF32) overscanVector->n) / ((psF32) n);
             newVec->data.F32[i] = psPolynomial1DEval(myPoly, x);
         }
     } else if (fit == PM_FIT_SPLINE) {
+        psS32 mustFreeSpline = 0;
+        // Fit a spline to the old overscan vector.
+        mySpline = (psSpline1D *) fitSpec;
+        // XXX: Does it make any sense to have a psSpline argument?
+        if (mySpline == NULL) {
+            mustFreeSpline = 1;
+        }
+
         //
         // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL,
@@ -380 +288 @@
         // properly when doing the spline eval.
         //
-        psSpline1D *mySpline = psVectorFitSpline1D(NULL, overscanVector);
+        //        mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL);
+        mySpline = psVectorFitSpline1D(NULL, overscanVector);
         if (mySpline == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to the psVector.\n");
+            psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n");
             return(NULL);
         }
+        //        PS_PRINT_SPLINE(mySpline);
 
         // For each element of the new vector, convert the x-ordinate to that
-        // of the old vector, use the fitted spline to determine the
+        // of the old vector, use the fitted polynomial to determine the
         // interpolated value at that point, and set the new vector.
-        for (psS32 i=0;i<n;i++) {
+        for (i=0;i<n;i++) {
             // Scale to [0 : overscanVector->n - 1]
             x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n);
             newVec->data.F32[i] = psSpline1DEval(mySpline, x);
         }
-
-        psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
-        if (ptrSpline != NULL) {
-            // Copy the resulting spline fit into ptrSpline.
-            PS_ASSERT_SPLINE(ptrSpline, NULL);
-            SplineCopy(ptrSpline, mySpline);
-        }
-        psFree(mySpline);
+        if (mustFreeSpline ==1) {
+            psFree(mySpline);
+        }
+        //        PS_VECTOR_PRINT_F32(newVec);
+
+
     } else {
         psError(PS_ERR_UNKNOWN, true, "unknown fit type.  Returning NULL.\n");
@@ -413 +321 @@
 }
 
+#endif
+
+// Produce an overscan vector from an array of pixels
+static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options
+                                const psArray *pixels, // Array of vectors containing the pixel values
+                                psStats *myStats // Statistic to use in reducing the overscan
+                               )
+{
+    // Reduce the overscans
+    psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row
+    psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate
+    psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting
+    for (int i = 0; i < pixels->n; i++) {
+        psVector *values = pixels->data[i]; // Vector with overscan values
+        if (values->n > 0) {
+            mask->data.U8[i] = 0;
+            ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1]
+            psVectorStats(myStats, values, NULL, NULL, 0);
+            double reducedVal = NAN; // Result of statistics
+            if (! p_psGetStatValue(myStats, &reducedVal)) {
+                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result "
+                        "of statistics on row %d.\n", i);
+                return NULL;
+            }
+            reduced->data.F32[i] = reducedVal;
+        } else if (overscanOpts->fitType == PM_FIT_NONE) {
+            psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the "
+                    "image, and no fit is requested.\n");
+            return NULL;
+        } else {
+            // We'll fit this one out
+            mask->data.U8[i] = 1;
+        }
+    }
+
+    // Fit the overscan, if required
+    switch (overscanOpts->fitType) {
+    case PM_FIT_NONE:
+        // No fitting --- that's easy.
+        break;
+    case PM_FIT_POLY_ORD:
+        if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
+                                   overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) {
+            psFree(overscanOpts->poly);
+            overscanOpts->poly = NULL;
+        }
+        if (! overscanOpts->poly) {
+            overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order);
+        }
+        psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
+        psFree(reduced);
+        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
+        break;
+    case PM_FIT_POLY_CHEBY:
+        if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
+                                   overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) {
+            psFree(overscanOpts->poly);
+            overscanOpts->poly = NULL;
+        }
+        if (! overscanOpts->poly) {
+            overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order);
+        }
+        psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
+        psFree(reduced);
+        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
+        break;
+    case PM_FIT_SPLINE:
+        // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an
+        // input spline
+        overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate);
+        psFree(reduced);
+        reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate);
+        break;
+    default:
+        psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType);
+        return NULL;
+        break;
+    }
+
+    psFree(ordinate);
+    psFree(mask);
+
+    return reduced;
+}
+
+
+
 /******************************************************************************
+XXX: The SDRS does not specify type support.  F32 is implemented here.
  *****************************************************************************/
-static psS32 GetOverscanSize(
-    psImage *inImg,
-    pmOverscanAxis overScanAxis)
-{
-    if (overScanAxis == PM_OVERSCAN_ROWS) {
-        return(inImg->numCols);
-    } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-        return(inImg->numRows);
-    } else if (overScanAxis == PM_OVERSCAN_ALL) {
-        return(1);
-    }
-    return(0);
-}
-
-/******************************************************************************
-GetOverscanAxis(in) this private routine determines the appropiate overscan
-axis from the parent cell metadata.
- 
-XXX: Verify the READDIR corresponds with my overscan axis.
- *****************************************************************************/
-static pmOverscanAxis GetOverscanAxis(pmReadout *in)
-{
-    psBool rc;
-    if ((in->parent != NULL) && (in->parent->concepts)) {
-        psS32 dir = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.READDIR");
-        if (rc == true) {
-            if (dir == 1) {
-                return(PM_OVERSCAN_ROWS);
-            } else if (dir == 2) {
-                return(PM_OVERSCAN_COLUMNS);
-            } else if (dir == 3) {
-                return(PM_OVERSCAN_ALL);
-            }
-        }
-    }
-
-    psLogMsg(__func__, PS_LOG_WARN,
-             "WARNING: pmSubtractBias.(): could not determine CELL.READDIR from in->parent metadata.  Setting overscan axis to PM_OVERSCAN_NONE.\n");
-    return(PM_OVERSCAN_NONE);
-}
-
-/******************************************************************************
-my_psListLength(list): determine the length of a psList.
- 
-XXX: Put this elsewhere.
- 
-XXX: psList.h now has a version of this function.  Use that instead.
- *****************************************************************************/
-
-static psS32 my_psListLength(
-    psList *list)
-{
-    psS32 length = 0;
-    psListElem *tmpElem = (psListElem *) list->head;
-    while (NULL != tmpElem) {
-        tmpElem = tmpElem->next;
-        length++;
-    }
-    return(length);
-}
-
-/******************************************************************************
-Note: this isn't needed anymore as of psModule SDRS 12-09.
- *****************************************************************************/
-static psBool OverscanReducePixel(
-    psImage *in,
-    psList *bias,
-    psStats *myStats)
-{
-    PS_ASSERT_PTR_NON_NULL(in, NULL);
-    PS_ASSERT_PTR_NON_NULL(bias, NULL);
-    PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
-    PS_ASSERT_PTR_NON_NULL(myStats, NULL);
-
-    // Allocate a psVector with one element per overscan image.
-    psS32 numOverscanImages = my_psListLength(bias);
-    psVector *statsAll = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
-    statsAll->n = statsAll->nalloc;
-    psListElem *tmpOverscan = (psListElem *) bias->head;
-    psS32 i = 0;
-    psF64 statValue;
-    //
-    // We loop through each overscan image, calculating the specified
-    // statistic on that image.
-    //
-    while (NULL != tmpOverscan) {
-        psImage *myOverscanImage = (psImage *) tmpOverscan->data;
-
-        PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
-        myStats = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff);
-        if (myStats == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
-            psFree(statsAll);
-            return(false);
-        }
-        if (false == p_psGetStatValue(myStats, &statValue)) {
-            psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-            psFree(statsAll);
-            return(false);
-        }
-        statsAll->data.F32[i] = statValue;
-        i++;
-        tmpOverscan = tmpOverscan->next;
-    }
-
-    //
-    // We reduce the individual stats for each overscan image to
-    // a single psF32.
-    //
-    myStats = psVectorStats(myStats, statsAll, NULL, NULL, 0);
-    if (myStats == NULL) {
-        psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
-        psFree(statsAll);
-        return(false);
-    }
-    if (false == p_psGetStatValue(myStats, &statValue)) {
-        psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-        psFree(statsAll);
-        return(false);
-    }
-
-    //
-    // Subtract the result and return.
-    //
-    ImageSubtractScalar(in, statValue);
-    psFree(statsAll);
-    return(in);
-}
-
-/******************************************************************************
-ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
-a single psImage to a column by combining all pixels from each row into a
-single pixel via requested statistic in myStats.
- *****************************************************************************/
-static psVector *ReduceOverscanImageToCol(
-    psImage *overscanImage,
-    psStats *myStats)
-{
-    psF64 statValue;
-    psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
-    psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
-    tmpRow->n = tmpRow->nalloc;
-    tmpCol->n = tmpCol->nalloc;
-
-    //
-    // For each row, we store all pixels in that row into a temporary psVector,
-    // then we run psVectorStats() on that vector.
-    //
-    for (psS32 i=0;i<overscanImage->numRows;i++) {
-        for (psS32 j=0;j<overscanImage->numCols;j++) {
-            tmpRow->data.F32[j] = overscanImage->data.F32[i][j];
-        }
-
-        psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
-        if (rc == NULL) {
-            psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
-            return(NULL);
-        }
-
-        if (false ==  p_psGetStatValue(rc, &statValue)) {
-            psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
-            return(NULL);
-        }
-
-        tmpCol->data.F32[i] = (psF32) statValue;
-    }
-    psFree(tmpRow);
-
-    return(tmpCol);
-}
-
-/******************************************************************************
-ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
-a single psImage to a row by combining all pixels from each column into a
-single pixel via requested statistic in myStats.
- *****************************************************************************/
-static psVector *ReduceOverscanImageToRow(
-    psImage *overscanImage,
-    psStats *myStats)
-{
-    psF64 statValue;
-    psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
-    psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
-    tmpRow->n = tmpRow->nalloc;
-    tmpCol->n = tmpCol->nalloc;
-
-    //
-    // For each column, we store all pixels in that column into a temporary psVector,
-    // then we run psVectorStats() on that vector.
-    //
-    for (psS32 i=0;i<overscanImage->numCols;i++) {
-        for (psS32 j=0;j<overscanImage->numRows;j++) {
-            tmpCol->data.F32[j] = overscanImage->data.F32[j][i];
-        }
-
-        psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
-        if (rc == NULL) {
-            psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
-            return(NULL);
-        }
-
-        if (false ==  p_psGetStatValue(rc, &statValue)) {
-            psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
-            return(NULL);
-        }
-
-        tmpRow->data.F32[i] = (psF32) statValue;
-    }
-    psFree(tmpCol);
-
-    return(tmpRow);
-}
-
-/******************************************************************************
-OverscanReduce(vecSize, bias, myStats): This private routine takes a psList of
-overscan images (in bias) and reduces them to a single psVector via the
-specified psStats struct.  The vector is then scaled to the length or the
-row/column in inImg.
- *****************************************************************************/
-static psVector* OverscanReduce(
-    psImage *inImg,
-    pmOverscanAxis overScanAxis,
-    psList *bias,
-    void *fitSpec,
-    pmFit fit,
-    psStats *myStats)
-{
-    if ((overScanAxis != PM_OVERSCAN_ROWS) && (overScanAxis != PM_OVERSCAN_COLUMNS)) {
-        psError(PS_ERR_UNKNOWN, true, "overScanAxis must be PM_OVERSCAN_ROWS or PM_OVERSCAN_COLUMNS\n");
-        return(NULL);
-    }
-    PS_ASSERT_PTR_NON_NULL(inImg, NULL);
-    PS_ASSERT_PTR_NON_NULL(bias, NULL);
-    PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
-    PS_ASSERT_PTR_NON_NULL(myStats, NULL);
-    //
-    // Allocate a psVector for the output of this routine.
-    //
-    psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
-    psVector *overscanVector = psVectorAlloc(vecSize, PS_TYPE_F32);
-    overscanVector->n = overscanVector->nalloc;
-
-    //
-    // Allocate an array of psVectors with one psVector per element of the
-    // final oversan column vector.  These psVectors will be used with
-    // psStats to reduce the multiple elements from each overscan column
-    // vector to a single final column vector.
-    //
-    psS32 numOverscanImages = my_psListLength(bias);
-    psVector **overscanVectors = (psVector **) psAlloc(numOverscanImages * sizeof(psVector *));
-    //    (*overscanVectors)->n = (*overscanVectors)->nalloc;
-    for (psS32 i = 0 ; i < numOverscanImages ; i++) {
-        overscanVectors[i] = NULL;
-    }
-
-    //
-    // We iterate through the list of overscan images.  For each image,
-    // we reduce it to a single column or row.  Save the overscan vector
-    // in overscanVectors[].
-    //
-    psListElem *tmpOverscan = (psListElem *) bias->head;
-    psS32 overscanID = 0;
-    while (tmpOverscan != NULL) {
-        psImage *tmpOverscanImage = (psImage *) tmpOverscan->data;
-        if (overScanAxis == PM_OVERSCAN_ROWS) {
-            overscanVectors[overscanID] = ReduceOverscanImageToRow(tmpOverscanImage, myStats);
-        } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-            overscanVectors[overscanID] = ReduceOverscanImageToCol(tmpOverscanImage, myStats);
-        }
-
-        tmpOverscan = tmpOverscan->next;
-        overscanID++;
-    }
-
-    //
-    // For each overscan vector, if necessary, we scale that column or
-    // row to vecSize.  Note: we should have already ensured that the
-    // fit is poly or spline.
-    //
-    for (psS32 i = 0 ; i < numOverscanImages ; i++) {
-        psVector *tmpOverscanVector = overscanVectors[i];
-
-        if (tmpOverscanVector->n != vecSize) {
-            overscanVectors[i] = ScaleOverscanVector(tmpOverscanVector, vecSize, fitSpec, fit);
-            if (overscanVectors[i] == NULL) {
-                psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.\n");
-                for (psS32 i = 0 ; i < numOverscanImages ; i++) {
-                    psFree(overscanVectors[i]);
-                }
-                psFree(overscanVectors);
-                psFree(tmpOverscanVector);
-                return(NULL);
-            }
-            psFree(tmpOverscanVector);
-        }
-    }
-
-    //
-    // We collect all elements in the overscan vectors for the various
-    // overscan images into a single psVector (tmpVec).  Then we call
-    // psStats on that vector to determine the final values for the
-    // overscan vector.
-    //
-    psVector *tmpVec = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
-    tmpVec->n = tmpVec->nalloc;
-    psF64 statValue;
-    for (psS32 i = 0 ; i < vecSize ; i++) {
-        // Collect the i-th elements from each overscan vector into a single vector.
-        for (psS32 j = 0 ; j < numOverscanImages ; j++) {
-            tmpVec->data.F32[j] = overscanVectors[j]->data.F32[i];
-        }
-
-        if (NULL == psVectorStats(myStats, tmpVec, NULL, NULL, 0)) {
-            psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
-            for (psS32 i = 0 ; i < numOverscanImages ; i++) {
-                psFree(overscanVectors[i]);
-            }
-            psFree(overscanVectors);
-            psFree(tmpVec);
-            return(NULL);
-        }
-        if (false == p_psGetStatValue(myStats, &statValue)) {
-            psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-            for (psS32 i = 0 ; i < numOverscanImages ; i++) {
-                psFree(overscanVectors[i]);
-            }
-            psFree(overscanVectors);
-            psFree(tmpVec);
-            return(NULL);
-        }
-
-        overscanVector->data.F32[i] = (psF32) statValue;
-    }
-
-    //
-    // We're done.  Free the intermediate overscan vectors.
-    //
-    psFree(tmpVec);
-    for (psS32 i = 0 ; i < numOverscanImages ; i++) {
-        psFree(overscanVectors[i]);
-    }
-    psFree(overscanVectors);
-
-    //
-    // Return the computed overscanVector
-    //
-    return(overscanVector);
-}
-
-/******************************************************************************
-RebinOverscanVector(overscanVector, nBinOrig, myStats): this private routine
-takes groups of nBinOrig elements in the input vector, combines them into a
-single pixel via myStats and psVectorStats(), and then outputs a vector of
-those pixels.
- *****************************************************************************/
-static psS32 RebinOverscanVector(
-    psVector *overscanVector,
-    psS32 nBinOrig,
-    psStats *myStats)
-{
-    psF64 statValue;
-    psS32 nBin;
-    if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
-        psS32 numBins = 1+((overscanVector->n)/nBinOrig);
-        psVector *myBin = psVectorAlloc(numBins, PS_TYPE_F32);
-        psVector *binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);
-        myBin->n = myBin->nalloc;
-        binVec->n = binVec->nalloc;
-
-        for (psS32 i=0;i<numBins;i++) {
-            for(psS32 j=0;j<nBinOrig;j++) {
-                if (overscanVector->n > ((i*nBinOrig)+j)) {
-                    binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j];
-                } else {
-                    // XXX: we get here if nBinOrig does not evenly divide
-                    // the overscanVector vector.  This is the last bin.  Should
-                    // we change the binVec->n to acknowledge that?
-                    binVec->n = j;
-                }
-            }
-            psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0);
-            if (rc == NULL) {
-                psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
-                return(-1);
-            }
-            if (false ==  p_psGetStatValue(rc, &statValue)) {
-                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                return(-1);
-            }
-            myBin->data.F32[i] = statValue;
-        }
-
-        // Change the effective size of overscanVector.
-        overscanVector->n = numBins;
-        for (psS32 i=0;i<numBins;i++) {
-            overscanVector->data.F32[i] = myBin->data.F32[i];
-        }
-        psFree(binVec);
-        psFree(myBin);
-        nBin = nBinOrig;
-    } else {
-        nBin = 1;
-    }
-
-    return(nBin);
-}
-
-/******************************************************************************
-FitOverscanVectorAndUnbin(inImg, overscanVector, overScanAxis, fitSpec, fit,
-nBin):  this private routine fits a psPolynomial or psSpline to the overscan
-vector.  It then creates a new vector, with a size determined by the input
-image, evaluates the psPolynomial or psSpline at each element in that vector,
-then returns that vector.
- *****************************************************************************/
-static psVector *FitOverscanVectorAndUnbin(
-    psImage *inImg,
-    psVector *overscanVector,
-    pmOverscanAxis overScanAxis,
-    void *fitSpec,
-    pmFit fit,
-    psS32 nBin)
-{
-    psPolynomial1D* myPoly = NULL;
-    psSpline1D *mySpline = NULL;
-    //
-    // Fit a polynomial or spline to the overscan vector.
-    //
-    if (fit == PM_FIT_POLYNOMIAL) {
-        myPoly = (psPolynomial1D *) fitSpec;
-        PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
-        PS_ASSERT_POLY1D(myPoly, NULL);
-        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
-        if (myPoly == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector.  Returning NULL.\n");
-            return(NULL);
-        }
-    } else if (fit == PM_FIT_SPLINE) {
-        mySpline = psVectorFitSpline1D(NULL, overscanVector);
-        if (mySpline == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning NULL.\n");
-            return(NULL);
-        }
-        if (fitSpec != NULL) {
-            // Copy the resulting spline fit into fitSpec.
-            psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
-            PS_ASSERT_SPLINE(ptrSpline, NULL);
-            SplineCopy(ptrSpline, mySpline);
-        }
-    }
-
-    //
-    // Evaluate the poly/spline at each pixel in the overscan row/column.
-    //
-    psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
-    psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
-    newVec->n = newVec->nalloc;
-    if ((nBin > 1) && (nBin < overscanVector->n)) {
-        for (psS32 i = 0 ; i < vecSize ; i++) {
-            if (fit == PM_FIT_POLYNOMIAL) {
-                newVec->data.F32[i] = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
-            } else if (fit == PM_FIT_SPLINE) {
-                newVec->data.F32[i] = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
-            }
-        }
-    } else {
-        for (psS32 i = 0 ; i < vecSize ; i++) {
-            if (fit == PM_FIT_POLYNOMIAL) {
-                newVec->data.F32[i] = psPolynomial1DEval(myPoly, (psF32) i);
-            } else if (fit == PM_FIT_SPLINE) {
-                newVec->data.F32[i] = psSpline1DEval(mySpline, (psF32) i);
-            }
-        }
-    }
-
-    psFree(mySpline);
-    psFree(overscanVector);
-    return(newVec);
-}
-
-
-
-/******************************************************************************
-UnbinOverscanVector(inImg, overscanVector, overScanAxis, nBin):  this private
-routine takes a psVector overscanVector that was previously binned by a factor
-of nBin, and then expands it to its original size, duplicated elements nBin
-times for each element in the input vector overscanVector.
- *****************************************************************************/
-static psVector *UnbinOverscanVector(
-    psImage *inImg,
-    psVector *overscanVector,
-    pmOverscanAxis overScanAxis,
-    psS32 nBin)
-{
-    psS32 vecSize = 0;
-
-    if (overScanAxis == PM_OVERSCAN_ROWS) {
-        vecSize = inImg->numCols;
-    } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-        vecSize = inImg->numRows;
-    }
-
-    psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
-    newVec->n = newVec->nalloc;
-    for (psS32 i = 0 ; i < vecSize ; i++) {
-        newVec->data.F32[i] = overscanVector->data.F32[i/nBin];
-    }
-
-    psFree(overscanVector);
-    return(newVec);
-}
-
-
-/******************************************************************************
-SubtractVectorFromImage(inImg, overscanVector, overScanAxis):  this private
-routine subtracts the overscanVector column-wise or row-wise from inImg.
- *****************************************************************************/
-static psImage *SubtractVectorFromImage(
-    psImage *inImg,
-    psVector *overscanVector,
-    pmOverscanAxis overScanAxis)
-{
-    //
-    // Subtract overscan vector row-wise from the image.
-    //
-    if (overScanAxis == PM_OVERSCAN_ROWS) {
-        for (psS32 i=0;i<inImg->numCols;i++) {
-            for (psS32 j=0;j<inImg->numRows;j++) {
-                inImg->data.F32[j][i]-= overscanVector->data.F32[i];
-            }
-        }
-    }
-
-    //
-    // Subtract overscan vector column-wise from the image.
-    //
-    if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-        for (psS32 i=0;i<inImg->numRows;i++) {
-            for (psS32 j=0;j<inImg->numCols;j++) {
-                inImg->data.F32[i][j]-= overscanVector->data.F32[i];
-            }
-        }
-    }
-
-    return(inImg);
-}
-
-
-
-typedef enum {
-    PM_ERROR_NO_SUBTRACTION,
-    PM_WARNING_NO_SUBTRACTION,
-    PM_ERROR_NO_BIAS_SUBTRACT,
-    PM_WARNING_NO_BIAS_SUBTRACT,
-    PM_ERROR_NO_DARK_SUBTRACT,
-    PM_WARNING_NO_DARK_SUBTRACT,
-    PM_OKAY
-} pmSubtractBiasAssertStatus;
-/******************************************************************************
-AssertCodeOverscan(....) this private routine verifies that the various input
-parameters to pmSubtractBias() are correct for overscan subtraction.
- *****************************************************************************/
-pmSubtractBiasAssertStatus AssertCodeOverscan(
-    pmReadout *in,
-    void *fitSpec,
-    pmFit fit,
-    bool overscan,
-    psStats *stat,
-    int nBinOrig,
-    const pmReadout *bias,
-    const pmReadout *dark)
-{
-
-    PS_ASSERT_READOUT_NON_NULL(in, PM_ERROR_NO_SUBTRACTION);
-    PS_ASSERT_READOUT_NON_EMPTY(in, PM_ERROR_NO_SUBTRACTION);
-    PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
-    PS_WARN_PTR_NON_NULL(in->parent);
-    if (in->parent != NULL) {
-        PS_WARN_PTR_NON_NULL(in->parent->concepts);
-    }
-
-    if (overscan == true) {
-        pmOverscanAxis overScanAxis = GetOverscanAxis(in);
-        PS_ASSERT_PTR_NON_NULL(stat, PM_ERROR_NO_SUBTRACTION);
-        PS_ASSERT_PTR_NON_NULL(in->bias, PM_ERROR_NO_SUBTRACTION);
-        PS_ASSERT_PTR_NON_NULL(in->bias->head, PM_ERROR_NO_SUBTRACTION);
-        //
-        // Check the type, size of each bias image.
-        //
-        psListElem *tmpOverscan = (psListElem *) in->bias->head;
-        psS32 numOverscans = 0;
-        while (NULL != tmpOverscan) {
-            numOverscans++;
-            psImage *myOverscanImage = (psImage *) tmpOverscan->data;
-            PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
-            // XXX: Get this right with the rows and columns.
-            if (overScanAxis == PM_OVERSCAN_ROWS) {
-                if (myOverscanImage->numRows != in->image->numRows) {
-                    psLogMsg(__func__, PS_LOG_WARN,
-                             "WARNING: pmSubtractBias.(): overscan image (# %d) has %d rows, input image has %d rows\n",
-                             numOverscans, myOverscanImage->numCols, in->image->numRows);
-                    if (fit == PM_FIT_NONE) {
-                        psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors.  Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n");
-                        return(PM_ERROR_NO_SUBTRACTION);
-                    }
-                }
-            } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-                if (myOverscanImage->numCols != in->image->numCols) {
-                    psLogMsg(__func__, PS_LOG_WARN,
-                             "WARNING: pmSubtractBias.(): overscan image (# %d) has %d columns, input image has %d columns\n",
-                             numOverscans, myOverscanImage->numCols, in->image->numCols);
-                    if (fit == PM_FIT_NONE) {
-                        psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors.  Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n");
-                        return(PM_ERROR_NO_SUBTRACTION);
-                    }
-                }
-            } else if (overScanAxis != PM_OVERSCAN_ALL) {
-                psError(PS_ERR_UNKNOWN, true, "Must specify and overscan axis.\n");
-                return(PM_ERROR_NO_SUBTRACTION);
-            }
-            tmpOverscan = tmpOverscan->next;
-        }
-    } else {
-        if (fit != PM_FIT_NONE) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmSubtractBias.(): overscan is FALSE and fit is not PM_FIT_NONE.\n");
-            return(PM_WARNING_NO_SUBTRACTION);
-        }
-    }
-
-    // XXX: I do not like the following spec since it's useless to specify
-    // a psSpline as the fitSpec.
-    if (0) {
-        if ((fitSpec == NULL) &&
-                ((fit != PM_FIT_NONE) || (overscan == true))) {
-            psError(PS_ERR_UNKNOWN, true, "fitSpec is NULL and fit is not PM_FIT_NONE or overscan is TRUE.\n");
-            return(PM_ERROR_NO_SUBTRACTION);
-        }
-    }
-
-    return(PM_OKAY);
-}
-
-/******************************************************************************
-AssertCodeBias(....) this private routine verifies that the various input
-parameters to pmSubtractBias() are correct for bias subtraction.
- *****************************************************************************/
-static pmSubtractBiasAssertStatus AssertCodeBias(
-    pmReadout *in,
-    void *fitSpec,
-    pmFit fit,
-    bool overscan,
-    psStats *stat,
-    int nBinOrig,
-    const pmReadout *bias,
-    const pmReadout *dark)
-{
-    if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) {
-        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows.  Returning in image\n");
-        return(PM_ERROR_NO_BIAS_SUBTRACT);
-    }
-    if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) {
-        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
-        return(PM_ERROR_NO_BIAS_SUBTRACT);
-    }
-
-    if (bias != NULL) {
-        PS_ASSERT_READOUT_NON_EMPTY(bias, PM_ERROR_NO_BIAS_SUBTRACT);
-        PS_ASSERT_READOUT_TYPE(bias, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
-    }
-    return(PM_OKAY);
-}
-
-/******************************************************************************
-AssertCodeDark(....) this private routine verifies that the various input
-parameters to pmSubtractBias() are correct for dark subtraction.
- *****************************************************************************/
-pmSubtractBiasAssertStatus AssertCodeDark(
-    pmReadout *in,
-    void *fitSpec,
-    pmFit fit,
-    bool overscan,
-    psStats *stat,
-    int nBinOrig,
-    const pmReadout *bias,
-    const pmReadout *dark)
-{
-    if ((in->image->numRows + in->row0 - dark->row0) > dark->image->numRows) {
-        psError(PS_ERR_UNKNOWN, true, "dark image does not have enough rows.  Returning in image\n");
-        return(PM_ERROR_NO_DARK_SUBTRACT);
-    }
-    if ((in->image->numCols + in->col0 - dark->col0) > dark->image->numCols) {
-        psError(PS_ERR_UNKNOWN, true, "dark image does not have enough columns.  Returning in image\n");
-        return(PM_ERROR_NO_DARK_SUBTRACT);
-    }
-
-    if (dark != NULL) {
-        PS_ASSERT_READOUT_NON_EMPTY(dark, PM_ERROR_NO_DARK_SUBTRACT);
-        PS_ASSERT_READOUT_TYPE(dark, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
-    }
-    return(PM_OKAY);
-}
-
-/******************************************************************************
-p_psDetermineTrimmedImage(): global routine: determines the region of the
-input pmReadout which will be operated on by the various detrend modules.  It
-does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout.
- 
-Use it this way:
-    PS_WARN_PTR_NON_NULL(in->parent);
-    if (in->parent != NULL) {
-        PS_WARN_PTR_NON_NULL(in->parent->concepts);
-    }
-    //
-    // Determine trimmed image from metadata.
-    //
-    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
- 
-XXX: Create a pmUtils.c file and put this routine there.
- *****************************************************************************/
-psImage *p_psDetermineTrimmedImage(pmReadout *in)
-{
-    if ((in->parent == NULL) || (in->parent->concepts == NULL)) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata (NULL).\n");
-        return(in->image);
-    }
-
-    psBool rc = false;
-    psImage *trimmedImg = NULL;
-    psRegion *trimRegion = psMetadataLookupPtr(&rc, in->parent->concepts,
-                           "CELL.TRIMSEC");
-    if (rc == false) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata.\n");
-        trimmedImg = in->image;
-    } else {
-        trimmedImg = psImageSubset(in->image, *trimRegion);
-    }
-
-    return(trimmedImg);
-}
-
-
-/******************************************************************************
-pmSubtractBias(....): see SDRS for complete specification.
- 
-XXX: Code and assert type support: U16, S32, F32.
-XXX: Add trace messages.
- *****************************************************************************/
-pmReadout *pmSubtractBias(
-    pmReadout *in,
-    void *fitSpec,
-    pmFit fit,
-    bool overscan,
-    psStats *stat,
-    int nBin,
-    const pmReadout *bias,
-    const pmReadout *dark)
+pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts,
+                          const pmReadout *bias, const pmReadout *dark)
 {
     psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
             "---- pmSubtractBias() begin ----\n");
-    //
-    // Check input parameters, generate warnings and errors.
-    //
-    if (PM_OKAY != AssertCodeOverscan(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
-        return(in);
-    }
-    //
-    // Determine trimmed image from metadata.
-    //
-    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
-
-    //
-    // Subtract overscan frames if necessary.
-    //
-    if (overscan == true) {
-        pmOverscanAxis overScanAxis = GetOverscanAxis(in);
-        //
-        //  Create a psStats data structure and determine the highest
-        //  priority stats option.
-        //
-        psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
-        if (stat != NULL) {
-            myStats->options = GenNewStatOptions(stat);
-        }
-
-        //
-        // Reduce overscan images to a single pixel, then subtract.
-        // This code is no longer required as of SDRS 12-09.
-        //
-        if (overScanAxis == PM_OVERSCAN_ALL) {
-            if (false == OverscanReducePixel(trimmedImg, in->bias, myStats)) {
+    PS_ASSERT_READOUT_NON_NULL(in, NULL);
+    PS_ASSERT_READOUT_NON_EMPTY(in, NULL);
+    PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL);
+
+    psImage *image = in->image;         // The input image
+
+    // Overscan processing
+    if (overscanOpts) {
+        // Check for an unallowable pmFit.
+        if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD &&
+                overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) {
+            psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d).  Returning original image.\n", overscanOpts->fitType);
+            return(in);
+        }
+
+        psList *overscans = in->bias; // List of the overscan images
+
+        psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // A new psStats, to avoid clobbering original
+        myStats->options = GenNewStatOptions(overscanOpts->stat);
+
+        // Reduce all overscan pixels to a single value
+        if (overscanOpts->single) {
+            psVector *pixels = psVectorAlloc(0, PS_TYPE_F32);
+            pixels->n = 0;
+            psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+            psImage *overscan = NULL;   // Overscan image from iterator
+            while ((overscan = psListGetAndIncrement(iter))) {
+                int index = pixels->n;  // Index
+                pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols);
+                // XXX Reimplement with memcpy
+                for (int i = 0; i < overscan->numRows; i++) {
+                    for (int j = 0; j < overscan->numCols; j++) {
+                        pixels->data.F32[index++] = overscan->data.F32[i][j];
+                    }
+                }
+
+            }
+            psFree(iter);
+
+            (void)psVectorStats(myStats, pixels, NULL, NULL, 0);
+            double reduced = NAN;     // Result of statistics
+            if (! p_psGetStatValue(myStats, &reduced)) {
+                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning input image.\n");
                 return(in);
             }
-            psFree(myStats);
+            (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32));
         } else {
-            //
-            // Reduce the overscan images to a single overscan vector.
-            //
-            psVector *overscanVector = OverscanReduce(in->image, overScanAxis,
-                                       in->bias, fitSpec,
-                                       fit, myStats);
-            if (overscanVector == NULL) {
-                psError(PS_ERR_UNKNOWN, false, "Could not reduce overscan images to a single overscan vector.  Returning in image\n");
-                psFree(myStats);
-                return(in);
+
+            // We do the regular overscan subtraction
+
+            bool readRows = psMetadataLookupBool(NULL, in->parent->concepts, "CELL.READDIR");// Read direction
+
+            if (readRows) {
+                // The read direction is rows
+                psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels
+                for (int i = 0; i < pixels->n; i++) {
+                    psVector *values = psVectorAlloc(0, PS_TYPE_F32);
+                    values->n = 0;
+                    pixels->data[i] = values;
+                }
+
+                // Pull the pixels out into the vectors
+                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                psImage *overscan = NULL; // Overscan image from iterator
+                while ((overscan = psListGetAndIncrement(iter))) {
+                    int diff = image->row0 - overscan->row0; // Offset between the two regions
+                    for (int i = MAX(0,diff); i < MIN(image->numRows, overscan->numRows + diff); i++) {
+                        // i is row on overscan
+                        // XXX Reimplement with memcpy
+                        psVector *values = pixels->data[i];
+                        int index = values->n; // Index in the vector
+                        values = psVectorRealloc(values, values->n + overscan->numCols);
+                        for (int j = 0; j < overscan->numCols; j++) {
+                            values->data.F32[index++] = overscan->data.F32[i][j];
+                        }
+                        values->n += overscan->numCols;
+                        pixels->data[i] = values; // Update the pointer in case it's moved
+                    }
+                }
+                psFree(iter);
+
+                // Reduce the overscans
+                psVector *reduced = overscanVector(overscanOpts, pixels, myStats);
+                psFree(pixels);
+                if (! reduced) {
+                    return in;
+                }
+
+                // Subtract row by row
+                for (int i = 0; i < image->numRows; i++) {
+                    for (int j = 0; j < image->numCols; j++) {
+                        image->data.F32[i][j] -= reduced->data.F32[i];
+                    }
+                }
+                psFree(reduced);
+
+            } else {
+                // The read direction is columns
+                psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels
+                for (int i = 0; i < pixels->n; i++) {
+                    psVector *values = psVectorAlloc(0, PS_TYPE_F32);
+                    values->n = 0;
+                    pixels->data[i] = values;
+                }
+
+                // Pull the pixels out into the vectors
+                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                psImage *overscan = NULL; // Overscan image from iterator
+                while ((overscan = psListGetAndIncrement(iter))) {
+                    int diff = image->col0 - overscan->col0; // Offset between the two regions
+                    for (int i = MAX(0,diff); i < MIN(image->numCols, overscan->numCols + diff); i++) {
+                        // i is column on overscan
+                        // XXX Reimplement with memcpy
+                        psVector *values = pixels->data[i];
+                        int index = values->n; // Index in the vector
+                        values = psVectorRealloc(values, values->n + overscan->numRows);
+                        for (int j = 0; j < overscan->numRows; j++) {
+                            values->data.F32[index++] = overscan->data.F32[i][j];
+                        }
+                        values->n += overscan->numRows;
+                        pixels->data[i] = values; // Update the pointer in case it's moved
+                    }
+                }
+                psFree(iter);
+
+                // Reduce the overscans
+                psVector *reduced = overscanVector(overscanOpts, pixels, myStats);
+                psFree(pixels);
+                if (! reduced) {
+                    return in;
+                }
+
+                // Subtract column by column
+                for (int i = 0; i < image->numCols; i++) {
+                    for (int j = 0; j < image->numRows; j++) {
+                        image->data.F32[j][i] -= reduced->data.F32[i];
+                    }
+                }
+                psFree(reduced);
             }
-
-            //
-            // Rebin the overscan vector if necessary.
-            //
-            psS32 newBin = RebinOverscanVector(overscanVector, nBin, myStats);
-            if (newBin < 0) {
-                psError(PS_ERR_UNKNOWN, false, "Could rebin the overscan vector.  Returning in image\n");
-                psFree(myStats);
-                return(in);
-            }
-
-            //
-            // If necessary, fit a psPolynomial or psSpline to the overscan vector.
-            // Then, unbin the overscan vector to appropriate length for the in image.
-            //
-            if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
-                overscanVector = FitOverscanVectorAndUnbin(trimmedImg, overscanVector, overScanAxis, fitSpec, fit, newBin);
-                if (overscanVector == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "Could not fit the polynomial or spline to the overscan vector.  Returning in image\n");
-                    psFree(myStats);
-                    return(in);
-                }
-            } else {
-                overscanVector = UnbinOverscanVector(trimmedImg, overscanVector, overScanAxis, newBin);
-            }
-
-            //
-            // Subtract the overscan vector from the input image.
-            //
-            SubtractVectorFromImage(trimmedImg, overscanVector, overScanAxis);
-            psFree(myStats);
-            psFree(overscanVector);
-        }
-    }
-
-    //
-    // Perform bias subtraction if necessary.
-    //
-    if (bias != NULL) {
-        if (PM_OKAY == AssertCodeBias(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
-            SubtractFrame(in, bias);
-        }
-    }
-
-    //
-    // Perform dark subtraction if necessary.
-    //
-    if (dark != NULL) {
-        if (PM_OKAY == AssertCodeDark(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
-            psBool rc;
-            psF32 scale = 0.0;
-            if (in->parent != NULL) {
-                scale = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.DARKTIME");
-                if (rc == false) {
-                    psLogMsg(__func__, PS_LOG_WARN,
-                             "WARNING: pmSubtractBias.(): could not determine CELL.FARKTIME from in->parent metadata.\n");
-                }
-            }
-            SubtractDarkFrame(in, dark, scale);
-        }
-    }
-
-    //
-    // All done.
-    //
-    psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
-            "---- pmSubtractBias() exit ----\n");
-    return(in);
-}
-
-
+        }
+        psFree(myStats);
+    } // End of overscan subtraction
+
+    // Bias frame subtraction
+    if (bias) {
+        SubtractFrame(in, bias, 1.0);
+    }
+
+    if (dark) {
+        // Get the scaling
+        float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME");
+        float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME");
+        SubtractFrame(in, dark, inTime/darkTime);
+    }
+
+    return in;
+}
+
+

trunk/psModules/src/imsubtract/pmSubtractBias.h

@@ -11 +11 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/imsubtract/pmSubtractSky.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/objects/pmModelGroup.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/objects/pmModelGroup.h

@@ -9 +9 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/objects/pmObjects.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-04 01:01:33 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -19 +19 @@
 #include "pmObjects.h"
 #include "pmModelGroup.h"
-/******************************************************************************
-pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members.
-*****************************************************************************/
-pmPeak *pmPeakAlloc(psS32 x,
-                    psS32 y,
-                    psF32 counts,
-                    pmPeakType class)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak));
-    tmp->x = x;
-    tmp->y = y;
-    tmp->counts = counts;
-    tmp->class = class;
 
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(tmp);
-}
-
-/******************************************************************************
-pmMomentsAlloc(): Allocate the pmMoments structure and initialize the members
-to zero.
-*****************************************************************************/
-pmMoments *pmMomentsAlloc()
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments));
-    tmp->x = 0.0;
-    tmp->y = 0.0;
-    tmp->Sx = 0.0;
-    tmp->Sy = 0.0;
-    tmp->Sxy = 0.0;
-    tmp->Sum = 0.0;
-    tmp->Peak = 0.0;
-    tmp->Sky = 0.0;
-    tmp->nPixels = 0;
-
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(tmp);
-}
-
-static void modelFree(pmModel *tmp)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    psFree(tmp->params);
-    psFree(tmp->dparams);
-    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-}
-
-static void sourceFree(pmSource *tmp)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    psFree(tmp->peak);
-    psFree(tmp->pixels);
-    psFree(tmp->weight);
-    psFree(tmp->mask);
-    psFree(tmp->moments);
-    psFree(tmp->modelPSF);
-    psFree(tmp->modelFLT);
-    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-}
-
-/******************************************************************************
-getRowVectorFromImage(): a private function which simply returns a
-psVector containing the specified row of data from the psImage.
- 
-XXX: Is there a better way to do this?
-XXX EAM: does this really need to alloc a new vector???
-*****************************************************************************/
-static psVector *getRowVectorFromImage(psImage *image,
-                                       psU32 row)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
-
-    psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32);
-    tmpVector->n = tmpVector->nalloc;
-    for (psU32 col = 0; col < image->numCols ; col++) {
-        tmpVector->data.F32[col] = image->data.F32[row][col];
-    }
-    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-    return(tmpVector);
-}
-
-/******************************************************************************
-myListAddPeak(): A private function which allocates a psArray, if the list
-argument is NULL, otherwise it adds the peak to that list.
-XXX EAM : changed the output to psArray
-XXX EAM : Switched row, col args
-XXX EAM : NOTE: this was changed in the call, so the new code is consistent
-*****************************************************************************/
-static psArray *myListAddPeak(psArray *list,
-                              psS32 row,
-                              psS32 col,
-                              psF32 counts,
-                              pmPeakType type)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);
-
-    if (list == NULL) {
-        list = psArrayAlloc(100);
-        list->n = 0;
-    }
-    psArrayAdd(list, 100, tmpPeak);
-    psFree (tmpPeak);
-    // XXX EAM : is this free appropriate?  (does psArrayAdd increment memory counter?)
-
-    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-    return(list);
-}
-
-
-/******************************************************************************
-bool checkRadius2(): private function which simply determines if the (x, y)
-point is within the radius of the specified peak.
- 
-XXX: macro this for performance.
-XXX: this is rather inefficient - at least compute and compare against radius^2
-*****************************************************************************/
-static bool checkRadius2(psF32 xCenter,
-                         psF32 yCenter,
-                         psF32 radius,
-                         psF32 x,
-                         psF32 y)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    /// XXX EAM should compare with hypot (x,y) for speed
-    if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) {
-        return(true);
-    }
-
-    psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
-    return(false);
-}
-
-// XXX: Macro this.
-static bool isItInThisRegion(const psRegion valid,
-                             psS32 x,
-                             psS32 y)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    if ((x >= valid.x0) &&
-            (x <= valid.x1) &&
-            (y >= valid.y0) &&
-            (y <= valid.y1)) {
-        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
-        return(true);
-    }
-    psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
-    return(false);
-}
-
-/******************************************************************************
-findValue(source, level, row, col, dir): a private function which determines
-the column coordinate of the model function which has the value "level".  If
-dir equals 0, then you loop leftwards from the peak pixel, otherwise,
-rightwards.
- 
-XXX: reverse order of row,col args?
- 
-XXX: Input row/col are in image coords.
- 
-XXX: The result is returned in image coords.
-*****************************************************************************/
-static psF32 findValue(pmSource *source,
-                       psF32 level,
-                       psU32 row,
-                       psU32 col,
-                       psU32 dir)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    //
-    // Convert coords to subImage space.
-    //
-    psU32 subRow = row - source->pixels->row0;
-    psU32 subCol = col - source->pixels->col0;
-
-    // Ensure that the starting column is allowable.
-    if (!((0 <= subCol) && (subCol < source->pixels->numCols))) {
-        psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range");
-        psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
-        return(NAN);
-    }
-    if (!((0 <= subRow) && (subRow < source->pixels->numRows))) {
-        psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
-        psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range");
-        return(NAN);
-    }
-
-    // XXX EAM : i changed this to match pmModelEval above, but see
-    // XXX EAM   the note below in pmSourceContour
-    psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
-    if (oldValue == level) {
-        psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-        return(((psF32) (subCol + source->pixels->col0)));
-    }
-
-    //
-    // We define variables incr and lastColumn so that we can use the same loop
-    // whether we are stepping leftwards, or rightwards.
-    //
-    psS32 incr;
-    psS32 lastColumn;
-    if (dir == 0) {
-        incr = -1;
-        lastColumn = -1;
-    } else {
-        incr = 1;
-        lastColumn = source->pixels->numCols;
-    }
-    subCol+=incr;
-
-    while (subCol != lastColumn) {
-        psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
-        if (oldValue == level) {
-            psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-            return((psF32) (subCol + source->pixels->col0));
-        }
-
-        if ((newValue <= level) && (level <= oldValue)) {
-            // This is simple linear interpolation.
-            psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-            return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) );
-        }
-
-        if ((oldValue <= level) && (level <= newValue)) {
-            // This is simple linear interpolation.
-            psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-            return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) );
-        }
-
-        subCol+=incr;
-    }
-
-    psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
-    return(NAN);
-}
-
-/******************************************************************************
-pmModelAlloc(): Allocate the pmModel structure, along with its parameters,
-and initialize the type member.  Initialize the params to 0.0.
-XXX EAM: simplifying code with pmModelParameterCount
-*****************************************************************************/
-pmModel *pmModelAlloc(pmModelType type)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel));
-
-    tmp->type = type;
-    tmp->chisq = 0.0;
-    tmp->nIter = 0;
-    tmp->radius = 0;
-    tmp->status = PM_MODEL_UNTRIED;
-
-    psS32 Nparams = pmModelParameterCount(type);
-    if (Nparams == 0) {
-        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
-        return(NULL);
-    }
-
-    tmp->params  = psVectorAlloc(Nparams, PS_TYPE_F32);
-    tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32);
-    tmp->params->n = tmp->params->nalloc;
-    tmp->dparams->n = tmp->dparams->nalloc;
-
-    for (psS32 i = 0; i < tmp->params->n; i++) {
-        tmp->params->data.F32[i] = 0.0;
-        tmp->dparams->data.F32[i] = 0.0;
-    }
-
-    psMemSetDeallocator(tmp, (psFreeFunc) modelFree);
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(tmp);
-}
-
-/******************************************************************************
-XXX EAM : we can now free these pixels - memory ref is incremented now
-*****************************************************************************/
-
-/******************************************************************************
-pmSourceAlloc(): Allocate the pmSource structure and initialize its members
-to NULL.
-*****************************************************************************/
-pmSource *pmSourceAlloc()
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource));
-    tmp->peak = NULL;
-    tmp->pixels = NULL;
-    tmp->weight = NULL;
-    tmp->mask = NULL;
-    tmp->moments = NULL;
-    tmp->modelPSF = NULL;
-    tmp->modelFLT = NULL;
-    tmp->type = 0;
-    psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
-
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(tmp);
-}
-
-/******************************************************************************
-pmFindVectorPeaks(vector, threshold): Find all local peaks in the given vector
-above the given threshold.  Returns a vector of type PS_TYPE_U32 containing
-the location (x value) of all peaks.
- 
-XXX: What types should be supported?  Only F32 is implemented.
- 
-XXX: We currently step through the input vector twice; once to determine the
-size of the output vector, then to set the values of the output vector.
-Depending upon actual use, this may need to be optimized.
-*****************************************************************************/
-psVector *pmFindVectorPeaks(const psVector *vector,
-                            psF32 threshold)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(vector, NULL);
-    PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL);
-    PS_ASSERT_VECTOR_TYPE(vector, PS_TYPE_F32, NULL);
-    int count = 0;
-    int n = vector->n;
-
-    //
-    // Special case: the input vector has a single element.
-    //
-    if (n == 1) {
-        psVector *tmpVector = NULL;
-        ;
-        if (vector->data.F32[0] > threshold) {
-            tmpVector = psVectorAlloc(1, PS_TYPE_U32);
-            tmpVector->n = 1;
-            tmpVector->data.U32[0] = 0;
-        } else {
-            tmpVector = psVectorAlloc(0, PS_TYPE_U32);
-        }
-        psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-        return(tmpVector);
-    }
-
-    //
-    // Determine if first pixel is a peak
-    //
-    if ((vector->data.F32[0] > vector->data.F32[1]) &&
-            (vector->data.F32[0] > threshold)) {
-        count++;
-    }
-
-    //
-    // Determine if interior pixels are peaks
-    //
-    for (psU32 i = 1; i < n-1 ; i++) {
-        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
-                (vector->data.F32[i] > vector->data.F32[i+1]) &&
-                (vector->data.F32[i] > threshold)) {
-            count++;
-        }
-    }
-
-    //
-    // Determine if last pixel is a peak
-    //
-    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
-            (vector->data.F32[n-1] > threshold)) {
-        count++;
-    }
-
-    //
-    // We know how many peaks exist, so we now allocate a psVector to store
-    // those peaks.
-    //
-    psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32);
-    tmpVector->n = tmpVector->nalloc;
-    count = 0;
-
-    //
-    // Determine if first pixel is a peak
-    //
-    if ((vector->data.F32[0] > vector->data.F32[1]) &&
-            (vector->data.F32[0] > threshold)) {
-        tmpVector->data.U32[count++] = 0;
-    }
-
-    //
-    // Determine if interior pixels are peaks
-    //
-    for (psU32 i = 1; i < (n-1) ; i++) {
-        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
-                (vector->data.F32[i] > vector->data.F32[i+1]) &&
-                (vector->data.F32[i] > threshold)) {
-            tmpVector->data.U32[count++] = i;
-        }
-    }
-
-    //
-    // Determine if last pixel is a peak
-    //
-    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
-            (vector->data.F32[n-1] > threshold)) {
-        tmpVector->data.U32[count++] = n-1;
-    }
-
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(tmpVector);
-}
-
-
-/******************************************************************************
-pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
-above the given threshold.  Returns a psArray containing location (x/y value)
-of all peaks.
- 
-XXX: I'm not convinced the peak type definition in the SDRS is mutually
-exclusive.  Some peaks can have multiple types.  Edges for sure.  Also, a
-digonal line with the same value at each point will have a peak for every
-point on that line.
- 
-XXX: This does not work if image has either a single row, or a single column.
- 
-XXX: In the output psArray elements, should we use the image row/column offsets?
-     Currently, we do not.
- 
-XXX: Merge with CVS 1.20.  This had the proper code for images with a single
-row or column.
-*****************************************************************************/
-psArray *pmFindImagePeaks(const psImage *image,
-                          psF32 threshold)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
-    if ((image->numRows == 1) || (image->numCols == 1)) {
-        psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns.");
-        psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
-        return(NULL);
-    }
-    psVector *tmpRow = NULL;
-    psU32 col = 0;
-    psU32 row = 0;
-    psArray *list = NULL;
-
-    //
-    // Find peaks in row 0 only.
-    //
-    row = 0;
-    tmpRow = getRowVectorFromImage((psImage *) image, row);
-    psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);
-
-    for (psU32 i = 0 ; i < row1->n ; i++ ) {
-        col = row1->data.U32[i];
-        //
-        // Determine if pixel (0,0) is a peak.
-        //
-        if (col == 0) {
-            if ( (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
-
-                if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
-                }
-            }
-        } else if (col < (image->numCols - 1)) {
-            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
-                if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
-                }
-            }
-
-        } else if (col == (image->numCols - 1)) {
-            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
-                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
-                if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
-                }
-            }
-
-        } else {
-            psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");
-        }
-    }
-    psFree (tmpRow);
-    psFree (row1);
-
-    //
-    // Exit if this image has a single row.
-    //
-    if (image->numRows == 1) {
-        psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-        return(list);
-    }
-
-    //
-    // Find peaks in interior rows only.
-    //
-    for (row = 1 ; row < (image->numRows - 1) ; row++) {
-        tmpRow = getRowVectorFromImage((psImage *) image, row);
-        row1 = pmFindVectorPeaks(tmpRow, threshold);
-
-        // Step through all local peaks in this row.
-        for (psU32 i = 0 ; i < row1->n ; i++ ) {
-            pmPeakType myType = PM_PEAK_UNDEF;
-            col = row1->data.U32[i];
-
-            if (col == 0) {
-                // If col==0, then we can not read col-1 pixels
-                if ((image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
-                    myType = PM_PEAK_EDGE;
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
-                }
-            } else if (col < (image->numCols - 1)) {
-                // This is an interior pixel
-                if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
-                        (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
-                        (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
-                    if (image->data.F32[row][col] > threshold) {
-                        if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) &&
-                                (image->data.F32[row][col] > image->data.F32[row-1][col]) &&
-                                (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&
-                                (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
-                                (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
-                                (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&
-                                (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
-                                (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {
-                            myType = PM_PEAK_LONE;
-                        }
-
-                        if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) ||
-                                (image->data.F32[row][col] == image->data.F32[row-1][col]) ||
-                                (image->data.F32[row][col] == image->data.F32[row-1][col+1]) ||
-                                (image->data.F32[row][col] == image->data.F32[row][col-1]) ||
-                                (image->data.F32[row][col] == image->data.F32[row][col+1]) ||
-                                (image->data.F32[row][col] == image->data.F32[row+1][col-1]) ||
-                                (image->data.F32[row][col] == image->data.F32[row+1][col]) ||
-                                (image->data.F32[row][col] == image->data.F32[row+1][col+1])) {
-                            myType = PM_PEAK_FLAT;
-                        }
-
-                        list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
-                    }
-                }
-            } else if (col == (image->numCols - 1)) {
-                // If col==numCols - 1, then we can not read col+1 pixels
-                if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
-                        (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                        (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col])) {
-                    myType = PM_PEAK_EDGE;
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
-                }
-            } else {
-                psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
-            }
-
-        }
-        psFree (tmpRow);
-        psFree (row1);
-    }
-
-    //
-    // Find peaks in the last row only.
-    //
-    row = image->numRows - 1;
-    tmpRow = getRowVectorFromImage((psImage *) image, row);
-    row1 = pmFindVectorPeaks(tmpRow, threshold);
-    for (psU32 i = 0 ; i < row1->n ; i++ ) {
-        col = row1->data.U32[i];
-        if (col == 0) {
-            if ( (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
-                if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
-                }
-            }
-        } else if (col < (image->numCols - 1)) {
-            if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row][col-1]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row][col+1])) {
-                if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
-                }
-            }
-
-        } else if (col == (image->numCols - 1)) {
-            if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row][col-1])) {
-                if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
-                }
-            }
-        } else {
-            psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
-        }
-    }
-    psFree (tmpRow);
-    psFree (row1);
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(list);
-}
-
-
-/******************************************************************************
-psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have
-a peak value above the given maximum, or fall outside the valid region.
- 
-XXX: Should the sky value be used when comparing the maximum?
- 
-XXX: warning message if valid is NULL?
- 
-XXX: changed API to create a NEW output psArray (should change name as well)
- 
-XXX: Do we free the psList elements of those culled peaks?
- 
-XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset?
-*****************************************************************************/
-psList *pmCullPeaks(psList *peaks,
-                    psF32 maxValue,
-                    const psRegion valid)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(peaks, NULL);
-
-    psListElem *tmpListElem = (psListElem *) peaks->head;
-    psS32 indexNum = 0;
-
-    //    printf("pmCullPeaks(): list size is %d\n", peaks->size);
-    while (tmpListElem != NULL) {
-        pmPeak *tmpPeak = (pmPeak *) tmpListElem->data;
-        if ((tmpPeak->counts > maxValue) ||
-                (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))) {
-            psListRemoveData(peaks, (psPtr) tmpPeak);
-        }
-
-        indexNum++;
-        tmpListElem = tmpListElem->next;
-    }
-
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(peaks);
-}
-
-// XXX EAM: I changed this to return a new, subset array
-//          rather than alter the existing one
-// XXX: Fix the *valid pointer.
-psArray *pmPeaksSubset(
-    psArray *peaks,
-    psF32 maxValue,
-    const psRegion valid)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(peaks, NULL);
-
-    psArray *output = psArrayAlloc (200);
-    output->n = 0;
-
-    psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n);
-
-    for (int i = 0; i < peaks->n; i++) {
-        pmPeak *tmpPeak = (pmPeak *) peaks->data[i];
-        if (tmpPeak->counts > maxValue)
-            continue;
-        if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))
-            continue;
-        psArrayAdd (output, 200, tmpPeak);
-    }
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(output);
-}
-
-/******************************************************************************
-pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this
-routine creates a new pmSource data structure and sets the following members:
-    ->pmPeak
-    ->pmMoments->sky
- 
-The sky value is set from the pixels in the square annulus surrounding the
-peak pixel.
- 
-We simply create a subSet image and mask the inner pixels, then call
-psImageStats on that subImage+mask.
- 
-XXX: The subImage has width of 1+2*outerRadius.  Verify with IfA.
- 
-XXX: Use static data structures for:
-     subImage
-     subImageMask
-     myStats
- 
-XXX: ensure that the inner and out radius fit in the actual image.  Should
-     we generate an error, or warning?  Currently an error.
- 
-XXX: Sync with IfA on whether the peak x/y coords are data structure coords,
-     or they use the image row/column offsets.
- 
-XXX: Should we simply set pmSource->peak = peak?  If so, should we increase
-the reference counter?  Or, should we copy the data structure?
- 
-XXX: Currently the subimage always has an even number of rows/columns.  Is
-     this correct?  Since there is a center pixel, maybe it should have an
-     odd number of rows/columns.
- 
-XXX: Use psTrace() for the print statements.
- 
-XXX: Don't use separate structs for the subimage and mask.  Use the source->
-     members.
-*****************************************************************************/
-
-bool pmSourceLocalSky(
-    pmSource *source,
-    psStatsOptions statsOptions,
-    psF32 Radius)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_IMAGE_NON_NULL(source->pixels, false);
-    PS_ASSERT_IMAGE_NON_NULL(source->mask, false);
-    PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_INT_POSITIVE(Radius, false);
-    PS_ASSERT_INT_NONNEGATIVE(Radius, false);
-
-    psImage *image = source->pixels;
-    psImage *mask  = source->mask;
-    pmPeak *peak  = source->peak;
-    psRegion srcRegion;
-
-    srcRegion = psRegionForSquare(peak->x, peak->y, Radius);
-    srcRegion = psRegionForImage(mask, srcRegion);
-
-    psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED);
-    psStats *myStats = psStatsAlloc(statsOptions);
-    myStats = psImageStats(myStats, image, mask, 0xff);
-    psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED);
-
-    psF64 tmpF64;
-    p_psGetStatValue(myStats, &tmpF64);
-    psFree(myStats);
-
-    if (isnan(tmpF64)) {
-        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
-        return(false);
-    }
-    source->moments = pmMomentsAlloc();
-    source->moments->Sky = (psF32) tmpF64;
-    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
-    return (true);
-}
-
-/******************************************************************************
-pmSourceMoments(source, radius): this function takes a subImage defined in the
-pmSource data structure, along with the peak location, and determines the
-various moments associated with that peak.
- 
-Requires the following to have been created:
-    pmSource
-    pmSource->peak
-    pmSource->pixels
-    pmSource->weight
-    pmSource->mask
- 
-XXX: The peak calculations are done in image coords, not subImage coords.
- 
-XXX EAM : this version clips input pixels on S/N
-XXX EAM : this version returns false for several reasons
-*****************************************************************************/
-# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
-
-bool pmSourceMoments(pmSource *source,
-                     psF32 radius)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
-    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
-
-    //
-    // XXX: Verify the setting for sky if source->moments == NULL.
-    //
-    psF32 sky = 0.0;
-    if (source->moments == NULL) {
-        source->moments = pmMomentsAlloc();
-    } else {
-        sky = source->moments->Sky;
-    }
-
-    //
-    // Sum = SUM (z - sky)
-    // X1  = SUM (x - xc)*(z - sky)
-    // X2  = SUM (x - xc)^2 * (z - sky)
-    // XY  = SUM (x - xc)*(y - yc)*(z - sky)
-    //
-    psF32 peakPixel = -PS_MAX_F32;
-    psS32 numPixels = 0;
-    psF32 Sum = 0.0;
-    psF32 X1 = 0.0;
-    psF32 Y1 = 0.0;
-    psF32 X2 = 0.0;
-    psF32 Y2 = 0.0;
-    psF32 XY = 0.0;
-    psF32 x  = 0;
-    psF32 y  = 0;
-    psF32 R2 = PS_SQR(radius);
-
-    psF32 xPeak = source->peak->x;
-    psF32 yPeak = source->peak->y;
-
-    // XXX why do I get different results for these two methods of finding Sx?
-    // XXX Sx, Sy would be better measured if we clip pixels close to sky
-    // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
-    // We loop through all pixels in this subimage (source->pixels), and for each
-    // pixel that is not masked, AND within the radius of the peak pixel, we
-    // proceed with the moments calculation.  need to do two loops for a
-    // numerically stable result.  first loop: get the sums.
-    // XXX EAM : mask == 0 is valid
-
-    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
-            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {
-                continue;
-            }
-
-            psF32 xDiff = col + source->pixels->col0 - xPeak;
-            psF32 yDiff = row + source->pixels->row0 - yPeak;
-
-            // XXX EAM : calculate xDiff, yDiff up front;
-            //           radius is just a function of (xDiff, yDiff)
-            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
-                continue;
-            }
-
-            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
-
-            // XXX EAM : check for valid S/N in pixel
-            // XXX EAM : should this limit be user-defined?
-            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {
-                continue;
-            }
-
-            Sum += pDiff;
-            X1  += xDiff * pDiff;
-            Y1  += yDiff * pDiff;
-            XY  += xDiff * yDiff * pDiff;
-
-            X2  += PS_SQR(xDiff) * pDiff;
-            Y2  += PS_SQR(yDiff) * pDiff;
-
-            peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel);
-            numPixels++;
-        }
-    }
-    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
-    if ((numPixels < 3) || (Sum <= 0)) {
-        psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
-        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
-        return (false);
-    }
-
-    psTrace (".psModules.pmSourceMoments", 5,
-             "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
-             sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
-
-    //
-    // first moment X  = X1/Sum + xc
-    // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
-    // Sxy             = XY / Sum
-    //
-    x = X1/Sum;
-    y = Y1/Sum;
-    if ((fabs(x) > radius) || (fabs(y) > radius)) {
-        psTrace (".psModules.pmSourceMoments", 5,
-                 "large centroid swing; invalid peak %d, %d\n",
-                 source->peak->x, source->peak->y);
-        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
-        return (false);
-    }
-
-    source->moments->x = x + xPeak;
-    source->moments->y = y + yPeak;
-
-    // XXX EAM : Sxy needs to have x*y subtracted
-    source->moments->Sxy = XY/Sum - x*y;
-    source->moments->Sum = Sum;
-    source->moments->Peak = peakPixel;
-    source->moments->nPixels = numPixels;
-
-    // XXX EAM : these values can be negative, so we need to limit the range
-    source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
-    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
-
-    psTrace (".psModules.pmSourceMoments", 4,
-             "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
-             sky, Sum, source->moments->x, source->moments->y,
-             source->moments->Sx, source->moments->Sy, source->moments->Sxy);
-
-    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
-    return(true);
-}
-
-// XXX EAM : I used
-int pmComparePeakAscend (const void **a, const void **b)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    pmPeak *A = *(pmPeak **)a;
-    pmPeak *B = *(pmPeak **)b;
-
-    psF32 diff;
-
-    diff = A->counts - B->counts;
-    if (diff < FLT_EPSILON) {
-        psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
-        return (-1);
-    } else if (diff > FLT_EPSILON) {
-        psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
-        return (+1);
-    }
-    psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);
-    return (0);
-}
-
-int pmComparePeakDescend (const void **a, const void **b)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    pmPeak *A = *(pmPeak **)a;
-    pmPeak *B = *(pmPeak **)b;
-
-    psF32 diff;
-
-    diff = A->counts - B->counts;
-    if (diff < FLT_EPSILON) {
-        psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
-        return (+1);
-    } else if (diff > FLT_EPSILON) {
-        psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
-        return (-1);
-    }
-    psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);
-    return (0);
-}
-
-/******************************************************************************
-pmSourcePSFClump(source, metadata): Find the likely PSF clump in the
-sigma-x, sigma-y plane. return 0,0 clump in case of error.
-*****************************************************************************/
-
-// XXX EAM include a S/N cutoff in selecting the sources?
-pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-
-    # define NPIX 10
-    # define SCALE 0.1
-
-    psArray *peaks  = NULL;
-    pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0};
-    pmPSFClump psfClump = emptyClump;
-
-    PS_ASSERT_PTR_NON_NULL(sources, emptyClump);
-    PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);
-
-    // find the sigmaX, sigmaY clump
-    {
-        psStats *stats  = NULL;
-        psImage *splane = NULL;
-        int binX, binY;
-        bool status;
-
-        psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
-        if (!status)
-            SX_MAX = 10.0;
-        psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
-        if (!status)
-            SY_MAX = 10.0;
-
-        // construct a sigma-plane image
-        // psImageAlloc does zero the data
-        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
-        for (int i = 0; i < splane->numRows; i++)
-        {
-            memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
-        }
-
-        // place the sources in the sigma-plane image (ignore 0,0 values?)
-        for (psS32 i = 0 ; i < sources->n ; i++)
-        {
-            pmSource *tmpSrc = (pmSource *) sources->data[i];
-            if (tmpSrc == NULL) {
-                continue;
-            }
-            if (tmpSrc->moments == NULL) {
-                continue;
-            }
-
-            // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
-            if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) {
-                continue;
-            }
-
-            // for the moment, force splane dimensions to be 10x10 image pix
-            binX = tmpSrc->moments->Sx/SCALE;
-            if (binX < 0)
-                continue;
-            if (binX >= splane->numCols)
-                continue;
-
-            binY = tmpSrc->moments->Sy/SCALE;
-            if (binY < 0)
-                continue;
-            if (binY >= splane->numRows)
-                continue;
-
-            splane->data.F32[binY][binX] += 1.0;
-        }
-
-        // find the peak in this image
-        stats = psStatsAlloc (PS_STAT_MAX);
-        stats = psImageStats (stats, splane, NULL, 0);
-        peaks = pmFindImagePeaks (splane, stats[0].max / 2);
-        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);
-
-        psFree (splane);
-        psFree (stats);
-
-    }
-    // XXX EAM : possible errors:
-    //           1) no peak in splane
-    //           2) no significant peak in splane
-
-    // measure statistics on Sx, Sy if Sx, Sy within range of clump
-    {
-        pmPeak *clump;
-        psF32 minSx, maxSx;
-        psF32 minSy, maxSy;
-        psVector *tmpSx = NULL;
-        psVector *tmpSy = NULL;
-        psStats *stats  = NULL;
-
-        // XXX EAM : this lets us takes the single highest peak
-        psArraySort (peaks, pmComparePeakDescend);
-        clump = peaks->data[0];
-        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);
-
-        // define section window for clump
-        minSx = clump->x * SCALE - 0.2;
-        maxSx = clump->x * SCALE + 0.2;
-        minSy = clump->y * SCALE - 0.2;
-        maxSy = clump->y * SCALE + 0.2;
-
-        tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32);
-        tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32);
-        tmpSx->n = 0;
-        tmpSy->n = 0;
-
-        // XXX clip sources based on flux?
-        // create vectors with Sx, Sy values in window
-        for (psS32 i = 0 ; i < sources->n ; i++)
-        {
-            pmSource *tmpSrc = (pmSource *) sources->data[i];
-
-            if (tmpSrc->moments->Sx < minSx)
-                continue;
-            if (tmpSrc->moments->Sx > maxSx)
-                continue;
-            if (tmpSrc->moments->Sy < minSy)
-                continue;
-            if (tmpSrc->moments->Sy > maxSy)
-                continue;
-            tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx;
-            tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy;
-            tmpSx->n++;
-            tmpSy->n++;
-            if (tmpSx->n == tmpSx->nalloc) {
-                psVectorRealloc (tmpSx, tmpSx->nalloc + 100);
-                psVectorRealloc (tmpSy, tmpSy->nalloc + 100);
-            }
-        }
-
-        // measures stats of Sx, Sy
-        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-
-        stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);
-        psfClump.X  = stats->clippedMean;
-        psfClump.dX = stats->clippedStdev;
-
-        stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);
-        psfClump.Y  = stats->clippedMean;
-        psfClump.dY = stats->clippedStdev;
-
-        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
-        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
-        // these values should be pushed on the metadata somewhere
-
-        psFree (stats);
-        psFree (peaks);
-        psFree (tmpSx);
-        psFree (tmpSy);
-    }
-
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return (psfClump);
-}
-
-/******************************************************************************
-pmSourceRoughClass(source, metadata): make a guess at the source
-classification.
- 
-XXX: push the clump info into the metadata?
- 
-XXX: How can this function ever return FALSE?
- 
-XXX EAM : add the saturated mask value to metadata
-*****************************************************************************/
-
-bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-
-    psBool rc = true;
-
-    int Nsat     = 0;
-    int Ngal     = 0;
-    int Nstar    = 0;
-    int Npsf     = 0;
-    int Ncr      = 0;
-    int Nsatstar = 0;
-    psRegion allArray = psRegionSet (0, 0, 0, 0);
-
-    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
-    starsn->n = 0;
-
-    // check return status value (do these exist?)
-    bool status;
-    psF32 RDNOISE    = psMetadataLookupF32 (&status, metadata, "RDNOISE");
-    psF32 GAIN       = psMetadataLookupF32 (&status, metadata, "GAIN");
-    psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
-    // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
-
-    // XXX allow clump size to be scaled relative to sigmas?
-    // make rough IDs based on clumpX,Y,DX,DY
-    for (psS32 i = 0 ; i < sources->n ; i++) {
-
-        pmSource *tmpSrc = (pmSource *) sources->data[i];
-
-        tmpSrc->peak->class = 0;
-
-        psF32 sigX = tmpSrc->moments->Sx;
-        psF32 sigY = tmpSrc->moments->Sy;
-
-        // calculate and save signal-to-noise estimates
-        psF32 S  = tmpSrc->moments->Sum;
-        psF32 A  = 4 * M_PI * sigX * sigY;
-        psF32 B  = tmpSrc->moments->Sky;
-        psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));
-        psF32 SN = (S * sqrt(GAIN) / RT);
-        tmpSrc->moments->SN = SN;
-
-        // XXX EAM : can we use the value of SATURATE if mask is NULL?
-        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED);
-
-        // saturated star (size consistent with PSF or larger)
-        // Nsigma should be user-configured parameter
-        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
-        if ((Nsatpix > 1) && big) {
-            tmpSrc->type = PM_SOURCE_SATSTAR;
-            Nsatstar ++;
-            continue;
-        }
-
-        // saturated object (not a star, eg bleed trails, hot pixels)
-        if (Nsatpix > 1) {
-            tmpSrc->type = PM_SOURCE_SATURATED;
-            Nsat ++;
-            continue;
-        }
-
-        // likely defect (too small to be stellar) (push out to 3 sigma)
-        // low S/N objects which are small are probably stellar
-        // only set candidate defects if
-        if ((sigX < 0.05) || (sigY < 0.05)) {
-            tmpSrc->type = PM_SOURCE_DEFECT;
-            Ncr ++;
-            continue;
-        }
-
-        // likely unsaturated galaxy (too large to be stellar)
-        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
-            tmpSrc->type = PM_SOURCE_GALAXY;
-            Ngal ++;
-            continue;
-        }
-
-        // the rest are probable stellar objects
-        starsn->data.F32[starsn->n] = SN;
-        starsn->n ++;
-        Nstar ++;
-
-        // PSF star (within 1.5 sigma of clump center
-        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
-        if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
-            tmpSrc->type = PM_SOURCE_PSFSTAR;
-            Npsf ++;
-            continue;
-        }
-
-        // random type of star
-        tmpSrc->type = PM_SOURCE_OTHER;
-    }
-
-    {
-        psStats *stats  = NULL;
-        stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-        stats = psVectorStats (stats, starsn, NULL, NULL, 0);
-        psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
-        psFree (stats);
-        psFree (starsn);
-    }
-
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:    %3d\n", Nstar);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:     %3d\n", Npsf);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:     %3d\n", Ngal);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:     %3d\n", Nsat);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:      %3d\n", Ncr);
-
-    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
-    return(rc);
-}
-
-/** pmSourceDefinePixels()
- *
- * Define psImage subarrays for the source located at coordinates x,y on the
- * image set defined by readout. The pixels defined by this operation consist of
- * a square window (of full width 2Radius+1) centered on the pixel which contains
- * the given coordinate, in the frame of the readout. The window is defined to
- * have limits which are valid within the boundary of the readout image, thus if
- * the radius would fall outside the image pixels, the subimage is truncated to
- * only consist of valid pixels. If readout->mask or readout->weight are not
- * NULL, matching subimages are defined for those images as well. This function
- * fails if no valid pixels can be defined (x or y less than Radius, for
- * example). This function should be used to define a region of interest around a
- * source, including both source and sky pixels.
- *
- * XXX: must code this.
- *
- */
-bool pmSourceDefinePixels(
-    pmSource *mySource,                 ///< Add comment.
-    pmReadout *readout,                 ///< Add comment.
-    psF32 x,                            ///< Add comment.
-    psF32 y,                            ///< Add comment.
-    psF32 Radius)                       ///< Add comment.
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented.  Returning FALSE.\n");
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(false);
-}
-
-/******************************************************************************
-pmSourceSetPixelsCircle(source, image, radius)
- 
-XXX: This was replaced by DefinePixels in SDRS.  Remove it.
-*****************************************************************************/
-bool pmSourceSetPixelsCircle(pmSource *source,
-                             const psImage *image,
-                             psF32 radius)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_IMAGE_NON_NULL(image, false);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->moments, false);
-    PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(radius, 0.0, false);
-
-    //
-    // We define variables for code readability.
-    //
-    // XXX: Since the peak->xy coords are in image, not subImage coords,
-    // these variables should be renamed for clarity (imageCenterRow, etc).
-    //
-    psS32 radiusS32 = (psS32) radius;
-    psS32 SubImageCenterRow = source->peak->y;
-    psS32 SubImageCenterCol = source->peak->x;
-    // XXX EAM : for the circle to stay on the image
-    // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
-    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - radiusS32);
-    psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1);
-    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
-    psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);
-
-    // XXX: Must recycle image.
-    // XXX EAM: this message reflects a programming error we know about.
-    //          i am setting it to a trace message which we can take out
-    if (source->pixels != NULL) {
-        psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4,
-                 "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL.  Freeing and reallocating.\n");
-        psFree(source->pixels);
-    }
-    source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
-                                   SubImageStartRow,
-                                   SubImageEndCol,
-                                   SubImageEndRow));
-
-    // XXX: Must recycle image.
-    if (source->mask != NULL) {
-        psFree(source->mask);
-    }
-    source->mask = psImageAlloc(source->pixels->numCols,
-                                source->pixels->numRows,
-                                PS_TYPE_U8); // XXX EAM : type was F32
-
-    //
-    // Loop through the subimage mask, initialize mask to 0 or 1.
-    // XXX EAM: valid pixels should have 0, not 1
-    for (psS32 row = 0 ; row < source->mask->numRows; row++) {
-        for (psS32 col = 0 ; col < source->mask->numCols; col++) {
-
-            if (checkRadius2((psF32) radiusS32,
-                             (psF32) radiusS32,
-                             radius,
-                             (psF32) col,
-                             (psF32) row)) {
-                source->mask->data.U8[row][col] = 0;
-            } else {
-                source->mask->data.U8[row][col] = 1;
-            }
-        }
-    }
-    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
-    return(true);
-}
-
-/******************************************************************************
-pmSourceModelGuess(source, model): This function allocates a new
-pmModel structure based on the given modelType specified in the argument list.
-The corresponding pmModelGuess function is returned, and used to
-supply the values of the params array in the pmModel structure.
- 
-XXX: Many parameters are based on the src->moments structure, which is in
-image, not subImage coords.  Therefore, the calls to the model evaluation
-functions will be in image, not subImage coords.  Remember this.
-*****************************************************************************/
-pmModel *pmSourceModelGuess(pmSource *source,
-                            pmModelType modelType)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source->moments, false);
-    PS_ASSERT_PTR_NON_NULL(source->peak, false);
-
-    pmModel *model = pmModelAlloc(modelType);
-
-    pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType);
-    modelGuessFunc(model, source);
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(model);
-}
-
-/******************************************************************************
-evalModel(source, level, row): a private function which evaluates the
-source->modelPSF function at the specified coords.  The coords are subImage, not
-image coords.
- 
-NOTE: The coords are in subImage source->pixel coords, not image coords.
- 
-XXX: reverse order of row,col args?
- 
-XXX: rename all coords in this file such that their name defines whether
-the coords is in subImage or image space.
- 
-XXX: This should probably be a public pmModules function.
- 
-XXX: Use static vectors for x.
- 
-XXX: Figure out if it's (row, col) or (col, row) for the model functions.
- 
-XXX: For a while, the first psVectorAlloc() was generating a seg fault during
-testing.  Try to reproduce that and debug.
-*****************************************************************************/
-
-// XXX EAM : I have made this a public function
-// XXX EAM : this now uses a pmModel as the input
-// XXX EAM : it was using src->type to find the model, not model->type
-psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(image, false);
-    PS_ASSERT_PTR_NON_NULL(model, false);
-    PS_ASSERT_PTR_NON_NULL(model->params, false);
-
-    // Allocate the x coordinate structure and convert row/col to image space.
-    //
-    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
-    x->n = x->nalloc;
-    x->data.F32[0] = (psF32) (col + image->col0);
-    x->data.F32[1] = (psF32) (row + image->row0);
-    psF32 tmpF;
-    pmModelFunc modelFunc;
-
-    modelFunc = pmModelFunc_GetFunction (model->type);
-    tmpF = modelFunc (NULL, model->params, x);
-    psFree(x);
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(tmpF);
-}
-
-/******************************************************************************
-pmSourceContour(src, img, level, mode): For an input subImage, and model, this
-routine returns a psArray of coordinates that evaluate to the specified level.
- 
-XXX: Probably should remove the "image" argument.
-XXX: What type should the output coordinate vectors consist of?  col,row?
-XXX: Why a pmArray output?
-XXX: doex x,y correspond with col,row or row/col?
-XXX: What is mode?
-XXX: The top, bottom of the contour is not correctly determined.
-XXX EAM : this function is using the model for the contour, but it should
-          be using only the image counts
-*****************************************************************************/
-psArray *pmSourceContour(pmSource *source,
-                         const psImage *image,
-                         psF32 level,
-                         pmContourType mode)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(image, false);
-    PS_ASSERT_PTR_NON_NULL(source->moments, false);
-    PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->modelFLT, false);
-    // XXX EAM : what is the purpose of modelPSF/modelFLT?
-
-    //
-    // Allocate data for x/y pairs.
-    //
-    psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);
-    psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);
-    xVec->n = xVec->nalloc;
-    yVec->n = yVec->nalloc;
-    //
-    // Start at the row with peak pixel, then decrement.
-    //
-    psS32 col = source->peak->x;
-    for (psS32 row = source->peak->y; row>= 0 ; row--) {
-        // XXX: yVec contain no real information.  Do we really need it?
-        yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);
-        yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);
-
-        // Starting at peak pixel, search leftwards for the column intercept.
-        psF32 leftIntercept = findValue(source, level, row, col, 0);
-        if (isnan(leftIntercept)) {
-            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
-            psFree(xVec);
-            psFree(yVec);
-            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
-            return(NULL);
-            //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
-        }
-        xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;
-
-        // Starting at peak pixel, search rightwards for the column intercept.
-
-        psF32 rightIntercept = findValue(source, level, row, col, 1);
-        if (isnan(rightIntercept)) {
-            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
-            psFree(xVec);
-            psFree(yVec);
-            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
-            return(NULL);
-            //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
-        }
-        psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);
-        xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;
-
-        // Set starting column for next row
-        col = (psS32) ((leftIntercept + rightIntercept) / 2.0);
-    }
-    //
-    // Start at the row (+1) with peak pixel, then increment.
-    //
-    col = source->peak->x;
-    for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) {
-        // XXX: yVec contain no real information.  Do we really need it?
-        yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);
-        yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);
-
-        // Starting at peak pixel, search leftwards for the column intercept.
-        psF32 leftIntercept = findValue(source, level, row, col, 0);
-        if (isnan(leftIntercept)) {
-            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
-            psFree(xVec);
-            psFree(yVec);
-            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
-            return(NULL);
-            //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
-        }
-        xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;
-
-        // Starting at peak pixel, search rightwards for the column intercept.
-        psF32 rightIntercept = findValue(source, level, row, col, 1);
-        if (isnan(rightIntercept)) {
-            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
-            psFree(xVec);
-            psFree(yVec);
-            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
-            return(NULL);
-            //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
-        }
-        xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;
-
-        // Set starting column for next row
-        col = (psS32) ((leftIntercept + rightIntercept) / 2.0);
-    }
-
-    //
-    // Allocate an array for result, store coord vectors there.
-    //
-    psArray *tmpArray = psArrayAlloc(2);
-    tmpArray->n = 2;
-    tmpArray->data[0] = (psPtr *) yVec;
-    tmpArray->data[1] = (psPtr *) xVec;
-    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
-    return(tmpArray);
-}
-
-// XXX EAM : these are better starting values, but should be available from metadata?
-#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
-#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
-/******************************************************************************
-pmSourceFitModel(source, model): must create the appropiate arguments to the
-LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
- 
-XXX: should there be a mask value?
-XXX EAM : fit the specified model (not necessarily the one in source)
-*****************************************************************************/
-bool pmSourceFitModel_v5(pmSource *source,
-                         pmModel *model,
-                         const bool PSF)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->moments, false);
-    PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
-    PS_ASSERT_PTR_NON_NULL(source->weight, false);
-    psBool fitStatus = true;
-    psBool onPic     = true;
-    psBool rc        = true;
-
-    // XXX EAM : is it necessary for the mask & weight to exist?  the
-    //           tests below could be conditions (!NULL)
-
-    psVector *params = model->params;
-    psVector *dparams = model->dparams;
-    psVector *paramMask = NULL;
-
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-
-    int nParams = PSF ? params->n - 4 : params->n;
-
-    // find the number of valid pixels
-    // XXX EAM : this loop and the loop below could just be one pass
-    //           using the psArrayAdd and psVectorExtend functions
-    psS32 count = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                count++;
-            }
-        }
-    }
-    if (count <  nParams + 1) {
-        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
-        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
-        model->status = PM_MODEL_BADARGS;
-        return(false);
-    }
-
-    // construct the coordinate and value entries
-    psArray *x = psArrayAlloc(count);
-    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
-    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
-    x->n = x->nalloc;
-    y->n = y->nalloc;
-    yErr->n = yErr->nalloc;
-    psS32 tmpCnt = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-                coord->n = 2;
-                // XXX: Convert i/j to image space:
-                // XXX EAM: coord order is (x,y) == (col,row)
-                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
-                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
-                x->data[tmpCnt] = (psPtr *) coord;
-                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
-                yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]);
-                // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
-                //           the minimization function calculates sq()
-                tmpCnt++;
-            }
-        }
-    }
-
-    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
-                            PM_SOURCE_FIT_MODEL_TOLERANCE);
-
-    // PSF model only fits first 4 parameters, FLT model fits all
-    if (PSF) {
-        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
-        paramMask->n = paramMask->nalloc;
-        for (int i = 0; i < 4; i++) {
-            paramMask->data.U8[i] = 0;
-        }
-        for (int i = 4; i < paramMask->n; i++) {
-            paramMask->data.U8[i] = 1;
-        }
-    }
-
-    // XXX EAM : covar must be F64?
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
-
-    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
-
-    psMinConstrain *constrain = psMinConstrainAlloc();
-    constrain->paramMask = paramMask;
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain,
-                                 x, y, yErr, modelFunc);
-    psFree(constrain);
-    for (int i = 0; i < dparams->n; i++) {
-        if ((paramMask != NULL) && paramMask->data.U8[i])
-            continue;
-        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
-    }
-
-    // save the resulting chisq, nDOF, nIter
-    model->chisq = myMin->value;
-    model->nIter = myMin->iter;
-    model->nDOF  = y->n - nParams;
-
-    // get the Gauss-Newton distance for fixed model parameters
-    if (paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
-        delta->n = delta->nalloc;
-        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
-        for (int i = 0; i < dparams->n; i++) {
-            if (!paramMask->data.U8[i])
-                continue;
-            dparams->data.F32[i] = delta->data.F64[i];
-        }
-        psFree (delta);
-    }
-
-    // set the model success or failure status
-    if (!fitStatus) {
-        model->status = PM_MODEL_NONCONVERGE;
-    } else {
-        model->status = PM_MODEL_SUCCESS;
-    }
-
-    // models can go insane: reject these
-    onPic &= (params->data.F32[2] >= source->pixels->col0);
-    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
-    onPic &= (params->data.F32[3] >= source->pixels->row0);
-    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
-    if (!onPic) {
-        model->status = PM_MODEL_OFFIMAGE;
-    }
-
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myMin);
-    psFree(covar);
-    psFree(paramMask);
-
-    rc = (onPic && fitStatus);
-    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
-    return(rc);
-}
-
-// XXX EAM : new version with parameter range limits and weight enhancement
-bool pmSourceFitModel (pmSource *source,
-                       pmModel *model,
-                       const bool PSF)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->moments, false);
-    PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
-    PS_ASSERT_PTR_NON_NULL(source->weight, false);
-
-    // XXX EAM : is it necessary for the mask & weight to exist?  the
-    //           tests below could be conditions (!NULL)
-
-    psBool fitStatus = true;
-    psBool onPic     = true;
-    psBool rc        = true;
-    psF32  Ro, ymodel;
-
-    psVector *params = model->params;
-    psVector *dparams = model->dparams;
-    psVector *paramMask = NULL;
-
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-
-    // XXX EAM : I need to use the sky value to constrain the weight model
-    int nParams = PSF ? params->n - 4 : params->n;
-    psF32 So = params->data.F32[0];
-
-    // find the number of valid pixels
-    // XXX EAM : this loop and the loop below could just be one pass
-    //           using the psArrayAdd and psVectorExtend functions
-    psS32 count = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                count++;
-            }
-        }
-    }
-    if (count <  nParams + 1) {
-        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
-        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
-        model->status = PM_MODEL_BADARGS;
-        return(false);
-    }
-
-    // construct the coordinate and value entries
-    psArray *x = psArrayAlloc(count);
-    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
-    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
-    x->n = x->nalloc;
-    y->n = y->nalloc;
-    yErr->n = yErr->nalloc;
-    psS32 tmpCnt = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-                coord->n = 2;
-                // XXX: Convert i/j to image space:
-                // XXX EAM: coord order is (x,y) == (col,row)
-                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
-                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
-                x->data[tmpCnt] = (psPtr *) coord;
-                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
-
-                // compare observed flux to model flux to adjust weight
-                ymodel = modelFunc (NULL, model->params, coord);
-
-                // this test enhances the weight based on deviation from the model flux
-                Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So)
-                        + PS_SQR(So));
-
-                // psMinimizeLMChi2_EAM takes wt = 1/dY^2
-                if (source->weight->data.F32[i][j] == 0) {
-                    yErr->data.F32[tmpCnt] = 0.0;
-                } else {
-                    yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
-                }
-                tmpCnt++;
-            }
-        }
-    }
-
-    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
-                            PM_SOURCE_FIT_MODEL_TOLERANCE);
-
-    // PSF model only fits first 4 parameters, FLT model fits all
-    if (PSF) {
-        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
-        paramMask->n = paramMask->nalloc;
-        for (int i = 0; i < 4; i++) {
-            paramMask->data.U8[i] = 0;
-        }
-        for (int i = 4; i < paramMask->n; i++) {
-            paramMask->data.U8[i] = 1;
-        }
-    }
-
-    // XXX EAM : I've added three types of parameter range checks
-    // XXX EAM : this requires my new psMinimization functions
-    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
-    psVector *beta_lim = NULL;
-    psVector *params_min = NULL;
-    psVector *params_max = NULL;
-
-    // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
-    //           in the SDRS, I define a new psMinimization which will take these in
-    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
-    modelLimits (&beta_lim, &params_min, &params_max);
-    for (int i = 0; i < params->n; i++) {
-        covar->data.F64[0][i] = beta_lim->data.F32[i];
-        covar->data.F64[1][i] = params_min->data.F32[i];
-        covar->data.F64[2][i] = params_max->data.F32[i];
-    }
-
-    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
-    psMinConstrain *constrain = psMinConstrainAlloc();
-    constrain->paramMask = paramMask;
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain,
-                                 x, y, yErr, modelFunc);
-    psFree(constrain);
-    for (int i = 0; i < dparams->n; i++) {
-        if ((paramMask != NULL) && paramMask->data.U8[i])
-            continue;
-        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
-    }
-
-    // XXX EAM: we need to do something (give an error?) if rc is false
-    // XXX EAM: psMinimizeLMChi2 does not check convergence
-
-    // XXX EAM: save the resulting chisq, nDOF, nIter
-    model->chisq = myMin->value;
-    model->nIter = myMin->iter;
-    model->nDOF  = y->n - nParams;
-
-    // XXX EAM get the Gauss-Newton distance for fixed model parameters
-    if (paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
-        delta->n = delta->nalloc;
-        psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);
-        for (int i = 0; i < dparams->n; i++) {
-            if (!paramMask->data.U8[i])
-                continue;
-            dparams->data.F32[i] = delta->data.F64[i];
-        }
-    }
-
-    // set the model success or failure status
-    if (!fitStatus) {
-        model->status = PM_MODEL_NONCONVERGE;
-    } else {
-        model->status = PM_MODEL_SUCCESS;
-    }
-
-    // XXX models can go insane: reject these
-    onPic &= (params->data.F32[2] >= source->pixels->col0);
-    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
-    onPic &= (params->data.F32[3] >= source->pixels->row0);
-    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
-    if (!onPic) {
-        model->status = PM_MODEL_OFFIMAGE;
-    }
-
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myMin);
-    psFree(covar);
-    psFree(paramMask);
-
-    rc = (onPic && fitStatus);
-    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
-    return(rc);
-}
-
-bool p_pmSourceAddOrSubModel(psImage *image,
-                             psImage *mask,
-                             pmModel *model,
-                             bool center,
-                             bool sky,
-                             bool add
-                                )
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-
-    PS_ASSERT_PTR_NON_NULL(model, false);
-    PS_ASSERT_IMAGE_NON_NULL(image, false);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
-
-    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
-    x->n = 2;
-    psVector *params = model->params;
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-    psS32 imageCol;
-    psS32 imageRow;
-    psF32 skyValue = params->data.F32[0];
-    psF32 pixelValue;
-
-    for (psS32 i = 0; i < image->numRows; i++) {
-        for (psS32 j = 0; j < image->numCols; j++) {
-            if ((mask != NULL) && mask->data.U8[i][j])
-                continue;
-
-            // XXX: Should you be adding the pixels for the entire subImage,
-            // or a radius of pixels around it?
-
-            // Convert i/j to imace coord space:
-            // XXX: Make sure you have col/row order correct.
-            // XXX EAM : 'center' option changes this
-            // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2]
-            if (center) {
-                imageCol = j - 0.5*image->numCols + model->params->data.F32[2];
-                imageRow = i - 0.5*image->numRows + model->params->data.F32[3];
-            } else {
-                imageCol = j + image->col0;
-                imageRow = i + image->row0;
-            }
-
-            x->data.F32[0] = (float) imageCol;
-            x->data.F32[1] = (float) imageRow;
-
-            // set the appropriate pixel value for this coordinate
-            if (sky) {
-                pixelValue = modelFunc (NULL, params, x);
-            } else {
-                pixelValue = modelFunc (NULL, params, x) - skyValue;
-            }
-
-
-            // add or subtract the value
-            if (add
-               ) {
-                image->data.F32[i][j] += pixelValue;
-            }
-            else {
-                image->data.F32[i][j] -= pixelValue;
-            }
-        }
-    }
-    psFree(x);
-    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
-    return(true);
-}
-
-
-
-/******************************************************************************
- *****************************************************************************/
-bool pmSourceAddModel(psImage *image,
-                      psImage *mask,
-                      pmModel *model,
-                      bool center,
-                      bool sky)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);
-    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
-    return(rc);
-}
-
-/******************************************************************************
- *****************************************************************************/
-bool pmSourceSubModel(psImage *image,
-                      psImage *mask,
-                      pmModel *model,
-                      bool center,
-                      bool sky)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);
-    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
-    return(rc);
-}
-
-bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask)
-{
-
-    float obsSum = 0;
-    float fitSum = 0;
-    float sky = model->params->data.F32[0];
-
-    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
-    fitSum = modelFluxFunc (model->params);
-
-    for (int ix = 0; ix < image->numCols; ix++) {
-        for (int iy = 0; iy < image->numRows; iy++) {
-            if (mask->data.U8[iy][ix])
-                continue;
-            obsSum += image->data.F32[iy][ix] - sky;
-        }
-    }
-    if (obsSum <= 0)
-        return false;
-    if (fitSum <= 0)
-        return false;
-
-    *fitMag = -2.5*log10(fitSum);
-    *obsMag = -2.5*log10(obsSum);
-    return (true);
-}
-

trunk/psModules/src/objects/pmPSF.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-04 01:01:33 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -18 +18 @@
 
 #include <pslib.h>
-#include "pmObjects.h"
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmModel.h"
+#include "pmSource.h"
+#include "pmGrowthCurve.h"
 #include "pmPSF.h"
 #include "pmModelGroup.h"
@@ -57 +63 @@
         return;
 
+    psFree (psf->ApTrend);
+    psFree (psf->growth);
     psFree (psf->params);
     return;
@@ -69 +77 @@
  X-center
  Y-center
- ???: Sky background value?
- ???: Zo?
+ Sky background value
+ Object Normalization
  *****************************************************************************/
-pmPSF *pmPSFAlloc (pmModelType type)
+pmPSF *pmPSFAlloc (pmModelType type, bool poissonErrors)
 {
     int Nparams;
@@ -83 +91 @@
     psf->dApResid = 0.0;
     psf->skyBias  = 0.0;
+    psf->skySat   = 0.0;
+    psf->poissonErrors = poissonErrors;
+
+    // the ApTrend components are (x, y, r2rflux, flux)
+    psf->ApTrend = psPolynomial4DAlloc (PS_POLYNOMIAL_ORD, 2, 2, 1, 1);
+    pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
+
+    if (psf->poissonErrors) {
+        psf->ChiTrend = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
+    } else {
+        psf->ChiTrend = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
+    }
+
+    // don't define a growth curve : user needs to choose radius bins
+    psf->growth = NULL;
 
     Nparams = pmModelParameterCount (type);
@@ -94 +117 @@
     for (int i = 0; i < psf->params->n; i++) {
         // XXX EAM : make this a user-defined value?
-        psf->params->data[i] = psPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
+        psf->params->data[i] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1);
     }
 
@@ -174 +197 @@
 
 /*****************************************************************************
-pmModelFromPSF (*modelFLT, *psf):  use the model position parameters to
+pmModelFromPSF (*modelEXT, *psf):  use the model position parameters to
 construct a realization of the PSF model at the object coordinates
  *****************************************************************************/
-pmModel *pmModelFromPSF (pmModel *modelFLT, pmPSF *psf)
-{
-
-    // need to define the relationship between the modelFLT and modelPSF ?
+pmModel *pmModelFromPSF (pmModel *modelEXT, pmPSF *psf)
+{
+
+    // need to define the relationship between the modelEXT and modelPSF ?
 
     // find function used to set the model parameters
@@ -189 +212 @@
 
     // set model parameters for this source based on PSF information
-    modelFromPSFFunc (modelPSF, modelFLT, psf);
+    modelFromPSFFunc (modelPSF, modelEXT, psf);
 
     return (modelPSF);
 }
+
+// zero and mask out all terms:
+static bool maskAllTerms (psPolynomial4D *trend)
+{
+
+    for (int i = 0; i < trend->nX + 1; i++) {
+        for (int j = 0; j < trend->nY + 1; j++) {
+            for (int k = 0; k < trend->nZ + 1; k++) {
+                for (int m = 0; m < trend->nT + 1; m++) {
+                    trend->mask[i][j][k][m] = 1;  // mask coeff
+                    trend->coeff[i][j][k][m] = 0;  // zero coeff
+                }
+            }
+        }
+    }
+    return true;
+}
+
+/***********************************************
+ * this function masks the psf.ApTrend polynomial 
+ * to enable the specific subset of the coefficients
+ **********************************************/
+bool pmPSF_MaskApTrend (pmPSF *psf, pmPSF_ApTrendOptions option)
+{
+
+    switch (option) {
+    case PM_PSF_NONE:
+        maskAllTerms (psf->ApTrend);
+        return true;
+
+    case PM_PSF_CONSTANT:
+        maskAllTerms (psf->ApTrend);
+        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
+        return true;
+
+    case PM_PSF_SKYBIAS:
+        maskAllTerms (psf->ApTrend);
+        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
+        psf->ApTrend->mask[0][0][1][0] = 0;  // unmask skybias
+        return true;
+
+    case PM_PSF_SKYSAT:
+        maskAllTerms (psf->ApTrend);
+        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
+        psf->ApTrend->mask[0][0][1][0] = 0;  // unmask skybias
+        psf->ApTrend->mask[0][0][0][1] = 0;  // unmask skybias
+        return true;
+
+    case PM_PSF_XY_LIN:
+        maskAllTerms (psf->ApTrend);
+        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
+        psf->ApTrend->mask[1][0][0][0] = 0;  // unmask x
+        psf->ApTrend->mask[0][1][0][0] = 0;  // unmask y
+        return true;
+
+    case PM_PSF_XY_QUAD:
+        maskAllTerms (psf->ApTrend);
+        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
+        psf->ApTrend->mask[1][0][0][0] = 0;  // unmask x
+        psf->ApTrend->mask[2][0][0][0] = 0;  // unmask x^2
+        psf->ApTrend->mask[1][1][0][0] = 0;  // unmask x y
+        psf->ApTrend->mask[0][1][0][0] = 0;  // unmask y
+        psf->ApTrend->mask[0][2][0][0] = 0;  // unmask y^2
+        return true;
+
+    case PM_PSF_SKY_XY_LIN:
+        maskAllTerms (psf->ApTrend);
+        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
+        psf->ApTrend->mask[1][0][0][0] = 0;  // unmask x
+        psf->ApTrend->mask[0][1][0][0] = 0;  // unmask y
+        psf->ApTrend->mask[0][0][1][0] = 0;  // unmask skybias
+        return true;
+
+    case PM_PSF_SKYSAT_XY_LIN:
+        maskAllTerms (psf->ApTrend);
+        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
+        psf->ApTrend->mask[1][0][0][0] = 0;  // unmask x
+        psf->ApTrend->mask[0][1][0][0] = 0;  // unmask y
+        psf->ApTrend->mask[0][0][1][0] = 0;  // unmask skybias
+        psf->ApTrend->mask[0][0][0][1] = 0;  // unmask skysat
+        return true;
+
+    case PM_PSF_ALL:
+    default:
+        maskAllTerms (psf->ApTrend);
+        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
+        psf->ApTrend->mask[0][0][1][0] = 0;  // unmask skybias
+        psf->ApTrend->mask[0][0][0][1] = 0;  // unmask skysat
+
+        psf->ApTrend->mask[1][0][0][0] = 0;  // unmask x
+        psf->ApTrend->mask[2][0][0][0] = 0;  // unmask x^2
+        psf->ApTrend->mask[1][1][0][0] = 0;  // unmask x y
+        psf->ApTrend->mask[0][1][0][0] = 0;  // unmask y
+        psf->ApTrend->mask[0][2][0][0] = 0;  // unmask y^2
+        return true;
+    }
+    return false;
+}

trunk/psModules/src/objects/pmPSF.h

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-04 01:01:33 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-17 18:10:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -13 +13 @@
 
 # include <pslib.h>
-# include "pmObjects.h"
-# include "pmPSF.h"
-# include "pmPSFtry.h"
-# include "pmModelGroup.h"
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmModel.h"
+#include "pmSource.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
+#include "pmPSFtry.h"
+#include "pmModelGroup.h"
+#include "pmSourceFitModel.h"
+#include "pmSourcePhotometry.h"
 
 // ********  pmPSFtry functions  **************************************************
@@ -33 +41 @@
     psFree (test->psf);
     psFree (test->sources);
-    psFree (test->modelFLT);
+    psFree (test->modelEXT);
     psFree (test->modelPSF);
     psFree (test->metric);
@@ -42 +50 @@
 
 // allocate a pmPSFtry based on the desired sources and the model (identified by name)
-pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName)
+pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors)
 {
 
@@ -51 +59 @@
     // XXX probably need to increment ref counter
     type           = pmModelSetType (modelName);
-    test->psf      = pmPSFAlloc (type);
+    test->psf      = pmPSFAlloc (type, poissonErrors);
     test->sources  = psMemIncrRefCounter(sources);
-    test->modelFLT = psArrayAlloc (sources->n);
+    test->modelEXT = psArrayAlloc (sources->n);
     test->modelPSF = psArrayAlloc (sources->n);
     test->metric   = psVectorAlloc (sources->n, PS_TYPE_F64);
@@ -59 +67 @@
     test->mask     = psVectorAlloc (sources->n, PS_TYPE_U8);
 
-    test->modelFLT->n = test->modelFLT->nalloc;
-    test->modelPSF->n = test->modelPSF->nalloc;
-    test->metric->n   = test->metric->nalloc;
-    test->fitMag->n   = test->fitMag->nalloc;
-    test->mask->n     = test->mask->nalloc;
-
-    for (int i = 0; i < test->modelFLT->n; i++) {
+    for (int i = 0; i < test->modelEXT->n; i++) {
         test->mask->data.U8[i]  = 0;
-        test->modelFLT->data[i] = NULL;
+        test->modelEXT->data[i] = NULL;
         test->modelPSF->data[i] = NULL;
         test->metric->data.F64[i] = 0;
@@ -86 +88 @@
 // mask values indicate the reason the source was rejected:
 
-pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS)
+pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors)
 {
     bool status;
@@ -93 +95 @@
     float x;
     float y;
-    int Nflt = 0;
+    int Next = 0;
     int Npsf = 0;
 
-    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName);
+    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors);
 
     // stage 1:  fit an independent model (freeModel) to all sources
@@ -108 +110 @@
 
         // set temporary object mask and fit object
-        // fit model as FLT, not PSF
-        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
-        status = pmSourceFitModel (source, model, false);
-        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
+        // fit model as EXT, not PSF
+
+        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED);
+        status = pmSourceFitModel (source, model, PM_SOURCE_FIT_EXT);
+        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PM_SOURCE_MASK_MARKED);
 
         // exclude the poor fits
         if (!status) {
-            psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
+            psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL;
             psFree (model);
             continue;
         }
-        psfTry->modelFLT->data[i] = model;
-        Nflt ++;
-    }
-    psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
-    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
-
-    // make this optional?
-    // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat");
+        psfTry->modelEXT->data[i] = model;
+        Next ++;
+    }
+    psLogMsg ("psphot.psftry", 4, "fit ext:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (EXT)\n", Next, sources->n);
 
     // stage 2: construct a psf (pmPSF) from this collection of model fits
-    pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask);
+    pmPSFFromModels (psfTry->psf, psfTry->modelEXT, psfTry->mask);
 
     // stage 3: refit with fixed shape parameters
@@ -139 +139 @@
 
         pmSource *source = psfTry->sources->data[i];
-        pmModel  *modelFLT = psfTry->modelFLT->data[i];
+        pmModel  *modelEXT = psfTry->modelEXT->data[i];
 
         // set shape for this model based on PSF
-        pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf);
+        pmModel *modelPSF = pmModelFromPSF (modelEXT, psfTry->psf);
         x = source->peak->x;
         y = source->peak->y;
 
-        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
-        status = pmSourceFitModel (source, modelPSF, true);
+        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED);
+        status = pmSourceFitModel (source, modelPSF, PM_SOURCE_FIT_PSF);
 
         // skip poor fits
@@ -162 +162 @@
         // XXX : use a different estimator for the local sky?
         // XXX : pass 'source' as input?
-        if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
+        if (!pmSourcePhotometryModel (&fitMag, modelPSF)) {
+            psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
+            goto next_source;
+        }
+        if (!pmSourcePhotometryAper (&obsMag, modelPSF, source->pixels, source->mask)) {
             psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
             goto next_source;
@@ -172 +176 @@
 
 next_source:
-        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
-
-    }
+        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PM_SOURCE_MASK_MARKED);
+
+    }
+    psfTry->psf->nPSFstars = Npsf;
+
     psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
 
-    // make this optional
-    // DumpModelFits (psfTry->modelPSF, "modelsPSF.dat");
+    // measure the chi-square trend as a function of flux (PAR[1])
+    // this should be linear for Poisson errors and quadratic for constant sky errors
+    psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    psVector *mask  = psVectorAlloc (psfTry->sources->n, PS_TYPE_MASK);
+
+    // write sources with models first
+    for (int i = 0; i < psfTry->sources->n; i++) {
+        pmModel *model = psfTry->modelPSF->data[i];
+        if (model == NULL) {
+            flux->data.F64[i] = 0.0;
+            chisq->data.F64[i] = 0.0;
+            mask->data.U8[i] = 1;
+        } else {
+            flux->data.F64[i] = model->params->data.F32[1];
+            chisq->data.F64[i] = model[0].chisq/model[0].nDOF;
+            mask->data.U8[i] = 0;
+        }
+    }
+    // use 3hi/3lo sigma clipping on the r2rflux vs metric fit
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->clipSigma = 3.0;
+    stats->clipIter = 3;
+    // linear clipped fit of ApResid to r2rflux
+    psVectorClipFitPolynomial1D (psfTry->psf->ChiTrend, stats, mask, 1, chisq, NULL, flux);
+    for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) {
+        psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i, psfTry->psf->ChiTrend->coeff[i]*pow(10000, i), psfTry->psf->ChiTrend->coeffErr[i]);
+    }
+    psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+    psFree (stats);
+    psFree (flux);
+    psFree (chisq);
 
     // XXX this function wants aperture radius for pmSourcePhotometry
-    pmPSFtryMetric_Alt (psfTry, RADIUS);
-    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
-              modelName,
-              psfTry->psf->ApResid,
-              psfTry->psf->dApResid,
-              psfTry->psf->skyBias);
+    pmPSFtryMetric (psfTry, RADIUS);
+    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
+              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
 
     return (psfTry);
 }
 
-
 bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
 {
-
-    float dBin;
-    int   nKeep, nSkip;
-
     // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
@@ -208 +236 @@
     //   we use an outlier rejection to avoid this bias
 
-    // rflux = ten(0.4*fitMag);
-    psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
-    rflux->n = rflux->nalloc;
+    // r2rflux = radius^2 * ten(0.4*fitMag);
+    psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     for (int i = 0; i < psfTry->sources->n; i++) {
         if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
             continue;
-        rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
-    }
-
-    // find min and max of (1/flux):
-    psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-    psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL);
-
-    // build binned versions of rflux, metric
-    dBin = (stats->max - stats->min) / 10.0;
-    psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64);
-    psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
-    psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
-    daBin->n = daBin->nalloc;
-    maskB->n = maskB->nalloc;
-    psFree (stats);
-
-    psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
-
-    // group data in daBin bins, measure lower 50% mean
-    for (int i = 0; i < daBin->n; i++) {
-
-        psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
-        tmp->n = 0;
-
-        // accumulate data within bin range
-        for (int j = 0; j < psfTry->sources->n; j++) {
-            // masked for: bad model fit, outlier in parameters
-            if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL)
-                continue;
-
-            // skip points with extreme dA values
-            if (fabs(psfTry->metric->data.F64[j]) > 0.5)
-                continue;
-
-            // skip points outside of this bin
-            if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
-                continue;
-            if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
-                continue;
-
-            tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j];
-            tmp->n ++;
-        }
-
-        // is this a valid point?
-        maskB->data.U8[i] = 0;
-        if (tmp->n < 2) {
-            maskB->data.U8[i] = 1;
-            psFree (tmp);
-            continue;
-        }
-
-        // dA values are contaminated with low outliers
-        // measure statistics only on upper 50% of points
-        // this would be easier if we could sort in reverse:
-        //
-        // psVectorSort (tmp, tmp);
-        // tmp->n = 0.5*tmp->n;
-        // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-        // psVectorStats (stats, tmp, NULL, NULL, 0);
-        // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
-
-        psVectorSort (tmp, tmp);
-        nKeep = 0.5*tmp->n;
-        nSkip = tmp->n - nKeep;
-
-        psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
-        tmp2->n = tmp2->nalloc;
-        for (int j = 0; j < tmp2->n; j++) {
-            tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
-        }
-
-        stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-        psVectorStats (stats, tmp2, NULL, NULL, 0);
-        psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
-
-        daBin->data.F64[i] = stats->sampleMedian;
-
-        psFree (stats);
-        psFree (tmp);
-        psFree (tmp2);
-    }
-
-    // linear fit to rfBin, daBin
-    psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
-    psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-    poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
-
-    // XXX EAM : this is the intended API (cycle 7? cycle 8?)
-    // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
-
-    // XXX EAM : replace this when the above version is implemented
-    // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL);
-
-    psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
-    psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
-
-    stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
-    stats = psVectorStats (stats, daResid, NULL, maskB, 1);
-
-    psfTry->psf->ApResid = poly->coeff[0];
-    psfTry->psf->dApResid = stats->clippedStdev;
-    psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
-
-    psFree (rflux);
-    psFree (rfBin);
-    psFree (daBin);
-    psFree (maskB);
-    psFree (daBinFit);
-    psFree (daResid);
-    psFree (poly);
-    psFree (stats);
-    psFree (fitstat);
-
-    return true;
-}
-
-bool pmPSFtryMetric_Alt (pmPSFtry *try
-                         , float RADIUS)
-{
-
-    // the measured (aperture - fit) magnitudes (dA == try->metric)
-    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
-    //     dA = dAo + dsky/flux
-    //   where flux is the flux of the star
-    // we fit this trend to find the infinite flux aperture correction (dAo),
-    //   the nominal sky bias (dsky), and the error on dAo
-    // the values of dA are contaminated by stars with close neighbors in the aperture
-    //   we use an outlier rejection to avoid this bias
-
-    // rflux = ten(0.4*fitMag);
-    psVector *rflux = psVectorAlloc (try
-                                     ->sources->n, PS_TYPE_F64);
-    rflux->n = rflux->nalloc;
-    for (int i = 0; i < try
-                ->sources->n; i++) {
-            if (try
-                    ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
-            rflux->data.F64[i] = pow(10.0, 0.4*try
-                                     ->fitMag->data.F64[i]);
-        }
-
-    // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
+        r2rflux->data.F64[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
+    }
+
+    // use 3hi/1lo sigma clipping on the r2rflux vs metric fit
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-
-    // XXX EAM
     stats->min = 1.0;
     stats->max = 3.0;
     stats->clipIter = 3;
 
-    // linear clipped fit to rfBin, daBin
-    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
-    poly = psVectorClipFitPolynomial1D (poly, stats, try
-                                            ->mask, PSFTRY_MASK_ALL, try
-                                                ->metric, NULL, rflux);
-    fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
-
-    try
-        ->psf->ApResid = poly->coeff[0];
-    try
-        ->psf->dApResid = stats->sampleStdev;
-    try
-        ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
-
-    fprintf (stderr, "*******************************************************************************\n");
-
-    FILE *f;
-    f = fopen ("apresid.dat", "w");
-    if (f == NULL)
-        psAbort ("pmPSFtry", "can't open output file");
-
-    for (int i = 0; i < try
-                ->sources->n; i++) {
-            fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try
-                         ->fitMag->data.F64[i], rflux->data.F64[i], try
-                             ->metric->data.F64[i], try
-                                 ->mask->data.U8[i]);
-        }
-    fclose (f);
-
-    psFree (rflux);
+    // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
+    // linear clipped fit of ApResid to r2rflux
+    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
+    poly = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, r2rflux);
+    psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+
+    pmPSF_MaskApTrend (psfTry->psf, PM_PSF_SKYBIAS);
+    psfTry->psf->ApTrend->coeff[0][0][0][0] = poly->coeff[0];
+    psfTry->psf->ApTrend->coeff[0][0][1][0] = 0;
+
+    psfTry->psf->skyBias  = poly->coeff[1];
+    psfTry->psf->ApResid  = poly->coeff[0];
+    psfTry->psf->dApResid = stats->sampleStdev;
+
+    psFree (r2rflux);
     psFree (poly);
     psFree (stats);
 
-    // psFree (daFit);
-    // psFree (daResid);
-
     return true;
 }
+
+/*
+  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
+  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
+  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
+*/
+

trunk/psModules/src/objects/pmPSFtry.h

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii