Index: /branches/pap_branch_080214/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- /branches/pap_branch_080214/psModules/src/imcombine/pmReadoutCombine.c (revision 16599) +++ /branches/pap_branch_080214/psModules/src/imcombine/pmReadoutCombine.c (revision 16600) @@ -113,21 +113,10 @@ } - pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true); + pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true, params->weights, + params->blank); psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); psStatsOptions combineStdev = 0; // Statistics option for weights if (params->weights) { - - if (!output->weight) { - output->weight = psImageAlloc(xSize, ySize, PS_TYPE_F32); - } - if (output->weight->numCols < xSize || output->weight->numRows < ySize) { - psImage *newWeight = psImageAlloc(xSize, ySize, PS_TYPE_F32); - psImageInit(newWeight, 0.0); - psImageOverlaySection(newWeight, output->weight, output->col0, output->row0, "="); - psFree(output->weight); - output->weight = newWeight; - } - if (first) { psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, Index: /branches/pap_branch_080214/psModules/src/imcombine/pmStack.c =================================================================== --- /branches/pap_branch_080214/psModules/src/imcombine/pmStack.c (revision 16599) +++ /branches/pap_branch_080214/psModules/src/imcombine/pmStack.c (revision 16600) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-02-13 02:55:33 $ + * @version $Revision: 1.16.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-22 19:18:11 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -23,4 +23,5 @@ #include "pmHDU.h" #include "pmFPA.h" +#include "pmReadoutStack.h" #include "pmConceptsAverage.h" @@ -42,4 +43,8 @@ psVector *weights; // Pixel weights psVector *sort; // Buffer for sorting (to get a robust estimator of the standard dev) +#if 0 + int x0, y0; // Offset from original image to combination region + int nx, ny; // Number of pixels to combine +#endif } combineBuffer; @@ -64,4 +69,11 @@ buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32); buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32); + +#if 0 + buffer->x0 = 0; + buffer->y0 = 0; + buffer->nx = 0; + buffer->ny = 0; +#endif return buffer; @@ -136,7 +148,20 @@ return false; } - *median = num % 2 ? (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0 : - sortBuffer->data.F32[num / 2]; - *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]); + *median = num % 2 ? sortBuffer->data.F32[num / 2] : + (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0 ; +#if 0 + if (num < NUM_DIRECT_STDEV) { +#endif + // If there are not many values, the direct standard deviation is better + double sum = 0.0; + for (int i = 0; i < num; i++) { + sum += PS_SQR(sortBuffer->data.F32[i] - *median); + } + *stdev = sqrt(sum / (float)(num - 1)); +#if 0 + } else { + *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]); + } +#endif return true; @@ -186,5 +211,5 @@ // Given a stack of images, combine with optional rejection. -// Pixels in the stack that are rejected are marked for subsequent +// Pixels in the stack that are rejected are marked for subsequent inspection static bool combinePixels(psImage *image, // Combined image, for output psImage *mask, // Combined mask, for output @@ -193,5 +218,5 @@ const psVector *weights, // Global (single value) weights for data, or NULL const psVector *reject, // Indices of pixels to reject, or NULL - int x, int y, // Coordinates of interest + int x, int y, // Coordinates of interest; frame of output image psMaskType maskVal, // Value to mask psMaskType bad, // Value to give bad pixels @@ -223,12 +248,13 @@ continue; } + int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout psImage *image = data->readout->image; // Image of interest psImage *weight = data->readout->weight; // Weight map of interest psImage *mask = data->readout->mask; // Mask of interest - pixelData->data.F32[i] = image->data.F32[y][x]; + pixelData->data.F32[i] = image->data.F32[yIn][xIn]; if (weight) { - pixelWeights->data.F32[i] = weight->data.F32[y][x]; - } - pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[y][x]; + pixelWeights->data.F32[i] = weight->data.F32[yIn][xIn]; + } + pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[yIn][xIn]; if (pixelMasks->data.PS_TYPE_MASK_DATA[i] & maskVal) { numBad++; @@ -567,4 +593,24 @@ } + // Get the sizes + psArray *stack = psArrayAlloc(input->n); + for (int i = 0; i < stack->n; i++) { + pmStackData *data = input->data[i]; // Stack data + stack->data[i] = psMemIncrRefCounter(data->readout); + } + 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, + stack)) { + psError(PS_ERR_UNKNOWN, false, "Input stack is not valid."); + psFree(stack); + return false; + } + psFree(stack); + pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, bad); + psTrace("psModules.imcombine", 1, "Combining [%d:%d,%d:%d] (%dx%d)\n", + minInputCols, maxInputCols, minInputRows, maxInputRows, xSize, ySize); + + // Buffer for combination combineBuffer *buffer = combineBufferAlloc(num, numIter == 0 ? PS_STAT_SAMPLE_MEAN : @@ -577,5 +623,5 @@ psImage *combinedWeight = combined->weight; // Combined mask - psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source + psArray *pixelMap = pixelMapGenerate(input, maxInputCols, maxInputRows); // Map of pixels to source psPixels *pixels = NULL; // Total list of pixels, with no duplicates for (int i = 0; i < num; i++) { @@ -588,5 +634,9 @@ } for (int i = 0; i < pixels->n; i++) { + // Pixel coordinates are in the frame of the original image int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest + if (x < minInputCols || x >= maxInputCols || y < minInputRows || y >= maxInputRows) { + continue; + } psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely combinePixels(combinedImage, combinedMask, combinedWeight, input, weights, reject, x, y, @@ -626,6 +676,6 @@ } - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { + for (int y = minInputRows; y < maxInputRows; y++) { + for (int x = minInputCols; x < maxInputCols; x++) { combinePixels(combinedImage, combinedMask, combinedWeights, input, weights, NULL, x, y, maskVal, bad, numIter, rej, buffer); Index: /branches/pap_branch_080214/psModules/src/imcombine/pmStackReject.c =================================================================== --- /branches/pap_branch_080214/psModules/src/imcombine/pmStackReject.c (revision 16599) +++ /branches/pap_branch_080214/psModules/src/imcombine/pmStackReject.c (revision 16600) @@ -24,5 +24,5 @@ int numRegions = subRegions->n; // Number of regions int numCols = 0, numRows = 0; // Size of original image - int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image --- should be 0,0 + int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image int size = 0; // Size of kernel for (int i = 0; i < numRegions; i++) { @@ -58,4 +58,5 @@ numRows = PS_MIN(valid->y1, numRows); } + psTrace("psModules.imcombine", 1, "Rejecting [%d:%d,%d:%d]\n", minCols, numCols, minRows, numRows); psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1), @@ -78,5 +79,5 @@ } pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest - if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, true)) { + if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, false, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i); psFree(convRO); @@ -127,8 +128,9 @@ // Threshold the convolved image psPixels *bad = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); // List of pixels that should be masked - for (int y = 0; y < convolved->numRows; y++) { - for (int x = 0; x < convolved->numCols; x++) { + for (int y = size; y < convolved->numRows - size; y++) { + for (int x = size; x < convolved->numCols - size; x++) { if (convolved->data.F32[y][x] > threshold) { - bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x, y); + // Pixel coordinates in "bad" correspond to the full image + bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x + minCols, y + minRows); } } @@ -137,13 +139,13 @@ // Now, grow the mask to include everything that touches a bad pixel in the convolution - mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0xff); - assert(mask->numCols == numCols && mask->numRows == numRows); + int x0 = minCols, y0 = minRows; // Offset for mask image + mask = psPixelsToMask(NULL, bad, psRegionSet(x0, numCols - 1, y0, numRows - 1), 0xff); for (int i = 0; i < bad->n; i++) { - int xPix = bad->data[i].x, yPix = bad->data[i].y; // Coordinates of interest + int xPix = bad->data[i].x - x0, yPix = bad->data[i].y - y0; // Coordinates in frame of mask image // Convolution limits int xMin = PS_MAX(xPix - size, 0); - int xMax = PS_MIN(xPix + size, numCols - 1); + int xMax = PS_MIN(xPix + size, mask->numCols - 1); int yMin = PS_MAX(yPix - size, 0); - int yMax = PS_MIN(yPix + size, numRows - 1); + int yMax = PS_MIN(yPix + size, mask->numRows - 1); for (int y = yMin; y <= yMax; y++) { for (int x = xMin; x <= xMax; x++) { @@ -153,8 +155,20 @@ } } - bad = psPixelsFromMask(bad, mask, 0xff); psFree(mask); + // Convert coordinates to frame of original image + for (int i = 0; i < bad->n; i++) { + int x = bad->data[i].x + x0; + int y = bad->data[i].y + y0; + if (x < 0 || x >= numCols || y < 0 || y >= numRows) { + psWarning("Bad pixel coordinate %d: %d,%d --- ignored.", + i, x, y); + continue; + } + bad->data[i].x = x; + bad->data[i].y = y; + } + return bad; } Index: /trunk/psModules/src/camera/pmReadoutStack.c =================================================================== --- /trunk/psModules/src/camera/pmReadoutStack.c (revision 16599) +++ /trunk/psModules/src/camera/pmReadoutStack.c (revision 16600) @@ -9,14 +9,15 @@ bool pmReadoutUpdateSize(pmReadout *readout, int minCols, int minRows, - int numCols, int numRows, bool mask) + int numCols, int numRows, bool mask, bool weight, + psMaskType blank) { PS_ASSERT_PTR_NON_NULL(readout, false); if (readout->image) { - *(psS32*) &(readout->col0) = PS_MIN(minCols, readout->col0); - *(psS32*) &(readout->row0) = PS_MIN(minRows, readout->row0); + readout->col0 = PS_MIN(minCols, readout->col0); + readout->row0 = PS_MIN(minRows, readout->row0); } else { - *(psS32*) &(readout->col0) = minCols; - *(psS32*) &(readout->row0) = minRows; + readout->col0 = minCols; + readout->row0 = minRows; } @@ -28,5 +29,5 @@ // 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, 0.0); + psImageInit(newImage, NAN); psImageOverlaySection(newImage, readout->image, readout->col0, readout->row0, "="); psFree(readout->image); @@ -40,8 +41,21 @@ if (readout->mask->numCols < numCols || readout->mask->numRows < numRows) { psImage *newMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); - psImageInit(newMask, 0); + 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); + } + 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; } } @@ -51,5 +65,5 @@ bool pmReadoutStackValidate(int *minInputColsPtr, int *maxInputColsPtr, int *minInputRowsPtr, - int *maxInputRowsPtr, int *numColsPtr, int *numRowsPtr, + int *maxInputRowsPtr, int *numColsPtr, int *numRowsPtr, const psArray *inputs) { Index: /trunk/psModules/src/camera/pmReadoutStack.h =================================================================== --- /trunk/psModules/src/camera/pmReadoutStack.h (revision 16599) +++ /trunk/psModules/src/camera/pmReadoutStack.h (revision 16600) @@ -9,5 +9,7 @@ int minCols, int minRows, ///< Minimum coordinates int numCols, int numRows, ///< Size of images - bool mask ///< Worry about the mask? + bool mask, ///< Worry about the mask? + bool weight, ///< Worry about the weight? + psMaskType blank ///< Mask value to give to blank pixels ); Index: /trunk/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 16599) +++ /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 16600) @@ -912,7 +912,7 @@ } - pmReadoutUpdateSize(shutter, minInputCols, minInputRows, xSize, ySize, false); + pmReadoutUpdateSize(shutter, minInputCols, minInputRows, xSize, ySize, false, false, maskVal); if (pattern) { - pmReadoutUpdateSize(pattern, minInputCols, minInputRows, xSize, ySize, false); + pmReadoutUpdateSize(pattern, minInputCols, minInputRows, xSize, ySize, false, false, maskVal); }