Index: /branches/eam_branch_20071023/psModules/src/camera/Makefile.am =================================================================== --- /branches/eam_branch_20071023/psModules/src/camera/Makefile.am (revision 15398) +++ /branches/eam_branch_20071023/psModules/src/camera/Makefile.am (revision 15399) @@ -27,5 +27,6 @@ pmFPAExtent.c \ pmCellSquish.c \ - pmReadoutStack.c + pmReadoutStack.c \ + pmReadoutFake.c pkginclude_HEADERS = \ @@ -53,5 +54,6 @@ pmFPAExtent.h \ pmCellSquish.h \ - pmReadoutStack.h + pmReadoutStack.h \ + pmReadoutFake.h CLEANFILES = *~ Index: /branches/eam_branch_20071023/psModules/src/camera/pmReadoutFake.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/camera/pmReadoutFake.c (revision 15399) +++ /branches/eam_branch_20071023/psModules/src/camera/pmReadoutFake.c (revision 15399) @@ -0,0 +1,111 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmModelClass.h" +#include "pmPeaks.h" +#include "pmSource.h" +#include "pmSourceUtils.h" + +#define MODEL_TYPE "PS_MODEL_RGAUSS" // Type of model to use + + +pmReadout *pmReadoutFakeFromSources(int numCols, int numRows, const psArray *sources, + float fwhm, float minFlux) +{ + PS_ASSERT_INT_LARGER_THAN(numCols, 0, NULL); + PS_ASSERT_INT_LARGER_THAN(numRows, 0, NULL); + PS_ASSERT_ARRAY_NON_NULL(sources, NULL); + PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, NULL); + PS_ASSERT_FLOAT_LARGER_THAN(minFlux, 0.0, NULL); + + pmReadout *readout = pmReadoutAlloc(NULL); // Output readout + readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + + int numSources = sources->n; // Number of stars + + pmModelType modelType = pmModelClassGetType(MODEL_TYPE); // Type of PSF model + assert(modelType >= 0); + pmModel *fakeModel = pmModelAlloc(modelType); + + float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma + + fakeModel->params->data.F32[PM_PAR_SKY] = 0.0; + fakeModel->params->data.F32[PM_PAR_I0] = 1.0; + fakeModel->params->data.F32[PM_PAR_XPOS] = NAN; + fakeModel->params->data.F32[PM_PAR_YPOS] = NAN; + fakeModel->params->data.F32[PM_PAR_SXX] = sigma * M_SQRT2; + fakeModel->params->data.F32[PM_PAR_SYY] = sigma * M_SQRT2; + fakeModel->params->data.F32[PM_PAR_SXY] = 0.0; + switch (modelType) { + case 0: // GAUSS + case 1: // PGAUSS + break; + case 2: // QGAUSS + fakeModel->params->data.F32[PM_PAR_7] = 1.0; + break; + case 3: // RGAUSS + fakeModel->params->data.F32[PM_PAR_7] = 2.0; + break; + default: + psAbort("Unsupported model type: %d", modelType); + } + + float flux0 = fakeModel->modelFlux(fakeModel->params); // Flux for central intensity of 1.0 + + + for (int i = 0; i < numSources; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!isfinite(source->psfMag)) { + continue; + } + float x, y; // Coordinates of source + if (source->modelPSF) { + x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + } else { + x = source->peak->xf; + y = source->peak->yf; + } + + fakeModel->params->data.F32[PM_PAR_XPOS] = x; + fakeModel->params->data.F32[PM_PAR_YPOS] = y; + fakeModel->params->data.F32[PM_PAR_I0] = powf(10.0, -0.4 * source->psfMag) / flux0; + + pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate + fakeSource->peak = pmPeakAlloc(x, y, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE); + float radius = fakeModel->modelRadius(fakeModel->params, minFlux); // Radius of interest for source + + if (!pmSourceDefinePixels(fakeSource, readout, x, y, radius)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source."); + psFree(readout); + psFree(fakeModel); + return NULL; + } + + if (!pmModelAdd(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add model of source to image."); + psFree(readout); + psFree(fakeModel); + return NULL; + } + psFree(fakeSource); + } + + psFree(fakeModel); + + return readout; +} Index: /branches/eam_branch_20071023/psModules/src/camera/pmReadoutFake.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/camera/pmReadoutFake.h (revision 15399) +++ /branches/eam_branch_20071023/psModules/src/camera/pmReadoutFake.h (revision 15399) @@ -0,0 +1,18 @@ +#ifndef PM_READOUT_FAKE_H +#define PM_READOUT_FAKE_H + +#include +#include +#include + +// Generate a fake readout from an array of sources +pmReadout *pmReadoutFakeFromSources(int numCols, int numRows, ///< Dimension of image + const psArray *sources, ///< Array of pmSource + float fwhm, ///< FWHM for sources + float minFlux ///< Minimum flux to bother about; for setting object size + ); + + + + +#endif Index: /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsAverage.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsAverage.c (revision 15398) +++ /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsAverage.c (revision 15399) @@ -97,4 +97,5 @@ int xBin = 0, yBin = 0; // Binning int x0 = 0, y0 = 0; // Offset + int xParity = 0, yParity = 0; // Parity int nCells = 0; // Number of cells; @@ -121,6 +122,8 @@ if (same) { // Only makes sense to update these if they are the same cell - x0 = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - y0 = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + x0 = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + y0 = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + xParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + yParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); } } else { @@ -154,4 +157,14 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Differing CELL.Y0 in use: %d vs %d\n", y0, psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0")); + success = false; + } + if (xParity != psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY")) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Differing CELL.XPARITY in use: %d vs %d\n", + xParity, psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY")); + success = false; + } + if (yParity != psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY")) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Differing CELL.YPARITY in use: %d vs %d\n", + yParity, psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY")); success = false; } @@ -189,4 +202,6 @@ MD_UPDATE(target->concepts, "CELL.X0", S32, x0); MD_UPDATE(target->concepts, "CELL.Y0", S32, y0); + MD_UPDATE(target->concepts, "CELL.XPARITY", S32, xParity); + MD_UPDATE(target->concepts, "CELL.YPARITY", S32, yParity); } Index: /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsAverage.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsAverage.h (revision 15398) +++ /branches/eam_branch_20071023/psModules/src/concepts/pmConceptsAverage.h (revision 15399) @@ -4,6 +4,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-11 01:05:59 $ + * @version $Revision: 1.9.10.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-29 01:40:54 $ * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ @@ -42,4 +42,5 @@ /// - CELL.TIMESYS /// - CELL.X0, CELL.Y0 +/// - CELL.XPARITY, CELL.YPARITY /// And for others, it takes the "worst" possible value: /// - CELL.SATURATION @@ -47,4 +48,5 @@ /// These concepts are only handled if the cells are all the same cell (mosaicking vs stacking): /// - CELL.X0, CELL.Y0 +/// - CELL.XPARITY, CELL.YPARITY bool pmConceptsAverageCells(pmCell *target,///< Target cell psList *sources, ///< List of source cells Index: /branches/eam_branch_20071023/psModules/src/config/pmConfig.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/config/pmConfig.c (revision 15398) +++ /branches/eam_branch_20071023/psModules/src/config/pmConfig.c (revision 15399) @@ -37,4 +37,7 @@ #define PS_DEFAULT_SITE ".ipprc" // Default site config file +#define DEFAULT_LOG STDERR_FILENO // Default file descriptor for log messages +#define DEFAULT_TRACE STDERR_FILENO // Default file descriptor for trace messages + static bool readCameraConfig = true; // Read the camera config on startup (with pmConfigRead)? static psArray *configPath = NULL; // Search path for configuration files @@ -59,4 +62,14 @@ psFree(config->arguments); psFree(config->database); + + // Close log and trace files + if (config->logFD != STDOUT_FILENO && config->logFD != STDERR_FILENO) { + close(config->logFD); + } + if (config->traceFD != STDOUT_FILENO && config->traceFD != STDERR_FILENO) { + close(config->traceFD); + } + + return; } @@ -79,4 +92,7 @@ config->defaultRecipe = NULL; + config->traceFD = DEFAULT_TRACE; + config->logFD = DEFAULT_LOG; + // the file structure is used to carry pmFPAfiles config->files = psMetadataAlloc (); @@ -124,4 +140,5 @@ } + void pmConfigSet(const char *path) { @@ -144,4 +161,6 @@ psFree(configPath); configPath = NULL; + + return; } @@ -379,6 +398,5 @@ psArgumentRemove(argNum, argc, argv); if (argNum >= *argc) { - psLogMsg("psModules.config", PS_LOG_WARN, - "-site command-line switch provided without the required filename --- ignored.\n"); + psWarning("-site command-line switch provided without the required filename --- ignored.\n"); } else { siteName = psStringCopy(argv[argNum]); @@ -446,45 +464,37 @@ argNum = psArgumentGet(*argc, argv, "-log"); - if (argNum > 0) - { + if (argNum > 0) { psArgumentRemove(argNum, argc, argv); if (argNum >= *argc) { - psLogMsg("psModules.config", PS_LOG_WARN, "-log command-line switch provided without the " - "required log destination --- ignored.\n"); + psWarning("-log command-line switch provided without the required log destination " + "--- ignored.\n"); } else { - if (!psLogSetDestination(psMessageDestination(argv[argNum]))) { - psLogMsg("psModules.config", PS_LOG_WARN, "Unable to set log destination to %s\n", - argv[argNum]); - } + config->logFD = psMessageDestination(argv[argNum]); psArgumentRemove(argNum, argc, argv); } - } else - { + } else { // Set logging destination psString logDest = psMetadataLookupStr(&mdok, config->site, "LOGDEST"); if (mdok && logDest && strlen(logDest) > 0) { - // 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 (!psLogSetDestination(psMessageDestination(logDest))) { - psLogMsg("psModules.config", PS_LOG_WARN, "Unable to set log destination to %s\n", - argv[argNum]); - } - } + config->logFD = psMessageDestination(logDest); + } + } + if (!psLogSetDestination(config->logFD)) { + psWarning("Unable to set log destination to file number %d --- ignored", config->logFD); } // Set trace levels psMetadata *trace = psMetadataLookupMetadata(&mdok, config->site, "TRACE"); - if (mdok && trace) - { + if (mdok && trace) { psMetadataIterator *traceIter = psMetadataIteratorAlloc(trace, PS_LIST_HEAD, NULL); // Iterator psMetadataItem *traceItem = NULL; // Item from MD iteration while ((traceItem = psMetadataGetAndIncrement(traceIter))) { if (traceItem->type != PS_DATA_S32) { - psLogMsg("psModules.config", PS_LOG_WARN, - "The level for trace component %s is not of type S32 (%x)\n", + psWarning("The level for trace component %s is not of type S32 (%x)\n", traceItem->name, traceItem->type); continue; } - psTrace("psModules.config", 7, "Setting trace level for %s to %d\n", traceItem->name, traceItem->data.S32); + psTrace("psModules.config", 7, "Setting trace level for %s to %d\n", + traceItem->name, traceItem->data.S32); (void)psTraceSetLevel(traceItem->name, traceItem->data.S32); } @@ -494,6 +504,5 @@ // Set trace formats psString traceFormat = psMetadataLookupStr(&mdok, config->site, "TRACEFORMAT"); - if (mdok && traceFormat) - { + if (mdok && traceFormat) { psTrace("psModules.config", 7, "Setting trace format to %s\n", traceFormat); (void)psTraceSetFormat(traceFormat); @@ -502,14 +511,23 @@ // Set trace destinations #ifndef PS_NO_TRACE - psString traceDest = psMetadataLookupStr(&mdok, config->site, "TRACEDEST"); - if (mdok && traceDest && strlen(traceDest) > 0) - { - psTrace("psModules.config", 7, "Setting trace destination to %s\n", traceDest); - // 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 (!psTraceSetDestination(psMessageDestination(traceDest))) { - psLogMsg("psModules.config", PS_LOG_WARN, "Unable to set trace destination to %s\n", traceDest); - - } + argNum = psArgumentGet(*argc, argv, "-tracedest"); + if (argNum > 0) { + psArgumentRemove(argNum, argc, argv); + if (argNum >= *argc) { + psWarning("-tracedest command-line switch provided without the required trace destination " + "--- ignored.\n"); + } else { + config->traceFD = psMessageDestination(argv[argNum]); + psArgumentRemove(argNum, argc, argv); + } + } else { + psString traceDest = psMetadataLookupStr(&mdok, config->site, "TRACEDEST"); + if (mdok && traceDest && strlen(traceDest) > 0) { + psTrace("psModules.config", 7, "Setting trace destination to %s\n", traceDest); + config->traceFD = psMessageDestination(traceDest); + } + } + if (!psTraceSetDestination(config->traceFD)) { + psWarning("Unable to set log destination to file %d\n", config->traceFD); } #endif @@ -741,11 +759,9 @@ psArgumentRemove(argNum, argc, argv); if (argNum >= *argc) { - psLogMsg("psModules.config", PS_LOG_WARN, - "-dbserver command-line switch provided without the required server name --- "); + psWarning("-dbserver command-line switch provided without the required server name --- "); } else { char *dbserver = argv[argNum]; // The camera configuration file to read if (!psMetadataAddStr(config->site, PS_LIST_TAIL, "DBSERVER", PS_META_REPLACE, NULL, dbserver)) { - psLogMsg("psModules.config", PS_LOG_WARN, - "failed to overwrite .ipprc DBSERVER value --- "); + psWarning("Failed to overwrite .ipprc DBSERVER value"); } @@ -758,11 +774,9 @@ psArgumentRemove(argNum, argc, argv); if (argNum >= *argc) { - psLogMsg("psModules.config", PS_LOG_WARN, - "-dbname command-line switch provided without the required database name --- "); + psWarning("-dbname command-line switch provided without the required database name"); } else { char *dbname = argv[argNum]; // The camera configuration file to read if (!psMetadataAddStr(config->site, PS_LIST_TAIL, "DBNAME", PS_META_REPLACE, NULL, dbname)) { - psLogMsg("psModules.config", PS_LOG_WARN, - "failed to overwrite .ipprc DBNAME value --- "); + psWarning("Failed to overwrite .ipprc DBNAME value"); } @@ -775,11 +789,9 @@ psArgumentRemove(argNum, argc, argv); if (argNum >= *argc) { - psLogMsg("psModules.config", PS_LOG_WARN, - "-dbuser command-line switch provided without the required database name --- "); + psWarning("-dbuser command-line switch provided without the required database name"); } else { char *dbuser = argv[argNum]; // The camera configuration file to read if (!psMetadataAddStr(config->site, PS_LIST_TAIL, "DBUSER", PS_META_REPLACE, NULL, dbuser)) { - psLogMsg("psModules.config", PS_LOG_WARN, - "failed to overwrite .ipprc DBUSER value --- "); + psWarning("Failed to overwrite .ipprc DBUSER value"); } @@ -792,11 +804,10 @@ psArgumentRemove(argNum, argc, argv); if (argNum >= *argc) { - psLogMsg("psModules.config", PS_LOG_WARN, - "-dbpassword command-line switch provided without the required password --- "); + psWarning("-dbpassword command-line switch provided without the required password"); } else { char *dbpassword = argv[argNum]; // The camera configuration file to read - if (!psMetadataAddStr(config->site, PS_LIST_TAIL, "DBPASSWORD", PS_META_REPLACE, NULL, dbpassword)) { - psLogMsg("psModules.config", PS_LOG_WARN, - "failed to overwrite .ipprc DBPASSWORD value --- "); + if (!psMetadataAddStr(config->site, PS_LIST_TAIL, "DBPASSWORD", PS_META_REPLACE, + NULL, dbpassword)) { + psWarning("Failed to overwrite .ipprc DBPASSWORD value"); } @@ -809,11 +820,10 @@ psArgumentRemove(argNum, argc, argv); if (argNum >= *argc) { - psLogMsg("psModules.config", PS_LOG_WARN, - "-dbpport command-line switch provided without the required port number --- "); + psWarning("-dbpport command-line switch provided without the required port number"); } else { char *dbport = argv[argNum]; // The camera configuration file to read - if (!psMetadataAddS32(config->site, PS_LIST_TAIL, "DBPORT", PS_META_REPLACE, NULL, (psS32)atoi(dbport))) { - psLogMsg("psModules.config", PS_LOG_WARN, - "failed to overwrite .ipprc DBPORT value --- "); + if (!psMetadataAddS32(config->site, PS_LIST_TAIL, "DBPORT", PS_META_REPLACE, NULL, + (psS32)atoi(dbport))) { + psWarning("Failed to overwrite .ipprc DBPORT value"); } @@ -963,6 +973,5 @@ result = true; } else { - psLogMsg("psModules.config", PS_LOG_WARN, - "Camera %s, format %s also matches header --- ignored.\n", + psWarning("Camera %s, format %s also matches header --- ignored.\n", cameraName, formatsItem->name); } @@ -1050,4 +1059,10 @@ } + // Now that we have the camera, we can apply recipes based on command-line arguments + if (readRecipes && !pmConfigReadRecipes(config, PM_RECIPE_SOURCE_CL)) { + psError(PS_ERR_IO, false, "Error reading recipes from camera config for %s", config->cameraName); + return NULL; + } + if (!config->format && !config->formatName) { config->formatName = name; @@ -1081,5 +1096,5 @@ if (!pmConfigFileRead(&camera, cameraPath, cameraName)) { - psLogMsg("psModules.config", PS_LOG_WARN, "Trouble reading reading camera configuration %s", cameraName); + psWarning("Trouble reading reading camera configuration %s", cameraName); psFree(camera); return NULL; @@ -1118,6 +1133,5 @@ } if (!(mdStatus01 && mdStatus02 && mdStatus03 && mdStatus04)) { - psLogMsg("psModules.config", PS_LOG_WARN, - "Could not determine database server, name, user, and password from site metadata.\n"); + psWarning("Could not determine database server, name, user, and password from site metadata.\n"); return NULL; } Index: /branches/eam_branch_20071023/psModules/src/config/pmConfig.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/config/pmConfig.h (revision 15398) +++ /branches/eam_branch_20071023/psModules/src/config/pmConfig.h (revision 15399) @@ -5,6 +5,7 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-04 20:15:48 $ + * @version $Revision: 1.31.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-29 01:40:54 $ + * * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ @@ -36,6 +37,5 @@ /// This structure stores the configuration information: the site, camera and recipe configuration, the /// command-line arguments, the pmFPAfiles used, and the database handle. -typedef struct -{ +typedef struct { psMetadata *site; ///< Site configuration psMetadata *camera; ///< Camera specification @@ -51,12 +51,7 @@ pmRecipeSource recipesRead; ///< Which recipe sources have been read psMetadata *recipeSymbols; ///< Where each recipe came from - - // dropping the argc, argv info from pmConfig - # if (0) - int *argc; ///< Number of command-line arguments - char **argv; ///< Command-line arguments (raw version) - # endif -} -pmConfig; + int traceFD; ///< File descriptor for trace messages + int logFD; ///< File descriptor for log messages +} pmConfig; /// Allocator for pmConfig Index: /branches/eam_branch_20071023/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/detrend/pmShutterCorrection.c (revision 15398) +++ /branches/eam_branch_20071023/psModules/src/detrend/pmShutterCorrection.c (revision 15399) @@ -832,11 +832,16 @@ for (int j = 0; j < MEASURE_SAMPLES; j++) { psRegion *region = data->regions->data[j]; // Region of interest - psImage *subImage = psImageSubset(readout->image, *region); // Sub-image + psRegion adjusted = *region; // Adjusted region, compensating for offsets + adjusted.x0 += readout->image->col0; + adjusted.x1 += readout->image->col0; + adjusted.y0 += readout->image->row0; + adjusted.y1 += readout->image->row0; + psImage *subImage = psImageSubset(readout->image, adjusted); // Sub-image psImage *subMask = NULL; // Sub-image of mask if (readout->mask) { - subMask = psImageSubset(readout->mask, *region); + subMask = psImageSubset(readout->mask, adjusted); } if (!psImageStats(stats, subImage, subMask, maskVal)) { - psString regionString = psRegionToString(*region); + psString regionString = psRegionToString(adjusted); psWarning("Unable to measure sample statistics at %s in image.\n", regionString); Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmStack.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmStack.c (revision 15398) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmStack.c (revision 15399) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-23 23:43:12 $ + * @version $Revision: 1.13.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-29 01:40:54 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -221,4 +221,7 @@ for (int i = 0; i < num; i++) { pmStackData *data = inputs->data[i]; // Stack data of interest + if (!data) { + continue; + } psImage *image = data->readout->image; // Image of interest psImage *weight = data->readout->weight; // Weight map of interest @@ -292,4 +295,7 @@ // Add the pixel as one to inspect pmStackData *data = inputs->data[j]; // Stack data of interest + if (!data) { + continue; + } data->pixels = psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y); // Mask it so it's not considered in other iterations within this function @@ -315,6 +321,8 @@ *num = input->n; - // The first is a template - pmStackData *data = input->data[0]; // First image off the rank + pmStackData *data = NULL; // First image off the rank, used as a template + for (int i = 0; !data && i < input->n; i++) { + data = input->data[i]; + } PS_ASSERT_PTR_NON_NULL(data, false); assert(psMemGetDeallocator(data) == (psFreeFunc)stackDataFree); // Ensure it's the right type @@ -338,4 +346,7 @@ for (int i = 1; i < *num; i++) { pmStackData *data = input->data[i]; // Stack data for this input + if (!data) { + continue; + } assert(psMemGetDeallocator(data) == (psFreeFunc)stackDataFree); // Ensure it's the right type if (!data->readout) { @@ -469,5 +480,8 @@ for (int i = 0; i < input->n; i++) { pmStackData *data = input->data[i]; - assert(data && data->pixels); + if (!data) { + continue; + } + assert(data->pixels); psPixels *pixels = data->pixels;// The pixels of interest for (int j = 0; j < pixels->n; j++) { @@ -544,4 +558,7 @@ for (int i = 0; i < num; i++) { pmStackData *data = input->data[i]; // Stack data for this input + if (!data) { + weights->data.F32[i] = 0.0; + } weights->data.F32[i] = data->weight; } @@ -561,4 +578,7 @@ for (int i = 0; i < num; i++) { pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels + if (!data) { + continue; + } pixels = psPixelsConcatenate(pixels, data->pixels); data->pixels = psPixelsRealloc(data->pixels, PIXEL_LIST_BUFFER); // Just in case more rejection @@ -596,4 +616,7 @@ for (int i = 0; i < num; i++) { pmStackData *data = input->data[i]; // Stack data for this input + if (!data) { + continue; + } data->pixels = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); } @@ -611,4 +634,7 @@ for (int i = 0; i < num; i++) { pmStackData *data = input->data[i]; // Stack data for this input + if (!data) { + continue; + } psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->pixels->n); } Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtraction.c (revision 15398) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtraction.c (revision 15399) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.65 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-17 02:45:40 $ + * @version $Revision: 1.65.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-29 01:40:54 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -907,4 +907,7 @@ for (int j = 0; j < numKernels; j++) { psKernel *convolution = convolutions->data[j]; // Convolution of reference + + // XXX Precalculate these values, so that we're not calculating the same thing for + // every pixel. double polynomial = 0.0; // Value of the polynomial for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) { @@ -928,5 +931,5 @@ deviations->data.F32[i] = devNorm * deviation; psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", - i, (int)stamp->x, (int)stamp->y, deviations->data.F32[i]); + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); totalSquareDev += PS_SQR(deviations->data.F32[i]); numStamps++; @@ -947,30 +950,47 @@ } + if (!isfinite(sigmaRej) || sigmaRej <= 0.0) { + psLogMsg("psModules.imcombine", PS_LOG_INFO, "RMS deviation: %f", sqrt(totalSquareDev / numStamps)); + psFree(deviations); + return 0; + } + + int numRejected = 0; // Number of stamps rejected + int numGood = 0; // Number of good stamps + double newSquareDev = 0.0; // New square deviation + float limit = sigmaRej * sqrt(totalSquareDev / numStamps); // Limit on maximum deviation psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit); - int numRejected = 0; for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > limit) { - // Mask out the stamp in the image so you it's not found again - psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, - (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); - numRejected++; - for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { - for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { - subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + if (deviations->data.F32[i] > limit) { + // Mask out the stamp in the image so you it's not found again + psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, + (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); + numRejected++; + for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { + for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { + subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; + } } - } - - // Set stamp for replacement - stamp->x = 0; - stamp->y = 0; - stamp->status = PM_SUBTRACTION_STAMP_REJECTED; - // Recalculate convolutions - psFree(stamp->convolutions); - stamp->convolutions = NULL; - } - } + + // Set stamp for replacement + stamp->x = 0; + stamp->y = 0; + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + // Recalculate convolutions + psFree(stamp->convolutions); + stamp->convolutions = NULL; + } else { + numGood++; + newSquareDev += PS_SQR(deviations->data.F32[i]); + } + } + } + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n", + numGood, numRejected, sqrt(totalSquareDev / numStamps), sqrt(newSquareDev / numGood)); psFree(deviations); Index: /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.c (revision 15398) +++ /branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionKernels.c (revision 15399) @@ -495,5 +495,5 @@ int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters int radius = inner; - while (radius < size) { + while (radius + fibIndex < size) { radius++; int fibNew = fibIndex + fibIndexMinus1; @@ -502,7 +502,5 @@ radius += fibIndex; fibNum++; - printf("%d ", radius); - } - printf("\n"); + } } Index: /branches/eam_branch_20071023/psModules/src/psmodules.h =================================================================== --- /branches/eam_branch_20071023/psModules/src/psmodules.h (revision 15398) +++ /branches/eam_branch_20071023/psModules/src/psmodules.h (revision 15399) @@ -106,3 +106,6 @@ #include +// The following headers are from random locations, here because they cross bounds +#include + #endif