Index: /branches/eam_branch_20071023/psModules/src/camera/pmFPAConstruct.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/camera/pmFPAConstruct.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/camera/pmFPAConstruct.c (revision 15516) @@ -45,5 +45,5 @@ psMetadata *cellData = psMetadataLookupMetadata(&status, cells, cellName); // The data for the particular cell if (!status || !cellData) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to find specs for cell %s: ignored\n", cellName); + psWarning("Unable to find specs for cell %s: ignored\n", cellName); } @@ -304,5 +304,5 @@ // Put in the cell data if (newCell->config) { - psLogMsg(__func__, PS_LOG_WARN, "Overwriting cell data in chip\n"); + psWarning("Overwriting cell data in chip\n"); psFree(newCell->config); // Make way! } @@ -1304,5 +1304,5 @@ const char *chipName = componentsItem->name; // Name of the chip if (componentsItem->type != PS_DATA_STRING) { - psLogMsg(__func__, PS_LOG_WARN, "Element %s in FPA within the camera configuration is not of " + psWarning("Element %s in FPA within the camera configuration is not of " "type STR (type=%x) --- ignored.\n", chipName, componentsItem->type); continue; Index: /branches/eam_branch_20071023/psModules/src/camera/pmFPAFlags.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/camera/pmFPAFlags.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/camera/pmFPAFlags.c (revision 15516) @@ -136,7 +136,7 @@ for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (chip == NULL) continue; - if (chip->data_exists) return true; + pmChip *chip = fpa->chips->data[i]; + if (chip == NULL) continue; + if (chip->data_exists) return true; } return false; @@ -170,33 +170,33 @@ if (view->chip == -1) { - bool exists = pmFPACheckDataStatus (fpa); - return exists; + bool exists = pmFPACheckDataStatus (fpa); + return exists; } if (view->chip >= fpa->chips->n) { - psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n); - return false; + psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n); + return false; } pmChip *chip = fpa->chips->data[view->chip]; if (view->cell == -1) { - bool exists = pmChipCheckDataStatus (chip); - return exists; + bool exists = pmChipCheckDataStatus (chip); + return exists; } if (view->cell >= chip->cells->n) { - psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n); - return false; + psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n); + return false; } pmCell *cell = chip->cells->data[view->cell]; if (view->readout == -1) { - bool exists = pmCellCheckDataStatus (cell); - return exists; + bool exists = pmCellCheckDataStatus (cell); + return exists; } if (view->readout >= cell->readouts->n) { - psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouds->n == %ld", view->readout, cell->readouts->n); - return false; + psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouds->n == %ld", view->readout, cell->readouts->n); + return false; } pmReadout *readout = cell->readouts->data[view->readout]; @@ -290,9 +290,9 @@ psArray *chips = fpa->chips; // Component chips if (chips == NULL) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: fpa->chips == NULL\n"); + psWarning("WARNING: fpa->chips == NULL\n"); return(0); } if ((chipNum >= chips->n) || (NULL == (pmChip *) chips->data[chipNum])) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: the specified chip (%d) does not exist.\n", chipNum); + psWarning("WARNING: the specified chip (%d) does not exist.\n", chipNum); return(0); } Index: /branches/eam_branch_20071023/psModules/src/camera/pmFPAMosaic.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/camera/pmFPAMosaic.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/camera/pmFPAMosaic.c (revision 15516) @@ -704,6 +704,6 @@ psArray *readouts = cell->readouts; // The array of readouts if (readouts->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Cell contains more than one readout (%ld) --- only the first will " - "be mosaicked.\n", readouts->n); + psWarning("Cell contains more than one readout (%ld) --- only the first will be mosaicked.\n", + readouts->n); } pmReadout *readout = readouts->data[0]; // The only readout we'll bother with @@ -754,15 +754,15 @@ int x0Target = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.X0"); if (!mdok) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.X0 is not set for the target cell; assuming 0.\n"); + psWarning("CELL.X0 is not set for the target cell; assuming 0.\n"); FIX_CONCEPT(targetCell->concepts, "CELL.X0", S32, 0); } int y0Target = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.Y0"); if (!mdok) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.Y0 is not set for the target cell; assuming 0.\n"); + psWarning("CELL.Y0 is not set for the target cell; assuming 0.\n"); FIX_CONCEPT(targetCell->concepts, "CELL.Y0", S32, 0); } int xParityTarget = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.XPARITY"); if (!mdok || (xParityTarget != -1 && xParityTarget != 1)) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.XPARITY is not set for the target cell; assuming 1.\n"); + psWarning("CELL.XPARITY is not set for the target cell; assuming 1.\n"); FIX_CONCEPT(targetCell->concepts, "CELL.XPARITY", S32, 1); xParityTarget = 1; @@ -770,5 +770,5 @@ int yParityTarget = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.YPARITY"); if (!mdok || (yParityTarget != -1 && yParityTarget != 1)) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.YPARITY is not set for the target cell; assuming 1.\n"); + psWarning("CELL.YPARITY is not set for the target cell; assuming 1.\n"); FIX_CONCEPT(targetCell->concepts, "CELL.YPARITY", S32, 1); yParityTarget = 1; @@ -863,25 +863,25 @@ int x0Target = psMetadataLookupS32(&mdok, targetChip->concepts, "CHIP.X0"); if (!mdok) { - psLogMsg(__func__, PS_LOG_WARN, "CHIP.X0 is not set for the target chip; assuming 0.\n"); + psWarning("CHIP.X0 is not set for the target chip; assuming 0.\n"); FIX_CONCEPT(targetChip->concepts, "CHIP.X0", S32, 0); } int y0Target = psMetadataLookupS32(&mdok, targetChip->concepts, "CHIP.Y0"); if (!mdok) { - psLogMsg(__func__, PS_LOG_WARN, "CHIP.Y0 is not set for the target chip; assuming 0.\n"); + psWarning("CHIP.Y0 is not set for the target chip; assuming 0.\n"); FIX_CONCEPT(targetChip->concepts, "CHIP.Y0", S32, 0); } x0Target += psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.X0"); if (!mdok) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.X0 is not set for the target cell; assuming 0.\n"); + psWarning("CELL.X0 is not set for the target cell; assuming 0.\n"); FIX_CONCEPT(targetCell->concepts, "CELL.X0", S32, 0); } y0Target += psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.Y0"); if (!mdok) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.Y0 is not set for the target cell; assuming 0.\n"); + psWarning("CELL.Y0 is not set for the target cell; assuming 0.\n"); FIX_CONCEPT(targetCell->concepts, "CELL.Y0", S32, 0); } int xParityChipTarget = psMetadataLookupS32(&mdok, targetChip->concepts, "CHIP.XPARITY"); if (!mdok || (xParityChipTarget != -1 && xParityChipTarget != 1)) { - psLogMsg(__func__, PS_LOG_WARN, "CHIP.XPARITY is not set for the target chip; assuming 1.\n"); + psWarning("CHIP.XPARITY is not set for the target chip; assuming 1.\n"); FIX_CONCEPT(targetChip->concepts, "CHIP.XPARITY", S32, 1); xParityChipTarget = 1; @@ -889,5 +889,5 @@ int yParityChipTarget = psMetadataLookupS32(&mdok, targetChip->concepts, "CHIP.YPARITY"); if (!mdok || (yParityChipTarget != -1 && yParityChipTarget != 1)) { - psLogMsg(__func__, PS_LOG_WARN, "CHIP.YPARITY is not set for the target chip; assuming 1.\n"); + psWarning("CHIP.YPARITY is not set for the target chip; assuming 1.\n"); FIX_CONCEPT(targetChip->concepts, "CHIP.YPARITY", S32, 1); yParityChipTarget = 1; @@ -895,5 +895,5 @@ int xParityCellTarget = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.XPARITY"); if (!mdok || (xParityCellTarget != -1 && xParityCellTarget != 1)) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.XPARITY is not set for the target cell; assuming 1.\n"); + psWarning("CELL.XPARITY is not set for the target cell; assuming 1.\n"); FIX_CONCEPT(targetCell->concepts, "CELL.XPARITY", S32, 1); xParityCellTarget = 1; @@ -901,5 +901,5 @@ int yParityCellTarget = psMetadataLookupS32(&mdok, targetCell->concepts, "CELL.YPARITY"); if (!mdok || (yParityCellTarget != -1 && yParityCellTarget != 1)) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.YPARITY is not set for the target cell; assuming 1.\n"); + psWarning("CELL.YPARITY is not set for the target cell; assuming 1.\n"); FIX_CONCEPT(targetCell->concepts, "CELL.YPARITY", S32, 1); yParityCellTarget = 1; @@ -1334,10 +1334,10 @@ pmHDU *sourceHDU = pmHDUGetHighest(source, firstSourceChip, firstSourceCell); // The HDU for the source if (!sourceHDU) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to find HDU in source FPA; unable to copy headers.\n"); + psWarning("Unable to find HDU in source FPA; unable to copy headers.\n"); return false; } pmHDU *targetHDU = pmHDUGetHighest(target, targetChip, targetCell); // The HDU for the target if (!targetHDU) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to find HDU in target FPA; unable to copy headers.\n"); + psWarning("Unable to find HDU in target FPA; unable to copy headers.\n"); return false; } Index: /branches/eam_branch_20071023/psModules/src/camera/pmHDU.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/camera/pmHDU.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/camera/pmHDU.c (revision 15516) @@ -119,5 +119,5 @@ if (*images) { - psLogMsg(__func__, PS_LOG_WARN, "HDU %s has already been read --- overwriting.\n", hdu->extname); + psWarning("HDU %s has already been read --- overwriting.\n", hdu->extname); psFree(*images); // Blow away anything existing } @@ -167,5 +167,5 @@ if (!images && !hdu->header) { - psLogMsg(__func__, PS_LOG_WARN, "Nothing to write for HDU %s\n", hdu->extname); + psWarning("Nothing to write for HDU %s\n", hdu->extname); return false; } Index: /branches/eam_branch_20071023/psModules/src/camera/pmHDUGenerate.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/camera/pmHDUGenerate.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/camera/pmHDUGenerate.c (revision 15516) @@ -134,6 +134,5 @@ psMetadataItem *biassecItem = psMetadataLookup(cell->concepts, "CELL.BIASSEC"); // Bias sections if (!biassecItem) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.BIASSEC has not been initialised in cell --- " - "ignored.\n"); + psWarning("CELL.BIASSEC has not been initialised in cell --- ignored.\n"); return false; } @@ -146,5 +145,5 @@ if (!mdok || (cellreaddir != 1 && cellreaddir != 2)) { // Probably unnecessary, but just in case.... - psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR is not set in cell --- ignored.\n"); + psWarning("CELL.READDIR is not set in cell --- ignored.\n"); return false; } @@ -152,5 +151,5 @@ *readdir = cellreaddir; } else if (*readdir != cellreaddir) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR does not match read direction for HDU --- ignored.\n"); + psWarning("CELL.READDIR does not match read direction for HDU --- ignored.\n"); return false; } @@ -255,6 +254,5 @@ psMetadataItem *trimsecItem = psMetadataLookup(cell->concepts, "CELL.TRIMSEC"); // Item with trimsec if (!trimsecItem || trimsecItem->type != PS_DATA_REGION) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.TRIMSEC has not been initialised in cell --- " - "ignored.\n"); + psWarning("CELL.TRIMSEC has not been initialised in cell --- ignored.\n"); continue; } @@ -264,5 +262,5 @@ if (!mdok || (cellreaddir != 1 && cellreaddir != 2)) { // Probably unnecessary, but just in case.... - psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR is not set in cell --- ignored.\n"); + psWarning("CELL.READDIR is not set in cell --- ignored.\n"); continue; } @@ -270,5 +268,5 @@ readdir = cellreaddir; } else if (readdir != cellreaddir) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.READDIR for cells within the HDU do not match!\n"); + psWarning("CELL.READDIR for cells within the HDU do not match!\n"); cellreaddir = readdir; } @@ -282,10 +280,10 @@ if (readout->mask && (readout->mask->numCols != image->numCols || readout->mask->numRows != image->numRows)) { - psLogMsg(__func__, PS_LOG_WARN, "Image and mask have different sizes (%dx%d vs %dx%d)!\n", + psWarning("Image and mask have different sizes (%dx%d vs %dx%d)!\n", image->numCols, image->numRows, readout->mask->numCols, readout->mask->numRows); } if (readout->weight && (readout->weight->numCols != image->numCols || readout->weight->numRows != image->numRows)) { - psLogMsg(__func__, PS_LOG_WARN, "Image and weight have different sizes (%dx%d vs %dx%d)!\n", + psWarning("Image and weight have different sizes (%dx%d vs %dx%d)!\n", image->numCols, image->numRows, readout->weight->numCols, readout->weight->numRows); } @@ -318,6 +316,6 @@ if (previous != current) { - psLogMsg(__func__, PS_LOG_WARN, "Images within the HDU are of different types " - "(%x vs %x) --- promoting\n", previous, current); + psWarning("Images within the HDU are of different types (%x vs %x) --- promoting\n", + previous, current); return PS_MAX(previous, current); } @@ -336,5 +334,5 @@ if (source->numCols != region->x1 - region->x0 || source->numRows != region->y1 - region->y0) { psString regionString = psRegionToString(*region); - psLogMsg(__func__, PS_LOG_WARN, "Image size (%dx%d) does not match region (%s).\n", + psWarning("Image size (%dx%d) does not match region (%s).\n", source->numCols, source->numRows, regionString); psFree(regionString); @@ -475,6 +473,6 @@ if (biassecs->n != readout->bias->n) { - psLogMsg(__func__, PS_LOG_WARN, "Number of bias sections (%ld) and number of biases (%ld)" - " do not match.\n", biassecs->n, readout->bias->n); + psWarning("Number of bias sections (%ld) and number of biases (%ld) do not match.\n", + biassecs->n, readout->bias->n); } psListIterator *biasIter = psListIteratorAlloc(readout->bias, PS_LIST_HEAD, false); // Iteratr Index: /branches/eam_branch_20071023/psModules/src/camera/pmReadoutFake.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/camera/pmReadoutFake.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/camera/pmReadoutFake.c (revision 15516) @@ -35,4 +35,5 @@ pmReadout *readout = pmReadoutAlloc(NULL); // Output readout readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->image, 0); int numSources = sources->n; // Number of stars @@ -86,4 +87,8 @@ fakeModel->params->data.F32[PM_PAR_I0] = powf(10.0, -0.4 * source->psfMag) / flux0; + psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n", + fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS], + fakeModel->params->data.F32[PM_PAR_I0]); + pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate fakeSource->peak = pmPeakAlloc(x, y, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE); Index: /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsAverage.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsAverage.h (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsAverage.h (revision 15516) @@ -4,6 +4,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.9.10.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-29 01:40:54 $ + * @version $Revision: 1.9.10.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-08 22:59:04 $ * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ Index: /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsStandard.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsStandard.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsStandard.c (revision 15516) @@ -375,5 +375,5 @@ assert(pattern); - psList *biassecs; // List of bias sections + psList *biassecs = NULL; // List of bias sections switch (concept->type) { Index: /branches/eam_branch_20071023/psModules/src/config/pmConfig.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/config/pmConfig.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/config/pmConfig.c (revision 15516) @@ -661,9 +661,15 @@ if (!cameraReadFormats(item->data.md, item->name)) { psWarning("Unable to read formats for camera %s: removed.\n", item->name); + psErrorStackPrint(stderr, "errors from read failure\n"); + psErrorClear(); psMetadataRemoveKey(cameras, item->name); + continue; } if (!cameraReadCalibrations(item->data.md, item->name)) { psWarning("Unable to read calibrations for camera %s: removed.\n", item->name); + psErrorStackPrint(stderr, "errors from read failure\n"); + psErrorClear(); psMetadataRemoveKey(cameras, item->name); + continue; } } Index: /branches/eam_branch_20071023/psModules/src/config/pmConfig.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/config/pmConfig.h (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/config/pmConfig.h (revision 15516) @@ -5,6 +5,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.31.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-29 01:40:54 $ + * @version $Revision: 1.31.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-08 22:59:04 $ * * Copyright 2005-2006 Institute for Astronomy, University of Hawaii Index: /branches/eam_branch_20071023/psModules/src/config/pmConfigCamera.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/config/pmConfigCamera.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/config/pmConfigCamera.c (revision 15516) @@ -464,5 +464,5 @@ while ((formatsItem = psMetadataGetAndIncrement(formatsIter))) { assert(formatsItem->type == PS_DATA_METADATA); // We should have read it by now! - psMetadata *format = formatsItem->data.V; // The camera format + psMetadata *format = formatsItem->data.md; // The camera format // Add a new RULE which uniquely describes the mosaicked format. this is needed so Index: /branches/eam_branch_20071023/psModules/src/config/pmConfigRecipes.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/config/pmConfigRecipes.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/config/pmConfigRecipes.c (revision 15516) @@ -5,5 +5,5 @@ #include #include -#include /* for strn?casecmp */ +#include /* for strn?casecmp */ #include #include @@ -33,6 +33,6 @@ options = psMetadataAlloc (); success = psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "OPTIONS", PS_DATA_METADATA, "", options); - assert (success); // type mismatch : OPTIONS already defined but wrong type - psFree (options); // drop extra reference + assert (success); // type mismatch : OPTIONS already defined but wrong type + psFree (options); // drop extra reference } @@ -44,6 +44,6 @@ recipe = psMetadataAlloc(); success = psMetadataAddPtr(options, PS_LIST_TAIL, recipeName, PS_DATA_METADATA, "", recipe); - assert (success); // type mismatch : OPTIONS already defined but wrong type - psFree (recipe); // drop extra reference + assert (success); // type mismatch : OPTIONS already defined but wrong type + psFree (recipe); // drop extra reference } return recipe; @@ -67,9 +67,9 @@ // master recipe files in the site:recipe location when they are built. if (config->site && (source & PM_RECIPE_SOURCE_SITE)) { - if (!loadRecipeSite(&status, config, config->site)) { + if (!loadRecipeSite(&status, config, config->site)) { psError(PS_ERR_IO, false, "Failed to read recipes from site config"); - return false; - } - psTrace ("psModules.config", 3, "read recipes from site config"); + return false; + } + psTrace ("psModules.config", 3, "read recipes from site config"); } @@ -81,11 +81,11 @@ if (!loadRecipeCamera(&status, config, config->camera)) { psError(PS_ERR_IO, false, "Failed to read recipes from camera config"); - return false; - } - if (status) { - psTrace ("psModules.config", 3, "read recipes from camera config"); - } else { - psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by camera config"); - } + return false; + } + if (status) { + psTrace ("psModules.config", 3, "read recipes from camera config"); + } else { + psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by camera config"); + } } @@ -94,25 +94,25 @@ if (!loadRecipeFromArguments(&status, config)) { psError(PS_ERR_IO, false, "Failed to read recipes from command-line arguments"); - return false; - } - if (status) { - psTrace ("psModules.config", 3, "read recipes from command-line arguments"); - } else { - psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied on command-line arguments"); - } + return false; + } + if (status) { + psTrace ("psModules.config", 3, "read recipes from command-line arguments"); + } else { + psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied on command-line arguments"); + } if (!loadRecipeSymbols(&status, config)) { psError(PS_ERR_IO, false, "Failed to read recipes from symbolic references"); - return false; - } - if (status) { - psTrace ("psModules.config", 3, "read recipes from symbolic references"); - } else { - psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by symbolic reference"); - } + return false; + } + if (status) { + psTrace ("psModules.config", 3, "read recipes from symbolic references"); + } else { + psLogMsg ("psModules.config", PS_LOG_DETAIL, "no recipe supplied by symbolic reference"); + } if (!loadRecipeOptions(&status, config)) { psError(PS_ERR_IO, false, "Failed to read recipes from symbolic references"); - return false; - } - if (status) { + return false; + } + if (status) { psTrace ("psModules.config", 3, "read recipes from command-line arguments"); } else { @@ -135,6 +135,6 @@ options = psMetadataAlloc(); success = psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "OPTIONS", PS_DATA_METADATA, "Command-line options specified with -D", options); - assert (success); // type mismatch : OPTIONS already defined but wrong type - psFree (options); // drop extra reference + assert (success); // type mismatch : OPTIONS already defined but wrong type + psFree (options); // drop extra reference } @@ -176,6 +176,6 @@ recipe = psMetadataAlloc(); success = psMetadataAddPtr(options, PS_LIST_TAIL, recipeName, PS_DATA_METADATA, "", recipe); - assert (success); // type mismatch : recipe already defined but wrong type - psFree (recipe); // drop extra reference + assert (success); // type mismatch : recipe already defined but wrong type + psFree (recipe); // drop extra reference } @@ -216,5 +216,5 @@ bool pmConfigLoadRecipeArguments (int *argc, char **argv, pmConfig *config) { - bool success; + bool success; psMetadata *recipes = psMetadataLookupMetadata(&success, config->arguments, "RECIPES"); @@ -222,6 +222,6 @@ recipes = psMetadataAlloc(); success = psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "RECIPES", PS_DATA_METADATA, "", recipes); - assert (success); - psFree (recipes); + assert (success); + psFree (recipes); } @@ -271,6 +271,6 @@ // Load the recipe files for SITE : REQUIRED static bool loadRecipeSite(bool *status, - pmConfig *config, // The configuration into which to read the recipes - psMetadata *source // The source configuration, from which to read the filenames + pmConfig *config, // The configuration into which to read the recipes + psMetadata *source // The source configuration, from which to read the filenames ) { @@ -300,17 +300,18 @@ if (fileItem->type != PS_DATA_STRING) { psError(PS_ERR_IO, true, "%s in site configuration RECIPES is not of type STR", fileItem->name); - return false; + return false; } // Read the recipe file. if we fail on a file, give a warning, but continue psMetadata *recipe = NULL; - if (!pmConfigFileRead(&recipe, fileItem->data.V, "recipe")) { - psError(PS_ERR_IO, false, "Failed to read recipe file %s listed in site configuration\n", (char *)fileItem->data.V); - return false; + if (!pmConfigFileRead(&recipe, fileItem->data.str, "recipe")) { + psError(PS_ERR_IO, false, "Failed to read recipe file %s listed in site configuration\n", + fileItem->data.str); + return false; } // if this named recipe exists, supplement it - psMetadataAdd(config->recipes, PS_LIST_TAIL, fileItem->name, PS_DATA_METADATA | PS_META_REPLACE, - fileItem->comment, recipe); + psMetadataAdd(config->recipes, PS_LIST_TAIL, fileItem->name, PS_DATA_METADATA | PS_META_REPLACE, + fileItem->comment, recipe); psFree(recipe); // Drop reference } @@ -326,6 +327,6 @@ // for sourceType == SITE | CAMERA, RECIPES contains a list of files to be read (pmConfigFileRead) static bool loadRecipeCamera(bool *status, // status variable - pmConfig *config, // The configuration into which to read the recipes - psMetadata *source // The source configuration, from which to read the filenames + pmConfig *config, // The configuration into which to read the recipes + psMetadata *source // The source configuration, from which to read the filenames ) { @@ -345,5 +346,5 @@ psMetadata *recipes = psMetadataLookupMetadata(&success, source, "RECIPES"); // The list of recipes if (!recipes) { - psTrace ("psModules.config", 3, "RECIPES not found in the camera configuration\n"); + psTrace ("psModules.config", 3, "RECIPES not found in the camera configuration\n"); return true; } @@ -355,8 +356,8 @@ psMetadataItem *fileItem = NULL; // MD item containing the filename, from recipe iteration while ((fileItem = psMetadataGetAndIncrement(recipesIter))) { - char *recipeName = fileItem->name; - char *recipeFile = fileItem->data.V; - - psTrace("psModules.config", 3, "Supplementing %s from %s within camera configuration.\n", recipeName, recipeFile); + char *recipeName = fileItem->name; + char *recipeFile = fileItem->data.str; + + psTrace("psModules.config", 3, "Supplementing %s from %s within camera configuration.\n", recipeName, recipeFile); // type mismatch is a serious error @@ -368,21 +369,22 @@ psMetadata *recipe = NULL; if (!pmConfigFileRead(&recipe, recipeFile, "recipe")) { - psError(PS_ERR_IO, false, "Failed to read recipe file %s listed in camera configuration\n", recipeFile); - return false; + psError(PS_ERR_IO, false, "Failed to read recipe file %s listed in camera configuration\n", + recipeFile); + return false; } // the named recipe must exist; supplement it - psMetadata *current = psMetadataLookupMetadata(NULL, config->recipes, recipeName); - if (!current) { - psError(PS_ERR_IO, false, "Failed to find recipe for %s in master recipe list", recipeName); - psFree(recipe); // Drop reference - return false; - } - - if (!psMetadataUpdate(current, recipe)) { - psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", recipeName); - psFree(recipe); // Drop reference - return false; - } + psMetadata *current = psMetadataLookupMetadata(NULL, config->recipes, recipeName); + if (!current) { + psError(PS_ERR_IO, false, "Failed to find recipe for %s in master recipe list", recipeName); + psFree(recipe); // Drop reference + return false; + } + + if (!psMetadataUpdate(current, recipe)) { + psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", recipeName); + psFree(recipe); // Drop reference + return false; + } psFree(recipe); // Drop reference } @@ -397,5 +399,5 @@ // entries for an existing recipe metadata static bool loadRecipeFromArguments(bool *status, - pmConfig *config // The configuration into which to read the recipes + pmConfig *config // The configuration into which to read the recipes ) { @@ -425,19 +427,19 @@ // increment the ref counter to protect the data - psMetadata *recipe = item->data.V; // Recipe of interest + psMetadata *recipe = item->data.md; // Recipe of interest // if this named recipe exists, supplement it psMetadata *current = psMetadataLookupMetadata(NULL, config->recipes, item->name); if (!current) { - psError(PS_ERR_IO, false, "Failed to find recipe for %s in master recipe list", item->name); - psFree(recipe); // Drop reference - return false; - } - psTrace("psModules.config", 3, "Supplementing %s from arguments.\n", item->name); + psError(PS_ERR_IO, false, "Failed to find recipe for %s in master recipe list", item->name); + psFree(recipe); // Drop reference + return false; + } + psTrace("psModules.config", 3, "Supplementing %s from arguments.\n", item->name); if (!psMetadataUpdate (current, recipe)) { - psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", item->name); - return false; - } + psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", item->name); + return false; + } } psFree(recipesIter); @@ -449,5 +451,5 @@ // entries for an existing recipe metadata static bool loadRecipeSymbols(bool *status, - pmConfig *config // The configuration into which to read the recipes + pmConfig *config // The configuration into which to read the recipes ) { @@ -465,30 +467,30 @@ while ((item = psMetadataGetAndIncrement(iter))) { assert(item->type == PS_DATA_STRING); // It should be this type: we put it in ourselves - const char *sourceName = item->data.V; // The name of the symbolic reference + const char *sourceName = item->data.str; // The name of the symbolic reference const char *targetName = item->name; - psTrace("psModules.config", 3, "Supplementing %s from %s.\n", targetName, sourceName); - - // the target recipe must exist; select it - psMetadata *targetMD = psMetadataLookupMetadata(&found, config->recipes, targetName); - if (!targetMD) { - psError(PS_ERR_IO, true, "Failed to find recipe for %s in master recipe list", targetName); - return false; - } + psTrace("psModules.config", 3, "Supplementing %s from %s.\n", targetName, sourceName); + + // the target recipe must exist; select it + psMetadata *targetMD = psMetadataLookupMetadata(&found, config->recipes, targetName); + if (!targetMD) { + psError(PS_ERR_IO, true, "Failed to find recipe for %s in master recipe list", targetName); + return false; + } // search for sourceName : it may be in config->recipes or target MD psMetadata *sourceMD = NULL; - sourceMD = psMetadataLookupMetadata(&found, config->recipes, sourceName); + sourceMD = psMetadataLookupMetadata(&found, config->recipes, sourceName); if (!sourceMD) { - sourceMD = psMetadataLookupMetadata(&found, targetMD, sourceName); - if (!sourceMD) { - psError(PS_ERR_IO, false, "Selected symbolic name %s does not exist in recipes", sourceName); - return false; - } - } - - if (!psMetadataUpdate(targetMD, sourceMD)) { - psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", targetName); - return false; - } + sourceMD = psMetadataLookupMetadata(&found, targetMD, sourceName); + if (!sourceMD) { + psError(PS_ERR_IO, false, "Selected symbolic name %s does not exist in recipes", sourceName); + return false; + } + } + + if (!psMetadataUpdate(targetMD, sourceMD)) { + psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", targetName); + return false; + } } psFree(iter); @@ -501,5 +503,5 @@ // entries for an existing recipe metadata static bool loadRecipeOptions(bool *status, - pmConfig *config // The configuration into which to read the recipes + pmConfig *config // The configuration into which to read the recipes ) { @@ -524,6 +526,6 @@ psMetadataItem *item = NULL; // MD item containing the filename, from recipe iteration while ((item = psMetadataGetAndIncrement(recipesIter))) { - char *recipeName = item->name; - psMetadata *recipe = item->data.V; + char *recipeName = item->name; + psMetadata *recipe = item->data.md; // type mismatch is a serious error @@ -535,11 +537,11 @@ psMetadata *current = psMetadataLookupMetadata(NULL, config->recipes, recipeName); if (!current) { - psError(PS_ERR_IO, false, "Selected recipe %s is not found in camera recipe", recipeName); - return false; - } + psError(PS_ERR_IO, false, "Selected recipe %s is not found in camera recipe", recipeName); + return false; + } if (!psMetadataUpdate (current, recipe)) { - psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", recipeName); - return false; - } + psError(PS_ERR_IO, false, "Failed to update recipe for %s from camera recipe", recipeName); + return false; + } } psFree(recipesIter); Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmReadoutCombine.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmReadoutCombine.c (revision 15516) @@ -315,24 +315,22 @@ // Combination if (!psVectorStats(stats, pixels, errors, mask, 1)) { - psError(PS_ERR_UNKNOWN, false, "Error in statistics."); - psFree(index); - psFree(pixels); - psFree(mask); - psFree(weights); - psFree(errors); - psFree(stats); - psFree(invScale); - return false; - } - - outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); - if (!isfinite(outputImage[yOut][xOut])) { + // Can't do much about it, but it's not worth worrying about + psErrorClear(); + outputImage[yOut][xOut] = NAN; outputMask[yOut][xOut] = params->blank; - } - if (params->weights) { - float stdev = psStatsGetValue(stats, combineStdev); - outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance - // XXXX this is not the correct formal error. - // also, the weighted mean is not obviously the correct thing here + if (params->weights) { + outputWeight[yOut][xOut] = NAN; + } + } else { + outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); + if (!isfinite(outputImage[yOut][xOut])) { + outputMask[yOut][xOut] = params->blank; + } + if (params->weights) { + float stdev = psStatsGetValue(stats, combineStdev); + outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance + // XXXX this is not the correct formal error. + // also, the weighted mean is not obviously the correct thing here + } } } Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmStack.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmStack.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmStack.c (revision 15516) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.13.8.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-29 01:40:54 $ + * @version $Revision: 1.13.8.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-08 22:59:04 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmStackReject.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmStackReject.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmStackReject.c (revision 15516) @@ -52,13 +52,15 @@ // Convolve the image with the kernel --- we're basically applying a matched filter and then thresholding - psImage *convolved = NULL; // Convolved image + pmReadout *convRO = pmReadoutAlloc(NULL); // Readout with convolved image + pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image + inRO->image = image; for (int i = 0; i < numRegions; i++) { psRegion *region = regions->data[i]; // Region of interest psVector *solution = solutions->data[i]; // Solution of interest - if (!pmSubtractionConvolve(&convolved, NULL, NULL, image, NULL, NULL, 0, - region, solution, kernels, true)) { + if (!pmSubtractionConvolve(convRO, inRO, NULL, NULL, 0, region, solution, kernels, + PM_SUBTRACTION_MODE_1, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i); - psFree(convolved); - psFree(image); + psFree(convRO); + psFree(inRO); return NULL; } @@ -73,6 +75,6 @@ if (!kernel) { psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); - psFree(convolved); - psFree(image); + psFree(convRO); + psFree(inRO); return NULL; } @@ -85,8 +87,10 @@ psFree(kernel); - psImage *subConv = psImageSubset(convolved, *region); // Sub-image of convolved image + psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32)); } - psFree(image); + psFree(inRO); + psImage *convolved = psMemIncrRefCounter(convRO->image); + psFree(convRO); // Threshold the convolved image Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtraction.h (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtraction.h (revision 15516) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-17 02:45:40 $ + * @version $Revision: 1.19.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-08 22:59:04 $ * * Copyright 2004-207 Institute for Astronomy, University of Hawaii @@ -16,18 +16,24 @@ #include -#include "pmSubtractionKernels.h" -#include "pmSubtractionStamps.h" + +#include +#include +#include +#include /// @addtogroup imcombine Image Combinations /// @{ +/// Mask values for the subtraction mask typedef enum { - PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking - PM_SUBTRACTION_MASK_REF = 0x01, // Reference image is bad - PM_SUBTRACTION_MASK_INPUT = 0x02, // Input image is bad - PM_SUBTRACTION_MASK_CONVOLVE = 0x04, // If convolved, would be bad - PM_SUBTRACTION_MASK_FOOTPRINT = 0x08, // Bad pixel within the stamp footprint - PM_SUBTRACTION_MASK_BORDER = 0x10, // Image border - PM_SUBTRACTION_MASK_REJ = 0x20, // Previously tried as a stamp, and rejected + PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking + PM_SUBTRACTION_MASK_BAD_1 = 0x01, // Image 1 is bad + PM_SUBTRACTION_MASK_BAD_2 = 0x02, // Image 2 is bad + PM_SUBTRACTION_MASK_CONVOLVE_1 = 0x04, // If image 1 is convolved, would be bad + PM_SUBTRACTION_MASK_CONVOLVE_2 = 0x08, // If image 2 is convolved, would be bad + PM_SUBTRACTION_MASK_FOOTPRINT_1 = 0x10, // Bad pixel within the stamp footprint of image 1 + PM_SUBTRACTION_MASK_FOOTPRINT_2 = 0x20, // Bad pixel within the stamp footprint of image 2 + PM_SUBTRACTION_MASK_BORDER = 0x40, // Image border + PM_SUBTRACTION_MASK_REJ = 0x80, // Previously tried as a stamp, and rejected } pmSubtractionMasks; @@ -45,5 +51,6 @@ bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, ///< Stamp to convolve const pmSubtractionKernels *kernels, ///< Kernel parameters - int footprint ///< Half-size of region over which to calculate equation + int footprint, ///< Half-size of region over which to calculate equation + pmSubtractionMode mode ///< Mode for subtraction (which to convolve) ); @@ -51,5 +58,6 @@ bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps const pmSubtractionKernels *kernels, ///< Kernel parameters - int footprint ///< Half-size of region over which to calculate equation + int footprint, ///< Half-size of region over which to calculate equation + pmSubtractionMode mode ///< Mode for subtraction (which to convolve) ); @@ -59,11 +67,18 @@ ); +/// Calculate deviations +psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, ///< Stamps + const psVector *solution, ///< Solution vector + int footprint, ///< Half-size of stamp + const pmSubtractionKernels *kernels, ///< Kernel parameters + pmSubtractionMode mode ///< Mode for subtraction + ); + /// Reject stamps int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, ///< Stamps + const psVector *deviations, ///< Deviations for each stamp psImage *subMask, ///< Subtraction mask - const psVector *solution, ///< Solution vector - int footprint, ///< Region to mask if stamp is bad float sigmaRej, ///< Number of RMS deviations above zero at which to reject - const pmSubtractionKernels *kernels ///< Kernel parameters + int footprint ///< Half-size of stamp ); @@ -81,9 +96,7 @@ /// Convolve image in preparation for subtraction -bool pmSubtractionConvolve(psImage **outImage, ///< Output image - psImage **outWeight, ///< Output weight map (or NULL) - psImage **outMask, ///< Output mask (or NULL) - const psImage *inImage, ///< Input image - const psImage *inWeight, ///< Input weight map (or NULL) +bool pmSubtractionConvolve(pmReadout *out, ///< Output image + const pmReadout *ro1, // Input image 1 + const pmReadout *ro2, // Input image 2 const psImage *subMask, ///< Subtraction mask (or NULL) psMaskType blank, ///< Mask value for blank regions @@ -91,4 +104,5 @@ const psVector *solution, ///< The solution vector const pmSubtractionKernels *kernels, ///< Kernel parameters + pmSubtractionMode mode, ///< Mode for subtraction bool useFFT ///< Use Fast Fourier Transform for the convolution? ); Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.c (revision 15516) @@ -4,5 +4,5 @@ #include -#include +#include #include @@ -488,5 +488,5 @@ PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); PS_ASSERT_INT_NONNEGATIVE(inner, NULL); - PS_ASSERT_INT_LESS_THAN(inner, size, NULL); + PS_ASSERT_INT_LESS_THAN_OR_EQUAL(inner, size, NULL); PS_ASSERT_INT_NONNEGATIVE(ringsOrder, NULL); @@ -572,5 +572,6 @@ for (int v = -size; v <= size; v++) { int v2 = PS_SQR(v); // Square of v - float vPoly = power(v, vOrder); // Value of v^vOrder +// float vPoly = power(v, vOrder); // Value of v^vOrder + float vPoly = powf(v/(float)size, vOrder); // Value of v^vOrder for (int u = -size; u <= size; u++) { @@ -578,5 +579,6 @@ int distance2 = u2 + v2; // Distance from the centre if (distance2 > lower2 && distance2 < upper2) { - float uPoly = power(u, uOrder); // Value of u^uOrder +// float uPoly = power(u, uOrder); // Value of u^uOrder + float uPoly = powf(u/(float)size, uOrder); // Value of u^uOrder float polyVal = uPoly * vPoly; // Value of polynomial Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.h (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.h (revision 15516) @@ -14,4 +14,14 @@ PM_SUBTRACTION_KERNEL_RINGS, ///< Rings Instead of the Normal Gaussian Subtraction } pmSubtractionKernelsType; + +/// Modes --- specifies which image to convolve +typedef enum { + PM_SUBTRACTION_MODE_ERR, // Error in the mode + PM_SUBTRACTION_MODE_TARGET, // Convolve image 1 to match target PSF + PM_SUBTRACTION_MODE_1, // Convolve image 1 + PM_SUBTRACTION_MODE_2, // Convolve image 2 + PM_SUBTRACTION_MODE_UNSURE, // Not sure yet which image to convolve so try to satisfy both + PM_SUBTRACTION_MODE_DUAL, // Dual convolution +} pmSubtractionMode; /// Kernels specification Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.c (revision 15516) @@ -50,6 +50,6 @@ static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read - const pmReadout *reference, // Reference readout - const pmReadout *input, // Input readout, or NULL to generate fake stamps + const pmReadout *ro1, // Readout 1 + const pmReadout *ro2, // Readout 2 const psImage *subMask, // Mask for subtraction, or NULL psImage *weight, // Weight map @@ -57,11 +57,11 @@ float threshold, // Threshold for stamp finding float stampSpacing, // Spacing between stamps - float targetWidth,// Target width for fake stamps int size, // Kernel half-size - int footprint // Convolution footprint for stamps + int footprint, // Convolution footprint for stamps + pmSubtractionMode mode // Mode for subtraction ) { psTrace("psModules.imcombine", 3, "Finding stamps...\n"); - *stamps = pmSubtractionStampsFind(*stamps, reference->image, subMask, region, threshold, stampSpacing); + *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, stampSpacing, mode); if (!*stamps) { psError(PS_ERR_UNKNOWN, false, "Unable to find stamps."); @@ -71,14 +71,6 @@ memCheck(" find stamps"); - if (!input && !pmSubtractionStampsGenerate(*stamps, targetWidth, footprint, size)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps."); - return false; - } - - memCheck(" generate stamps"); - psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionStampsExtract(*stamps, reference->image, input ? input->image : NULL, - weight, footprint, size)) { + if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, footprint, size)) { psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); return false; @@ -94,40 +86,40 @@ -bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *reference, const pmReadout *input, +bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *ro1, const pmReadout *ro2, int footprint, float regionSize, float stampSpacing, float threshold, - const psArray *sources, const char *stampsName, float targetWidth, + const psArray *sources, const char *stampsName, pmSubtractionKernelsType type, int size, int spatialOrder, const psVector *isisWidths, const psVector *isisOrders, int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad, - psMaskType maskBlank, float badFrac) + psMaskType maskBlank, float badFrac, pmSubtractionMode mode) { PS_ASSERT_PTR_NON_NULL(convolved, false); - PS_ASSERT_PTR_NON_NULL(reference, false); - PS_ASSERT_IMAGE_NON_NULL(reference->image, false); - PS_ASSERT_IMAGE_TYPE(reference->image, PS_TYPE_F32, false); - if (reference->mask) { - PS_ASSERT_IMAGE_NON_NULL(reference->mask, false); - PS_ASSERT_IMAGE_TYPE(reference->mask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(reference->mask, reference->image, false); - } - if (reference->weight) { - PS_ASSERT_IMAGE_NON_NULL(reference->weight, false); - PS_ASSERT_IMAGE_TYPE(reference->weight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(reference->weight, reference->image, false); - } - if (input) { - PS_ASSERT_IMAGE_NON_NULL(input->image, false); - PS_ASSERT_IMAGE_TYPE(input->image, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input->image, reference->image, false); - if (input->mask) { - PS_ASSERT_IMAGE_NON_NULL(input->mask, false); - PS_ASSERT_IMAGE_TYPE(input->mask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input->mask, reference->image, false); - } - if (input->weight) { - PS_ASSERT_IMAGE_NON_NULL(input->weight, false); - PS_ASSERT_IMAGE_TYPE(input->weight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input->weight, reference->image, false); + PS_ASSERT_PTR_NON_NULL(ro1, false); + PS_ASSERT_IMAGE_NON_NULL(ro1->image, false); + PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false); + if (ro1->mask) { + PS_ASSERT_IMAGE_NON_NULL(ro1->mask, false); + PS_ASSERT_IMAGE_TYPE(ro1->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->mask, ro1->image, false); + } + if (ro1->weight) { + PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false); + PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false); + } + if (ro2) { + PS_ASSERT_IMAGE_NON_NULL(ro2->image, false); + PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->image, ro1->image, false); + if (ro2->mask) { + PS_ASSERT_IMAGE_NON_NULL(ro2->mask, false); + PS_ASSERT_IMAGE_TYPE(ro2->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->mask, ro1->image, false); + } + if (ro2->weight) { + PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false); + PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false); } } else if (!stampsName && !sources) { @@ -144,5 +136,4 @@ } // stampsName may be anything - // targetWidth can be just about anything (except maybe negative, but it can be NAN) // We'll check kernel type when we allocate the kernels PS_ASSERT_INT_POSITIVE(size, false); @@ -196,20 +187,18 @@ } - psImage *inImage = NULL, *inMask = NULL; // Input image, mask, weight - if (input) { - inImage = input->image; - inMask = input->mask; - } - // Where does our weight map come from? psImage *weight = NULL; // Weight image to use - if (input && input->weight) { - weight = input->weight; - } else if (reference->weight) { - weight = reference->weight; - } else if (input) { - weight = input->image; + if (ro1->weight && ro2 && ro2->weight) { + weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight); + } else if (ro1->weight) { + weight = psMemIncrRefCounter(ro1->weight); + } else if (ro2) { + if (ro2->weight) { + weight = psMemIncrRefCounter(ro2->weight); + } else { + weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); + } } else { - weight = reference->image; + weight = psMemIncrRefCounter(ro1->image); } @@ -221,12 +210,13 @@ psVector *solution = NULL; // Solution to match PSF pmSubtractionKernels *kernels = NULL; // Kernel basis functions - - int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions + int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions memCheck("start"); - subMask = pmSubtractionMask(reference->mask, inMask, maskBad, size, footprint, badFrac, useFFT); + subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskBad, size, footprint, + badFrac, useFFT); if (!subMask) { psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask."); + psFree(weight); return false; } @@ -260,39 +250,41 @@ if (sources) { - stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, - input ? 0 : 2 * footprint); + stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, mode); } else if (stampsName && strlen(stampsName) > 0) { - // Read stamps from file - psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName); - const char *stampFormat = input ? "%f %f" : "%f %f %f"; // Format for reading stamp file - psArray *stampsData = stampsData = psVectorsReadFromFile(stampsName, stampFormat); - psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes - if (!stampsData) { - psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName); - goto ERROR; - } - xStamp = stampsData->data[0]; - yStamp = stampsData->data[1]; - if (!input) { - fluxStamp = stampsData->data[2]; - } - - // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions - psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - - stamps = pmSubtractionStampsSet(xStamp, yStamp, fluxStamp, reference->image, subMask, - region, stampSpacing, input ? 0 : footprint + size); - psFree(stampsData); + stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, stampSpacing, mode); + } + + // We get the stamps here; we will also attempt to get stamps at the first iteration, but it + // doesn't matter. + if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing, + size, footprint, mode)) { + goto MATCH_ERROR; + } + + if (mode == PM_SUBTRACTION_MODE_UNSURE || mode == PM_SUBTRACTION_MODE_TARGET) { + pmSubtractionMode newMode = pmSubtractionOrder(stamps, footprint); // Subtraction mode + switch (newMode) { + case PM_SUBTRACTION_MODE_1: + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2."); + break; + case PM_SUBTRACTION_MODE_2: + if (mode == PM_SUBTRACTION_MODE_TARGET) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Input PSF is larger than target PSF --- can't match image."); + goto MATCH_ERROR; + } + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1."); + break; + default: + psError(PS_ERR_UNKNOWN, false, "Unable to determine subtraction order."); + goto MATCH_ERROR; + } + mode = newMode; } // Define kernel basis functions if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { - if (!getStamps(&stamps, reference, input, subMask, weight, NULL, - threshold, stampSpacing, targetWidth, size, footprint)) { - goto ERROR; - } kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, - stamps, footprint, optThreshold); + stamps, footprint, optThreshold, mode); if (!kernels) { psErrorClear(); @@ -314,13 +306,13 @@ psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); - if (!getStamps(&stamps, reference, input, subMask, weight, region, - threshold, stampSpacing, targetWidth, size, footprint)) { - goto ERROR; + if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing, + size, footprint, mode)) { + goto MATCH_ERROR; } psTrace("psModules.imcombine", 3, "Calculating equation...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, footprint)) { + if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; + goto MATCH_ERROR; } @@ -331,15 +323,27 @@ if (!solution) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; + goto MATCH_ERROR; } memCheck(" solve equation"); + psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, + mode); // Deviations for each stamp + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + goto MATCH_ERROR; + } + + memCheck(" calculate deviations"); + psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); - numRejected = pmSubtractionRejectStamps(stamps, subMask, solution, footprint, rej, kernels); + numRejected = pmSubtractionRejectStamps(stamps, deviations, subMask, rej, footprint); if (numRejected < 0) { psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); - goto ERROR; - } + psFree(deviations); + goto MATCH_ERROR; + } + psFree(deviations); + memCheck(" reject stamps"); } @@ -350,7 +354,14 @@ if (!solution) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; - } - (void)pmSubtractionRejectStamps(stamps, subMask, solution, footprint, NAN, kernels); + goto MATCH_ERROR; + } + psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, + mode); // Deviations for each stamp + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + goto MATCH_ERROR; + } + pmSubtractionRejectStamps(stamps, deviations, subMask, NAN, footprint); + psFree(deviations); } psFree(stamps); @@ -372,5 +383,5 @@ psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); psFree(convKernels); - goto ERROR; + goto MATCH_ERROR; } @@ -380,5 +391,5 @@ psFree(kernel); psFree(convKernels); - goto ERROR; + goto MATCH_ERROR; } psFree(kernel); @@ -416,9 +427,8 @@ psTrace("psModules.imcombine", 2, "Convolving...\n"); - if (!pmSubtractionConvolve(&convolved->image, &convolved->weight, &convolved->mask, - reference->image, reference->weight, subMask, maskBlank, region, - solution, kernels, useFFT)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image."); - goto ERROR; + if (!pmSubtractionConvolve(convolved, ro1, ro2, subMask, maskBlank, region, + solution, kernels, mode, useFFT)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); + goto MATCH_ERROR; } psFree(kernels); @@ -461,8 +471,10 @@ psFree(subMask); subMask = NULL; + psFree(weight); + weight = NULL; if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) { psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image."); - goto ERROR; + goto MATCH_ERROR; } @@ -472,5 +484,5 @@ return true; -ERROR: +MATCH_ERROR: psFree(region); psFree(regionString); @@ -479,4 +491,131 @@ psFree(stamps); psFree(solution); + psFree(weight); return false; } + +// Calculate the second order moments for an image +static float subtractionOrderMoment(const psKernel *kernel, // Image for which to measure moments + int radius // Maximum radius + ) +{ + assert(kernel && kernel->kernel); + + int xMin = PS_MAX(kernel->xMin, -radius), xMax = PS_MIN(kernel->xMax, radius); // Bounds in x + int yMin = PS_MAX(kernel->yMin, -radius), yMax = PS_MIN(kernel->yMax, radius); // Bounds in y + + float xCentroid = 0.0, yCentroid = 0.0; // Centroid (first moment) + float sum = 0.0; // Sum (zero-th moment) + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + xCentroid += kernel->kernel[y][x] * x; + yCentroid += kernel->kernel[y][x] * y; + sum += kernel->kernel[y][x]; + } + } + xCentroid /= sum; + yCentroid /= sum; + + float eta20 = 0.0, eta02 = 0.0; // Second moments + for (int y = yMin; y <= yMax; y++) { + float yDiff = y - yCentroid; + for (int x = xMin; x <= xMax; x++) { + float xDiff = x - xCentroid; + eta20 += PS_SQR(xDiff) * kernel->kernel[y][x]; + eta02 += PS_SQR(yDiff) * kernel->kernel[y][x]; + } + } + + // Normalise to calculate the scale-invariant + float sum2 = PS_SQR(sum); + eta20 /= sum2; + eta02 /= sum2; + // eta11 /= sum2; + + return eta20 + eta02; +} + +#if 0 +// Calculate the deviations for a particular subtraction order +static psVector *subtractionOrderDeviation(float *sumKernel, // Sum of the kernel + pmSubtractionStampList *stamps, // Stamps to convolve + const pmSubtractionKernels *kernels, // Kernel basis functions + int footprint, // Stamp footprint + pmSubtractionMode mode // Mode of subtraction + ) +{ + assert(stamps); + assert(footprint >= 0); + assert(mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_2); + + if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return NULL; + } + + psVector *solution = pmSubtractionSolveEquation(NULL, stamps); + if (!solution) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return NULL; + } + + if (sumKernel) { + float sum = 0.0; // Sum of the kernel + psImage *image = pmSubtractionKernelImage(solution, kernels, 0.0, 0.0); // Image of kernel + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + sum += image->data.F32[y][x]; + } + } + psFree(image); + *sumKernel = sum; + } + + psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, mode); + psFree(solution); + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + return NULL; + } + + return deviations; +} +#endif + +pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, int radius) +{ + PS_ASSERT_INT_POSITIVE(radius, PM_SUBTRACTION_MODE_ERR); + + psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask for stamps + psVector *moments = psVectorAlloc(stamps->num, PS_TYPE_F32); // Moments + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) { + mask->data.PS_TYPE_MASK_DATA[i] = 0xff; + continue; + } + mask->data.PS_TYPE_MASK_DATA[i] = 0; + moments->data.F32[i] = subtractionOrderMoment(stamp->image1, radius) / + subtractionOrderMoment(stamp->image2, radius); + } + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + if (!psVectorStats(stats, moments, NULL, mask, 0xff)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio."); + psFree(mask); + psFree(moments); + psFree(stats); + return PM_SUBTRACTION_MODE_ERR; + } + psFree(moments); + psFree(mask); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median ratio of second moments: %lf", stats->robustMedian); + pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2); + psFree(stats); + + return mode; +} + + + Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.h (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionMatch.h (revision 15516) @@ -4,12 +4,13 @@ #include -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmSubtractionKernels.h" +#include +#include +#include +#include - +/// Match two images bool pmSubtractionMatch(pmReadout *convolved, ///< Output convolved data - const pmReadout *reference, ///< Reference data - const pmReadout *input, ///< Input data + const pmReadout *ro1, ///< Image 1 + const pmReadout *ro2, ///< Image 2 // Stamp parameters int footprint, ///< Stamp half-size @@ -19,5 +20,4 @@ const psArray *sources, ///< Sources for stamps const char *stampsName, ///< Filename for stamps - float targetWidth, ///< Width of PSF for simulated target // Kernel parameters pmSubtractionKernelsType type, ///< Kernel type @@ -38,6 +38,13 @@ psMaskType maskBad, ///< Value to mask psMaskType maskBlank, ///< Mask for blank region - float badFrac ///< Maximum fraction of bad input pixels to accept + float badFrac, ///< Maximum fraction of bad input pixels to accept + pmSubtractionMode mode ///< Mode of subtraction ); +/// Determine which image to convolve +pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, ///< Stamps that have been extracted + int footprint ///< Stamp half-size + ); + + #endif Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionParams.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionParams.c (revision 15516) @@ -50,4 +50,21 @@ #endif +/// Select the appropriate convolution, given the kernel basis function and subtraction mode +static inline psKernel *selectConvolution(const pmSubtractionStamp *stamp, // Stamp + int kernelIndex, // Index for kernel component + pmSubtractionMode mode // Mode of subtraction + ) +{ + switch (mode) { + case PM_SUBTRACTION_MODE_1: + return stamp->convolutions1->data[kernelIndex]; + case PM_SUBTRACTION_MODE_2: + return stamp->convolutions2->data[kernelIndex]; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + return NULL; // Unreached +} + // Accumulate cross-term sums for a stamp static void accumulateCross(double *sumI, // Sum of I(x)/sigma(x)^2 @@ -55,14 +72,15 @@ double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2 const pmSubtractionStamp *stamp, // Stamp with weight - const psKernel *input, // Input image, I(x) + const psKernel *target, // Target stamp int kernelIndex, // Index for kernel component - int footprint // Size of region of interest + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction ) { psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution @@ -82,9 +100,10 @@ const pmSubtractionStamp *stamp, // Stamp with input and weight int kernelIndex, // Index for kernel component - int footprint // Size of region of interest + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction ) { psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { @@ -100,18 +119,19 @@ } -static double accumulateChi2(psKernel *input, // Input stamp +static double accumulateChi2(const psKernel *target, // Target stamp pmSubtractionStamp *stamp, // Stamp with weight int kernelIndex, // Index for kernel component double coeff, // Coefficient of convolution double bg, // Background term - int footprint // Size of region of interest + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction ) { double chi2 = 0.0; psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution @@ -125,17 +145,28 @@ // Return the initial value of chi^2 -static double initialChi2(psKernel *input, // Input stamp +static double initialChi2(const psKernel *target, // Target stamp const pmSubtractionStamp *stamp, // Stamp with weight - int footprint // Size of convolution + int footprint, // Size of convolution + pmSubtractionMode mode // Mode of subtraction ) { psKernel *weight = stamp->weight; // Weight map - psKernel *reference = stamp->reference; // Reference stamp + psKernel *source; // Source stamp + switch (mode) { + case PM_SUBTRACTION_MODE_1: + source = stamp->image1; + break; + case PM_SUBTRACTION_MODE_2: + source = stamp->image2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } double chi2 = 0.0; // Chi^2 for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight - psF32 *ref = &reference->kernel[y][-footprint]; // Derference reference + psF32 *ref = &source->kernel[y][-footprint]; // Derference reference for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { float diff = *in - *ref; // Temporary value @@ -148,16 +179,16 @@ // Subtract a convolution from the input -static void subtractConvolution(psKernel *input, // Input stamp +static void subtractConvolution(psKernel *target, // Target stamp const pmSubtractionStamp *stamp, // Stamp with weight int kernelIndex, // Index for kernel component float coeff, // Coefficient of subtraction float bg, // Background term - int footprint // Size of region of interest - ) -{ - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest - + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction + ) +{ + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, in++, conv++) { @@ -173,5 +204,5 @@ int spatialOrder, const psVector *fwhms, int maxOrder, const pmSubtractionStampList *stamps, int footprint, - float tolerance) + float tolerance, pmSubtractionMode mode) { if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { @@ -210,5 +241,5 @@ // Need to save the stamp inputs --- we're changing the values! int numStamps = stamps->num; // Number of stamps - psArray *inputs = psArrayAlloc(numStamps); // Deep copies of the inputs + psArray *targets = psArrayAlloc(numStamps); // Deep copies of the targets psVector *badStamps = psVectorAlloc(numStamps, PS_TYPE_U8); // Mark the bad stamps psVectorInit(badStamps, 0); @@ -219,7 +250,17 @@ continue; } - psKernel *input = stamp->input; // Input image of interest - psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image - inputs->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint); + psKernel *target; // Target image of interest + switch (mode) { + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + psImage *copy = psImageCopy(NULL, target->image, PS_TYPE_F32); // Copy of the image + targets->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint); psFree(copy); // Drop reference } @@ -238,7 +279,7 @@ continue; } - if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { + if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i); - psFree(inputs); + psFree(targets); psFree(kernels); psFree(badStamps); @@ -258,5 +299,5 @@ "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n", i, (int)stamp->x, (int)stamp->y); - psFree(inputs); + psFree(targets); psFree(kernels); psFree(badStamps); @@ -265,8 +306,8 @@ for (int j = 0; j < numKernels; j++) { - accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint); - } - - lastChi2 += initialChi2(inputs->data[i], stamp, footprint); + accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint, mode); + } + + lastChi2 += initialChi2(targets->data[i], stamp, footprint, mode); numPixels += PS_SQR(2 * footprint + 1); } @@ -297,5 +338,5 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest - accumulateCross(&sumI, &sumII, &sumIC, stamp, inputs->data[j], i, footprint); + accumulateCross(&sumI, &sumII, &sumIC, stamp, targets->data[j], i, footprint, mode); } @@ -310,5 +351,5 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest - chi2 += accumulateChi2(inputs->data[j], stamp, i, coeff, bg, footprint); + chi2 += accumulateChi2(targets->data[j], stamp, i, coeff, bg, footprint, mode); } @@ -328,5 +369,5 @@ if (bestIndex == -1) { psError(PS_ERR_UNKNOWN, false, "Unable to find best kernel component in round %d.", iter); - psFree(inputs); + psFree(targets); psFree(sumC); psFree(sumCC); @@ -345,5 +386,5 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest - subtractConvolution(inputs->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint); + subtractConvolution(targets->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint, mode); } @@ -361,5 +402,5 @@ lastChi2 = bestChi2; } - psFree(inputs); + psFree(targets); psFree(sumC); psFree(sumCC); @@ -401,5 +442,5 @@ pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest psArray *convolutions = convNew->data[j]; // Convolutions for this stamp - convolutions->data[rank] = psMemIncrRefCounter(stamp->convolutions->data[i]); + convolutions->data[rank] = psMemIncrRefCounter(selectConvolution(stamp, i, mode)); } } @@ -420,6 +461,18 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - psFree(stamp->convolutions); - stamp->convolutions = convNew->data[i]; + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + switch (mode) { + case PM_SUBTRACTION_MODE_1: + stamp->convolutions1 = convNew->data[i]; + stamp->convolutions2 = NULL; + break; + case PM_SUBTRACTION_MODE_2: + stamp->convolutions1 = NULL; + stamp->convolutions2 = convNew->data[i]; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } } Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionParams.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionParams.h (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionParams.h (revision 15516) @@ -15,5 +15,6 @@ const pmSubtractionStampList *stamps, ///< Stamps int footprint, ///< Convolution footprint for stamps - float tolerance ///< Maximum difference in chi^2 + float tolerance, ///< Maximum difference in chi^2 + pmSubtractionMode mode // Mode for subtraction ); Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionStamps.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionStamps.c (revision 15516) @@ -4,4 +4,5 @@ #include +#include #include @@ -45,8 +46,9 @@ psFree(stamp->matrix); psFree(stamp->vector); - psFree(stamp->reference); - psFree(stamp->input); + psFree(stamp->image1); + psFree(stamp->image2); psFree(stamp->weight); - psFree(stamp->convolutions); + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); } @@ -65,5 +67,6 @@ // Is this position unmasked? static bool checkStampMask(int x, int y, // Coordinates of stamp - const psImage *mask + const psImage *mask, // Mask + pmSubtractionMode mode // Mode for subtraction ) { @@ -74,7 +77,23 @@ return false; } - return (mask->data.PS_TYPE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER | - PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_REJ)) ? - false : true; + + psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ; // Mask value + switch (mode) { + case PM_SUBTRACTION_MODE_1: + case PM_SUBTRACTION_MODE_TARGET: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1; + break; + case PM_SUBTRACTION_MODE_2: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_2; + break; + case PM_SUBTRACTION_MODE_UNSURE: + case PM_SUBTRACTION_MODE_DUAL: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1 | PM_SUBTRACTION_MASK_FOOTPRINT_2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + + return (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? false : true; } @@ -140,8 +159,9 @@ stamp->status = PM_SUBTRACTION_STAMP_INIT; - stamp->reference = NULL; - stamp->input = NULL; + stamp->image1 = NULL; + stamp->image2 = NULL; stamp->weight = NULL; - stamp->convolutions = NULL; + stamp->convolutions1 = NULL; + stamp->convolutions2 = NULL; return stamp; @@ -151,5 +171,5 @@ pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image, const psImage *subMask, const psRegion *region, - float threshold, float spacing) + float threshold, float spacing, pmSubtractionMode mode) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -226,5 +246,5 @@ for (int y = subRegion->y0; y <= subRegion->y1 ; y++) { for (int x = subRegion->x0; x <= subRegion->y1 ; x++) { - if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxStamp) { + if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) { fluxStamp = image->data.F32[y][x]; xStamp = x; @@ -242,10 +262,13 @@ // Reset the postage stamps since we're making a new stamp - psFree(stamp->reference); - psFree(stamp->input); + psFree(stamp->image1); + psFree(stamp->image2); psFree(stamp->weight); - stamp->reference = stamp->input = stamp->weight = NULL; - - stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + stamp->image1 = stamp->image2 = stamp->weight = NULL; + stamp->convolutions1 = stamp->convolutions2 = NULL; + + stamp->status = PM_SUBTRACTION_STAMP_FOUND; numFound++; psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n", @@ -266,5 +289,6 @@ pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux, const psImage *image, const psImage *subMask, - const psRegion *region, float spacing, int exclusionZone) + const psRegion *region, float spacing, pmSubtractionMode mode) + { PS_ASSERT_VECTOR_NON_NULL(x, NULL); @@ -289,46 +313,6 @@ } PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); - PS_ASSERT_INT_NONNEGATIVE(exclusionZone, NULL); int numStars = x->n; // Number of stars - psVector *exclude = psVectorAlloc(numStars, PS_TYPE_U8); // Exclude a star? - psVectorInit(exclude, 0); - - // Apply exclusion zone around each star --- important when we're convolving to a specified PSF; less so - // when we're trying to get a good subtraction. - if (exclusionZone > 0) { - psTrace("psModules.imcombine", 2, "Applying exclusion zone of %d pixels for stamps\n", exclusionZone); - // We use something resembling a binary search --- coordinates are sorted in the x dimension, and then - // to exclude stars within a nominated distance, we need only examine (i.e., calculate the - // 2-dimensional distance to) other stars in the list that are within that distance of the x - // coordinate of the star of interest. By marking both stars when we find a violation of the - // exclusion zone, we need only search upwards from the x coordinate of the star of interest. - - int minDistance2 = PS_SQR(exclusionZone); // Minimum square distance for other stars - psVector *indexes = psVectorSortIndex(NULL, x); // Indices for sorting in x - for (int i = 0; i < numStars - 1; i++) { - int iIndex = indexes->data.S32[i]; // Index for star i - if (exclude->data.U8[iIndex]) { - continue; - } - float ix = x->data.F32[iIndex], iy = y->data.F32[iIndex]; // Coordinates for star i - float jx, jy; // Coordinates for star j - for (int j = i + 1, jIndex = indexes->data.S32[j]; - j < numStars && (jx = x->data.F32[jIndex]) < ix + exclusionZone; - j++, jIndex = indexes->data.S32[j]) { - jy = y->data.F32[jIndex]; - if (PS_SQR(ix - jx) + PS_SQR(iy - jy) < minDistance2) { - exclude->data.U8[iIndex] = 0xff; - exclude->data.U8[jIndex] = 0xff; - psTrace("psModules.imcombine", 5, "Excluding stamps %d,%d and %d,%d\n", - (int)x->data.F32[iIndex], (int)y->data.F32[iIndex], - (int)x->data.F32[jIndex], (int)y->data.F32[jIndex]); - // Would 'break' here, but there might be more stamps within the exclusion zone. - } - } - } - psFree(indexes); - } - pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, region, spacing); // Stamp list @@ -347,8 +331,4 @@ // Put the stars into their appropriate subregions for (int i = 0; i < numStars; i++) { - if (exclude->data.U8[i]) { - // Star has been excluded - continue; - } float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp @@ -357,5 +337,5 @@ continue; } - if (!checkStampMask(xPix, yPix, subMask)) { + if (!checkStampMask(xPix, yPix, subMask, mode)) { // Not a good stamp continue; @@ -386,5 +366,4 @@ } } - psFree(exclude); // Sort the list by flux, with the brightest last @@ -421,35 +400,27 @@ -bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *reference, psImage *input, +bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2, psImage *weight, int footprint, int kernelSize) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); - PS_ASSERT_IMAGE_NON_NULL(reference, false); - PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false); - if (input) { - PS_ASSERT_IMAGE_NON_NULL(input, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false); - PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false); - } - if (weight) { - PS_ASSERT_IMAGE_NON_NULL(weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(weight, reference, false); - PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); - } + PS_ASSERT_IMAGE_NON_NULL(image1, false); + PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, false); + if (image2) { + PS_ASSERT_IMAGE_NON_NULL(image2, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(image2, image1, false); + PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false); + } + PS_ASSERT_IMAGE_NON_NULL(weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false); + PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); PS_ASSERT_INT_NONNEGATIVE(footprint, false); PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); - int numCols = reference->numCols, numRows = reference->numRows; // Size of images + int numCols = image1->numCols, numRows = image1->numRows; // Size of images int size = kernelSize + footprint; // Size of postage stamps - - if (!weight) { - // Use the input (or if I must, the reference) as a rough approximation to the variance map, and HOPE - // that it's positive. - weight = input ? input : reference; - } for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { + if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_FOUND) { continue; } @@ -468,19 +439,26 @@ } + // Catch memory leaks --- these should have been freed and NULLed before + assert(stamp->image1 == NULL); + assert(stamp->image2 == NULL); + assert(stamp->weight == NULL); + psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest - psImage *refSub = psImageSubset(reference, region); // Subimage with stamp - stamp->reference = psKernelAllocFromImage(refSub, size, size); - psFree(refSub); // Drop reference - - if (input) { - psImage *inSub = psImageSubset(input, region); // Subimage with stamp - stamp->input = psKernelAllocFromImage(inSub, size, size); - psFree(inSub); // Drop reference + psImage *sub1 = psImageSubset(image1, region); // Subimage with stamp + stamp->image1 = psKernelAllocFromImage(sub1, size, size); + psFree(sub1); // Drop reference + + if (image2) { + psImage *sub2 = psImageSubset(image2, region); // Subimage with stamp + stamp->image2 = psKernelAllocFromImage(sub2, size, size); + psFree(sub2); // Drop reference } psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp stamp->weight = psKernelAllocFromImage(wtSub, size, size); - psFree(wtSub); // Drop reference + psFree(wtSub); // Drop reference + + stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; } @@ -488,5 +466,5 @@ } - +#if 0 bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize) { @@ -517,7 +495,7 @@ float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame - psFree(stamp->input); - stamp->input = psKernelAlloc(-size, size, -size, size); - psKernel *input = stamp->input; // Target stamp + psFree(stamp->image2); + stamp->image2 = psKernelAlloc(-size, size, -size, size); + psKernel *target = stamp->image2; // Target stamp // Put in a Waussian, just for fun! @@ -525,5 +503,5 @@ for (int u = -size; u <= size; u++) { float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma)); - input->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z)); + target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z)); } } @@ -533,9 +511,9 @@ return true; } - +#endif pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask, const psRegion *region, float spacing, - int exclusionZone) + pmSubtractionMode mode) { PS_ASSERT_ARRAY_NON_NULL(sources, NULL); @@ -561,5 +539,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, - spacing, exclusionZone); // Stamps to return + spacing, mode); // Stamps to return psFree(x); psFree(y); @@ -572,2 +550,31 @@ return stamps; } + + +pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *subMask, + const psRegion *region, float spacing, + pmSubtractionMode mode) +{ + PS_ASSERT_STRING_NON_EMPTY(filename, NULL); + // Let pmSubtractionStampsSet take care of the rest of the assertions + + const char *format = (mode == PM_SUBTRACTION_MODE_TARGET ? "%f %f %f" : "%f %f"); // Format of file + psArray *data = psVectorsReadFromFile(filename, format); + if (!data) { + psError(PS_ERR_IO, false, "Unable to read stamps file %s", filename); + return NULL; + } + psVector *x = data->data[0], *y = data->data[1]; // Stamp positions + psVector *flux = (mode == PM_SUBTRACTION_MODE_TARGET ? data->data[2] : NULL); // Stamp fluxes + + // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions + psBinaryOp(x, x, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + + pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, spacing, + mode); + psFree(data); + + return stamps; + +} Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionStamps.h (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionStamps.h (revision 15516) @@ -9,7 +9,8 @@ typedef enum { PM_SUBTRACTION_STAMP_INIT, ///< Initial state + PM_SUBTRACTION_STAMP_FOUND, ///< Found a suitable source for this stamp + PM_SUBTRACTION_STAMP_CALCULATE, ///< Calculate matrix and vector values for this stamp PM_SUBTRACTION_STAMP_USED, ///< Use this stamp PM_SUBTRACTION_STAMP_REJECTED, ///< This stamp has been rejected - PM_SUBTRACTION_STAMP_CALCULATE, ///< Calculate matrix and vector values for this stamp PM_SUBTRACTION_STAMP_NONE ///< No stamp in this region } pmSubtractionStampStatus; @@ -54,8 +55,9 @@ float flux; ///< Flux float xNorm, yNorm; ///< Normalised position - psKernel *reference; ///< Reference image postage stamp - psKernel *input; ///< Input image postage stamp + psKernel *image1; ///< Reference image postage stamp + psKernel *image2; ///< Input image postage stamp psKernel *weight; ///< Weight image postage stamp - psArray *convolutions; ///< Convolutions of the reference for each kernel component + psArray *convolutions1; ///< Convolutions of image 1 for each kernel component + psArray *convolutions2; ///< Convolutions of image 2 for each kernel component psImage *matrix; ///< Associated matrix psVector *vector; ///< Assoicated vector @@ -72,11 +74,9 @@ const psRegion *region, ///< Region to search, or NULL float threshold, ///< Threshold for stamps in the image - float spacing ///< Rough spacing for stamps + float spacing, ///< Rough spacing for stamps + pmSubtractionMode mode ///< Mode for subtraction ); /// Set stamps based on a list of x,y -/// -/// We may optionally apply an exclusion zone around each star --- this is important when we're convolving to -/// a specified PSF; less so when we're only trying to get a good subtraction. pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp const psVector *y, ///< y coordinates for each stamp @@ -86,14 +86,23 @@ const psRegion *region, ///< Region to search, or NULL float spacing, ///< Rough spacing for stamps - int exclusionZone ///< Exclusion zone around each star + pmSubtractionMode mode ///< Mode for subtraction ); -/// +/// Set stamps based on a list of sources pmSubtractionStampList *pmSubtractionStampsSetFromSources( - const psArray *sources, ///< Sources for each stamp - const psImage *subMask, ///< Mask, or NULL - const psRegion *region, ///< Region to search, or NULL - float spacing, ///< Rough spacing for stamps - int exclusionZone ///< Exclusion zone around each star + const psArray *sources, ///< Sources for each stamp + const psImage *subMask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + float spacing, ///< Rough spacing for stamps + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Set stamps based on values in a file +pmSubtractionStampList *pmSubtractionStampsSetFromFile( + const char *filename, ///< Filename of file containing x,y (or x,y,flux) on each line + const psImage *subMask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + float spacing, ///< Rough spacing for stamps + pmSubtractionMode mode ///< Mode for subtraction ); @@ -107,10 +116,3 @@ ); -/// Generate target for stamps based on list -bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, ///< Stamps - float fwhm, ///< FWHM for each stamp - int footprint, ///< Stamp footprint size - int size ///< Kernel half-size - ); - #endif Index: /branches/eam_branch_20071023/psModules/src/objects/pmTrend2D.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/objects/pmTrend2D.c (revision 15515) +++ /branches/eam_branch_20071023/psModules/src/objects/pmTrend2D.c (revision 15516) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.4.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-29 22:42:34 $ + * @version $Revision: 1.4.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-08 22:59:04 $ * * Copyright 2004 Institute for Astronomy, University of Hawaii @@ -11,4 +11,8 @@ // XXXX: Ignore (2) +// XXXX: ignore + +// XXXX: ignore + #ifdef HAVE_CONFIG_H #include @@ -32,5 +36,8 @@ pmTrend2D *pmTrend2DAlloc (pmTrend2DMode mode, psImage *image, int nXtrend, int nYtrend, psStats *stats) { - assert (image); + PS_ASSERT_PTR_NON_NULL(stats, NULL); + if (mode == PM_TREND_MAP) { + PS_ASSERT_PTR_NON_NULL(image, NULL); + } pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D)); @@ -73,9 +80,15 @@ break; } - + // XXX: Put a more graceful error here. default: psAbort ("error"); } return (trend); +} + +bool psMemCheckTrend2D(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmTrend2DFree); } @@ -123,4 +136,5 @@ pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats) { + PS_ASSERT_PTR_NON_NULL(stats, NULL); pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D)); psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree); @@ -164,4 +178,5 @@ default: + // XXX: Put a more graceful error here. psAbort ("error"); } @@ -169,5 +184,8 @@ } -bool pmTrend2DFit (pmTrend2D *trend, psVector *mask, psMaskType maskVal, psVector *x, psVector *y, psVector *f, psVector *df) { +bool pmTrend2DFit (pmTrend2D *trend, psVector *mask, psMaskType maskVal, psVector *x, + psVector *y, psVector *f, psVector *df) +{ + PS_ASSERT_PTR_NON_NULL(trend, false); bool status; @@ -199,10 +217,9 @@ } -double pmTrend2DEval (pmTrend2D *trend, float x, float y) { +double pmTrend2DEval (pmTrend2D *trend, float x, float y) +{ + if (!trend) return 0.0; double result; - - assert (trend); - switch (trend->mode) { case PM_TREND_POLY_ORD: @@ -221,6 +238,7 @@ } -psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y) { - +psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y) +{ + PS_ASSERT_PTR_NON_NULL(trend, NULL); psVector *result;