Index: trunk/psModules/src/camera/pmFPA.c =================================================================== --- trunk/psModules/src/camera/pmFPA.c (revision 18829) +++ trunk/psModules/src/camera/pmFPA.c (revision 18830) @@ -176,7 +176,12 @@ readout->col0 = 0; readout->row0 = 0; - readout->imageScan = 0; - readout->maskScan = 0; - readout->weightScan = 0; + + readout->thisImageScan = 0; + readout->thisMaskScan = 0; + readout->thisWeightScan = 0; + + readout->lastImageScan = 0; + readout->lastMaskScan = 0; + readout->lastWeightScan = 0; } @@ -269,7 +274,11 @@ tmpReadout->col0 = 0; - tmpReadout->imageScan = 0; - tmpReadout->maskScan = 0; - tmpReadout->weightScan = 0; + tmpReadout->thisImageScan = 0; + tmpReadout->thisMaskScan = 0; + tmpReadout->thisWeightScan = 0; + + tmpReadout->lastImageScan = 0; + tmpReadout->lastMaskScan = 0; + tmpReadout->lastWeightScan = 0; return(tmpReadout); Index: trunk/psModules/src/camera/pmFPA.h =================================================================== --- trunk/psModules/src/camera/pmFPA.h (revision 18829) +++ trunk/psModules/src/camera/pmFPA.h (revision 18830) @@ -6,6 +6,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-06-05 01:31:33 $ + * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-01 00:01:25 $ * Copyright 2005-2006 Institute for Astronomy, University of Hawaii */ @@ -127,5 +127,10 @@ bool file_exists; ///< Does the file for this readout exist (read case only)? bool data_exists; ///< Does the data for this readout exist (read case only)? - int imageScan, maskScan, weightScan;///< Separate tracking numbers for reading images incrementally + int thisImageScan; ///< start scan for next/current read + int lastImageScan; ///< start scan of the last read + int thisMaskScan; ///< start scan for next/current read + int lastMaskScan; ///< start scan of the last read + int thisWeightScan; ///< start scan for next/current read + int lastWeightScan; ///< start scan of the last read } pmReadout; Index: trunk/psModules/src/camera/pmFPARead.c =================================================================== --- trunk/psModules/src/camera/pmFPARead.c (revision 18829) +++ trunk/psModules/src/camera/pmFPARead.c (revision 18830) @@ -45,18 +45,40 @@ // Return pointer to appropriate value for tracking scans -static int *readoutScansByType(pmReadout *readout, // Readout of interest - fpaReadType type // Type of image +static int readoutGetThisScan(pmReadout *readout, // Readout of interest + fpaReadType type // Type of image ) { switch (type) { case FPA_READ_TYPE_IMAGE: - return &readout->imageScan; + return readout->thisImageScan; case FPA_READ_TYPE_MASK: - return &readout->maskScan; + return readout->thisMaskScan; case FPA_READ_TYPE_WEIGHT: - return &readout->weightScan; + return readout->thisWeightScan; default: psAbort("Unknown read type: %x\n", type); } +} + +// Return pointer to appropriate value for tracking scans +static int readoutSetLastScan(pmReadout *readout, // Readout of interest + fpaReadType type, // Type of image + int numScans // requested number of scans + ) +{ + switch (type) { + case FPA_READ_TYPE_IMAGE: + readout->lastImageScan = readout->thisImageScan + numScans; + return readout->lastImageScan; + case FPA_READ_TYPE_MASK: + readout->lastMaskScan = readout->thisMaskScan + numScans; + return readout->lastMaskScan; + case FPA_READ_TYPE_WEIGHT: + readout->lastWeightScan = readout->thisWeightScan + numScans; + return readout->lastWeightScan; + default: + psAbort("Unknown read type: %x\n", type); + } + return false; } @@ -134,21 +156,18 @@ } - -// Determine readout scan properties: the next and last scans -// Requires that cellNumReadouts() has been called before (for header and concepts to have been read) -// In the process, adjusts the TRIMSEC -static bool readoutScanProperties(int *next, // Index of next scan - int *last, // Index of last scan - pmReadout *readout, // Readout of interest - int numScans, // Number of scans to read at a time - fpaReadType type, // Type of image - pmConfig *config // Configuration +// Does the current readout, with scans set for a new read, represent any real data, or is it +// beyond the end? Requires that cellNumReadouts() has been called before (for header and +// concepts to have been read) In the process, adjusts the TRIMSEC +static bool readoutHaveMoreScans(bool *result, // true : more data to read + pmReadout *readout, // Readout of interest + int numScans, // Number of scans to read at a time + fpaReadType type, // Type of image + pmConfig *config // Configuration ) { - assert(next); - assert(last); + assert(result); assert(readout); - psImage *image = *readoutImageByType(readout, type); // Appropriate image from readout + *result = false; if (!pmConceptsReadCell(readout->parent, PM_CONCEPT_SOURCE_DEFAULTS | PM_CONCEPT_SOURCE_DATABASE, @@ -163,4 +182,5 @@ PS_ASSERT_PTR_NON_NULL(cell, false); pmHDU *hdu = pmHDUFromCell(cell); // HDU for data + bool mdok = true; // Status of MD lookup psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim sections @@ -195,11 +215,22 @@ // Calculate the segment offset and upper limit if (numScans == 0) { - *next = (readdir == 1) ? trimsec->y0 : trimsec->x0; - } else { - *next = image ? *readoutScansByType(readout, type) + numScans : 0; - } - *last = (readdir == 1) ? trimsec->y1 : trimsec->x1; - + // in the case of numScans == 0, we never call this funtion unless the data has not yet + // been read. thus, only if the delta is should we return false (ie, trimsec defines an empty region) + int thisScan = (readdir == 1) ? trimsec->y0 : trimsec->x0; + int lastScan = (readdir == 1) ? trimsec->y1 : trimsec->x1; + *result = (lastScan > thisScan); + return true; + } + + // allow multiple threads to read different segments into different readouts + // this thread may not have read a segment yet; but others have. trust readout->imageScan + + // is the start scan of the read less than the last possible scan? + int thisScan = readoutGetThisScan(readout, type); + int lastScan = (readdir == 1) ? trimsec->y1 : trimsec->x1; + + *result = (lastScan > thisScan); return true; + // XXX fold this and the above case together } @@ -217,8 +248,13 @@ psImage *image = *readoutImageByType(readout, type); + // XXX this may not be the valid test in a multithread environment. consider a fileGroup of + // N readouts, but numScans set to 0. only the first should report that it requires data, + // even if all readouts lack the image pointer. if (!image) { return true; - } else if (numScans == 0) { - // Can only read the entire image once + } + + // If we have already read an image, this result implies we are done (no more to read) + if (numScans == 0) { return false; } @@ -235,12 +271,11 @@ } - int next; // Next position - int last; // Last position - if (!readoutScanProperties(&next, &last, readout, numScans, type, config)) { + bool haveData; + if (!readoutHaveMoreScans(&haveData, readout, numScans, type, config)) { psError(PS_ERR_UNKNOWN, false, "Unable to determine readout properties."); return false; } - return (next < last); + return haveData; } @@ -423,14 +458,14 @@ } - int next; // Next position - int last; // Last position - if (!readoutScanProperties(&next, &last, readout, numScans, type, config)) { + bool haveData; + if (!readoutHaveMoreScans (&haveData, readout, type, numScans, config)) { psError(PS_ERR_UNKNOWN, false, "Unable to determine readout properties."); return false; } - if (next >= last) { + if (!haveData) { psError(PS_ERR_IO, true, "No more of the readout to read."); return false; } + pmHDU *hdu = pmHDUFromCell(cell); // The HDU @@ -459,6 +494,4 @@ } - - // Check the third dimension int naxis = psMetadataLookupS32(&mdok, hdu->header, "NAXIS"); // The number of axes @@ -478,11 +511,11 @@ } - // Update number of scans read - int *scansTracker = readoutScansByType(readout, type); // Tracker for how many scans have been read - *scansTracker = next; - if (next == 0) { + // Determine the numbe of scans to read + int thisScan = readoutGetThisScan(readout, type); + int lastScan = readoutSetLastScan (readout, type, numScans); + if (thisScan == 0) { overlap = 0; } - next -= overlap; + thisScan -= overlap; // Calculate limits, adjust readout->row0,col0 @@ -491,5 +524,5 @@ if (readdir == 1) { // Reading rows - readout->row0 = next; + readout->row0 = thisScan; readout->col0 = trimsec->x0; if (numScans == 0) { @@ -498,5 +531,5 @@ } else { // Reading cols - readout->col0 = next; + readout->col0 = thisScan; readout->row0 = trimsec->y0; if (numScans == 0) { @@ -504,5 +537,4 @@ } } - int upper = next + numScans + overlap; // Upper limit to next section // Blow away existing data. @@ -510,5 +542,5 @@ psFree(*image); *image = NULL; - *image = readoutReadComponent(*image, fits, trimsec, readdir, next, upper, z, bad, pixelTypes[type]); + *image = readoutReadComponent(*image, fits, trimsec, readdir, thisScan, lastScan, z, bad, pixelTypes[type]); // Read overscans only for "image" type --- weights and masks shouldn't record overscans @@ -528,6 +560,5 @@ psRegion *biassec = NULL; // Bias section from iteration while ((biassec = psListGetAndIncrement(biassecsIter))) { - psImage *bias = readoutReadComponent(NULL, fits, biassec, readdir, next, upper, z, - bad, pixelTypes[type]); // The bias + psImage *bias = readoutReadComponent(NULL, fits, biassec, readdir, thisScan, lastScan, z, bad, pixelTypes[type]); // The bias psListAdd(readout->bias, PS_LIST_TAIL, bias); psFree(bias); // Drop reference Index: trunk/psModules/src/camera/pmReadoutStack.c =================================================================== --- trunk/psModules/src/camera/pmReadoutStack.c (revision 18829) +++ trunk/psModules/src/camera/pmReadoutStack.c (revision 18830) @@ -8,4 +8,41 @@ #include "pmReadoutStack.h" + +// generate the specified image +// XXX should it be an error for the image to exist? +psImage *pmReadoutSetAnalysisImage(pmReadout *readout, // Readout containing image + const char *name, // Name of image in analysis metadata + int numCols, int numRows, // Expected size of image + psElemType type, // Expected type of image + double init // Initial value + ) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_STRING_NON_EMPTY(name, false); + + psImage *image = psImageAlloc(numCols, numRows, type); + + if (!psMetadataAddImage(readout->analysis, PS_LIST_TAIL, name, 0, "Analysis image from " __FILE__, image)) { + psAbort ("analysis image already exists"); + } + psImageInit(image, init); + + psFree (image); // we still have a view on readout->analysis + return image; +} + +// retrieve the specified image +// XXX not sure why this should call psMemIncrRefCounter +psImage *pmReadoutGetAnalysisImage(pmReadout *readout, // Readout containing image + const char *name // Name of image in analysis metadata + ) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_STRING_NON_EMPTY(name, false); + + bool mdok; // Status of MD lookup + psImage *image = psMetadataLookupPtr(&mdok, readout->analysis, name); + return image; +} psImage *pmReadoutAnalysisImage(pmReadout *readout, // Readout containing image @@ -40,58 +77,30 @@ } -bool pmReadoutUpdateSize(pmReadout *readout, int minCols, int minRows, - int numCols, int numRows, bool mask, bool weight, - psMaskType blank) -{ - PS_ASSERT_PTR_NON_NULL(readout, false); - - if (readout->image) { - readout->col0 = PS_MIN(minCols, readout->col0); - readout->row0 = PS_MIN(minRows, readout->row0); - } else { - readout->col0 = minCols; - readout->row0 = minRows; - } - - // (reAllocate the images - if (!readout->image) { - readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - psImageInit(readout->image, NAN); - } - if (readout->image->numCols < numCols || readout->image->numRows < numRows) { - // Generate the new output image by extending the current one, or making a whole new one - psImage *newImage = psImageAlloc(numCols, numRows, PS_TYPE_F32); - psImageInit(newImage, NAN); - psImageOverlaySection(newImage, readout->image, readout->col0, readout->row0, "="); - psFree(readout->image); - readout->image = newImage; - } +// XXX for the moment, use col0, row0, numCols, numRows supplied from the outside +bool pmReadoutStackDefineOutput(pmReadout *readout, int col0, int row0, int numCols, int numRows, bool mask, bool weight, psMaskType blank) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + + // XXX is this an error? + if (readout->image) return false; + readout->col0 = col0; + readout->row0 = row0; + + // allocate the images + readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->image, NAN); if (mask) { - if (!readout->mask) { - readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); - psImageInit(readout->mask, blank); - } - if (readout->mask->numCols < numCols || readout->mask->numRows < numRows) { - psImage *newMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); - psImageInit(newMask, blank); - psImageOverlaySection(newMask, readout->mask, readout->col0, readout->row0, "="); - psFree(readout->mask); - readout->mask = newMask; - } + // XXX is this an error? + if (readout->mask) return false; + readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + psImageInit(readout->mask, blank); } if (weight) { - if (!readout->weight) { - readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - psImageInit(readout->weight, NAN); - } - if (readout->weight->numCols < numCols || readout->weight->numRows < numRows) { - psImage *newWeight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - psImageInit(newWeight, NAN); - psImageOverlaySection(newWeight, readout->weight, readout->col0, readout->row0, "="); - psFree(readout->weight); - readout->weight = newWeight; - } + // XXX is this an error? + if (readout->weight) return false; + readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->weight, NAN); } @@ -99,23 +108,14 @@ } -bool pmReadoutStackValidate(int *minInputColsPtr, int *maxInputColsPtr, int *minInputRowsPtr, - int *maxInputRowsPtr, int *numColsPtr, int *numRowsPtr, - const psArray *inputs) +bool pmReadoutStackSetOutputSize(int *col0, int *row0, int *numCols, int *numRows, const psArray *inputs) { PS_ASSERT_ARRAY_NON_NULL(inputs, false); - PS_ASSERT_PTR_NON_NULL(minInputColsPtr, false); - PS_ASSERT_PTR_NON_NULL(maxInputColsPtr, false); - PS_ASSERT_PTR_NON_NULL(minInputRowsPtr, false); - PS_ASSERT_PTR_NON_NULL(maxInputRowsPtr, false); - PS_ASSERT_PTR_NON_NULL(numColsPtr, false); - PS_ASSERT_PTR_NON_NULL(numRowsPtr, false); - - // Step through each readout in the input image list to determine how big of an output image is needed to - // combine these input images. - int maxInputCols = 0; // The largest input column value - int maxInputRows = 0; // The largest input row value - int minInputCols = INT_MAX; // The smallest input column value - int minInputRows = INT_MAX; // The smallest input row value - int xSize = 0, ySize = 0; // The size of the output image + PS_ASSERT_PTR_NON_NULL(col0, false); + PS_ASSERT_PTR_NON_NULL(row0, false); + PS_ASSERT_PTR_NON_NULL(numCols, false); + PS_ASSERT_PTR_NON_NULL(numRows, false); + + // Step through each readout in the input image list to determine how big of an output + // image is needed to combine these input images. int xMin = INT_MAX; @@ -123,4 +123,6 @@ int xMax = 0; int yMax = 0; + int xSize = 0; + int ySize = 0; // The size of the output image bool valid = false; // Do we have a single valid input? @@ -128,11 +130,5 @@ pmReadout *readout = inputs->data[i]; // Readout of interest - if (!readout) { - continue; - } - if (!readout->image) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Input readout %ld has NULL image.\n", i); - return false; - } + if (!readout) continue; // use the trimsec to define the max full range of the output pixels @@ -150,4 +146,129 @@ yMax = PS_MAX(yMax, trimsec->y1); } + valid = true; + psTrace("psModules.camera", 7, "Readout %ld: trimsec: %f,%f - %f,%f\n", i, trimsec->x0, trimsec->y0, trimsec->x1, trimsec->y1); + } + + *col0 = xMin; + *row0 = yMin; + *numCols = xSize; + *numRows = ySize; + + if (!valid) { + psError(PS_ERR_UNKNOWN, false, "No valid input readouts."); + } + return valid; +} + +bool pmReadoutUpdateSize(pmReadout *readout, int minCols, int minRows, + int numCols, int numRows, bool mask, bool weight, + psMaskType blank) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + + if (readout->image) { + readout->col0 = PS_MIN(minCols, readout->col0); + readout->row0 = PS_MIN(minRows, readout->row0); + } else { + readout->col0 = minCols; + readout->row0 = minRows; + } + + // (reAllocate the images + if (!readout->image) { + readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->image, NAN); + } + if (readout->image->numCols < numCols || readout->image->numRows < numRows) { + // Generate the new output image by extending the current one, or making a whole new one + psImage *newImage = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(newImage, NAN); + psImageOverlaySection(newImage, readout->image, readout->col0, readout->row0, "="); + psFree(readout->image); + readout->image = newImage; + } + + if (mask) { + if (!readout->mask) { + readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + psImageInit(readout->mask, blank); + } + if (readout->mask->numCols < numCols || readout->mask->numRows < numRows) { + psImage *newMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + psImageInit(newMask, blank); + psImageOverlaySection(newMask, readout->mask, readout->col0, readout->row0, "="); + psFree(readout->mask); + readout->mask = newMask; + } + } + + if (weight) { + if (!readout->weight) { + readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->weight, NAN); + } + if (readout->weight->numCols < numCols || readout->weight->numRows < numRows) { + psImage *newWeight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(newWeight, NAN); + psImageOverlaySection(newWeight, readout->weight, readout->col0, readout->row0, "="); + psFree(readout->weight); + readout->weight = newWeight; + } + } + + return true; +} + +bool pmReadoutStackValidate(int *minInputColsPtr, int *maxInputColsPtr, int *minInputRowsPtr, + int *maxInputRowsPtr, int *numColsPtr, int *numRowsPtr, + const psArray *inputs) +{ + PS_ASSERT_ARRAY_NON_NULL(inputs, false); + PS_ASSERT_PTR_NON_NULL(minInputColsPtr, false); + PS_ASSERT_PTR_NON_NULL(maxInputColsPtr, false); + PS_ASSERT_PTR_NON_NULL(minInputRowsPtr, false); + PS_ASSERT_PTR_NON_NULL(maxInputRowsPtr, false); + PS_ASSERT_PTR_NON_NULL(numColsPtr, false); + PS_ASSERT_PTR_NON_NULL(numRowsPtr, false); + + // Step through each readout in the input image list to determine how big of an output image is needed to + // combine these input images. + int maxInputCols = 0; // The largest input column value + int maxInputRows = 0; // The largest input row value + int minInputCols = INT_MAX; // The smallest input column value + int minInputRows = INT_MAX; // The smallest input row value + int xSize = 0, ySize = 0; // The size of the output image + + int xMin = INT_MAX; + int yMin = INT_MAX; + int xMax = 0; + int yMax = 0; + + bool valid = false; // Do we have a single valid input? + for (long i = 0; i < inputs->n; i++) { + pmReadout *readout = inputs->data[i]; // Readout of interest + + if (!readout) { + continue; + } + if (!readout->image) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Input readout %ld has NULL image.\n", i); + return false; + } + + // use the trimsec to define the max full range of the output pixels + pmCell *cell = readout->parent; // The parent cell + bool mdok = true; // Status of MD lookup + psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section + if (!mdok || !trimsec || psRegionIsNaN(*trimsec)) { + psWarning("CELL.TRIMSEC is not set for readout %ld --- ignored.\n", i); + } else { + xSize = PS_MAX(xSize, trimsec->x1 - trimsec->x0); + ySize = PS_MAX(ySize, trimsec->y1 - trimsec->y0); + xMin = PS_MIN(xMin, trimsec->x0); + xMax = PS_MAX(xMax, trimsec->x1); + yMin = PS_MIN(yMin, trimsec->y0); + yMax = PS_MAX(yMax, trimsec->y1); + } valid = true; @@ -158,6 +279,6 @@ minInputRows = PS_MAX(yMin, PS_MIN(minInputRows, readout->row0)); maxInputRows = PS_MIN(yMax, PS_MAX(maxInputRows, readout->row0 + readout->image->numRows)); - psTrace("psModules.camera", 7, "Readout %ld: offset %d,%d; size %dx%d\n", i, - readout->col0, readout->row0, readout->image->numCols, readout->image->numRows); + + psTrace("psModules.camera", 7, "Readout %ld: offset %d,%d; size %dx%d\n", i, readout->col0, readout->row0, readout->image->numCols, readout->image->numRows); } Index: trunk/psModules/src/camera/pmReadoutStack.h =================================================================== --- trunk/psModules/src/camera/pmReadoutStack.h (revision 18829) +++ trunk/psModules/src/camera/pmReadoutStack.h (revision 18830) @@ -24,4 +24,18 @@ ); +psImage *pmReadoutSetAnalysisImage(pmReadout *readout, // Readout containing image + const char *name, // Name of image in analysis metadata + int numCols, int numRows, // Expected size of image + psElemType type, // Expected type of image + double init // Initial value + ); + +// retrieve the specified image +// XXX not sure why this should call psMemIncrRefCounter +psImage *pmReadoutGetAnalysisImage(pmReadout *readout, // Readout containing image + const char *name // Name of image in analysis metadata + ); + + /// Return an image from analysis metadata, produced while stacking psImage *pmReadoutAnalysisImage(pmReadout *readout, // Readout containing image @@ -32,3 +46,8 @@ ); +// XXX for the moment, use col0, row0, numCols, numRows supplied from the outside +bool pmReadoutStackDefineOutput(pmReadout *readout, int col0, int row0, int numCols, int numRows, bool mask, bool weight, psMaskType blank); + +bool pmReadoutStackSetOutputSize(int *col0, int *row0, int *numCols, int *numRows, const psArray *inputs); + #endif Index: trunk/psModules/src/config/pmConfig.c =================================================================== --- trunk/psModules/src/config/pmConfig.c (revision 18829) +++ trunk/psModules/src/config/pmConfig.c (revision 18830) @@ -1627,7 +1627,6 @@ if (trunc) { - if(truncate(filename, 0) != 0) { - psError(PS_ERR_IO, true, "Failed to truncate file, %s\n", filename); - nebServerFree(server); + if(truncate(path, 0) != 0) { + psError(PS_ERR_IO, true, "Failed to truncate Nebulous file %s (real name %s)\n", filename, path); return NULL; } Index: trunk/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 18829) +++ trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 18830) @@ -42,4 +42,122 @@ } +// check the input parameters and set up the output images +bool pmReadoutCombinePrepare(pmReadout *output, const psArray *inputs, const pmCombineParams *params) +{ + // Check inputs + PS_ASSERT_PTR_NON_NULL(output, false); + PS_ASSERT_ARRAY_NON_NULL(inputs, false); + PS_ASSERT_PTR_NON_NULL(params, false); + PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracLow, 0.0, 1.0, false); + PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracHigh, 0.0, 1.0, false); + + // valid combintion statistic? + bool valid = false; + valid |= (params->combine == PS_STAT_SAMPLE_MEAN); + valid |= (params->combine == PS_STAT_SAMPLE_MEDIAN); + valid |= (params->combine == PS_STAT_ROBUST_MEDIAN); + valid |= (params->combine == PS_STAT_FITTED_MEAN); + valid |= (params->combine == PS_STAT_CLIPPED_MEAN); + if (!valid) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Combination method is not SAMPLE_MEAN, SAMPLE_MEDIAN, " + "ROBUST_MEDIAN, FITTED_MEAN or CLIPPED_MEAN.\n"); + return false; + } + + // weights exist if weights desired? + for (int i = 0; i < inputs->n; i++) { + pmReadout *readout = inputs->data[i]; // Readout of interest + if (params->weights && !readout->weight) { + psError(PS_ERR_UNEXPECTED_NULL, true, + "Rejection based on weights requested, but no weights supplied for image %d.\n", i); + return false; + } + } + + pmHDU *hdu = pmHDUFromReadout(output); // Output HDU + if (!hdu) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find HDU for readout.\n"); + return false; + } + + // set the output header metadata + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Combining using statistic: %x", params->combine); + if (!hdu->header) { + hdu->header = psMetadataAlloc(); + } + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + + // note the clipping parameters, if used + if (params->combine == PS_STAT_CLIPPED_MEAN) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", params->iter, params->rej); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + // note the use of weights + if (params->weights) { + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, + "Using input weights to combine images", ""); + } + + // note the rejection fraction + float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep + if (keepFrac != 1.0) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Min/max rejection: %f high, %f low, keep %d", + params->fracHigh, params->fracLow, params->nKeep); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + // note the mask value actually used + psMaskType maskVal = params->maskVal; // The mask value + if (maskVal) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Mask for combination: %x", maskVal); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + // determine the output image size based on the input images + int row0, col0, numCols, numRows; + if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, inputs)) { + psError(PS_ERR_UNKNOWN, false, "problem setting output readout size."); + return false; + } + + // generate the required output images based on the specified sizes + pmReadoutStackDefineOutput(output, col0, row0, numCols, numRows, true, params->weights, params->blank); + psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); + + // these calls allocate and save the requested images on the output analysis metadata + psImage *counts = pmReadoutSetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT, numCols, numRows, PS_TYPE_U16, 0); + if (!counts) { + return false; + } + psImage *sigma = pmReadoutSetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA, numCols, numRows, PS_TYPE_F32, NAN); + if (!sigma) { + return false; + } + + // Update the "concepts" + psList *inputCells = psListAlloc(NULL); // List of cells + for (long i = 0; i < inputs->n; i++) { + pmReadout *readout = inputs->data[i]; // Readout of interest + psListAdd(inputCells, PS_LIST_TAIL, readout->parent); + } + bool success = pmConceptsAverageCells(output->parent, inputCells, NULL, NULL, true); + psFree(inputCells); + + // set these even though the values are not yet set + output->data_exists = true; + output->parent->data_exists = true; + output->parent->parent->data_exists = true; + + return success; +} // XXX: Maybe add support for S16 and S32 types. Currently, only F32 supported. @@ -61,21 +179,4 @@ PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracLow, 0.0, 1.0, false); PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracHigh, 0.0, 1.0, false); - if (params->combine != PS_STAT_SAMPLE_MEAN && params->combine != PS_STAT_SAMPLE_MEDIAN && - params->combine != PS_STAT_ROBUST_MEDIAN && params->combine != PS_STAT_FITTED_MEAN && - params->combine != PS_STAT_CLIPPED_MEAN) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Combination method is not SAMPLE_MEAN, SAMPLE_MEDIAN, " - "ROBUST_MEDIAN, FITTED_MEAN or CLIPPED_MEAN.\n"); - return false; - } - for (int i = 0; i < inputs->n; i++) { - pmReadout *readout = inputs->data[i]; // Readout of interest - if (params->weights && !readout->weight) { - psError(PS_ERR_UNEXPECTED_NULL, true, - "Rejection based on weights requested, but no weights supplied for image %d.\n", i); - return false; - } - } - - bool first = !output->image; // First pass through? pmHDU *hdu = pmHDUFromReadout(output); // Output HDU @@ -85,13 +186,8 @@ } - if (first) { - psString comment = NULL; // Comment to add to header - psStringAppend(&comment, "Combining using statistic: %x", params->combine); - if (!hdu->header) { - hdu->header = psMetadataAlloc(); - } - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - } + pthread_t id = pthread_self(); + char name[64]; + sprintf (name, "%x", (unsigned int) id); + psTimerStart (name); psStatsOptions combineStdev = 0; // Statistics option for weights @@ -118,43 +214,23 @@ stats->clipSigma = params->rej; stats->clipIter = params->iter; - - if (first) { - psString comment = NULL; // Comment to add to header - psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", - params->iter, params->rej); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - } - } - if (params->weights && first) { - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, - "Using input weights to combine images", ""); - } + } + + psImage *counts = pmReadoutGetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT); + if (!counts) { + return false; + } + psImage *sigma = pmReadoutGetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA); + if (!sigma) { + return false; + } + + stats->options |= combineStdev; int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine int xSize, ySize; // Size of the output image - if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, - inputs)) { + if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, inputs)) { psError(PS_ERR_UNKNOWN, false, "No valid input readouts."); return false; } - - pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true, params->weights, - params->blank); - psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); - - psImage *counts = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT, xSize, ySize, - PS_TYPE_U16, 0); - if (!counts) { - return false; - } - psImage *sigma = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA, xSize, ySize, - PS_TYPE_F32, NAN); - if (!sigma) { - psFree(counts); - return false; - } - - stats->options |= combineStdev; // We loop through each pixel in the output image. We loop through each input readout. We determine if @@ -180,19 +256,5 @@ float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep - if (keepFrac != 1.0 && first) { - psString comment = NULL; // Comment to add to header - psStringAppend(&comment, "Min/max rejection: %f high, %f low, keep %d", - params->fracHigh, params->fracLow, params->nKeep); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - } - psMaskType maskVal = params->maskVal; // The mask value - if (maskVal && first) { - psString comment = NULL; // Comment to add to header - psStringAppend(&comment, "Mask for combination: %x", maskVal); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - } #ifndef PS_NO_TRACE @@ -226,5 +288,5 @@ int yOut = i - output->row0; // y position on output readout #ifdef SHOW_BUSY - + if (psTraceGetLevel("psModules.imcombine") > 9) { printf("Processing row %d\r", i); @@ -242,11 +304,4 @@ int xIn = j - readout->col0; // x position on input readout psImage *image = readout->image; // The readout image - - #if 0 // This should have been taken care of already: - // Check bounds - if (xIn < 0 || xIn >= image->numCols || yIn < 0 || yIn >= image->numRows) { - continue; - } - #endif pixelsData[r] = image->data.F32[yIn][xIn]; @@ -348,4 +403,5 @@ } #endif + psFree(index); psFree(pixels); @@ -355,21 +411,7 @@ psFree(stats); psFree(invScale); - psFree(counts); - psFree(sigma); - - // Update the "concepts" - psList *inputCells = psListAlloc(NULL); // List of cells - for (long i = 0; i < inputs->n; i++) { - pmReadout *readout = inputs->data[i]; // Readout of interest - psListAdd(inputCells, PS_LIST_TAIL, readout->parent); - } - bool success = pmConceptsAverageCells(output->parent, inputCells, NULL, NULL, true); - psFree(inputCells); - - output->data_exists = true; - output->parent->data_exists = true; - output->parent->parent->data_exists = true; - - return success; + + // fprintf (stderr, "done with combine %x : %f sec\n", (unsigned int) id, psTimerMark (name)); + return true; } Index: trunk/psModules/src/imcombine/pmReadoutCombine.h =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.h (revision 18829) +++ trunk/psModules/src/imcombine/pmReadoutCombine.h (revision 18830) @@ -5,6 +5,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-29 03:10:17 $ + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-01 00:01:26 $ * Copyright 2004-2006 Institute for Astronomy, University of Hawaii */ @@ -36,4 +36,7 @@ ); +// check the input parameters and set up the output images +bool pmReadoutCombinePrepare(pmReadout *output, const psArray *inputs, const pmCombineParams *params); + /// Combine multiple readouts, applying zero and scale, with optional minmax clipping bool pmReadoutCombine(pmReadout *output,///< Output readout; altered and returned Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 18829) +++ trunk/psModules/src/objects/Makefile.am (revision 18830) @@ -4,4 +4,13 @@ libpsmodulesobjects_la_LDFLAGS = -release $(PACKAGE_VERSION) libpsmodulesobjects_la_SOURCES = \ + pmDetections.c \ + pmSpan.c \ + pmFootprint.c \ + pmFootprintArrayGrow.c \ + pmFootprintArraysMerge.c \ + pmFootprintAssignPeaks.c \ + pmFootprintFind.c \ + pmFootprintFindAtPoint.c \ + pmFootprintIDs.c \ pmPeaks.c \ pmMoments.c \ @@ -45,4 +54,7 @@ pkginclude_HEADERS = \ + pmDetections.h \ + pmSpan.h \ + pmFootprint.h \ pmPeaks.h \ pmMoments.h \