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

    r25120 r26035  
    5454    psFree(stamp->image1);
    5555    psFree(stamp->image2);
    56     psFree(stamp->variance);
     56    psFree(stamp->weight);
    5757    psFree(stamp->convolutions1);
    5858    psFree(stamp->convolutions2);
    … …  
    180180
    181181pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region,
    182                                                     int footprint, float spacing)
     182                                                    int footprint, float spacing, float sysErr)
    183183{
    184184    pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return
    … …  
    221221    list->flux = NULL;
    222222    list->footprint = footprint;
     223    list->sysErr = sysErr;
    223224
    224225    return list;
    … …  
    256257        outStamp->image1 = inStamp->image1 ? psKernelCopy(inStamp->image1) : NULL;
    257258        outStamp->image2 = inStamp->image2 ? psKernelCopy(inStamp->image2) : NULL;
    258         outStamp->variance = inStamp->variance ? psKernelCopy(inStamp->variance) : NULL;
     259        outStamp->weight = inStamp->weight ? psKernelCopy(inStamp->weight) : NULL;
    259260
    260261        if (inStamp->convolutions1) {
    … …  
    305306    stamp->image1 = NULL;
    306307    stamp->image2 = NULL;
    307     stamp->variance = NULL;
     308    stamp->weight = NULL;
    308309    stamp->convolutions1 = NULL;
    309310    stamp->convolutions2 = NULL;
    … …  
    322323                                                const psImage *image2, const psImage *subMask,
    323324                                                const psRegion *region, float thresh1, float thresh2,
    324                                                 int size, int footprint, float spacing,
     325                                                int size, int footprint, float spacing, float sysErr,
    325326                                                pmSubtractionMode mode)
    326327{
    … …  
    379380
    380381    if (!stamps) {
    381         stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing);
     382        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr);
    382383    }
    383384
    … …  
    452453                psFree(stamp->image1);
    453454                psFree(stamp->image2);
    454                 psFree(stamp->variance);
     455                psFree(stamp->weight);
    455456                psFree(stamp->convolutions1);
    456457                psFree(stamp->convolutions2);
    457                 stamp->image1 = stamp->image2 = stamp->variance = NULL;
     458                stamp->image1 = stamp->image2 = stamp->weight = NULL;
    458459                stamp->convolutions1 = stamp->convolutions2 = NULL;
    459460
    … …  
    487488                                               const psImage *image, const psImage *subMask,
    488489                                               const psRegion *region, int size, int footprint,
    489                                                float spacing, pmSubtractionMode mode)
     490                                               float spacing, float sysErr, pmSubtractionMode mode)
    490491
    491492{
    … …  
    507508    int numStars = x->n;                // Number of stars
    508509    pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
    509                                                                  region, footprint, spacing); // Stamp list
     510                                                                 region, footprint, sysErr,
     511                                                                 spacing); // Stamp list
    510512    int numStamps = stamps->num;        // Number of stamps
    511513
    … …  
    616618        PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false);
    617619    }
    618     PS_ASSERT_IMAGE_NON_NULL(variance, false);
    619     PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false);
    620     PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false);
    621     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
     620    if (variance) {
     621        PS_ASSERT_IMAGE_NON_NULL(variance, false);
     622        PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false);
     623        PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false);
     624        PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
     625    }
    622626
    623627    int numCols = image1->numCols, numRows = image1->numRows; // Size of images
    … …  
    646650        assert(stamp->image1 == NULL);
    647651        assert(stamp->image2 == NULL);
    648         assert(stamp->variance == NULL);
     652        assert(stamp->weight == NULL);
    649653
    650654        psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest
    … …  
    660664        }
    661665
    662         psImage *wtSub = psImageSubset(variance, region); // Subimage with stamp
    663         stamp->variance = psKernelAllocFromImage(wtSub, size, size);
    664         psFree(wtSub);                  // Drop reference
     666        psKernel *weight = stamp->weight = psKernelAlloc(-size, size, -size, size); // Weight = 1/variance
     667        if (variance) {
     668            psImage *varSub = psImageSubset(variance, region); // Subimage with stamp
     669            psKernel *var = psKernelAllocFromImage(varSub, size, size); // Variance postage stamp
     670            if (isfinite(stamps->sysErr) && stamps->sysErr > 0) {
     671                float sysErr = 0.5 * stamps->sysErr; // Systematic error
     672                psKernel *image1 = stamp->image1, *image2 = stamp->image2; // Input images
     673                for (int y = -size; y <= size; y++) {
     674                    for (int x = -size; x <= size; x++) {
     675                        weight->kernel[y][x] = 1.0 / (var->kernel[y][x] +
     676                                                      sysErr * (image1->kernel[y][x] + image2->kernel[y][x]));
     677                    }
     678                }
     679            } else {
     680                for (int y = -size; y <= size; y++) {
     681                    for (int x = -size; x <= size; x++) {
     682                        weight->kernel[y][x] = 1.0 / var->kernel[y][x];
     683                    }
     684                }
     685            }
     686            psFree(var);
     687            psFree(varSub);
     688        } else {
     689            psImageInit(weight->image, 1.0);
     690        }
    665691
    666692        stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
    … …  
    672698pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
    673699                                                          const psImage *subMask, const psRegion *region,
    674                                                           int size, int footprint, float spacing,
     700                                                          int size, int footprint, float spacing, float sysErr,
    675701                                                          pmSubtractionMode mode)
    676702{
    … …  
    702728
    703729    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size,
    704                                                             footprint, spacing, mode); // Stamps to return
     730                                                            footprint, spacing, sysErr, mode); // Stamps
    705731    psFree(x);
    706732    psFree(y);
    … …  
    716742pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image,
    717743                                                       const psImage *subMask, const psRegion *region,
    718                                                        int size, int footprint, float spacing,
     744                                                       int size, int footprint, float spacing, float sysErr,
    719745                                                       pmSubtractionMode mode)
    720746{
    … …  
    734760
    735761    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint,
    736                                                             spacing, mode);
     762                                                            spacing, sysErr, mode);
    737763    psFree(data);
    738764
Note: See TracChangeset for help on using the changeset viewer.