IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 4, 2009, 4:52:05 PM (17 years ago)
Author:
Paul Price
Message:

Soften errors, both in the least-squares matrix/vector construction and in the rejection stage. Softening in the rejection stage seems to be very important so that bright stamps aren't rejected (after which the solution can go crazy). Now carrying around a weight image = 1/variance instead of the variance image. This can be softened, saving many divisions in the matrix/vector construction.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/psModules/src/imcombine/pmSubtractionParams.c

    r21363 r26035  
    7171                            double *sumII, // Sum of I(x)^2/sigma(x)^2
    7272                            double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2
    73                             const pmSubtractionStamp *stamp, // Stamp with variance
     73                            const pmSubtractionStamp *stamp, // Stamp
    7474                            const psKernel *target, // Target stamp
    7575                            int kernelIndex, // Index for kernel component
     
    7878    )
    7979{
    80     psKernel *variance = stamp->variance;   // Variance, sigma(x)^2
     80    psKernel *weight = stamp->weight;   // Weight image
    8181    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    8282
    8383    for (int y = -footprint; y <= footprint; y++) {
    8484        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    85         psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
     85        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    8686        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    8787        for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) {
    88             double temp = *in / *wt; // Temporary product
     88            double temp = *in * *wt; // Temporary product
    8989            *sumI += temp;
    9090            *sumII += *in * temp;
     
    9898static void accumulateConvolutions(double *sumC, // Sum of conv(x)/sigma(x)^2
    9999                                   double *sumCC, // Sum of conv(x)^2/sigma(x)^2
    100                                    const pmSubtractionStamp *stamp, // Stamp with input and variance
     100                                   const pmSubtractionStamp *stamp, // Stamp with input and weight
    101101                                   int kernelIndex, // Index for kernel component
    102102                                   int footprint, // Size of region of interest
     
    104104    )
    105105{
    106     psKernel *variance = stamp->variance;   // Variance, sigma(x)^2
     106    psKernel *weight = stamp->weight;   // Weight image
    107107    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    108108
    109109    for (int y = -footprint; y <= footprint; y++) {
    110         psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
     110        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    111111        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    112112        for (int x = -footprint; x <= footprint; x++, wt++, conv++) {
    113             double convNoise = *conv / *wt; // Temporary product
     113            double convNoise = *conv * *wt; // Temporary product
    114114            *sumC += convNoise;
    115115            *sumCC += *conv * convNoise;
     
    120120
    121121static double accumulateChi2(const psKernel *target, // Target stamp
    122                              pmSubtractionStamp *stamp, // Stamp with variance
     122                             pmSubtractionStamp *stamp, // Stamp with weight
    123123                             int kernelIndex, // Index for kernel component
    124124                             double coeff, // Coefficient of convolution
     
    129129{
    130130    double chi2 = 0.0;
    131     psKernel *variance = stamp->variance;   // Variance, sigma(x)^2
     131    psKernel *weight = stamp->weight;   // Weight image
    132132    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    133133
    134134    for (int y = -footprint; y <= footprint; y++) {
    135135        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    136         psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
     136        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    137137        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    138138        for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) {
    139             chi2 += PS_SQR(*in - bg - coeff * *conv) / *wt;
     139            chi2 += PS_SQR(*in - bg - coeff * *conv) * *wt;
    140140        }
    141141    }
     
    146146// Return the initial value of chi^2
    147147static double initialChi2(const psKernel *target, // Target stamp
    148                           const pmSubtractionStamp *stamp, // Stamp with variance
     148                          const pmSubtractionStamp *stamp, // Stamp
    149149                          int footprint, // Size of convolution
    150150                          pmSubtractionMode mode // Mode of subtraction
    151151    )
    152152{
    153     psKernel *variance = stamp->variance;   // Variance map
     153    psKernel *weight = stamp->weight;   // Weight image
    154154    psKernel *source;                   // Source stamp
    155155    switch (mode) {
     
    167167    for (int y = -footprint; y <= footprint; y++) {
    168168        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    169         psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
     169        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    170170        psF32 *ref = &source->kernel[y][-footprint]; // Derference reference
    171171        for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) {
    172172            float diff = *in - *ref;    // Temporary value
    173             chi2 += PS_SQR(diff) / *wt;
     173            chi2 += PS_SQR(diff) * *wt;
    174174        }
    175175    }
     
    180180// Subtract a convolution from the input
    181181static void subtractConvolution(psKernel *target, // Target stamp
    182                                 const pmSubtractionStamp *stamp, // Stamp with variance
     182                                const pmSubtractionStamp *stamp, // Stamp
    183183                                int kernelIndex, // Index for kernel component
    184184                                float coeff, // Coefficient of subtraction
     
    288288
    289289        // This sum is invariant to the kernel
    290         psKernel *variance = stamp->variance; // Variance map for stamp
     290        psKernel *weight = stamp->weight; // Weight image
    291291        for (int v = -footprint; v <= footprint; v++) {
    292             psF32 *wt = &variance->kernel[v][-footprint]; // Dereference variance map
     292            psF32 *wt = &weight->kernel[v][-footprint]; // Dereference weight
    293293            for (int u = -footprint; u <= footprint; u++, wt++) {
    294                 sum1 += 1.0 / *wt;
     294                sum1 += 1.0 * *wt;
    295295            }
    296296        }
Note: See TracChangeset for help on using the changeset viewer.