Changeset 16628 for trunk/psModules/src/imcombine/pmStack.c

Timestamp:

Feb 22, 2008, 5:01:45 PM (18 years ago)

Author:

Paul Price

Message:

Removing old code that was #if-ed out.

File:

• trunk/psModules/src/imcombine/pmStack.c (modified) (3 diffs)

trunk/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @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 +479 @@
 
 
-#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 +696 @@
 }
 
-#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