Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 16626) +++ trunk/psModules/src/imcombine/pmStack.c (revision 16628) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-02-23 03:00:30 $ + * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-23 03:01:45 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii * @@ -479,93 +479,4 @@ -#if 0 -// Examine a pixel carefully to see if it should be rejected in the stack by convolving all the input images -// to the same seeing, and clipping there. -static bool convolveTest(psVector *reject, // Image indices to reject for this pixel (output) - psVector *inspect, // Images indices to inspect for this pixel (input) - psArray *input, // The input pmStackData - int xPix, int yPix, // Pixel coordinates of interest - psVector *seeing, // Seeing convolution for each image - psMaskType maskVal, // Value to mask - float extent, // Extent of convolution (Gaussian sigmas) - float threshold// Rejection threshold (standard deviations) - ) -{ - assert(reject && reject->type.type == PS_TYPE_U16); - assert(input); - assert(seeing); - assert(seeing->n == input->n); - - int num = input->n; // Number of input images - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); // Statistics - psVector *values = psVectorAlloc(num, PS_TYPE_F32); // Values from convolving the image - psVector *valuesMask = psVectorAlloc(num, PS_TYPE_MASK); // Mask for convolution values - - for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stacking data - int radius = extent * seeing->data.F32[i]; // How much to convolve - psImage *image = data->readout->image; // Image to convolve - psImage *mask = data->readout->mask; // Image mask - - int xMin = PS_MAX(xPix - radius, 0); - int xMax = PS_MIN(xPix + radius, image->numCols - 1); - int yMin = PS_MAX(yPix - radius, 0); - int yMax = PS_MIN(yPix + radius, image->numRows - 1); - - float sum = 0.0; // Sum of pixels in convolution - float sumKernel = 0.0; // Sum of kernel - for (int y = yMin; y <= yMax; y++) { - int v2 = PS_SQR(y - yPix); - for (int x = xMin; x <= xMax; x++) { - int u2 = PS_SQR(x - xPix); - float kernel = expf(-(u2 + v2) / seeing->data.F32[i]); - if (!(mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal)) { - sum += image->data.F32[y][x] * kernel; - sumKernel += kernel; - } - } - } - - if (sumKernel != 0.0) { - values->data.F32[i] = sum / sumKernel; - valuesMask->data.PS_TYPE_MASK_DATA[i] = 0; - } else { - // Everything's masked out; don't include this pixel in the analysis - values->data.F32[i] = NAN; - valuesMask->data.PS_TYPE_MASK_DATA[i] = MASK_BAD; - } - } - - // Mask the ones we're inspecting - for (int i = 0; i < inspect->n; i++) { - valuesMask->data.PS_TYPE_MASK_DATA[inspect->data.U16[i]] |= MASK_SUSPECT; - } - - if (!psVectorStats(stats, values, NULL, valuesMask, MASK_BAD | MASK_SUSPECT)) { - psFree(stats); - psFree(values); - psFree(valuesMask); - return false; - } - - // Inspect the pixels of interest - float mean = stats->sampleMean; // Mean value - float limit = threshold * stats->sampleStdev; // Rejection threshold - for (int j = 0; j < inspect->n; j++) { - int index = inspect->data.U16[j]; // Index of interest - if (fabsf(values->data.F32[index] - mean) > limit) { - reject->data.U16[reject->n++] = index; - } - } - - psFree(stats); - psFree(values); - psFree(valuesMask); - - return true; -} -#endif - - // Generate a "pixel map". // @@ -785,86 +696,2 @@ } -#if 0 -bool pmStackReject(psArray *input, psMaskType maskVal, float extent, float threshold) -{ - bool haveWeights; // Do we have the sky weight maps? - bool havePixels; // Do we have lists of pixels? - int num; // Number of inputs - int numCols, numRows; // Size of (sky) images - if (!validateInputData(&haveWeights, &havePixels, &num, &numCols, &numRows, input)) { - return false; - } - if (!havePixels) { - psError(PS_ERR_UNEXPECTED_NULL, true, "No pixels have been flagged for possible rejection."); - return false; - } - PS_ASSERT_FLOAT_LARGER_THAN(extent, 0.0, false); - - // Want to concatenate the list of pixels, but also want to know where each came from - psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source - psPixels *pixels = NULL; // Total list of pixels, with no duplicates - for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels - pixels = psPixelsConcatenate(pixels, data->pixels); - // Throw out the old in preparation for the new - data->pixels = psPixelsRealloc(data->pixels, 0); - } - - // Get the convolution widths - psVector *seeing = psVectorAlloc(num, PS_TYPE_F32); // Seeing for each image - const float fwhm2sigma = 1.0 / (2.0*sqrt((double)2.0*(log((double)2.0)))); // Conversion FWHM --> sigma - float seeingMax = -INFINITY; // Maximum FWHM - for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stacking data - seeing->data.F32[i] = data->seeing * fwhm2sigma; - if (data->seeing > seeingMax) { - seeingMax = data->seeing; - } - } - (void)psBinaryOp(seeing, seeing, "*", seeing); - (void)psBinaryOp(seeing, psScalarAlloc(PS_SQR(seeingMax), PS_TYPE_F32), "-", seeing); - (void)psUnaryOp(seeing, seeing, "sqrt"); - - - psVector *reject = psVectorAllocEmpty(num, PS_TYPE_U16); // Image indices to reject for a given pixel - - for (int i = 0; i < pixels->n; i++) { - int x = pixels->data[i].x; // x coordinate of interest - int y = pixels->data[i].y; // y coordinate of interest - psVector *inspect = pixelMapQuery(pixelMap, x, y); // Inspect these images closely - - // Can daisy-chain multiple tests here - reject->n = 0; - - if (!convolveTest(reject, inspect, input, x, y, seeing, maskVal, extent, threshold)) { - // Some sort of error message here; NAN and mask all the pixels - for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stacking data - psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y); - } - continue; - } - - // Add the rejected pixels to the official rejection list - for (int i = 0; i < reject->n; i++) { - pmStackData *data = input->data[reject->data.U16[i]]; // Stacking data - psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y); - } - } - psFree(reject); - psFree(seeing); - psFree(pixels); - psFree(pixelMap); - -#ifndef PS_NO_TRACE - if (psTraceGetLevel("psModules.imcombine") >= 5) { - for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stack data for this input - psTrace("psModules.imcombine", 5, "Image %d: %ld pixels rejected.\n", i, data->pixels->n); - } - } -#endif - - return true; -} -#endif