IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 19164


Ignore:
Timestamp:
Aug 22, 2008, 12:45:36 PM (18 years ago)
Author:
Paul Price
Message:

Reworking the masking so as to avoid exploding the number of bad pixels in the convolved image. Introduced the bad/poor distinction. Bad pixels that are convolved now have a large 'poor' halo in the mask, with a smaller 'bad' halo. Got rid of the 'FOOTPRINT' mask value (internal subtraction mask): there's no need to convolve the mask in order to find out if a sprinkling of footprints are bad (not sure this is so efficient when a list of sources isn't supplied), and I needed the mask bits to signify the difference between 'bad' and 'poor'. Tested on a single subtraction, which seems to work.

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

Legend:

Unmodified
Added
Removed

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

    r18508 r19164  
    8787        }
    8888        pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
    89         if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, false, true)) {
     89        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, 0, 1.0, region, kernel, false, true)) {
    9090            psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
    9191            psFree(convRO);

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

    r19148 r19164  
    249249    psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, size, -size)); // Cut off the edges
    250250
    251     if (threaded) {
    252         psMutexLock(target);
    253     }
    254     psImage *subTarget = psImageSubset(target, region); // Target for this subregion
    255     if (threaded) {
    256         psMutexUnlock(target);
    257     }
     251    int xMin = region.x0, xMax = region.x1, yMin = region.y0, yMax = region.y1; // Bounds of region
    258252    if (background != 0.0) {
    259         psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32));
     253        for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
     254            for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) {
     255                target->data.F32[yTarget][xTarget] = convolved->data.F32[ySource][xSource] + background;
     256            }
     257        }
    260258    } else {
    261         int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
    262         for (int y = 0; y < subTarget->numRows; y++) {
    263             memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes);
     259        int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
     260        for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
     261            memcpy(&target->data.F32[yTarget][xMin], &convolved->data.F32[ySource][size], numBytes);
    264262        }
    265263    }
    266264    psFree(subConv);
    267265    psFree(convolved);
    268     if (threaded) {
    269         psMutexLock(target);
    270     }
    271     psFree(subTarget);
    272     if (threaded) {
    273         psMutexUnlock(target);
    274     }
    275266
    276267    return;
    … …  
    291282            for (int v = -size; v <= size; v++) {
    292283                for (int u = -size; u <= size; u++) {
    293                     target->data.F32[y][x] += kernel->kernel[v][u] *
    294                         image->data.F32[y - v][x - u];
     284                    target->data.F32[y][x] += kernel->kernel[v][u] * image->data.F32[y - v][x - u];
    295285                }
    296286            }
    … …  
    303293static inline void convolveRegion(psImage *convImage, // Convolved image (output)
    304294                                  psImage *convWeight, // Convolved weight map (output), or NULL
     295                                  psImage *convMask, // Convolve mask (output), or NULL
    305296                                  psKernel **kernelImage, // Convolution kernel for the image
    306297                                  psKernel **kernelWeight, // Convolution kernel for the weight map, or NULL
    … …  
    308299                                  psImage *weight, // Weight map to convolve, or NULL
    309300                                  psImage *subMask, // Subtraction mask
    310                                   psMaskType maskVal, // Mask value to apply in convolution
    311301                                  const pmSubtractionKernels *kernels, // Kernels
    312302                                  psImage *polyValues, // Polynomial values
    313303                                  float background, // Background value to apply
    314304                                  psRegion region, // Region to convolve
     305                                  float poorFrac, // Fraction for "poor"
    315306                                  bool useFFT,  // Use FFT to convolve?
    316307                                  bool wantDual // Want the dual convolution?
    … …  
    318309{
    319310    *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual);
    320     if (weight) {
     311    if (weight || subMask) {
    321312        *kernelWeight = varianceKernel(*kernelWeight, *kernelImage);
    322313    }
    323314
     315    psMaskType maskSource;              // Mask these values when
     316    psMaskType maskTarget;              // Set this mask value when convolving
     317    if (kernels->mode == PM_SUBTRACTION_MODE_1 || (kernels->mode == PM_SUBTRACTION_MODE_DUAL && !wantDual)) {
     318        maskSource = PM_SUBTRACTION_MASK_BAD_1;
     319        maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_BAD_1;
     320    } else {
     321        maskSource = PM_SUBTRACTION_MASK_BAD_2;
     322        maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_BAD_2;
     323    }
     324
     325    // Convolve the image and weight
    324326    if (useFFT) {
    325327        // Use Fast Fourier Transform to do the convolution
    326328        // This provides a big speed-up for large kernels
    327 
    328         convolveFFT(convImage, image, subMask, maskVal, *kernelImage, region, background, kernels->size);
     329        convolveFFT(convImage, image, subMask, maskSource, *kernelImage, region, background, kernels->size);
    329330        if (weight) {
    330             convolveFFT(convWeight, weight, subMask, maskVal, *kernelWeight, region, 0.0, kernels->size);
     331            convolveFFT(convWeight, weight, subMask, maskSource, *kernelWeight, region, 0.0, kernels->size);
    331332        }
    332333    } else {
    … …  
    334335        if (weight) {
    335336            convolveDirect(convWeight, weight, *kernelWeight, region, 0.0, kernels->size);
     337        }
     338    }
     339
     340    // Convolve the mask for bad pixels
     341    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
     370        if (box > 0) {
     371            bool threaded = pmSubtractionThreaded(); // Are we running threaded?
     372            if (threaded) {
     373                psMutexLock(subMask);
     374            }
     375            psImage *image = psImageSubset(subMask, psRegionSet(region.x0 + xMin, region.x1 + xMax,
     376                                                                region.y0 + yMin, region.y1 + yMax));
     377            if (threaded) {
     378                psMutexUnlock(subMask);
     379            }
     380
     381            psImage *convolved = psImageConvolveMask(NULL, image, maskSource, maskTarget,
     382                                                     -box, box, -box, box); // Convolved subtraction mask
     383
     384            if (threaded) {
     385                psMutexLock(subMask);
     386            }
     387            psFree(image);
     388            if (threaded) {
     389                psMutexUnlock(subMask);
     390            }
     391
     392            int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds
     393            int numBytes = (colMax - colMin) * PSELEMTYPE_SIZEOF(PS_TYPE_MASK); // Number of bytes to copy
     394
     395            // Since we're copying back into the source, we need to lock
     396            if (threaded) {
     397                psMutexLock(subMask);
     398            }
     399            for (int yTarget = rowMin, ySource = -yMin; yTarget < rowMax; yTarget++, ySource++) {
     400                memcpy(&subMask->data.PS_TYPE_MASK_DATA[yTarget][colMin],
     401                       &convolved->data.PS_TYPE_MASK_DATA[ySource][-xMin], numBytes);
     402            }
     403            if (threaded) {
     404                psMutexUnlock(subMask);
     405            }
     406
     407            // No need to lock: we own this
     408            psFree(convolved);
    336409        }
    337410    }
    … …  
    753826
    754827// XXX Put kernelImage, kernelWeight and polyValues on thread-dependent data
    755 static bool subtractionConvolvePatch(int numCols, int numRows, int x0, int y0,
    756                                      pmReadout *out1, pmReadout *out2, psImage *convMask,
    757                                      const pmReadout *ro1, const pmReadout *ro2, psImage *subMask,
    758                                      psMaskType blank, psMaskType maskSource, psMaskType maskTarget,
    759                                      const psRegion *region, const pmSubtractionKernels *kernels,
    760                                      bool doBG, bool useFFT)
     828static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image
     829                                     int x0, int y0, // Offsets for image
     830                                     pmReadout *out1, pmReadout *out2, // Output readouts
     831                                     psImage *convMask, // Output convolved mask
     832                                     const pmReadout *ro1, const pmReadout *ro2, // Input readouts
     833                                     psImage *subMask, // Input subtraction mask
     834                                     psMaskType maskBad, // Mask value to give bad pixels
     835                                     psMaskType maskPoor, // Mask value to give poor pixels
     836                                     float poorFrac, // Fraction for "poor"
     837                                     const psRegion *region, // Patch to convolve
     838                                     const pmSubtractionKernels *kernels, // Kernels
     839                                     bool doBG, // Add in background when convolving?
     840                                     bool useFFT // Use FFT to do the convolution?
     841    )
    761842{
    762843    int size = kernels->size;           // Half-size of kernel
    … …  
    782863
    783864    if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    784         convolveRegion(out1->image, out1->weight, &kernelImage, &kernelWeight, ro1->image, ro1->weight,
    785                        subMask, maskSource, kernels, polyValues, background, *region, useFFT,
    786                        false);
     865        convolveRegion(out1->image, out1->weight, convMask, &kernelImage, &kernelWeight,
     866                       ro1->image, ro1->weight, subMask, kernels, polyValues, background,
     867                       *region, poorFrac, useFFT, false);
    787868    }
    788869    if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    789         convolveRegion(out2->image, out2->weight, &kernelImage, &kernelWeight, ro2->image, ro2->weight,
    790                        subMask, maskSource, kernels, polyValues, background, *region, useFFT,
    791                        kernels->mode == PM_SUBTRACTION_MODE_DUAL);
     870        convolveRegion(out2->image, out2->weight, convMask, &kernelImage, &kernelWeight,
     871                       ro2->image, ro2->weight, subMask, kernels, polyValues, background,
     872                       *region, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL);
    792873    }
    793874
    … …  
    796877    psFree(polyValues);
    797878
    798     // Propagate the mask
     879    // Propagate the subtraction mask
    799880    if (subMask) {
    800881        psMaskType **target = convMask->data.PS_TYPE_MASK_DATA; // Target mask
    801882        psMaskType **source = subMask->data.PS_TYPE_MASK_DATA; // Source mask
    802883
     884        psMaskType poor, bad;           // Mask value for "poor" and "bad" pixels in subtraction mask
     885        switch (kernels->mode) {
     886          case PM_SUBTRACTION_MODE_1:
     887            poor = PM_SUBTRACTION_MASK_CONVOLVE_1;
     888            bad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 | PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
     889            break;
     890          case PM_SUBTRACTION_MODE_2:
     891            poor = PM_SUBTRACTION_MASK_CONVOLVE_2;
     892            bad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 | PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
     893            break;
     894          case PM_SUBTRACTION_MODE_DUAL:
     895            poor = PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_1;
     896            bad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 | PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 |
     897                PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
     898            break;
     899          default:
     900            psAbort("Unrecognised subtraction mode: %x", kernels->mode);
     901        }
     902
    803903        for (int y = yMin; y < yMax; y++) {
    804904            for (int x = xMin; x < xMax; x++) {
    805                 if (source[y][x] & maskTarget) {
    806                     target[y][x] |= blank;
     905                // Pixels marked as "bad" shouldn't also be marked as "poor"
     906                if (source[y][x] & bad) {
     907                    target[y][x] |= maskBad;
     908                } else if (source[y][x] & poor){
     909                    target[y][x] |= maskPoor;
    807910                }
    808911            }
    … …  
    848951    const pmReadout *ro2 = args->data[8]; // Input readout 2
    849952    psImage *subMask = args->data[9]; // Subtraction mask
    850     psMaskType blank = PS_SCALAR_VALUE(args->data[10], U8); // Output mask value
    851     psMaskType maskSource = PS_SCALAR_VALUE(args->data[11], U8); // Mask value corresponding to source image
    852     psMaskType maskTarget = PS_SCALAR_VALUE(args->data[12], U8); // Mask value corresponding to target image
     953    psMaskType maskBad = PS_SCALAR_VALUE(args->data[10], U8); // Output mask value for bad pixels
     954    psMaskType maskPoor = PS_SCALAR_VALUE(args->data[11], U8); // Output mask value for poor pixels
     955    float poorFrac = PS_SCALAR_VALUE(args->data[12], F32); // Fraction for "poor"
    853956    const psRegion *region = args->data[13]; // Region to convolve
    854957    const pmSubtractionKernels *kernels = args->data[14]; // Kernels
    … …  
    856959    bool useFFT = PS_SCALAR_VALUE(args->data[16], U8); // Use FFT for convolution?
    857960
    858     return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask,
    859                                     ro1, ro2, subMask, blank, maskSource, maskTarget,
    860                                     region, kernels, doBG, useFFT);
     961    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, subMask,
     962                                    maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT);
    861963}
    862964
    863965bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2,
    864                            psImage *subMask, psMaskType blank, const psRegion *region,
    865                            const pmSubtractionKernels *kernels, bool doBG, bool useFFT)
     966                           psImage *subMask, psMaskType maskBad, psMaskType maskPoor, float poorFrac,
     967                           const psRegion *region, const pmSubtractionKernels *kernels,
     968                           bool doBG, bool useFFT)
    866969{
    867970    int numCols = 0, numRows = 0;       // Image dimensions
    … …  
    9981101        yMin = PS_MAX(region->y0, yMin);
    9991102        yMax = PS_MIN(region->y1, yMax);
    1000     }
    1001 
    1002     psMaskType maskSource;              // Mask these pixels when convolving
    1003     psMaskType maskTarget;              // Mark these pixels as bad when propagating the subtractionMask
    1004     switch (kernels->mode) {
    1005       case PM_SUBTRACTION_MODE_1:
    1006         maskSource = PM_SUBTRACTION_MASK_BAD_1;
    1007         maskTarget = PM_SUBTRACTION_MASK_BAD_2 | PM_SUBTRACTION_MASK_CONVOLVE_1;
    1008         break;
    1009       case PM_SUBTRACTION_MODE_2:
    1010         maskSource = PM_SUBTRACTION_MASK_BAD_2;
    1011         maskTarget = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_CONVOLVE_2;
    1012         break;
    1013       case PM_SUBTRACTION_MODE_DUAL:
    1014         maskSource = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
    1015         maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2;
    1016         break;
    1017       default:
    1018         psAbort("Unsupported subtraction mode: %x", kernels->mode);
    10191103    }
    10201104
    … …  
    10451129                psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const
    10461130                psArrayAdd(args, 1, (psImage*)subMask); // Casting away const
    1047                 PS_ARRAY_ADD_SCALAR(args, blank, PS_TYPE_U8);
    1048                 PS_ARRAY_ADD_SCALAR(args, maskSource, PS_TYPE_U8);
    1049                 PS_ARRAY_ADD_SCALAR(args, maskTarget, PS_TYPE_U8);
     1131                PS_ARRAY_ADD_SCALAR(args, maskBad, PS_TYPE_U8);
     1132                PS_ARRAY_ADD_SCALAR(args, maskPoor, PS_TYPE_U8);
     1133                PS_ARRAY_ADD_SCALAR(args, poorFrac, PS_TYPE_F32);
    10501134                psArrayAdd(args, 1, subRegion);
    10511135                psArrayAdd(args, 1, (pmSubtractionKernels*)kernels); // Casting away const
    … …  
    10601144            } else {
    10611145                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, subMask,
    1062                                          blank, maskSource, maskTarget, subRegion, kernels, doBG, useFFT);
     1146                                         maskBad, maskPoor, poorFrac, subRegion, kernels, doBG, useFFT);
    10631147            }
    10641148            psFree(subRegion);

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

    r19122 r19164  
    66 * @author GLG, MHPCC
    77 *
    8  * @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
    9  * @date $Date: 2008-08-19 00:44:42 $
     8 * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
     9 * @date $Date: 2008-08-22 22:45:36 $
    1010 * Copyright 2004-207 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    2626/// Mask values for the subtraction mask
    2727typedef enum {
    28     PM_SUBTRACTION_MASK_CLEAR       = 0x00, // No masking
    29     PM_SUBTRACTION_MASK_BAD_1       = 0x01, // Image 1 is bad
    30     PM_SUBTRACTION_MASK_BAD_2       = 0x02, // Image 2 is bad
    31     PM_SUBTRACTION_MASK_CONVOLVE_1  = 0x04, // If image 1 is convolved, would be bad
    32     PM_SUBTRACTION_MASK_CONVOLVE_2  = 0x08, // If image 2 is convolved, would be bad
    33     PM_SUBTRACTION_MASK_FOOTPRINT_1 = 0x10, // Bad pixel within the stamp footprint of image 1
    34     PM_SUBTRACTION_MASK_FOOTPRINT_2 = 0x20, // Bad pixel within the stamp footprint of image 2
    35     PM_SUBTRACTION_MASK_BORDER      = 0x40, // Image border
    36     PM_SUBTRACTION_MASK_REJ         = 0x80, // Previously tried as a stamp, and rejected
     28    PM_SUBTRACTION_MASK_CLEAR          = 0x00, // No masking
     29    PM_SUBTRACTION_MASK_BAD_1          = 0x01, // Image 1 is bad
     30    PM_SUBTRACTION_MASK_BAD_2          = 0x02, // Image 2 is bad
     31    PM_SUBTRACTION_MASK_CONVOLVE_1     = 0x04, // If image 1 is convolved, would be poor or bad
     32    PM_SUBTRACTION_MASK_CONVOLVE_2     = 0x08, // If image 2 is convolved, would be poor or bad
     33    PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 = 0x10, // If image 1 is convolved, would be bad
     34    PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 = 0x20, // If image 2 is convolved, would be bad
     35    PM_SUBTRACTION_MASK_BORDER         = 0x40, // Image border
     36    PM_SUBTRACTION_MASK_REJ            = 0x80, // Previously tried as a stamp, and rejected
    3737} pmSubtractionMasks;
    3838
    … …  
    104104                           const pmReadout *ro2, // Input image 2
    105105                           psImage *subMask, ///< Subtraction mask (or NULL)
    106                            psMaskType blank, ///< Mask value for blank regions
     106                           psMaskType maskBad, ///< Mask value to give bad pixels
     107                           psMaskType maskPoor, ///< Mask value to give poor pixels
     108                           float poorFrac, ///< Fraction for "poor"
    107109                           const psRegion *region, ///< Region to convolve (or NULL)
    108110                           const pmSubtractionKernels *kernels, ///< Kernel parameters

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

    r17729 r19164  
    138138
    139139    if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_1,
    140                              PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1,
     140                             PM_SUBTRACTION_MASK_CONVOLVE_1,
    141141                             -size, size, -size, size)) {
    142142        psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1.");
    … …  
    145145    }
    146146    if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_2,
    147                              PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2,
     147                             PM_SUBTRACTION_MASK_CONVOLVE_2,
    148148                             -size, size, -size, size)) {
    149149        psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2.");
    150         psFree(mask);
    151         return NULL;
    152     }
    153     if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1,
    154                              PM_SUBTRACTION_MASK_FOOTPRINT_1,
    155                              -footprint, footprint, -footprint, footprint)) {
    156         psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1.");
    157         psFree(mask);
    158         return NULL;
    159     }
    160     if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2,
    161                              PM_SUBTRACTION_MASK_FOOTPRINT_2,
    162                              -footprint, footprint, -footprint, footprint)) {
    163         psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2.");
    164150        psFree(mask);
    165151        return NULL;

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

    r18962 r19164  
    121121                        int inner, int ringsOrder, int binning, float penalty,
    122122                        bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
    123                         int iter, float rej, psMaskType maskBad, psMaskType maskBlank,
    124                         float badFrac, pmSubtractionMode mode)
     123                        int iter, float rej, psMaskType maskVal, psMaskType maskBad, psMaskType maskPoor,
     124                        float poorFrac, float badFrac, pmSubtractionMode subMode)
    125125{
    126     if (mode != PM_SUBTRACTION_MODE_2) {
     126    if (subMode != PM_SUBTRACTION_MODE_2) {
    127127        PM_ASSERT_READOUT_NON_NULL(conv1, false);
    128128        if (conv1->image) {
    … …  
    139139        }
    140140    }
    141     if (mode != PM_SUBTRACTION_MODE_1) {
     141    if (subMode != PM_SUBTRACTION_MODE_1) {
    142142        PM_ASSERT_READOUT_NON_NULL(conv2, false);
    143143        if (conv2->image) {
    … …  
    190190    PS_ASSERT_INT_POSITIVE(iter, false);
    191191    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
     192    // Don't care about maskVal
    192193    // Don't care about maskBad
    193     // Don't care about maskBlank
     194    // Don't care about maskPoor
     195    PS_ASSERT_FLOAT_LARGER_THAN(poorFrac, 0.0, NULL);
     196    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(poorFrac, 1.0, NULL);
    194197    if (isfinite(badFrac)) {
    195198        PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL);
    … …  
    230233    memCheck("start");
    231234
    232     subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskBad, size, footprint,
     235    subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint,
    233236                                badFrac, useFFT);
    234237    if (!subMask) {
    … …  
    267270            if (sources) {
    268271                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, footprint,
    269                                                            stampSpacing, mode);
     272                                                           stampSpacing, subMode);
    270273            } else if (stampsName && strlen(stampsName) > 0) {
    271274                stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, footprint,
    272                                                         stampSpacing, mode);
     275                                                        stampSpacing, subMode);
    273276            }
    274277
    … …  
    276279            // doesn't matter.
    277280            if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing,
    278                            size, footprint, mode)) {
     281                           size, footprint, subMode)) {
    279282                goto MATCH_ERROR;
    280283            }
    281284
    282             if (mode == PM_SUBTRACTION_MODE_UNSURE) {
     285            if (subMode == PM_SUBTRACTION_MODE_UNSURE) {
    283286                // Get backgrounds
    284287                psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background
    285288                psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
    286289                psVector *buffer = NULL;// Buffer for stats
    287                 if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskBad, rng)) {
     290                if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskVal, rng)) {
    288291                    psError(PS_ERR_UNKNOWN, false, "Unable to measure background of image 1.");
    289292                    psFree(bgStats);
    … …  
    293296                }
    294297                float bg1 = psStatsGetValue(bgStats, BG_STAT); // Background for image 1
    295                 if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskBad, rng)) {
     298                if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskVal, rng)) {
    296299                    psError(PS_ERR_UNKNOWN, false, "Unable to measure background of image 2.");
    297300                    psFree(bgStats);
    … …  
    317320                    goto MATCH_ERROR;
    318321                }
    319                 mode = newMode;
     322                subMode = newMode;
    320323            }
    321324
    … …  
    323326            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    324327                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
    325                                                           stamps, footprint, optThreshold, penalty, mode);
     328                                                          stamps, footprint, optThreshold, penalty, subMode);
    326329                if (!kernels) {
    327330                    psErrorClear();
    … …  
    332335                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
    333336                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
    334                                                        inner, binning, ringsOrder, penalty, mode);
     337                                                       inner, binning, ringsOrder, penalty, subMode);
    335338            }
    336339
    … …  
    344347                }
    345348
    346                 if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) {
     349                if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    347350                    psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
    348351                                     PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
    349352                    psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
    350                                      PS_META_DUPLICATE_OK, "Subtraction kernels", mode);
     353                                     PS_META_DUPLICATE_OK, "Subtraction kernels", subMode);
    351354                    psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
    352355                                     PS_DATA_REGION | PS_META_DUPLICATE_OK,
    353356                                     "Region over which subtraction was performed", subRegion);
    354357                }
    355                 if (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) {
     358                if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    356359                    psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
    357360                                     PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
    358361                    psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
    359                                      PS_META_DUPLICATE_OK, "Subtraction kernels", mode);
     362                                     PS_META_DUPLICATE_OK, "Subtraction kernels", subMode);
    360363                    psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
    361364                                     PS_DATA_REGION | PS_META_DUPLICATE_OK,
    … …  
    374377
    375378                if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing,
    376                                size, footprint, mode)) {
     379                               size, footprint, subMode)) {
    377380                    goto MATCH_ERROR;
    378381                }
    … …  
    434437            stamps = NULL;
    435438
    436             if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) {
     439            if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    437440                psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_NUM,
    438441                                 PS_META_DUPLICATE_OK, "Number of good stamps", numStamps);
    … …  
    440443                                 PS_META_DUPLICATE_OK, "RMS deviation of stamps", rmsStamps);
    441444            }
    442             if (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) {
     445            if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    443446                psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_NUM,
    444447                                 PS_META_DUPLICATE_OK, "Number of good stamps", numStamps);
    … …  
    454457                int fullSize = 2 * size + 1 + 1;    // Full size of kernel
    455458                int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize;
    456                 psImage *convKernels = psImageAlloc((mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) *
     459                psImage *convKernels = psImageAlloc((subMode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) *
    457460                                                    imageSize - 1 +
    458                                                     (mode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0),
     461                                                    (subMode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0),
    459462                                                    imageSize - 1, PS_TYPE_F32);
    460463                psImageInit(convKernels, NAN);
    … …  
    479482                        psFree(kernel);
    480483
    481                         if (mode == PM_SUBTRACTION_MODE_DUAL) {
     484                        if (subMode == PM_SUBTRACTION_MODE_DUAL) {
    482485                            kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
    483486                                                              (float)j / (float)KERNEL_MOSAIC,
    … …  
    534537
    535538            psTrace("psModules.imcombine", 2, "Convolving...\n");
    536             if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBlank, region, kernels,
    537                                        true, useFFT)) {
     539            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBad, maskPoor, poorFrac,
     540                                       region, kernels, true, useFFT)) {
    538541                psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
    539542                goto MATCH_ERROR;
    … …  
    541544
    542545            // Set the variance factors
    543             switch (mode) {
     546            switch (subMode) {
    544547              case PM_SUBTRACTION_MODE_1: {
    545548                  recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false),
    … …  
    565568
    566569            // There is data in the readout now
    567             if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) {
     570            if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    568571                conv1->data_exists = true;
    569572                if (conv1->parent) {
    … …  
    572575                }
    573576            }
    574             if (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) {
     577            if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    575578                conv2->data_exists = true;
    576579                if (conv2->parent) {
    … …  
    590593    weight = NULL;
    591594
    592     if (!pmSubtractionBorder(conv1->image, conv1->weight, conv1->mask, size, maskBlank)) {
     595    if (!pmSubtractionBorder(conv1->image, conv1->weight, conv1->mask, size, maskBad)) {
    593596        psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
    594597        goto MATCH_ERROR;
    … …  
    600603#ifdef TESTING
    601604    {
    602         if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) {
     605        if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    603606            psFits *fits = psFitsOpen("convolved1.fits", "w");
    604607            psFitsWriteImage(fits, NULL, conv1->image, 0, NULL);
    … …  
    606609        }
    607610
    608         if (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) {
     611        if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    609612            psFits *fits = psFitsOpen("convolved2.fits", "w");
    610613            psFitsWriteImage(fits, NULL, conv2->image, 0, NULL);

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

    r18962 r19164  
    4646                        int iter,       ///< Rejection iterations
    4747                        float rej,      ///< Rejection threshold
    48                         psMaskType maskBad, ///< Value to mask
    49                         psMaskType maskBlank, ///< Mask for blank region
     48                        psMaskType maskVal, ///< Value to mask for input
     49                        psMaskType maskBad, ///< Mask for output bad pixels
     50                        psMaskType maskPoor, ///< Mask for output poor pixels
     51                        float poorFrac, ///< Fraction for "poor"
    5052                        float badFrac,   ///< Maximum fraction of bad input pixels to accept
    5153                        pmSubtractionMode mode ///< Mode of subtraction; may be modified

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

    r17297 r19164  
    7373static bool checkStampMask(int x, int y, // Coordinates of stamp
    7474                           const psImage *mask, // Mask
    75                            pmSubtractionMode mode // Mode for subtraction
     75                           pmSubtractionMode mode, // Mode for subtraction
     76                           int footprint // Footprint to check for Bad Stuff
    7677                           )
    7778{
    … …  
    7980        return true;
    8081    }
    81     if (x < 0 || x >= mask->numCols || y < 0 || y >= mask->numRows) {
    82         return false;
    83     }
    84 
    85     psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ; // Mask value
     82
     83    bool clean = true;                  // Is the footprint clean?
     84    int numCols = mask->numCols, numRows = mask->numRows; // Size of image
     85
     86    // Check the footprint bounds
     87    if (x < footprint || x >= numCols - footprint || y < footprint || y >= numRows - footprint) {
     88        clean = false;
     89    }
     90
     91    // Determine mask value
     92    psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
    8693    switch (mode) {
    8794      case PM_SUBTRACTION_MODE_1:
    88         maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1;
     95        maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_1;
    8996        break;
    9097      case PM_SUBTRACTION_MODE_2:
    91         maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_2;
     98        maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_2;
    9299        break;
    93100      case PM_SUBTRACTION_MODE_UNSURE:
    94101      case PM_SUBTRACTION_MODE_DUAL:
    95         maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1 | PM_SUBTRACTION_MASK_FOOTPRINT_2;
     102        maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2;
    96103        break;
    97104      default:
    … …  
    99106    }
    100107
    101     return (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? false : true;
     108    // Check the immediate pixel
     109    if (clean && (mask->data.PS_TYPE_MASK_DATA[y][x] & (maskVal | PM_SUBTRACTION_MASK_REJ))) {
     110        clean = false;
     111    }
     112
     113    int xMin = PS_MAX(x - footprint, 0), xMax = PS_MIN(x + footprint, numCols - 1); // Bounds in x
     114    int yMin = PS_MAX(y - footprint, 0), yMax = PS_MIN(y + footprint, numRows - 1); // Bounds in y
     115
     116    // Check the footprint
     117    if (clean) {
     118        for (int j = yMin; j <= yMax; j++) {
     119            for (int i = xMin; i <= xMax; i++) {
     120                if (mask->data.PS_TYPE_MASK_DATA[j][i] & maskVal) {
     121                    clean = false;
     122                    goto CHECK_STAMP_MASK_DONE;
     123                }
     124            }
     125        }
     126    }
     127CHECK_STAMP_MASK_DONE:
     128
     129    if (!clean) {
     130        // Mask the footprint, so we don't select something near it again
     131        for (int j = yMin; j <= yMax; j++) {
     132            for (int i = xMin; i <= xMax; i++) {
     133                mask->data.PS_TYPE_MASK_DATA[j][i] |= PM_SUBTRACTION_MASK_REJ;
     134            }
     135        }
     136    }
     137
     138    return clean;
    102139}
    103140
    … …  
    259296            for (int y = subRegion->y0; y <= subRegion->y1; y++) {
    260297                for (int x = subRegion->x0; x <= subRegion->x1; x++) {
    261                     if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) {
     298                    if (checkStampMask(x, y, subMask, mode, footprint) && image->data.F32[y][x] > fluxStamp) {
    262299                        fluxStamp = image->data.F32[y][x];
    263300                        xStamp = x;
    … …  
    353390            continue;
    354391        }
    355         if (!checkStampMask(xPix, yPix, subMask, mode)) {
     392        if (!checkStampMask(xPix, yPix, subMask, mode, footprint)) {
    356393            // Not a good stamp
    357394            psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask",
Note: See TracChangeset for help on using the changeset viewer.