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

    r15286 r15443  
    44
    55#include <stdio.h>
     6#include <string.h>
    67#include <pslib.h>
    78
    … …  
    4546    psFree(stamp->matrix);
    4647    psFree(stamp->vector);
    47     psFree(stamp->reference);
    48     psFree(stamp->input);
     48    psFree(stamp->image1);
     49    psFree(stamp->image2);
    4950    psFree(stamp->weight);
    50     psFree(stamp->convolutions);
     51    psFree(stamp->convolutions1);
     52    psFree(stamp->convolutions2);
    5153}
    5254
    … …  
    6567// Is this position unmasked?
    6668static bool checkStampMask(int x, int y, // Coordinates of stamp
    67                            const psImage *mask
     69                           const psImage *mask, // Mask
     70                           pmSubtractionMode mode // Mode for subtraction
    6871                           )
    6972{
    … …  
    7477        return false;
    7578    }
    76     return (mask->data.PS_TYPE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER |
    77                                                   PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_REJ)) ?
    78         false : true;
     79
     80    psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ; // Mask value
     81    switch (mode) {
     82      case PM_SUBTRACTION_MODE_1:
     83      case PM_SUBTRACTION_MODE_TARGET:
     84        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1;
     85        break;
     86      case PM_SUBTRACTION_MODE_2:
     87        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_2;
     88        break;
     89      case PM_SUBTRACTION_MODE_UNSURE:
     90      case PM_SUBTRACTION_MODE_DUAL:
     91        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1 | PM_SUBTRACTION_MASK_FOOTPRINT_2;
     92        break;
     93      default:
     94        psAbort("Unsupported subtraction mode: %x", mode);
     95    }
     96
     97    return (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? false : true;
    7998}
    8099
    … …  
    140159    stamp->status = PM_SUBTRACTION_STAMP_INIT;
    141160
    142     stamp->reference = NULL;
    143     stamp->input = NULL;
     161    stamp->image1 = NULL;
     162    stamp->image2 = NULL;
    144163    stamp->weight = NULL;
    145     stamp->convolutions = NULL;
     164    stamp->convolutions1 = NULL;
     165    stamp->convolutions2 = NULL;
    146166
    147167    return stamp;
    … …  
    151171pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image,
    152172                                                const psImage *subMask, const psRegion *region,
    153                                                 float threshold, float spacing)
     173                                                float threshold, float spacing, pmSubtractionMode mode)
    154174{
    155175    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    … …  
    226246            for (int y = subRegion->y0; y <= subRegion->y1 ; y++) {
    227247                for (int x = subRegion->x0; x <= subRegion->y1 ; x++) {
    228                     if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxStamp) {
     248                    if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) {
    229249                        fluxStamp = image->data.F32[y][x];
    230250                        xStamp = x;
    … …  
    242262
    243263            // Reset the postage stamps since we're making a new stamp
    244             psFree(stamp->reference);
    245             psFree(stamp->input);
     264            psFree(stamp->image1);
     265            psFree(stamp->image2);
    246266            psFree(stamp->weight);
    247             stamp->reference = stamp->input = stamp->weight = NULL;
    248 
    249             stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
     267            psFree(stamp->convolutions1);
     268            psFree(stamp->convolutions2);
     269            stamp->image1 = stamp->image2 = stamp->weight = NULL;
     270            stamp->convolutions1 = stamp->convolutions2 = NULL;
     271
     272            stamp->status = PM_SUBTRACTION_STAMP_FOUND;
    250273            numFound++;
    251274            psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n",
    … …  
    266289pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux,
    267290                                               const psImage *image, const psImage *subMask,
    268                                                const psRegion *region, float spacing, int exclusionZone)
     291                                               const psRegion *region, float spacing, pmSubtractionMode mode)
     292
    269293{
    270294    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
    … …  
    289313    }
    290314    PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL);
    291     PS_ASSERT_INT_NONNEGATIVE(exclusionZone, NULL);
    292315
    293316    int numStars = x->n;                // Number of stars
    294     psVector *exclude = psVectorAlloc(numStars, PS_TYPE_U8); // Exclude a star?
    295     psVectorInit(exclude, 0);
    296 
    297     // Apply exclusion zone around each star --- important when we're convolving to a specified PSF; less so
    298     // when we're trying to get a good subtraction.
    299     if (exclusionZone > 0) {
    300         psTrace("psModules.imcombine", 2, "Applying exclusion zone of %d pixels for stamps\n", exclusionZone);
    301         // We use something resembling a binary search --- coordinates are sorted in the x dimension, and then
    302         // to exclude stars within a nominated distance, we need only examine (i.e., calculate the
    303         // 2-dimensional distance to) other stars in the list that are within that distance of the x
    304         // coordinate of the star of interest.  By marking both stars when we find a violation of the
    305         // exclusion zone, we need only search upwards from the x coordinate of the star of interest.
    306 
    307         int minDistance2 = PS_SQR(exclusionZone); // Minimum square distance for other stars
    308         psVector *indexes = psVectorSortIndex(NULL, x); // Indices for sorting in x
    309         for (int i = 0; i < numStars - 1; i++) {
    310             int iIndex = indexes->data.S32[i]; // Index for star i
    311             if (exclude->data.U8[iIndex]) {
    312                 continue;
    313             }
    314             float ix = x->data.F32[iIndex], iy = y->data.F32[iIndex]; // Coordinates for star i
    315             float jx, jy;                   // Coordinates for star j
    316             for (int j = i + 1, jIndex = indexes->data.S32[j];
    317                  j < numStars && (jx = x->data.F32[jIndex]) < ix + exclusionZone;
    318                  j++, jIndex = indexes->data.S32[j]) {
    319                 jy = y->data.F32[jIndex];
    320                 if (PS_SQR(ix - jx) + PS_SQR(iy - jy) < minDistance2) {
    321                     exclude->data.U8[iIndex] = 0xff;
    322                     exclude->data.U8[jIndex] = 0xff;
    323                     psTrace("psModules.imcombine", 5, "Excluding stamps %d,%d and %d,%d\n",
    324                             (int)x->data.F32[iIndex], (int)y->data.F32[iIndex],
    325                             (int)x->data.F32[jIndex], (int)y->data.F32[jIndex]);
    326                     // Would 'break' here, but there might be more stamps within the exclusion zone.
    327                 }
    328             }
    329         }
    330         psFree(indexes);
    331     }
    332 
    333317    pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
    334318                                                                 region, spacing); // Stamp list
    … …  
    347331    // Put the stars into their appropriate subregions
    348332    for (int i = 0; i < numStars; i++) {
    349         if (exclude->data.U8[i]) {
    350             // Star has been excluded
    351             continue;
    352         }
    353333        float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp
    354334        int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp
    … …  
    357337            continue;
    358338        }
    359         if (!checkStampMask(xPix, yPix, subMask)) {
     339        if (!checkStampMask(xPix, yPix, subMask, mode)) {
    360340            // Not a good stamp
    361341            continue;
    … …  
    386366        }
    387367    }
    388     psFree(exclude);
    389368
    390369    // Sort the list by flux, with the brightest last
    … …  
    421400
    422401
    423 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *reference, psImage *input,
     402bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
    424403                                psImage *weight, int footprint, int kernelSize)
    425404{
    426405    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    427     PS_ASSERT_IMAGE_NON_NULL(reference, false);
    428     PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false);
    429     if (input) {
    430         PS_ASSERT_IMAGE_NON_NULL(input, false);
    431         PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);
    432         PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false);
    433     }
    434     if (weight) {
    435         PS_ASSERT_IMAGE_NON_NULL(weight, false);
    436         PS_ASSERT_IMAGES_SIZE_EQUAL(weight, reference, false);
    437         PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
    438     }
     406    PS_ASSERT_IMAGE_NON_NULL(image1, false);
     407    PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, false);
     408    if (image2) {
     409        PS_ASSERT_IMAGE_NON_NULL(image2, false);
     410        PS_ASSERT_IMAGES_SIZE_EQUAL(image2, image1, false);
     411        PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false);
     412    }
     413    PS_ASSERT_IMAGE_NON_NULL(weight, false);
     414    PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false);
     415    PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
    439416    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
    440417    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
    441418
    442     int numCols = reference->numCols, numRows = reference->numRows; // Size of images
     419    int numCols = image1->numCols, numRows = image1->numRows; // Size of images
    443420    int size = kernelSize + footprint; // Size of postage stamps
    444 
    445     if (!weight) {
    446         // Use the input (or if I must, the reference) as a rough approximation to the variance map, and HOPE
    447         // that it's positive.
    448         weight = input ? input : reference;
    449     }
    450421
    451422    for (int i = 0; i < stamps->num; i++) {
    452423        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    453         if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
     424        if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_FOUND) {
    454425            continue;
    455426        }
    … …  
    468439        }
    469440
     441        // Catch memory leaks --- these should have been freed and NULLed before
     442        assert(stamp->image1 == NULL);
     443        assert(stamp->image2 == NULL);
     444        assert(stamp->weight == NULL);
     445
    470446        psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest
    471447
    472         psImage *refSub = psImageSubset(reference, region); // Subimage with stamp
    473         stamp->reference = psKernelAllocFromImage(refSub, size, size);
    474         psFree(refSub);                 // Drop reference
    475 
    476         if (input) {
    477             psImage *inSub = psImageSubset(input, region); // Subimage with stamp
    478             stamp->input = psKernelAllocFromImage(inSub, size, size);
    479             psFree(inSub);              // Drop reference
     448        psImage *sub1 = psImageSubset(image1, region); // Subimage with stamp
     449        stamp->image1 = psKernelAllocFromImage(sub1, size, size);
     450        psFree(sub1);                   // Drop reference
     451
     452        if (image2) {
     453            psImage *sub2 = psImageSubset(image2, region); // Subimage with stamp
     454            stamp->image2 = psKernelAllocFromImage(sub2, size, size);
     455            psFree(sub2);               // Drop reference
    480456        }
    481457
    482458        psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp
    483459        stamp->weight = psKernelAllocFromImage(wtSub, size, size);
    484         psFree(wtSub);                 // Drop reference
     460        psFree(wtSub);                  // Drop reference
     461
     462        stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
    485463    }
    486464
    … …  
    488466}
    489467
    490 
     468#if 0
    491469bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize)
    492470{
    … …  
    517495        float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame
    518496
    519         psFree(stamp->input);
    520         stamp->input = psKernelAlloc(-size, size, -size, size);
    521         psKernel *input = stamp->input; // Target stamp
     497        psFree(stamp->image2);
     498        stamp->image2 = psKernelAlloc(-size, size, -size, size);
     499        psKernel *target = stamp->image2; // Target stamp
    522500
    523501        // Put in a Waussian, just for fun!
    … …  
    525503            for (int u = -size; u <= size; u++) {
    526504                float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma));
    527                 input->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));
     505                target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));
    528506            }
    529507        }
    … …  
    533511    return true;
    534512}
    535 
     513#endif
    536514
    537515pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask,
    538516                                                          const psRegion *region, float spacing,
    539                                                           int exclusionZone)
     517                                                          pmSubtractionMode mode)
    540518{
    541519    PS_ASSERT_ARRAY_NON_NULL(sources, NULL);
    … …  
    561539
    562540    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region,
    563                                                             spacing, exclusionZone); // Stamps to return
     541                                                            spacing, mode); // Stamps to return
    564542    psFree(x);
    565543    psFree(y);
    … …  
    572550    return stamps;
    573551}
     552
     553
     554pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *subMask,
     555                                                       const psRegion *region, float spacing,
     556                                                       pmSubtractionMode mode)
     557{
     558    PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
     559    // Let pmSubtractionStampsSet take care of the rest of the assertions
     560
     561    const char *format = (mode == PM_SUBTRACTION_MODE_TARGET ? "%f %f %f" : "%f %f"); // Format of file
     562    psArray *data = psVectorsReadFromFile(filename, format);
     563    if (!data) {
     564        psError(PS_ERR_IO, false, "Unable to read stamps file %s", filename);
     565        return NULL;
     566    }
     567    psVector *x = data->data[0], *y = data->data[1]; // Stamp positions
     568    psVector *flux = (mode == PM_SUBTRACTION_MODE_TARGET ? data->data[2] : NULL); // Stamp fluxes
     569
     570    // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
     571    psBinaryOp(x, x, "-", psScalarAlloc(1.0, PS_TYPE_F32));
     572    psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32));
     573
     574    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, spacing,
     575                                                            mode);
     576    psFree(data);
     577
     578    return stamps;
     579
     580}
Note: See TracChangeset for help on using the changeset viewer.