Index: /trunk/psModules/src/astrom/pmAstrometry.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometry.c (revision 6872) +++ /trunk/psModules/src/astrom/pmAstrometry.c (revision 6873) @@ -13,6 +13,6 @@ * 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 Index: /trunk/psModules/src/astrom/pmAstrometry.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometry.h (revision 6872) +++ /trunk/psModules/src/astrom/pmAstrometry.h (revision 6873) @@ -8,6 +8,6 @@ * @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 Index: /trunk/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 6872) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 6873) @@ -8,6 +8,6 @@ * @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 Index: /trunk/psModules/src/config/pmConfig.c =================================================================== --- /trunk/psModules/src/config/pmConfig.c (revision 6872) +++ /trunk/psModules/src/config/pmConfig.c (revision 6873) @@ -3,6 +3,6 @@ * @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,9 +11,49 @@ #include #include +#include +#include +#include +#include +#include #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,24 +63,74 @@ * */ -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,23 +147,13 @@ 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,5 +177,5 @@ "-site command-line switch provided without the required filename --- ignored.\n"); } else { - siteName = argv[argNum]; + siteName = psStringCopy(argv[argNum]); psArgumentRemove(argNum, argc, argv); } @@ -106,4 +186,7 @@ if (!siteName) { siteName = getenv(PS_SITE); + if (siteName) { + siteName = psStringCopy (siteName); + } } @@ -111,46 +194,33 @@ // 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,9 +234,11 @@ } 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,5 +255,5 @@ // 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,5 +267,5 @@ // 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,5 +278,5 @@ // 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,16 +290,22 @@ // 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,5 +314,5 @@ // 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,14 +342,16 @@ 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,10 +362,10 @@ 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,5 +377,6 @@ if (strncmp(ruleItem->data.V, headerItem->data.V, strlen(ruleItem->data.V)) != 0) { - match = false; + psFree(ruleIter); + return false; } break; @@ -305,5 +386,6 @@ ruleItem->data.S32, headerItem->data.S32); if (ruleItem->data.S32 != headerItem->data.S32) { - match = false; + psFree(ruleIter); + return false; } break; @@ -312,5 +394,6 @@ ruleItem->data.F32, headerItem->data.F32); if (ruleItem->data.F32 != headerItem->data.F32) { - match = false; + psFree(ruleIter); + return false; } break; @@ -319,5 +402,6 @@ ruleItem->data.F64, headerItem->data.F64); if (ruleItem->data.F64 != headerItem->data.F64) { - match = false; + psFree(ruleIter); + return false; } break; @@ -329,93 +413,196 @@ 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,15 +610,9 @@ 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,27 +623,111 @@ 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; +} Index: /trunk/psModules/src/config/pmConfig.h =================================================================== --- /trunk/psModules/src/config/pmConfig.h (revision 6872) +++ /trunk/psModules/src/config/pmConfig.h (revision 6873) @@ -3,6 +3,6 @@ * @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 Index: /trunk/psModules/src/detrend/pmFlatField.c =================================================================== --- /trunk/psModules/src/detrend/pmFlatField.c (revision 6872) +++ /trunk/psModules/src/detrend/pmFlatField.c (revision 6873) @@ -24,6 +24,6 @@ * @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 Index: /trunk/psModules/src/detrend/pmMaskBadPixels.c =================================================================== --- /trunk/psModules/src/detrend/pmMaskBadPixels.c (revision 6872) +++ /trunk/psModules/src/detrend/pmMaskBadPixels.c (revision 6873) @@ -24,6 +24,6 @@ * @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 Index: /trunk/psModules/src/detrend/pmMaskBadPixelsErrors.h =================================================================== --- /trunk/psModules/src/detrend/pmMaskBadPixelsErrors.h (revision 6872) +++ /trunk/psModules/src/detrend/pmMaskBadPixelsErrors.h (revision 6873) @@ -12,6 +12,6 @@ * @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 Index: /trunk/psModules/src/detrend/pmNonLinear.c =================================================================== --- /trunk/psModules/src/detrend/pmNonLinear.c (revision 6872) +++ /trunk/psModules/src/detrend/pmNonLinear.c (revision 6873) @@ -10,6 +10,6 @@ * @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 Index: /trunk/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- /trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 6872) +++ /trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 6873) @@ -5,6 +5,6 @@ * @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,15 +36,287 @@ } +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,6 +690,6 @@ 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,6 +726,6 @@ } - 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,11 +736,9 @@ 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); } - - /** Index: /trunk/psModules/src/imcombine/pmReadoutCombine.h =================================================================== --- /trunk/psModules/src/imcombine/pmReadoutCombine.h (revision 6872) +++ /trunk/psModules/src/imcombine/pmReadoutCombine.h (revision 6873) @@ -5,6 +5,6 @@ * @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 Index: /trunk/psModules/src/imsubtract/pmSubtractBias.c =================================================================== --- /trunk/psModules/src/imsubtract/pmSubtractBias.c (revision 6872) +++ /trunk/psModules/src/imsubtract/pmSubtractBias.c (revision 6873) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// 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,28 +11,28 @@ * @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 -#include -#include -#include "pslib.h" + #if HAVE_CONFIG_H #include #endif + +#include #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,105 +62,98 @@ }\ -/*****************************************************************************/ -/* 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;iimage->numRows;i++) { - for (psS32 j=0;jimage->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;iimage->numRows;i++) { - for (psS32 j=0;jimage->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;iimage->numRows;i++) { - for (psS32 j=0;jimage->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;inumRows;i++) { @@ -166,4 +164,5 @@ return(image); } +#endif /****************************************************************************** @@ -179,12 +178,4 @@ 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,11 +185,4 @@ } - 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,5 +206,5 @@ 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,97 +215,7 @@ } -/****************************************************************************** -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,14 +224,16 @@ 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,21 +243,25 @@ // 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,9 +270,17 @@ // 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;in) / ((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,26 +288,26 @@ // 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;in - 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,882 +321,257 @@ } +#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;inumRows;i++) { - for (psS32 j=0;jnumCols;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;inumCols;i++) { - for (psS32 j=0;jnumRows;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;in > ((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;idata.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;inumCols;i++) { - for (psS32 j=0;jnumRows;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;inumRows;i++) { - for (psS32 j=0;jnumCols;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; +} + + Index: /trunk/psModules/src/imsubtract/pmSubtractBias.h =================================================================== --- /trunk/psModules/src/imsubtract/pmSubtractBias.h (revision 6872) +++ /trunk/psModules/src/imsubtract/pmSubtractBias.h (revision 6873) @@ -11,6 +11,6 @@ * @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 Index: /trunk/psModules/src/imsubtract/pmSubtractSky.h =================================================================== --- /trunk/psModules/src/imsubtract/pmSubtractSky.h (revision 6872) +++ /trunk/psModules/src/imsubtract/pmSubtractSky.h (revision 6873) @@ -6,6 +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 Index: /trunk/psModules/src/objects/pmModelGroup.c =================================================================== --- /trunk/psModules/src/objects/pmModelGroup.c (revision 6872) +++ /trunk/psModules/src/objects/pmModelGroup.c (revision 6873) @@ -6,6 +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 Index: /trunk/psModules/src/objects/pmModelGroup.h =================================================================== --- /trunk/psModules/src/objects/pmModelGroup.h (revision 6872) +++ /trunk/psModules/src/objects/pmModelGroup.h (revision 6873) @@ -9,6 +9,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 Index: /trunk/psModules/src/objects/pmObjects.c =================================================================== --- /trunk/psModules/src/objects/pmObjects.c (revision 6872) +++ /trunk/psModules/src/objects/pmObjects.c (revision 6873) @@ -6,6 +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,2012 +19,3 @@ #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, ¶ms_min, ¶ms_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); -} - Index: /trunk/psModules/src/objects/pmPSF.c =================================================================== --- /trunk/psModules/src/objects/pmPSF.c (revision 6872) +++ /trunk/psModules/src/objects/pmPSF.c (revision 6873) @@ -6,6 +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,5 +18,11 @@ #include -#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,4 +63,6 @@ return; + psFree (psf->ApTrend); + psFree (psf->growth); psFree (psf->params); return; @@ -69,8 +77,8 @@ 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,4 +91,19 @@ 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,5 +117,5 @@ 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,11 +197,11 @@ /***************************************************************************** -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,6 +212,104 @@ // 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; +} Index: /trunk/psModules/src/objects/pmPSF.h =================================================================== --- /trunk/psModules/src/objects/pmPSF.h (revision 6872) +++ /trunk/psModules/src/objects/pmPSF.h (revision 6873) @@ -6,6 +6,4 @@ * @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 Index: /trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtry.c (revision 6872) +++ /trunk/psModules/src/objects/pmPSFtry.c (revision 6873) @@ -5,6 +5,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 @@ -13,8 +13,16 @@ # include -# 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,5 +41,5 @@ psFree (test->psf); psFree (test->sources); - psFree (test->modelFLT); + psFree (test->modelEXT); psFree (test->modelPSF); psFree (test->metric); @@ -42,5 +50,5 @@ // 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,7 +59,7 @@ // 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,13 +67,7 @@ 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,5 +88,5 @@ // 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,8 +95,8 @@ 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,26 +110,24 @@ // 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,13 +139,13 @@ 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,5 +162,9 @@ // 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,31 +176,55 @@ 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,197 +236,43 @@ // 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) +*/ + Index: /trunk/psModules/src/objects/pmPSFtry.h =================================================================== --- /trunk/psModules/src/objects/pmPSFtry.h (revision 6872) +++ /trunk/psModules/src/objects/pmPSFtry.h (revision 6873) @@ -6,6 +6,4 @@ * @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