Index: trunk/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 16600) +++ trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 16604) @@ -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: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 16600) +++ trunk/psModules/src/imcombine/pmStack.c (revision 16604) @@ -8,6 +8,4 @@ * @author GLG, MHPCC * - * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-02-14 23:33:09 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -23,4 +21,5 @@ #include "pmHDU.h" #include "pmFPA.h" +#include "pmReadoutStack.h" #include "pmConceptsAverage.h" @@ -42,4 +41,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 +67,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 +146,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 +209,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 +216,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 +246,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++; @@ -563,4 +587,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 : @@ -573,5 +617,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++) { @@ -584,5 +628,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, @@ -622,6 +670,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: trunk/psModules/src/imcombine/pmStackReject.c =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.c (revision 16600) +++ trunk/psModules/src/imcombine/pmStackReject.c (revision 16604) @@ -25,5 +25,6 @@ int numRegions = regions->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++) { psRegion *region = regions->data[i]; // Region of interest @@ -46,4 +47,5 @@ return NULL; } + psTrace("psModules.imcombine", 1, "Rejecting [%d:%d,%d:%d]\n", minCols, numCols, minRows, numRows); psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols, 0, numRows), 0x01); // Mask image @@ -56,6 +58,12 @@ inRO->image = image; for (int i = 0; i < numRegions; i++) { - psRegion *region = regions->data[i]; // Region of interest - if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernels, true)) { + psRegion *region = subRegions->data[i]; // Region of interest + if (valid && (region->x0 > valid->x1 || region->x1 < valid->x0 || + region->y0 > valid->y1 || region->y1 < valid->y0)) { + // Region is outside of our sub-image + continue; + } + pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest + 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); @@ -94,8 +102,9 @@ // Threshold the convolved image psPixels *bad = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); // List of pixels that should be masked - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < 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); } } @@ -103,14 +112,14 @@ psFree(convolved); - // Now, we want to convolve the original pixels properly - mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols, 0, numRows), 0xff); - int size = kernels->size; // Size of kernels + // Now, grow the mask to include everything that touches a bad pixel in the convolution + 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++) { @@ -119,8 +128,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/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 16600) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 16604) @@ -3,7 +3,4 @@ * @author Paul Price, IfA * @author GLG, MHPCC - * - * @version $Revision: 1.80 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-02-14 23:33:09 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -645,5 +642,5 @@ bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2, const psImage *subMask, psMaskType blank, const psRegion *region, - const pmSubtractionKernels *kernels, bool useFFT) + const pmSubtractionKernels *kernels, bool doBG, bool useFFT) { PS_ASSERT_PTR_NON_NULL(out1, false); @@ -707,21 +704,20 @@ } - psImage *image, *weight; // Image and weight map to convolve + const pmReadout *source; // Source for image parameters switch (kernels->mode) { case PM_SUBTRACTION_MODE_TARGET: case PM_SUBTRACTION_MODE_1: case PM_SUBTRACTION_MODE_DUAL: - image = ro1->image; - weight = ro1->weight; + source = ro1; break; case PM_SUBTRACTION_MODE_2: - image = ro2->image; - weight = ro2->weight; + source = ro2; break; default: psAbort("Unsupported subtraction mode: %x", kernels->mode); } - - int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions + psImage *image = source->image, *weight = source->weight; // Image and weight map to convolve + int numCols = image->numCols, numRows = image->numRows; // Image dimensions + int x0 = source->col0, y0 = source->row0; // Image offset psImage *convImage1 = out1->image; // Convolved image @@ -815,6 +811,6 @@ // (with the stamps) that it does not vary rapidly on this scale. polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, xNorm, yNorm); - float background = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background term - + float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : + 0.0; // Background term psRegion subRegion = psRegionSet(i, xSubMax, j, ySubMax); // Sub-region to convolve Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 16600) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 16604) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-12-07 01:57:15 $ + * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-22 19:24:42 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -92,4 +92,5 @@ const psRegion *region, ///< Region to convolve (or NULL) const pmSubtractionKernels *kernels, ///< Kernel parameters + bool doBG, ///< Apply background term? bool useFFT ///< Use Fast Fourier Transform for the convolution? ); Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 16600) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 16604) @@ -434,5 +434,6 @@ psTrace("psModules.imcombine", 2, "Convolving...\n"); - if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBlank, region, kernels, useFFT)) { + if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBlank, region, kernels, + true, useFFT)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); goto MATCH_ERROR;