IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 2, 2007, 4:28:24 PM (19 years ago)
Author:
Paul Price
Message:

Adding function pmSubtractionOrder (choice of name is probably not the best) to determine which of the two images under consideration should be convolved to match the other. Originally was doing this by solving for a RING kernel of width 1, going both ways, and comparing the deviations for each. However, when doing stacks this didn't work (convolving the wider fake Gaussian image gave smaller deviations than convolving the narrower input image), so have settled on measuring the second moments for each stamp, and using the ratio of moments between the two images to determine which is wider; this seems to work for both subtracting and stacking images. In the process, plugged a few memory leaks, and added code to support convolving either way (or both ways; this will be useful when it comes time to code the dual convolution algorithm).

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15247 r15443  
    5050#endif
    5151
     52/// Select the appropriate convolution, given the kernel basis function and subtraction mode
     53static inline psKernel *selectConvolution(const pmSubtractionStamp *stamp, // Stamp
     54                                          int kernelIndex, // Index for kernel component
     55                                          pmSubtractionMode mode // Mode of subtraction
     56    )
     57{
     58    switch (mode) {
     59      case PM_SUBTRACTION_MODE_1:
     60        return stamp->convolutions1->data[kernelIndex];
     61      case PM_SUBTRACTION_MODE_2:
     62        return stamp->convolutions2->data[kernelIndex];
     63      default:
     64        psAbort("Unsupported subtraction mode: %x", mode);
     65    }
     66    return NULL;                        // Unreached
     67}
     68
    5269// Accumulate cross-term sums for a stamp
    5370static void accumulateCross(double *sumI, // Sum of I(x)/sigma(x)^2
    … …  
    5572                            double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2
    5673                            const pmSubtractionStamp *stamp, // Stamp with weight
    57                             const psKernel *input, // Input image, I(x)
     74                            const psKernel *target, // Target stamp
    5875                            int kernelIndex, // Index for kernel component
    59                             int footprint // Size of region of interest
     76                            int footprint, // Size of region of interest
     77                            pmSubtractionMode mode // Mode of subtraction
    6078    )
    6179{
    6280    psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
    63     psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
     81    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    6482
    6583    for (int y = -footprint; y <= footprint; y++) {
    66         psF32 *in = &input->kernel[y][-footprint]; // Dereference input
     84        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    6785        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    6886        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    … …  
    82100                                   const pmSubtractionStamp *stamp, // Stamp with input and weight
    83101                                   int kernelIndex, // Index for kernel component
    84                                    int footprint // Size of region of interest
     102                                   int footprint, // Size of region of interest
     103                                   pmSubtractionMode mode // Mode of subtraction
    85104    )
    86105{
    87106    psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
    88     psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
     107    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    89108
    90109    for (int y = -footprint; y <= footprint; y++) {
    … …  
    100119}
    101120
    102 static double accumulateChi2(psKernel *input, // Input stamp
     121static double accumulateChi2(const psKernel *target, // Target stamp
    103122                             pmSubtractionStamp *stamp, // Stamp with weight
    104123                             int kernelIndex, // Index for kernel component
    105124                             double coeff, // Coefficient of convolution
    106125                             double bg,  // Background term
    107                              int footprint // Size of region of interest
     126                             int footprint, // Size of region of interest
     127                             pmSubtractionMode mode // Mode of subtraction
    108128    )
    109129{
    110130    double chi2 = 0.0;
    111131    psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
    112     psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
     132    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    113133
    114134    for (int y = -footprint; y <= footprint; y++) {
    115         psF32 *in = &input->kernel[y][-footprint]; // Dereference input
     135        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    116136        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    117137        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    … …  
    125145
    126146// Return the initial value of chi^2
    127 static double initialChi2(psKernel *input, // Input stamp
     147static double initialChi2(const psKernel *target, // Target stamp
    128148                          const pmSubtractionStamp *stamp, // Stamp with weight
    129                           int footprint // Size of convolution
     149                          int footprint, // Size of convolution
     150                          pmSubtractionMode mode // Mode of subtraction
    130151    )
    131152{
    132153    psKernel *weight = stamp->weight;   // Weight map
    133     psKernel *reference = stamp->reference; // Reference stamp
     154    psKernel *source;                   // Source stamp
     155    switch (mode) {
     156      case PM_SUBTRACTION_MODE_1:
     157        source = stamp->image1;
     158        break;
     159      case PM_SUBTRACTION_MODE_2:
     160        source = stamp->image2;
     161        break;
     162      default:
     163        psAbort("Unsupported subtraction mode: %x", mode);
     164    }
    134165
    135166    double chi2 = 0.0;                  // Chi^2
    136167    for (int y = -footprint; y <= footprint; y++) {
    137         psF32 *in = &input->kernel[y][-footprint]; // Dereference input
     168        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    138169        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    139         psF32 *ref = &reference->kernel[y][-footprint]; // Derference reference
     170        psF32 *ref = &source->kernel[y][-footprint]; // Derference reference
    140171        for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) {
    141172            float diff = *in - *ref;    // Temporary value
    … …  
    148179
    149180// Subtract a convolution from the input
    150 static void subtractConvolution(psKernel *input, // Input stamp
     181static void subtractConvolution(psKernel *target, // Target stamp
    151182                                const pmSubtractionStamp *stamp, // Stamp with weight
    152183                                int kernelIndex, // Index for kernel component
    153184                                float coeff, // Coefficient of subtraction
    154185                                float bg, // Background term
    155                                 int footprint // Size of region of interest
    156     )
    157 {
    158     psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
    159 
     186                                int footprint, // Size of region of interest
     187                                pmSubtractionMode mode // Mode of subtraction
     188    )
     189{
     190    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    160191    for (int y = -footprint; y <= footprint; y++) {
    161         psF32 *in = &input->kernel[y][-footprint]; // Dereference input
     192        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    162193        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    163194        for (int x = -footprint; x <= footprint; x++, in++, conv++) {
    … …  
    173204                                                      int spatialOrder, const psVector *fwhms, int maxOrder,
    174205                                                      const pmSubtractionStampList *stamps, int footprint,
    175                                                       float tolerance)
     206                                                      float tolerance, pmSubtractionMode mode)
    176207{
    177208    if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) {
    … …  
    210241    // Need to save the stamp inputs --- we're changing the values!
    211242    int numStamps = stamps->num;        // Number of stamps
    212     psArray *inputs = psArrayAlloc(numStamps); // Deep copies of the inputs
     243    psArray *targets = psArrayAlloc(numStamps); // Deep copies of the targets
    213244    psVector *badStamps = psVectorAlloc(numStamps, PS_TYPE_U8); // Mark the bad stamps
    214245    psVectorInit(badStamps, 0);
    … …  
    219250            continue;
    220251        }
    221         psKernel *input = stamp->input; // Input image of interest
    222         psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image
    223         inputs->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint);
     252        psKernel *target;               // Target image of interest
     253        switch (mode) {
     254          case PM_SUBTRACTION_MODE_1:
     255            target = stamp->image2;
     256            break;
     257          case PM_SUBTRACTION_MODE_2:
     258            target = stamp->image1;
     259            break;
     260          default:
     261            psAbort("Unsupported subtraction mode: %x", mode);
     262        }
     263        psImage *copy = psImageCopy(NULL, target->image, PS_TYPE_F32); // Copy of the image
     264        targets->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint);
    224265        psFree(copy);                   // Drop reference
    225266    }
    … …  
    238279            continue;
    239280        }
    240         if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
     281        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) {
    241282            psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
    242             psFree(inputs);
     283            psFree(targets);
    243284            psFree(kernels);
    244285            psFree(badStamps);
    … …  
    258299                    "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n",
    259300                    i, (int)stamp->x, (int)stamp->y);
    260             psFree(inputs);
     301            psFree(targets);
    261302            psFree(kernels);
    262303            psFree(badStamps);
    … …  
    265306
    266307        for (int j = 0; j < numKernels; j++) {
    267             accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint);
    268         }
    269 
    270         lastChi2 += initialChi2(inputs->data[i], stamp, footprint);
     308            accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint, mode);
     309        }
     310
     311        lastChi2 += initialChi2(targets->data[i], stamp, footprint, mode);
    271312        numPixels += PS_SQR(2 * footprint + 1);
    272313    }
    … …  
    297338                }
    298339                pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
    299                 accumulateCross(&sumI, &sumII, &sumIC, stamp, inputs->data[j], i, footprint);
     340                accumulateCross(&sumI, &sumII, &sumIC, stamp, targets->data[j], i, footprint, mode);
    300341            }
    301342
    … …  
    310351                }
    311352                pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
    312                 chi2 += accumulateChi2(inputs->data[j], stamp, i, coeff, bg, footprint);
     353                chi2 += accumulateChi2(targets->data[j], stamp, i, coeff, bg, footprint, mode);
    313354            }
    314355
    … …  
    328369        if (bestIndex == -1) {
    329370            psError(PS_ERR_UNKNOWN, false, "Unable to find best kernel component in round %d.", iter);
    330             psFree(inputs);
     371            psFree(targets);
    331372            psFree(sumC);
    332373            psFree(sumCC);
    … …  
    345386            }
    346387            pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
    347             subtractConvolution(inputs->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint);
     388            subtractConvolution(targets->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint, mode);
    348389        }
    349390
    … …  
    361402        lastChi2 = bestChi2;
    362403    }
    363     psFree(inputs);
     404    psFree(targets);
    364405    psFree(sumC);
    365406    psFree(sumCC);
    … …  
    401442                pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
    402443                psArray *convolutions = convNew->data[j]; // Convolutions for this stamp
    403                 convolutions->data[rank] = psMemIncrRefCounter(stamp->convolutions->data[i]);
     444                convolutions->data[rank] = psMemIncrRefCounter(selectConvolution(stamp, i, mode));
    404445            }
    405446        }
    … …  
    420461        }
    421462        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    422         psFree(stamp->convolutions);
    423         stamp->convolutions = convNew->data[i];
     463        psFree(stamp->convolutions1);
     464        psFree(stamp->convolutions2);
     465        switch (mode) {
     466          case PM_SUBTRACTION_MODE_1:
     467            stamp->convolutions1 = convNew->data[i];
     468            stamp->convolutions2 = NULL;
     469            break;
     470          case PM_SUBTRACTION_MODE_2:
     471            stamp->convolutions1 = NULL;
     472            stamp->convolutions2 = convNew->data[i];
     473            break;
     474          default:
     475            psAbort("Unsupported subtraction mode: %x", mode);
     476        }
    424477    }
    425478
Note: See TracChangeset for help on using the changeset viewer.