IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 19282


Ignore:
Timestamp:
Aug 28, 2008, 6:04:19 PM (18 years ago)
Author:
Paul Price
Message:

We want a bad/poor mask distinction in the mask convolution for stacking. Again, otherwise the number of bad pixels just explodes. Since convolving a mask with the PSF-matching kernel is something we want to do in both stack rejection and normal convolution, pulled the code that determines the size of the 'bad' region into a separate function in pmSubtraction.c. Then pmStackReject basically does the same thing that's happening in pmSubtractionConvolve/convolveRegion: go through the image by footprints, convolve the footprint and plug the result back in to the original image. Small API change for pmStackReject. This is very effective in reducing the amount of bad pixels in the stacked image.

Location:
trunk/psModules/src/imcombine
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmStackReject.c

    r19164 r19282  
    44
    55#include <stdio.h>
     6#include <string.h>
    67#include <pslib.h>
    78
    … …  
    1415
    1516
    16 psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold,
    17                         const psArray *subRegions, const psArray *kernels)
     17psPixels *pmStackReject(const psPixels *in, int numCols, int numRows, float threshold, float poorFrac,
     18                        const psArray *subRegions, const psArray *subKernels)
    1819{
    1920    PS_ASSERT_PIXELS_NON_NULL(in, NULL);
    … …  
    2122    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(threshold, 1.0, NULL);
    2223    PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL);
    23     PS_ASSERT_ARRAY_NON_NULL(kernels, NULL);
    24     PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, NULL);
     24    PS_ASSERT_ARRAY_NON_NULL(subKernels, NULL);
     25    PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, subKernels, NULL);
    2526
    2627    // Trivial case
    … …  
    2930    }
    3031
    31     // Get the original image size
    32     int numRegions = subRegions->n;        // Number of regions
    33     int numCols = 0, numRows = 0;       // Size of original image
    34     int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image
     32    // Check consistency of kernels
     33    int numRegions = subRegions->n;     // Number of regions
    3534    int size = 0;                       // Size of kernel
    3635    for (int i = 0; i < numRegions; i++) {
    37         psRegion *subRegion = subRegions->data[i]; // Region of interest
    38         if (subRegion->x0 < minCols) {
    39             minCols = subRegion->x0;
    40         }
    41         if (subRegion->y0 < minRows) {
    42             minRows = subRegion->y0;
    43         }
    44         if (subRegion->x1 > numCols) {
    45             numCols = subRegion->x1;
    46         }
    47         if (subRegion->y1 > numRows) {
    48             numRows = subRegion->y1;
    49         }
    50 
    51         pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
     36        pmSubtractionKernels *kernels = subKernels->data[i]; // Kernel of interest
    5237        if (size == 0) {
    53             size = kernel->size;
    54         } else if (kernel->size != size) {
     38            size = kernels->size;
     39        } else if (kernels->size != size) {
    5540            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Kernel sizes are not identical: %d vs %d",
    56                     size, kernel->size);
     41                    size, kernels->size);
    5742            return NULL;
    5843        }
    5944    }
    6045
    61     // Adjust the size for the size of the subimage
    62     if (valid) {
    63         minCols = PS_MAX(valid->x0, minCols);
    64         minRows = PS_MAX(valid->y0, minRows);
    65         numCols = PS_MIN(valid->x1, numCols);
    66         numRows = PS_MIN(valid->y1, numRows);
    67     }
    68     psTrace("psModules.imcombine", 1, "Rejecting [%d:%d,%d:%d]\n", minCols, numCols, minRows, numRows);
    69 
    70     psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1),
    71                                    0x01); // Mask
     46    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols - 1, 0, numRows - 1), 0x01); // Mask
    7247    psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve
    7348    psFree(mask);
    … …  
    7752    pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
    7853    inRO->image = image;
    79     inRO->col0 = minCols;
    80     inRO->row0 = minRows;
    8154    for (int i = 0; i < numRegions; i++) {
    8255        psRegion *region = subRegions->data[i]; // Region of interest
    83         if (valid && (region->x0 > valid->x1 || region->x1 < valid->x0 ||
    84                       region->y0 > valid->y1 || region->y1 < valid->y0)) {
    85             // Region is outside of our sub-image
    86             continue;
    87         }
    88         pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
    89         if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, 0, 1.0, region, kernel, false, true)) {
     56        pmSubtractionKernels *kernels = subKernels->data[i]; // Kernel of interest
     57        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, 0, 1.0, region, kernels, false, true)) {
    9058            psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
    9159            psFree(convRO);
    … …  
    9866
    9967        // Image of the kernel at the centre of the region
    100         float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernel->numCols/2.0) /
    101             (float)kernel->numCols;
    102         float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernel->numRows/2.0) /
    103             (float)kernel->numRows;
    104         psImage *image = pmSubtractionKernelImage(kernel, xNorm, yNorm, false);
     68        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernels->numCols/2.0) /
     69            (float)kernels->numCols;
     70        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernels->numRows/2.0) /
     71            (float)kernels->numRows;
     72        psImage *image = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
    10573        if (!image) {
    10674            psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
    … …  
    11886
    11987        // Range for normalisation
    120         int yMin = PS_MAX(minRows, region->y0) - inRO->row0;
    121         int yMax = PS_MIN(numRows - 1, region->y1) - inRO->row0;
    122         int xMin = PS_MAX(minCols, region->x0) - inRO->col0;
    123         int xMax = PS_MIN(numCols - 1, region->x1) - inRO->col0;
     88        int xMin = PS_MAX(0, region->x0), xMax = PS_MIN(numCols - 1, region->x1);
     89        int yMin = PS_MAX(0, region->y0), yMax = PS_MIN(numRows - 1, region->y1);
    12490        psTrace("psModules.imcombine", 2, "Normalising convolved mask image by %f over %d:%d,%d:%d\n",
    12591                sum, xMin, xMax, yMin, yMax);
     92        sum = 1.0 / sum;
    12693        for (int y = yMin; y <= yMax; y++) {
    12794            for (int x = xMin; x <= xMax; x++) {
    128                 convRO->image->data.F32[y][x] /= sum;
     95                convRO->image->data.F32[y][x] *= sum;
    12996            }
    13097        }
    … …  
    152119        for (int x = size; x < convolved->numCols - size; x++) {
    153120            if (convolved->data.F32[y][x] > threshold) {
    154                 // Pixel coordinates in "bad" correspond to the full image
    155                 bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x + minCols, y + minRows);
     121                bad = psPixelsAdd(bad, bad->nalloc, x, y);
    156122            }
    157123        }
    … …  
    160126
    161127    // Now, grow the mask to include everything that touches a bad pixel in the convolution
    162     int x0 = minCols, y0 = minRows;     // Offset for mask image
    163     mask = psPixelsToMask(NULL, bad, psRegionSet(x0, numCols - 1, y0, numRows - 1), 0xff);
    164     for (int i = 0; i < bad->n; i++) {
    165         int xPix = bad->data[i].x - x0, yPix = bad->data[i].y - y0; // Coordinates in frame of mask image
    166         // Convolution limits
    167         int xMin = PS_MAX(xPix - size, 0);
    168         int xMax = PS_MIN(xPix + size, mask->numCols - 1);
    169         int yMin = PS_MAX(yPix - size, 0);
    170         int yMax = PS_MIN(yPix + size, mask->numRows - 1);
    171         for (int y = yMin; y <= yMax; y++) {
    172             for (int x = xMin; x <= xMax; x++) {
    173                 assert(x < mask->numCols && y < mask->numRows);
    174                 mask->data.PS_TYPE_MASK_DATA[y][x] = 0xff;
     128    // Mask values:
     129    // 0x0f: we think this is bad
     130    // 0xf0: this is within the convolution kernel of a bad pixel
     131    mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0x0f);
     132    for (int i = 0; i < subRegions->n; i++) {
     133        psRegion *region = subRegions->data[i]; // Subtraction region
     134        pmSubtractionKernels *kernels = subKernels->data[i]; // Subtraction kernel
     135
     136        int size = kernels->size;           // Half-size of kernel
     137        int fullSize = 2 * size + 1;        // Full size of kernel
     138
     139        // Get region for convolution: [xMin:xMax,yMin:yMax]
     140        int xMin = PS_MAX(region->x0, size), xMax = PS_MIN(region->x1, numCols - size);
     141        int yMin = PS_MAX(region->y0, size), yMax = PS_MIN(region->y1, numRows - size);
     142
     143        psImage *polyValues = NULL;     // Pre-calculated polynomial values
     144        for (int j = yMin; j < yMax; j += fullSize) {
     145            int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
     146            for (int i = xMin; i < xMax; i += fullSize) {
     147                int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest
     148
     149                polyValues = p_pmSubtractionPolynomialFromCoords(polyValues, kernels, numCols, numRows,
     150                                                                 i + size + 1, j + size + 1);
     151                int box = p_pmSubtractionBadRadius(NULL, kernels, polyValues,
     152                                                   false, poorFrac); // Radius of bad box
     153                if (box > 0) {
     154                    // Convolve a subimage, then stick it in the original
     155                    psImage *subMask = psImageSubset(mask, psRegionSet(i - box, xSubMax + box,
     156                                                                       j - box, ySubMax + box)); // Subimage
     157                    psImage *convolved = psImageConvolveMask(NULL, subMask, 0x0f, 0xf0,
     158                                                             -box, box, -box, box); // Convolved mask
     159                    psFree(subMask);
     160
     161                    int numBytes = (xSubMax - i) * PSELEMTYPE_SIZEOF(PS_TYPE_MASK); // Number of bytes to copy
     162                    psAssert(convolved->numCols - 2 * box == xSubMax - i, "Bad number of columns");
     163                    psAssert(convolved->numRows - 2 * box == ySubMax - j, "Bad number of rows");
     164
     165                    for (int yTarget = j, ySource = box; yTarget < ySubMax; yTarget++, ySource++) {
     166                        memcpy(&mask->data.PS_TYPE_MASK_DATA[yTarget][i],
     167                               &convolved->data.PS_TYPE_MASK_DATA[ySource][box], numBytes);
     168                    }
     169                    psFree(convolved);
     170                }
    175171            }
    176172        }
     173        psFree(polyValues);
    177174    }
    178175    bad = psPixelsFromMask(bad, mask, 0xff);
    179     psFree(mask);
    180 
    181     // Convert coordinates to frame of original image
    182     for (int i = 0; i < bad->n; i++) {
    183         int x = bad->data[i].x + x0;
    184         int y = bad->data[i].y + y0;
    185         if (x < 0 || x >= numCols || y < 0 || y >= numRows) {
    186             psWarning("Bad pixel coordinate %d: %d,%d --- ignored.",
    187                       i, x, y);
    188             continue;
    189         }
    190         bad->data[i].x = x;
    191         bad->data[i].y = y;
    192     }
    193176
    194177    return bad;

  • trunk/psModules/src/imcombine/pmStackReject.h

    r16607 r19282  
    1010/// We apply a matched filter to the corresponding mask image, and threshold to find the original pixels
    1111psPixels *pmStackReject(const psPixels *in, ///< List of pixels in the convolved image
    12                         const psRegion *valid, ///< Valid region to consider
     12                        int numCols, int numRows, ///< Size of image of interest
    1313                        float threshold, ///< Threshold on convolved image, 0..1
     14                        float poorFrac, ///< Fraction for "poor"
    1415                        const psArray *regions, ///< Array of image regions for image
    1516                        const psArray *kernels ///< Array of kernel parameters for each region

  • trunk/psModules/src/imcombine/pmSubtraction.c

    r19271 r19282  
    300300                                  psImage *subMask, // Subtraction mask
    301301                                  const pmSubtractionKernels *kernels, // Kernels
    302                                   psImage *polyValues, // Polynomial values
     302                                  const psImage *polyValues, // Polynomial values
    303303                                  float background, // Background value to apply
    304304                                  psRegion region, // Region to convolve
    … …  
    340340    // Convolve the mask for bad pixels
    341341    if (subMask && convMask) {
    342         psKernel *kernel = *kernelWeight; // Kernel of interest
    343         int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds
    344 
    345         // Determine the thresholds
    346         double sumKernel2 = 0.0;        // Sum of the kernel-squared
    347         for (int y = yMin; y <= yMax; y++) {
    348             for (int x = xMin; x <= xMax; x++) {
    349                 sumKernel2 += kernel->kernel[y][x];
    350             }
    351         }
    352         float threshold = sumKernel2 * poorFrac; // Threshold between poor and bad
    353 
    354         // Get bounds of threshold region
    355         // Start with the entire kernel, and keep reducing the size of the box until it drops below threshold
    356         int box = kernels->size;                    // Size of box with bad pixels
    357         for (double sumBox = sumKernel2; box > 0; box--) {
    358             for (int x = -box; x <= box; x++) {
    359                 sumBox -= kernel->kernel[-box][x] + kernel->kernel[box][x];
    360             }
    361             for (int y = -box + 1; y <= box - 1; y++) {
    362                 // Note: not doing corners
    363                 sumBox -= kernel->kernel[y][-box] + kernel->kernel[y][box];
    364             }
    365             if (sumBox < threshold) {
    366                 break;
    367             }
    368         }
    369 
     342        int box = p_pmSubtractionBadRadius(*kernelImage, kernels, polyValues,
     343                                           wantDual, poorFrac); // Size of bad box
    370344        if (box > 0) {
    371345            int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds
    … …  
    415389    return;
    416390}
     391
     392
    417393
    418394//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    430406    psFree(conv);
    431407    return convolved;
     408}
     409
     410int p_pmSubtractionBadRadius(psKernel *preKernel, const pmSubtractionKernels *kernels,
     411                             const psImage *polyValues, bool wantDual, float poorFrac)
     412{
     413    psKernel *kernel;                   // Kernel to use
     414    if (!preKernel) {
     415        kernel = solvedKernel(NULL, kernels, polyValues, wantDual);
     416    } else {
     417        kernel = psMemIncrRefCounter(preKernel);
     418    }
     419    PS_ASSERT_IMAGE_NON_NULL(polyValues, -1);
     420
     421    int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds
     422
     423    // Determine the threshold between bad and poor
     424    double sumKernel2 = 0.0;            // Sum of the kernel-squared
     425    for (int y = yMin; y <= yMax; y++) {
     426        for (int x = xMin; x <= xMax; x++) {
     427            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
     428        }
     429    }
     430    float threshold = sumKernel2 * poorFrac; // Threshold between poor and bad
     431
     432    // Get bounds of threshold region
     433    // Start with the entire kernel, and keep reducing the size of the box until it drops below threshold
     434    int box = kernels->size;                    // Size of box with bad pixels
     435    for (double sumBox = sumKernel2; box > 0; box--) {
     436        for (int x = -box; x <= box; x++) {
     437            sumBox -= PS_SQR(kernel->kernel[-box][x]) + PS_SQR(kernel->kernel[box][x]);
     438        }
     439        for (int y = -box + 1; y <= box - 1; y++) {
     440            // Note: not doing corners
     441            sumBox -= PS_SQR(kernel->kernel[y][-box]) + PS_SQR(kernel->kernel[y][box]);
     442        }
     443        if (sumBox < threshold) {
     444            break;
     445        }
     446    }
     447
     448    psFree(kernel);
     449
     450    return box;
     451}
     452
     453psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels,
     454                                             int numCols, int numRows, int x, int y)
     455{
     456    assert(kernels);
     457    assert(numCols > 0 && numRows > 0);
     458
     459    // Size to use when calculating normalised coordinates (different from actual size when convolving
     460    // subimage)
     461    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
     462    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
     463
     464    // Normalised coordinates
     465    float yNorm = 2.0 * (float)(y - yNormSize/2.0) / (float)yNormSize;
     466    float xNorm = 2.0 * (float)(x - xNormSize/2.0) / (float)xNormSize;
     467
     468    return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm);
    432469}
    433470
    … …  
    847884    int xMin = region->x0, xMax = region->x1, yMin = region->y0, yMax = region->y1; // Bounds of patch
    848885
    849     // Size to use when calculating normalised coordinates (different from actual size when convolving
    850     // subimage)
    851     int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
    852     int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
    853 
    854     // Normalised coordinates
    855     float yNorm = 2.0 * (float)(yMin + y0 + size + 1 - yNormSize/2.0) / (float)yNormSize; // Normalised coord
    856     float xNorm = 2.0 * (float)(xMin + x0 + size + 1 - xNormSize/2.0) / (float)xNormSize; // Normalised coord
    857 
    858886    psKernel *kernelImage = NULL;       // Kernel for the images
    859887    psKernel *kernelWeight = NULL;      // Kernel for the weight maps
    … …  
    861889    // Only generate polynomial values every kernel footprint, since we have already assumed
    862890    // (with the stamps) that it does not vary rapidly on this scale.
    863     psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder,
    864                                                     xNorm, yNorm); // Pre-calculated polynomial values
     891    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
     892                                                              xMin + x0 + size + 1,
     893                                                              yMin + y0 + size + 1);
    865894    float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term
    866895

  • trunk/psModules/src/imcombine/pmSubtraction.h

    r19164 r19282  
    66 * @author GLG, MHPCC
    77 *
    8  * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
    9  * @date $Date: 2008-08-22 22:45:36 $
     8 * @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
     9 * @date $Date: 2008-08-29 04:04:19 $
    1010 * Copyright 2004-207 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    126126    );
    127127
     128/// Given pixel coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j
     129///
     130/// Same as p_pmSubtractionPolynomial except that the normalisation is applied
     131psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, ///< Output matrix, or NULL
     132                                             const pmSubtractionKernels *kernels, ///< Kernel parameters
     133                                             int numCols, int numRows, ///< Size of image of interest
     134                                             int x, int y ///< Position of interest
     135    );
     136
     137/// Return the radius from the centre of the convolution kernel that distinguishes "bad" and "poor" pixels
     138int p_pmSubtractionBadRadius(psKernel *preKernel, ///< Pre-calculated convolution kernel
     139                             const pmSubtractionKernels *kernels, ///< Kernel parameters
     140                             const psImage *polyValues, ///< Polynomial values
     141                             bool wantDual, ///< Calculate for the dual kernel?
     142                             float poorFrac ///< Fraction for "poor"
     143    );
     144
    128145/// @}
    129146#endif
Note: See TracChangeset for help on using the changeset viewer.