IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 10, 2007, 10:10:05 AM (19 years ago)
Author:
Paul Price
Message:

Preparing to take a (long) list of sources as candidate stamps. Changed the list of stamps from a simple array to a more complicated structure so that I can carry around a list of candidate stamps for each region of interest; when one is rejected, it is replaced by the next brightest within the same region. Optionally apply an exclusion zone around stamps, which is important when we're convolving to a specified PSF (as opposed to convolving to match an image).

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r14701 r14801  
    99#include "pmSubtractionStamps.h"
    1010
     11#define STAMP_LIST_BUFFER 20            // Number of stamps to add to list at a time
     12
    1113//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    1214// Private (file-static) functions
    1315//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    1416
     17// Free function for pmSubtractionStampList
     18static void subtractionStampListFree(pmSubtractionStampList *list // Stamp list to free
     19                                     )
     20{
     21    psFree(list->stamps);
     22    psFree(list->regions);
     23    psFree(list->x);
     24    psFree(list->y);
     25    psFree(list->flux);
     26}
     27
    1528// Free function for pmSubtractionStamp
    16 void subtractionStampFree(pmSubtractionStamp *stamp // Stamp to free
    17     )
     29static void subtractionStampFree(pmSubtractionStamp *stamp // Stamp to free
     30                                 )
    1831{
    1932    psFree(stamp->matrix);
    … …  
    2841static bool checkStampRegion(int x, int y, // Coordinates of stamp
    2942                             const psRegion *region // Region of interest
    30     )
     43                             )
    3144{
    3245    if (!region) {
    … …  
    4053static bool checkStampMask(int x, int y, // Coordinates of stamp
    4154                           const psImage *mask
    42     )
     55                           )
    4356{
    4457    if (!mask) {
    … …  
    5467//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    5568
    56 pmSubtractionStamp *pmSubtractionStampAlloc(pmSubtractionStampStatus status)
     69pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region,
     70                                                    float spacing)
     71{
     72    pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return
     73    psMemSetDeallocator(list, (psFreeFunc)subtractionStampListFree);
     74
     75    // Get region in which to find stamps: [xMin:xMax,yMin:yMax]
     76    int xMin = 0, xMax = numCols, yMin = 0, yMax = numRows;
     77    if (region) {
     78        xMin = PS_MAX(region->x0, xMin);
     79        xMax = PS_MIN(region->x1, xMax);
     80        yMin = PS_MAX(region->y0, yMin);
     81        yMax = PS_MIN(region->y1, yMax);
     82    }
     83    int xSize = xMax - xMin, ySize = yMax - yMin; // Size of region of interest
     84    int xStamps = xSize / spacing + 1, yStamps = ySize / spacing + 1; // Number of stamps in x and y
     85
     86    list->num = xStamps * yStamps;
     87    list->stamps = psArrayAlloc(list->num);
     88    list->regions = psArrayAlloc(list->num);
     89
     90    for (int y = 0, index = 0; y < yStamps; y++) {
     91        int yStart = yMin + y * ySize / (yStamps + 1); // Subregion starts here
     92        int yStop = yMin + (y + 1) * ySize / (yStamps + 1) - 1; // Subregion stops here
     93        assert(yStart >= yMin && yStop < yMax);
     94
     95        for (int x = 0; x < xStamps; x++, index++) {
     96            int xStart = xMin + x * xSize / (xStamps + 1); // Subregion starts here
     97            int xStop = xMin + (x + 1) * xSize / (xStamps + 1) - 1; // Subregion stops here
     98            assert(xStart >= xMin && xStop < xMax);
     99
     100            list->stamps->data[index] = pmSubtractionStampAlloc();
     101            list->regions->data[index] = psRegionAlloc(xStart, xStop, yStart, yStop);
     102        }
     103    }
     104
     105    list->x = NULL;
     106    list->y = NULL;
     107    list->flux = NULL;
     108
     109    return list;
     110}
     111
     112pmSubtractionStamp *pmSubtractionStampAlloc(void)
    57113{
    58114    pmSubtractionStamp *stamp = psAlloc(sizeof(pmSubtractionStamp)); // Stamp to return
    … …  
    66122    stamp->matrix = NULL;
    67123    stamp->vector = NULL;
    68     stamp->status = status;
     124    stamp->status = PM_SUBTRACTION_STAMP_INIT;
    69125
    70126    stamp->reference = NULL;
    … …  
    77133
    78134
    79 psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *subMask,
    80                                  const psRegion *region, float threshold, float spacing)
     135pmSubtractionStampList *pmSubtractionFindStamps(pmSubtractionStampList *stamps, const psImage *image,
     136                                                const psImage *subMask, const psRegion *region,
     137                                                float threshold, float spacing)
    81138{
    82139    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    … …  
    107164    int numRows = image->numRows, numCols = image->numCols; // Size of image
    108165
    109     // Get region in which to find stamps: [xMin:xMax,yMin:yMax]
    110     int xMin = 0, xMax = numCols, yMin = 0, yMax = numRows;
    111     if (region) {
    112         xMin = PS_MAX(region->x0, xMin);
    113         xMax = PS_MIN(region->x1, xMax);
    114         yMin = PS_MAX(region->y0, yMin);
    115         yMax = PS_MIN(region->y1, yMax);
    116     }
    117     int xSize = xMax - xMin, ySize = yMax - yMin; // Size of region of interest
    118 
    119     // Number of stamps
    120     int xNumStamps = xSize / spacing + 1;
    121     int yNumStamps = ySize / spacing + 1;
    122 
     166    if (!stamps) {
     167        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, spacing);
     168    }
     169
     170    int numStamps = stamps->num;        // Number of stamp regions
    123171    int numFound = 0;                   // Number of stamps found
    124 
    125     if (stamps) {
    126         PS_ASSERT_INT_EQUAL(stamps->n, xNumStamps * yNumStamps, NULL);
    127         // Just in case, check for NULL stamps
    128         for (int i = 0; i < xNumStamps * yNumStamps; i++) {
    129             if (!stamps->data[i]) {
    130                 stamps->data[i] = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_REJECTED);
     172    for (int i = 0; i < numStamps; i++) {
     173        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     174
     175        // Only find a new stamp if we need to
     176        if (stamp->status != PM_SUBTRACTION_STAMP_REJECTED && stamp->status != PM_SUBTRACTION_STAMP_INIT) {
     177            continue;
     178        }
     179
     180        float xStamp = 0, yStamp = 0;   // Coordinates of stamp
     181        float fluxStamp = NAN;          // Flux of stamp
     182        bool goodStamp = false;         // Found a good stamp?
     183
     184        // A couple different ways of finding stamps:
     185        if (stamps->x && stamps->y) {
     186            // Get the next stamp from the list
     187            psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Coordinate lists
     188            psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes
     189
     190            // Take stamp off the top of the (sorted) list
     191            if (xList->n > 0) {
     192                int index = xList->n - 1; // Index of new stamp
     193                xStamp = xList->data.F32[index];
     194                yStamp = yList->data.F32[index];
     195                fluxStamp = fluxList->data.F32[index];
     196
     197                // Chop off the top of the list
     198                xList->n = index;
     199                yList->n = index;
     200                fluxList->n = index;
     201
     202                goodStamp = true;
    131203            }
    132         }
    133     } else {
    134         stamps = psArrayAlloc(xNumStamps * yNumStamps);
    135         for (int i = 0; i < xNumStamps * yNumStamps ; i++) {
    136             stamps->data[i] = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_INIT);
    137         }
    138     }
    139 
    140     for (int j = 0, index = 0; j < yNumStamps; j++) {
    141         for (int i = 0; i < xNumStamps; i++, index++) {
    142             pmSubtractionStamp *stamp = stamps->data[index]; // Stamp of interest
    143 
    144             // Only find a new stamp if we need to
    145             if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED ||
    146                 stamp->status == PM_SUBTRACTION_STAMP_INIT) {
    147 
    148                 // Find maximum non-masked value in the image section,
    149                 // but don't include a footprint around the edge
    150                 float fluxBest = threshold; // Flux of best stamp pixel
    151                 int xBest = 0, yBest = 0; // Coordinates of best stamp
    152 
    153                 // Bounds of region to search for stamp
    154                 int yStart = yMin + j * ySize / (yNumStamps + 1);
    155                 int yStop = yMin + (j + 1) * ySize / (yNumStamps + 1) - 1;
    156                 int xStart = xMin + i * xSize / (xNumStamps + 1);
    157                 int xStop = xMin + (i + 1) * xSize / (xNumStamps + 1) - 1;
    158                 assert(yStop < numRows && xStop < numCols && yStart >= 0 && xStart >= 0);
    159 
    160                 for (int y = yStart; y <= yStop ; y++) {
    161                     for (int x = xStart; x <= xStop ; x++) {
    162                         if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxBest) {
    163                             fluxBest = image->data.F32[y][x];
    164                             xBest = x;
    165                             yBest = y;
    166                         }
     204        } else {
     205            // Use a simple method of automatically finding stamps --- take the highest pixel in the subregion
     206            fluxStamp = threshold;
     207            psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest
     208            for (int y = subRegion->y0; y <= subRegion->y1 ; y++) {
     209                for (int x = subRegion->x0; x <= subRegion->y1 ; x++) {
     210                    if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxStamp) {
     211                        fluxStamp = image->data.F32[y][x];
     212                        xStamp = x;
     213                        yStamp = y;
     214                        goodStamp = true;
    167215                    }
    168216                }
    169 
    170                 stamp->x = xBest;
    171                 stamp->y = yBest;
    172                 stamp->flux = fluxBest;
    173 
    174                 // Reset the postage stamps since we're making a new stamp
    175                 psFree(stamp->reference);
    176                 psFree(stamp->input);
    177                 psFree(stamp->weight);
    178                 stamp->reference = stamp->input = stamp->weight = NULL;
    179 
    180                 if (fluxBest > threshold) {
    181                     stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
    182                     numFound++;
    183                     psTrace("psModules.imcombine", 5, "Found stamp in region (%d,%d): %d,%d\n",
    184                             i, j, (int)stamp->x, (int)stamp->y);
    185                 } else {
    186                     stamp->status = PM_SUBTRACTION_STAMP_NONE;
    187                 }
    188217            }
    189218        }
    190     }
    191 
    192     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps in %dx%d regions",
    193              numFound, xNumStamps, yNumStamps);
     219
     220        if (goodStamp) {
     221            stamp->x = xStamp;
     222            stamp->y = yStamp;
     223            stamp->flux = fluxStamp;
     224
     225            // Reset the postage stamps since we're making a new stamp
     226            psFree(stamp->reference);
     227            psFree(stamp->input);
     228            psFree(stamp->weight);
     229            stamp->reference = stamp->input = stamp->weight = NULL;
     230
     231            stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
     232            numFound++;
     233            psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n",
     234                    i, (int)stamp->x, (int)stamp->y);
     235        } else {
     236            stamp->status = PM_SUBTRACTION_STAMP_NONE;
     237        }
     238    }
     239
     240    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps", numFound);
    194241
    195242    return stamps;
    … …  
    197244
    198245
    199 psArray *pmSubtractionSetStamps(const psVector *x, const psVector *y, const psVector *flux,
    200                                 const psImage *subMask, const psRegion *region)
     246pmSubtractionStampList *pmSubtractionSetStamps(const psVector *x, const psVector *y, const psVector *flux,
     247                                               const psImage *image, const psImage *subMask,
     248                                               const psRegion *region, float spacing, int exclusionZone)
    201249{
    202250    PS_ASSERT_VECTOR_NON_NULL(x, false);
    … …  
    209257        PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, false);
    210258        PS_ASSERT_VECTORS_SIZE_EQUAL(flux, x, false);
     259    } else {
     260        PS_ASSERT_IMAGE_NON_NULL(image, false);
    211261    }
    212262    if (subMask) {
    213263        PS_ASSERT_IMAGE_NON_NULL(subMask, false);
    214264        PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false);
    215     }
    216 
    217     int num = x->n;                     // Number of stamps
    218     psArray *stamps = psArrayAllocEmpty(num); // Stamps, to return
    219 
    220     for (int i = 0; i < num; i++) {
    221         int xStamp = x->data.F32[i] + 0.5, yStamp = y->data.F32[i] + 0.5; // Pixel coords of stamp
    222         if (!checkStampRegion(xStamp, yStamp, region)) {
    223             continue;
    224         }
    225         if (!checkStampMask(xStamp, yStamp, subMask)) {
    226             continue;
    227         }
    228 
    229         pmSubtractionStamp *stamp = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_CALCULATE); // Stamp
    230         stamp->x = x->data.F32[i];
    231         stamp->y = y->data.F32[i];
    232         if (flux) {
    233             stamp->flux = flux->data.F32[i];
    234         }
    235         stamps->data[stamps->n++] = stamp;
    236     }
     265        if (image) {
     266            PS_ASSERT_IMAGE_NON_NULL(image, false);
     267            PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, false);
     268        }
     269    }
     270
     271    int numStars = x->n;                // Number of stars
     272    psVector *exclude = psVectorAlloc(numStars, PS_TYPE_U8); // Exclude a star?
     273    psVectorInit(exclude, 0);
     274
     275    // Apply exclusion zone around each star --- important when we're convolving to a specified PSF; less so
     276    // when we're trying to get a good subtraction.
     277    if (exclusionZone > 0) {
     278        psTrace("psModules.imcombine", 2, "Applying exclusion zone of %d pixels for stamps\n", exclusionZone);
     279        // We use something resembling a binary search --- coordinates are sorted in the x dimension, and then
     280        // to exclude stars within a nominated distance, we need only examine (i.e., calculate the
     281        // 2-dimensional distance to) other stars in the list that are within that distance of the x
     282        // coordinate of the star of interest.  By marking both stars when we find a violation of the
     283        // exclusion zone, we need only search upwards from the x coordinate of the star of interest.
     284
     285        int minDistance2 = PS_SQR(exclusionZone); // Minimum square distance for other stars
     286        psVector *indexes = psVectorSortIndex(NULL, x); // Indices for sorting in x
     287        for (int i = 0; i < numStars; i++) {
     288            int iIndex = indexes->data.S32[i]; // Index for star i
     289            if (exclude->data.U8[iIndex]) {
     290                continue;
     291            }
     292            float ix = x->data.F32[iIndex], iy = y->data.F32[iIndex]; // Coordinates for star i
     293            float jx, jy;                   // Coordinates for star j
     294            for (int j = i + 1, jIndex = indexes->data.S32[j];
     295                 (jx = x->data.F32[jIndex]) < ix + exclusionZone && j < numStars;
     296                 j++, jIndex = indexes->data.S32[j]) {
     297                jy = y->data.F32[jIndex];
     298                if (PS_SQR(ix - jx) + PS_SQR(iy - jy) < minDistance2) {
     299                    exclude->data.U8[iIndex] = 0xff;
     300                    exclude->data.U8[jIndex] = 0xff;
     301                    psTrace("psModules.imcombine", 5, "Excluding stamps %d,%d and %d,%d\n",
     302                            (int)x->data.F32[iIndex], (int)y->data.F32[iIndex],
     303                            (int)x->data.F32[jIndex], (int)y->data.F32[jIndex]);
     304                    // Would 'break' here, but there might be more stamps within the exclusion zone.
     305                }
     306            }
     307        }
     308        psFree(indexes);
     309    }
     310
     311    pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
     312                                                                 region, spacing); // Stamp list
     313    int numStamps = stamps->num;        // Number of stamps
     314
     315    // Initialise the lists
     316    stamps->x = psArrayAlloc(numStamps);
     317    stamps->y = psArrayAlloc(numStamps);
     318    stamps->flux = psArrayAlloc(numStamps);
     319    for (int i = 0; i < numStamps; i++) {
     320        stamps->x->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32);
     321        stamps->y->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32);
     322        stamps->flux->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32);
     323    }
     324
     325    // Put the stars into their appropriate subregions
     326    for (int i = 0; i < numStars; i++) {
     327        if (exclude->data.U8[i]) {
     328            // Star has been excluded
     329            continue;
     330        }
     331        float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp
     332        int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp
     333        if (!checkStampRegion(xPix, yPix, region)) {
     334            // It's not in the big region
     335            continue;
     336        }
     337        if (!checkStampMask(xPix, yPix, subMask)) {
     338            // Not a good stamp
     339            continue;
     340        }
     341
     342        for (int j = 0; j < numStamps; j++) {
     343            psRegion *subRegion = stamps->regions->data[j]; // Subregion of interest
     344            if (checkStampRegion(xPix, yPix, subRegion)) {
     345                psVector *xList = stamps->x->data[j], *yList = stamps->y->data[j]; // Pixel lists
     346                psVector *fluxList = stamps->flux->data[j]; // Flux list
     347
     348                int index = xList->n;   // Index of new stamp candidate
     349                xList->data.F32[index] = xStamp;
     350                yList->data.F32[index] = yStamp;
     351
     352                if (flux) {
     353                    fluxList->data.F32[index] = flux->data.F32[i];
     354                } else {
     355                    fluxList->data.F32[index] = image->data.F32[yPix][xPix];
     356                }
     357
     358                psVectorExtend(xList, STAMP_LIST_BUFFER, 1);
     359                psVectorExtend(yList, STAMP_LIST_BUFFER, 1);
     360                psVectorExtend(fluxList, STAMP_LIST_BUFFER, 1);
     361
     362                break;
     363            }
     364        }
     365    }
     366    psFree(exclude);
     367
     368    // Sort the list by flux, with the brightest last
     369    for (int i = 0; i < numStamps; i++) {
     370        psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Pixel lists
     371        psVector *fluxList = stamps->flux->data[i]; // Flux list
     372
     373        psVector *indexes = psVectorSortIndex(NULL, fluxList); // Indices to sort flux
     374        int num = indexes->n;           // Number of candidate stamps in this subregion
     375
     376        psVector *xSorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of x list
     377        psVector *ySorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of y list
     378        psVector *fluxSorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of flux list
     379        for (int j = 0; j < num; j++) {
     380            int k = indexes->data.S32[j]; // Sorted index
     381            xSorted->data.F32[j] = xList->data.F32[k];
     382            ySorted->data.F32[j] = yList->data.F32[k];
     383            fluxSorted->data.F32[j] = fluxList->data.F32[k];
     384        }
     385        psFree(indexes);
     386
     387        psFree(stamps->x->data[i]);
     388        psFree(stamps->y->data[i]);
     389        psFree(stamps->flux->data[i]);
     390
     391        stamps->x->data[i] = xSorted;
     392        stamps->y->data[i] = ySorted;
     393        stamps->flux->data[i] = fluxSorted;
     394    }
     395
    237396
    238397    return stamps;
    … …  
    240399
    241400
    242 bool pmSubtractionExtractStamps(psArray *stamps, psImage *reference, psImage *input, psImage *weight,
    243                                 int footprint, int kernelSize)
    244 {
    245     PS_ASSERT_ARRAY_NON_NULL(stamps, false);
     401bool pmSubtractionExtractStamps(pmSubtractionStampList *stamps, psImage *reference, psImage *input,
     402                                psImage *weight, int footprint, int kernelSize)
     403{
     404    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    246405    PS_ASSERT_IMAGE_NON_NULL(reference, false);
    247406    PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false);
    … …  
    263422
    264423    if (!weight) {
    265         // Use the input as a rough approximation to the variance map, and HOPE that it's positive.
     424        // Use the input (or if I must, the reference) as a rough approximation to the variance map, and HOPE
     425        // that it's positive.
    266426        weight = input ? input : reference;
    267427    }
    268428
    269     for (int i = 0; i < stamps->n; i++) {
    270         pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
     429    for (int i = 0; i < stamps->num; i++) {
     430        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    271431        if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
    272432            continue;
    … …  
    307467
    308468
    309 bool pmSubtractionGenerateStamps(psArray *stamps, float fwhm, int footprint, int kernelSize)
    310 {
    311     PS_ASSERT_ARRAY_NON_NULL(stamps, false);
     469bool pmSubtractionGenerateStamps(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize)
     470{
     471    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    312472    PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, false);
    313473    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
     474    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
    314475
    315476    int size = kernelSize + footprint; // Size of postage stamps
    316     int num = stamps->n;                // Number of stamps
     477    int num = stamps->num;              // Number of stamps
    317478    float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma
    318479
    319480    for (int i = 0; i < num; i++) {
    320         pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
     481        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    321482        if (!(stamp->status & PM_SUBTRACTION_STAMP_CALCULATE)) {
    322483            continue;
    … …  
    347508    }
    348509
    349     return stamps;
    350 }
     510    return true;
     511}
Note: See TracChangeset for help on using the changeset viewer.