IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16600


Ignore:
Timestamp:
Feb 22, 2008, 9:18:11 AM (18 years ago)
Author:
Paul Price
Message:

ppStack is now working with the incremental reads. Needed to modify a few things in pmStack*.c so that it would combine (and reject) subimages into a large image. Added weights to pmReadoutStack.

Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/pap_branch_080214/psModules/src/imcombine/pmReadoutCombine.c

    r15433 r16600  
    113113    }
    114114
    115     pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true);
     115    pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true, params->weights,
     116                        params->blank);
    116117    psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
    117118
    118119    psStatsOptions combineStdev = 0; // Statistics option for weights
    119120    if (params->weights) {
    120 
    121         if (!output->weight) {
    122             output->weight = psImageAlloc(xSize, ySize, PS_TYPE_F32);
    123         }
    124         if (output->weight->numCols < xSize || output->weight->numRows < ySize) {
    125             psImage *newWeight = psImageAlloc(xSize, ySize, PS_TYPE_F32);
    126             psImageInit(newWeight, 0.0);
    127             psImageOverlaySection(newWeight, output->weight, output->col0, output->row0, "=");
    128             psFree(output->weight);
    129             output->weight = newWeight;
    130         }
    131 
    132121        if (first) {
    133122            psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,

  • branches/pap_branch_080214/psModules/src/imcombine/pmStack.c

    r16420 r16600  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2008-02-13 02:55:33 $
     10 *  @version $Revision: 1.16.2.1 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2008-02-22 19:18:11 $
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    1313 *
    … …  
    2323#include "pmHDU.h"
    2424#include "pmFPA.h"
     25#include "pmReadoutStack.h"
    2526#include "pmConceptsAverage.h"
    2627
    … …  
    4243    psVector *weights;                  // Pixel weights
    4344    psVector *sort;                     // Buffer for sorting (to get a robust estimator of the standard dev)
     45#if 0
     46    int x0, y0;                         // Offset from original image to combination region
     47    int nx, ny;                         // Number of pixels to combine
     48#endif
    4449} combineBuffer;
    4550
    … …  
    6469    buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
    6570    buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32);
     71
     72#if 0
     73    buffer->x0 = 0;
     74    buffer->y0 = 0;
     75    buffer->nx = 0;
     76    buffer->ny = 0;
     77#endif
    6678
    6779    return buffer;
    … …  
    136148        return false;
    137149    }
    138     *median = num % 2 ? (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0 :
    139         sortBuffer->data.F32[num / 2];
    140     *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]);
     150    *median = num % 2 ? sortBuffer->data.F32[num / 2] :
     151        (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0 ;
     152#if 0
     153    if (num < NUM_DIRECT_STDEV) {
     154#endif
     155        // If there are not many values, the direct standard deviation is better
     156        double sum = 0.0;
     157        for (int i = 0; i < num; i++) {
     158            sum += PS_SQR(sortBuffer->data.F32[i] - *median);
     159        }
     160        *stdev = sqrt(sum / (float)(num - 1));
     161#if 0
     162    } else {
     163        *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]);
     164    }
     165#endif
    141166
    142167    return true;
    … …  
    186211
    187212// Given a stack of images, combine with optional rejection.
    188 // Pixels in the stack that are rejected are marked for subsequent
     213// Pixels in the stack that are rejected are marked for subsequent inspection
    189214static bool combinePixels(psImage *image, // Combined image, for output
    190215                          psImage *mask, // Combined mask, for output
    … …  
    193218                          const psVector *weights, // Global (single value) weights for data, or NULL
    194219                          const psVector *reject, // Indices of pixels to reject, or NULL
    195                           int x, int y, // Coordinates of interest
     220                          int x, int y, // Coordinates of interest; frame of output image
    196221                          psMaskType maskVal, // Value to mask
    197222                          psMaskType bad, // Value to give bad pixels
    … …  
    223248            continue;
    224249        }
     250        int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
    225251        psImage *image = data->readout->image; // Image of interest
    226252        psImage *weight = data->readout->weight; // Weight map of interest
    227253        psImage *mask = data->readout->mask; // Mask of interest
    228         pixelData->data.F32[i] = image->data.F32[y][x];
     254        pixelData->data.F32[i] = image->data.F32[yIn][xIn];
    229255        if (weight) {
    230             pixelWeights->data.F32[i] = weight->data.F32[y][x];
    231         }
    232         pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[y][x];
     256            pixelWeights->data.F32[i] = weight->data.F32[yIn][xIn];
     257        }
     258        pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[yIn][xIn];
    233259        if (pixelMasks->data.PS_TYPE_MASK_DATA[i] & maskVal) {
    234260            numBad++;
    … …  
    567593    }
    568594
     595    // Get the sizes
     596    psArray *stack = psArrayAlloc(input->n);
     597    for (int i = 0; i < stack->n; i++) {
     598        pmStackData *data = input->data[i]; // Stack data
     599        stack->data[i] = psMemIncrRefCounter(data->readout);
     600    }
     601    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
     602    int xSize, ySize;                   // Size of the output image
     603    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
     604                                stack)) {
     605        psError(PS_ERR_UNKNOWN, false, "Input stack is not valid.");
     606        psFree(stack);
     607        return false;
     608    }
     609    psFree(stack);
     610    pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, bad);
     611    psTrace("psModules.imcombine", 1, "Combining [%d:%d,%d:%d] (%dx%d)\n",
     612            minInputCols, maxInputCols, minInputRows, maxInputRows, xSize, ySize);
     613
     614
    569615    // Buffer for combination
    570616    combineBuffer *buffer = combineBufferAlloc(num, numIter == 0 ? PS_STAT_SAMPLE_MEAN :
    … …  
    577623        psImage *combinedWeight = combined->weight; // Combined mask
    578624
    579         psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source
     625        psArray *pixelMap = pixelMapGenerate(input, maxInputCols, maxInputRows); // Map of pixels to source
    580626        psPixels *pixels = NULL;            // Total list of pixels, with no duplicates
    581627        for (int i = 0; i < num; i++) {
    … …  
    588634        }
    589635        for (int i = 0; i < pixels->n; i++) {
     636            // Pixel coordinates are in the frame of the original image
    590637            int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest
     638            if (x < minInputCols || x >= maxInputCols || y < minInputRows || y >= maxInputRows) {
     639                continue;
     640            }
    591641            psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely
    592642            combinePixels(combinedImage, combinedMask, combinedWeight, input, weights, reject, x, y,
    … …  
    626676        }
    627677
    628         for (int y = 0; y < numRows; y++) {
    629             for (int x = 0; x < numCols; x++) {
     678        for (int y = minInputRows; y < maxInputRows; y++) {
     679            for (int x = minInputCols; x < maxInputCols; x++) {
    630680                combinePixels(combinedImage, combinedMask, combinedWeights, input, weights, NULL, x, y,
    631681                              maskVal, bad, numIter, rej, buffer);

  • branches/pap_branch_080214/psModules/src/imcombine/pmStackReject.c

    r16476 r16600  
    2424    int numRegions = subRegions->n;        // Number of regions
    2525    int numCols = 0, numRows = 0;       // Size of original image
    26     int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image --- should be 0,0
     26    int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image
    2727    int size = 0;                       // Size of kernel
    2828    for (int i = 0; i < numRegions; i++) {
    … …  
    5858        numRows = PS_MIN(valid->y1, numRows);
    5959    }
     60    psTrace("psModules.imcombine", 1, "Rejecting [%d:%d,%d:%d]\n", minCols, numCols, minRows, numRows);
    6061
    6162    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1),
    … …  
    7879        }
    7980        pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
    80         if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, true)) {
     81        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, false, true)) {
    8182            psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
    8283            psFree(convRO);
    … …  
    127128    // Threshold the convolved image
    128129    psPixels *bad = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); // List of pixels that should be masked
    129     for (int y = 0; y < convolved->numRows; y++) {
    130         for (int x = 0; x < convolved->numCols; x++) {
     130    for (int y = size; y < convolved->numRows - size; y++) {
     131        for (int x = size; x < convolved->numCols - size; x++) {
    131132            if (convolved->data.F32[y][x] > threshold) {
    132                 bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x, y);
     133                // Pixel coordinates in "bad" correspond to the full image
     134                bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x + minCols, y + minRows);
    133135            }
    134136        }
    … …  
    137139
    138140    // Now, grow the mask to include everything that touches a bad pixel in the convolution
    139     mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0xff);
    140     assert(mask->numCols == numCols && mask->numRows == numRows);
     141    int x0 = minCols, y0 = minRows;     // Offset for mask image
     142    mask = psPixelsToMask(NULL, bad, psRegionSet(x0, numCols - 1, y0, numRows - 1), 0xff);
    141143    for (int i = 0; i < bad->n; i++) {
    142         int xPix = bad->data[i].x, yPix = bad->data[i].y; // Coordinates of interest
     144        int xPix = bad->data[i].x - x0, yPix = bad->data[i].y - y0; // Coordinates in frame of mask image
    143145        // Convolution limits
    144146        int xMin = PS_MAX(xPix - size, 0);
    145         int xMax = PS_MIN(xPix + size, numCols - 1);
     147        int xMax = PS_MIN(xPix + size, mask->numCols - 1);
    146148        int yMin = PS_MAX(yPix - size, 0);
    147         int yMax = PS_MIN(yPix + size, numRows - 1);
     149        int yMax = PS_MIN(yPix + size, mask->numRows - 1);
    148150        for (int y = yMin; y <= yMax; y++) {
    149151            for (int x = xMin; x <= xMax; x++) {
    … …  
    153155        }
    154156    }
    155 
    156157    bad = psPixelsFromMask(bad, mask, 0xff);
    157158    psFree(mask);
    158159
     160    // Convert coordinates to frame of original image
     161    for (int i = 0; i < bad->n; i++) {
     162        int x = bad->data[i].x + x0;
     163        int y = bad->data[i].y + y0;
     164        if (x < 0 || x >= numCols || y < 0 || y >= numRows) {
     165            psWarning("Bad pixel coordinate %d: %d,%d --- ignored.",
     166                      i, x, y);
     167            continue;
     168        }
     169        bad->data[i].x = x;
     170        bad->data[i].y = y;
     171    }
     172
    159173    return bad;
    160174}

  • trunk/psModules/src/camera/pmReadoutStack.c

    r15974 r16600  
    99
    1010bool pmReadoutUpdateSize(pmReadout *readout, int minCols, int minRows,
    11                          int numCols, int numRows, bool mask)
     11                         int numCols, int numRows, bool mask, bool weight,
     12                         psMaskType blank)
    1213{
    1314    PS_ASSERT_PTR_NON_NULL(readout, false);
    1415
    1516    if (readout->image) {
    16         *(psS32*) &(readout->col0) = PS_MIN(minCols, readout->col0);
    17         *(psS32*) &(readout->row0) = PS_MIN(minRows, readout->row0);
     17        readout->col0 = PS_MIN(minCols, readout->col0);
     18        readout->row0 = PS_MIN(minRows, readout->row0);
    1819    } else {
    19         *(psS32*) &(readout->col0) = minCols;
    20         *(psS32*) &(readout->row0) = minRows;
     20        readout->col0 = minCols;
     21        readout->row0 = minRows;
    2122    }
    2223
    … …  
    2829        // Generate the new output image by extending the current one, or making a whole new one
    2930        psImage *newImage = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    30         psImageInit(newImage, 0.0);
     31        psImageInit(newImage, NAN);
    3132        psImageOverlaySection(newImage, readout->image, readout->col0, readout->row0, "=");
    3233        psFree(readout->image);
    … …  
    4041        if (readout->mask->numCols < numCols || readout->mask->numRows < numRows) {
    4142            psImage *newMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
    42             psImageInit(newMask, 0);
     43            psImageInit(newMask, blank);
    4344            psImageOverlaySection(newMask, readout->mask, readout->col0, readout->row0, "=");
    4445            psFree(readout->mask);
    4546            readout->mask = newMask;
     47        }
     48    }
     49
     50    if (weight) {
     51        if (!readout->weight) {
     52            readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     53        }
     54        if (readout->weight->numCols < numCols || readout->weight->numRows < numRows) {
     55            psImage *newWeight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     56            psImageInit(newWeight, NAN);
     57            psImageOverlaySection(newWeight, readout->weight, readout->col0, readout->row0, "=");
     58            psFree(readout->weight);
     59            readout->weight = newWeight;
    4660        }
    4761    }
    … …  
    5165
    5266bool pmReadoutStackValidate(int *minInputColsPtr, int *maxInputColsPtr, int *minInputRowsPtr,
    53                             int *maxInputRowsPtr, int *numColsPtr, int *numRowsPtr, 
     67                            int *maxInputRowsPtr, int *numColsPtr, int *numRowsPtr,
    5468                            const psArray *inputs)
    5569{

  • trunk/psModules/src/camera/pmReadoutStack.h

    r13872 r16600  
    99                         int minCols, int minRows, ///< Minimum coordinates
    1010                         int numCols, int numRows, ///< Size of images
    11                          bool mask      ///< Worry about the mask?
     11                         bool mask,     ///< Worry about the mask?
     12                         bool weight,   ///< Worry about the weight?
     13                         psMaskType blank ///< Mask value to give to blank pixels
    1214    );
    1315

  • trunk/psModules/src/detrend/pmShutterCorrection.c

    r15382 r16600  
    912912    }
    913913
    914     pmReadoutUpdateSize(shutter, minInputCols, minInputRows, xSize, ySize, false);
     914    pmReadoutUpdateSize(shutter, minInputCols, minInputRows, xSize, ySize, false, false, maskVal);
    915915    if (pattern) {
    916         pmReadoutUpdateSize(pattern, minInputCols, minInputRows, xSize, ySize, false);
     916        pmReadoutUpdateSize(pattern, minInputCols, minInputRows, xSize, ySize, false, false, maskVal);
    917917    }
    918918
Note: See TracChangeset for help on using the changeset viewer.