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/pmSubtractionEquation.c

    r26031 r26035  
    1717//#define TESTING                         // TESTING output for debugging; may not work with threads!
    1818
    19 #define USE_VARIANCE                    // Include variance in equation?
     19#define USE_WEIGHT                      // Include weight (1/variance) in equation?
    2020
    2121
    … …  
    2929                                  const psKernel *input, // Input image (target)
    3030                                  const psKernel *reference, // Reference image (convolution source)
    31                                   const psKernel *variance,  // Variance image
     31                                  const psKernel *weight,  // Weight image
    3232                                  const psArray *convolutions,         // Convolutions for each kernel
    3333                                  const pmSubtractionKernels *kernels, // Kernels
    … …  
    7474                for (int x = - footprint; x <= footprint; x++) {
    7575                    double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
    76 #ifdef USE_VARIANCE
    77                     cc /= variance->kernel[y][x];
     76#ifdef USE_WEIGHT
     77                    cc *= weight->kernel[y][x];
    7878#endif
    7979                    sumCC += cc;
    … …  
    102102                double rc = ref * conv;
    103103                double c = conv;
    104 #ifdef USE_VARIANCE
    105                 float varVal = 1.0 / variance->kernel[y][x];
    106                 ic *= varVal;
    107                 rc *= varVal;
    108                 c *= varVal;
     104#ifdef USE_WEIGHT
     105                float wtVal = weight->kernel[y][x];
     106                ic *= wtVal;
     107                rc *= wtVal;
     108                c *= wtVal;
    109109#endif
    110110                sumIC += ic;
    … …  
    137137            double rr = PS_SQR(ref);
    138138            double one = 1.0;
    139 #ifdef USE_VARIANCE
    140             float varVal = 1.0 / variance->kernel[y][x];
    141             rr *= varVal;
    142             ir *= varVal;
    143             in *= varVal;
    144             ref *= varVal;
    145             one *= varVal;
     139#ifdef USE_WEIGHT
     140            float wtVal = weight->kernel[y][x];
     141            rr *= wtVal;
     142            ir *= wtVal;
     143            in *= wtVal;
     144            ref *= wtVal;
     145            one *= wtVal;
    146146#endif
    147147            sumRR += rr;
    … …  
    169169                                      const psKernel *image1, // Image 1
    170170                                      const psKernel *image2, // Image 2
    171                                       const psKernel *variance,  // Variance image
     171                                      const psKernel *weight,  // Weight image
    172172                                      const psArray *convolutions1, // Convolutions of image 1 for each kernel
    173173                                      const psArray *convolutions2, // Convolutions of image 2 for each kernel
    … …  
    227227                    double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
    228228                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
    229 #ifdef USE_VARIANCE
    230                     aa /= variance->kernel[y][x];
    231                     bb /= variance->kernel[y][x];
    232                     ab /= variance->kernel[y][x];
     229#ifdef USE_WEIGHT
     230                    float wtVal = weight->kernel[y][x];
     231                    aa *= wtVal;
     232                    bb *= wtVal;
     233                    ab *= wtVal;
    233234#endif
    234235                    sumAA += aa;
    … …  
    258259                for (int x = - footprint; x <= footprint; x++) {
    259260                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
    260 #ifdef USE_VARIANCE
    261                     ab /= variance->kernel[y][x];
     261#ifdef USE_WEIGHT
     262                    ab *= weight->kernel[y][x];
    262263#endif
    263264                    sumAB += ab;
    … …  
    295296                double i1i2 = i1 * i2;
    296297
    297 #ifdef USE_VARIANCE
    298                 float varVal = 1.0 / variance->kernel[y][x];
    299                 ai2 *= varVal;
    300                 bi2 *= varVal;
    301                 ai1 *= varVal;
    302                 bi1 *= varVal;
    303                 i1i2 *= varVal;
    304                 a *= varVal;
    305                 b *= varVal;
    306                 i2 *= varVal;
     298#ifdef USE_WEIGHT
     299                float wtVal = weight->kernel[y][x];
     300                ai2 *= wtVal;
     301                bi2 *= wtVal;
     302                ai1 *= wtVal;
     303                bi1 *= wtVal;
     304                i1i2 *= wtVal;
     305                a *= wtVal;
     306                b *= wtVal;
     307                i2 *= wtVal;
    307308#endif
    308309
    … …  
    351352            double i1i2 = i1 * i2;
    352353
    353 #ifdef USE_VARIANCE
    354             float varVal = 1.0 / variance->kernel[y][x];
    355             i1 *= varVal;
    356             i1i1 *= varVal;
    357             one *= varVal;
    358             i2 *= varVal;
    359             i1i2 *= varVal;
     354#ifdef USE_WEIGHT
     355            float wtVal = weight->kernel[y][x];
     356            i1 *= wtVal;
     357            i1i1 *= wtVal;
     358            one *= wtVal;
     359            i2 *= wtVal;
     360            i1i2 *= wtVal;
    360361#endif
    361362
    … …  
    619620      case PM_SUBTRACTION_MODE_1:
    620621        status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
    621                                        stamp->variance, stamp->convolutions1, kernels, polyValues,
     622                                       stamp->weight, stamp->convolutions1, kernels, polyValues,
    622623                                       footprint);
    623624        break;
    624625      case PM_SUBTRACTION_MODE_2:
    625626        status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
    626                                        stamp->variance, stamp->convolutions2, kernels, polyValues,
     627                                       stamp->weight, stamp->convolutions2, kernels, polyValues,
    627628                                       footprint);
    628629        break;
    … …  
    639640#endif
    640641        status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2,
    641                                            stamp->matrixX, stamp->image1, stamp->image2, stamp->variance,
     642                                           stamp->matrixX, stamp->image1, stamp->image2, stamp->weight,
    642643                                           stamp->convolutions1, stamp->convolutions2, kernels, polyValues,
    643644                                           footprint);
    … …  
    11631164
    11641165        // Calculate residuals
    1165         psKernel *variance = stamp->variance; // Variance postage stamp
     1166        psKernel *weight = stamp->weight; // Weight postage stamp
    11661167        psImageInit(residual->image, 0.0);
    11671168        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
    … …  
    12471248        for (int y = - footprint; y <= footprint; y++) {
    12481249            for (int x = - footprint; x <= footprint; x++) {
    1249                 double dev = PS_SQR(residual->kernel[y][x]) / variance->kernel[y][x];
     1250                double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
    12501251                deviation += dev;
    12511252#ifdef TESTING
    … …  
    12901291            psFitsClose(fits);
    12911292        }
    1292         if (stamp->variance) {
     1293        if (stamp->weight) {
    12931294            psString filename = NULL;
    1294             psStringAppend(&filename, "stamp_variance_%03d.fits", i);
     1295            psStringAppend(&filename, "stamp_weight_%03d.fits", i);
    12951296            psFits *fits = psFitsOpen(filename, "w");
    12961297            psFree(filename);
    1297             psFitsWriteImage(fits, NULL, stamp->variance->image, 0, NULL);
     1298            psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
    12981299            psFitsClose(fits);
    12991300        }
Note: See TracChangeset for help on using the changeset viewer.