IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 16, 2007, 1:51:51 PM (19 years ago)
Author:
Paul Price
Message:

Re-instituting pmImageCombine.[ch], to get ppStac working. Cleaned up and simplified pmCombineImages(); now using a function that does the combination for each pixel. Other functions are not yet tested!

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r10551 r11115  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2006-12-08 11:39:30 $
     10 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2007-01-16 23:51:51 $
    1212 *
    1313 *  XXX: pmRejectPixels() has a known bug with the pmImageTransform() call.
    … …  
    2626#include "pslib.h"
    2727
    28 //
    29 // The following macros define how big the initial pixel list will be, and
    30 // how much it should be incremented when realloc'ed.
    31 //
    32 #define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH 100
    33 #define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC 100
    34 #define PS_COMBINE_IMAGE_MAX_QUESTIONABLE_PIXELS 1000
    35 /******************************************************************************
    36 pmCombineImages(combine, questionablePixels, images, errors, masks, maskVal,
    37                 pixels, numIter, sigmaClip, stats)
    38  
    39 XXX: Allocate a dummy psStats structure so that we don't destroy away its data.
    40  *****************************************************************************/
     28#define PIXEL_LIST_BUFFER 100           // Size of the pixel list buffer
     29
     30// Data structure for use as a buffer in combining pixels
     31typedef struct
     32{
     33    psVector *pixels;                   // Pixel values
     34    psVector *masks;                    // Pixel masks
     35    psVector *errors;                   // Pixel errors
     36    psStats *stats;                     // Statistics to use with combination
     37}
     38combineBuffer;
     39
     40void combineBufferFree(combineBuffer *buffer)
     41{
     42    psFree(buffer->pixels);
     43    psFree(buffer->masks);
     44    psFree(buffer->errors);
     45    psFree(buffer->stats);
     46}
     47
     48combineBuffer *combineBufferAlloc(long numImages // Number of images that will be combined
     49                                 )
     50{
     51    combineBuffer *buffer = psAlloc(sizeof(combineBuffer));
     52    psMemSetDeallocator(buffer, (psFreeFunc)combineBufferFree);
     53
     54    buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32);
     55    buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK);
     56    buffer->errors = psVectorAlloc(numImages, PS_TYPE_F32);
     57    buffer->stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     58
     59    return buffer;
     60}
     61
     62
     63static bool combinePixels(psImage *combine, // Combined image, for output
     64                          psArray *questionablePixels, // Array of rejection masks
     65                          int x, int y, // Position in the images
     66                          const psArray *images, // Array of input images
     67                          const psArray *errors, // Array of input error images
     68                          const psArray *masks, // Array of input masks
     69                          psU32 maskVal, // Mask value
     70                          psS32 numIter, // Number of rejection iterations
     71                          psF32 sigmaClip, // Number of standard deviations at which to reject
     72                          combineBuffer *buffer // Buffer for combination; to avoid multiple allocations
     73                         )
     74{
     75    assert(combine);
     76    assert(x >= 0 && x < combine->numCols);
     77    assert(y >= 0 && y < combine->numRows);
     78    assert(images);
     79    int numImages = images->n;          // Number of images to combine
     80    if (masks) {
     81        assert(masks->n == numImages);
     82    }
     83    if (errors) {
     84        assert(errors->n == numImages);
     85    }
     86    assert(numIter >= 0);
     87    assert(sigmaClip > 0);
     88    assert(!questionablePixels || (questionablePixels && questionablePixels->n == numImages));
     89
     90    if (buffer) {
     91        psMemIncrRefCounter(buffer);
     92    } else {
     93        buffer = combineBufferAlloc(numImages);
     94    }
     95
     96    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
     97    psVector *pixelErrors = buffer->errors; // Errors for the pixel of interest
     98    psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest
     99    psStats *stats = buffer->stats;     // Statistics for combination
     100
     101    //
     102    // Loop through each image, extract the pixel/mask/error data into psVectors.
     103    //
     104    if (!masks) {
     105        psVectorInit(pixelMasks, 0);
     106    }
     107    if (!errors) {
     108        pixelErrors = NULL;
     109    }
     110    for (int i = 0; i < numImages; i++) {
     111        // Set the pixel data
     112        psImage *image = images->data[i]; // Image of interest
     113        pixelData->data.F32[i] = image->data.F32[y][x];
     114        // Set the pixel mask data, if necessary
     115        if (masks) {
     116            psImage *mask = masks->data[i]; // Mask of interest
     117            pixelMasks->data.U8[i] = mask->data.U8[y][x];
     118        }        // Set the pixel error data, if necessary
     119        if (errors) {
     120            psImage *error = errors->data[i]; // Error image of interest
     121            pixelErrors->data.F32[i] = error->data.F32[y][x];
     122        }
     123    }
     124
     125    //
     126    // Iterate on the pixels, rejecting outliers
     127    //
     128    for (int iter = 0; iter < numIter; iter++) {
     129        // Combine all the pixels, using the specified stat.
     130        if (!psVectorStats(stats, pixelData, pixelErrors, pixelMasks, maskVal)) {
     131            combine->data.F32[y][x] = NAN;
     132            psFree(buffer);
     133            return false;
     134        }
     135        float combinedPixel = stats->sampleMean; // Value of the combination
     136
     137        if (iter == 0) {
     138            // Save the value produced with no rejection, since it may be useful later
     139            // (if the rejection turns out to be unnecessary)
     140            combine->data.F32[y][x] = combinedPixel;
     141        }
     142
     143        //
     144        // Reject all pixels that lie more that sigmaClip standard deviations from
     145        // the combined pixel value.
     146        //
     147        int numRejects = 0;     // Number of rejections
     148        float stdev = stats->sampleStdev;
     149        for (int i = 0; i < numImages; i++) {
     150            if (!(pixelMasks->data.U8[i] & maskVal) &&
     151                    fabs(pixelData->data.F32[i] - combinedPixel) > sigmaClip * stdev) {
     152                // Reject pixel as questionable
     153                numRejects++;
     154                pixelMasks->data.U8[i] = maskVal;
     155                if (questionablePixels) {
     156                    // Mark the pixel as questionable
     157                    psPixels *qp = questionablePixels->data[i]; // Questionable pixels for this image
     158                    int qpNum = qp->n; // Number of QPs in the image of interest
     159                    if (qpNum >= qp->nalloc) {
     160                        // Grow dynamically, if required
     161                        qp = psPixelsRealloc(qp, qp->nalloc + PIXEL_LIST_BUFFER);
     162                        questionablePixels->data[i] = qp;
     163                    }
     164                    qp->data[qpNum].x = x;
     165                    qp->data[qpNum].y = y;
     166                    qp->n++;
     167                }
     168            }
     169        }
     170
     171        //
     172        // If the number of rejected pixels is zero, then there's no point continuing the loop.
     173        //
     174        if (numRejects == 0) {
     175            break;
     176        }
     177
     178        //
     179        // XXX: Is it possible to have all pixels rejected?  If so, we should exit the loop.
     180        //
     181    }
     182
     183    psFree(buffer);
     184    return true;
     185}
     186
    41187psImage *pmCombineImages(
    42188    psImage *combine,                   ///< Combined image (output)
    … …  
    48194    const psPixels *pixels,             ///< Pixels to combine
    49195    psS32 numIter,                      ///< Number of rejection iterations
    50     psF32 sigmaClip,                    ///< Number of standard deviations at which to reject
    51     psStats *stats)               ///< Statistics to use in the combination
    52 {
    53 
    54     PS_ASSERT_PTR_NON_NULL(images, combine);
    55     psU32 numImages = images->n;
     196    psF32 sigmaClip                     ///< Number of standard deviations at which to reject
     197)
     198{
    56199    psTrace("psModules.imcombine", 3, "Calling pmCombineImages(%ld)\n", images->n);
    57200
    58     if (errors != NULL) {
    59         if (images->n != errors->n) {
    60             psError(PS_ERR_UNKNOWN, true, "images and errors args must have same length (%ld != %ld)\n",
    61                     images->n, errors->n);
    62             return(combine);
    63         }
    64     }
    65     if (masks != NULL) {
    66         if (images->n != masks->n) {
    67             psError(PS_ERR_UNKNOWN, true, "images and masks args must have same length (%ld != %ld)\n",
    68                     images->n, masks->n);
    69             return(combine);
    70         }
    71     }
    72 
    73     psImage *tmpImg = (psImage *) images->data[0];
    74     psU32 numRows = tmpImg->numRows;
    75     psU32 numCols = tmpImg->numCols;
    76 
    77     //
    78     // Check that all images have the appropriate size and type.
    79     //
    80     for (psS32 i = 0 ; i < numImages ; i++) {
    81         psImage *tmpDataImg = (psImage *) images->data[i];
    82         PS_ASSERT_IMAGE_NON_NULL(tmpDataImg, combine);
    83         PS_ASSERT_IMAGE_TYPE(tmpDataImg, PS_TYPE_F32, combine);
    84         if ((tmpDataImg->numRows != numRows) || (tmpDataImg->numCols != numCols)) {
    85             psError(PS_ERR_UNKNOWN, true, "image %d has size (%d, %d); should be (%d, %d).\n",
    86                     i, tmpDataImg->numRows, tmpDataImg->numCols, numRows, numCols);
    87         }
    88 
    89         if (errors != NULL) {
    90             psImage *tmpErrorImg = (psImage *) errors->data[i];
    91             PS_ASSERT_IMAGE_NON_NULL(tmpErrorImg, combine);
    92             PS_ASSERT_IMAGE_TYPE(tmpErrorImg, PS_TYPE_F32, combine);
    93             PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpErrorImg, NULL);
    94         }
    95 
    96         if (masks != NULL) {
    97             psImage *tmpMaskImg = (psImage *) masks->data[i];
    98             PS_ASSERT_IMAGE_NON_NULL(tmpMaskImg, combine);
    99             PS_ASSERT_IMAGE_TYPE(tmpMaskImg, PS_TYPE_U8, combine);
    100             PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpMaskImg, NULL);
    101         }
    102     }
    103     PS_ASSERT_PTR_NON_NULL(stats, combine);
    104     psStatsOptions statistic = psStatsSingleOption(stats->options);
    105     if (statistic == 0) {
    106         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %x\n", statistic);
    107         return combine;
    108     }
     201    PS_ASSERT_ARRAY_NON_NULL(images, NULL);
     202    PS_ASSERT_INT_POSITIVE(images->n, NULL);
     203    long numImages = images->n;          // Number of images
     204    int numCols = ((psImage*)images->data[0])->numCols; // Size in x
     205    int numRows = ((psImage*)images->data[0])->numRows; // Size in y
     206
     207    if (combine) {
     208        PS_ASSERT_IMAGE_NON_NULL(combine, NULL);
     209        PS_ASSERT_IMAGE_SIZE(combine, numCols, numRows, NULL);
     210    }
     211    if (questionablePixels && !*questionablePixels) {
     212        PS_ASSERT_ARRAY_NON_NULL(*questionablePixels, NULL);
     213        PS_ASSERT_ARRAY_SIZE(*questionablePixels, numImages, 0);
     214    }
     215    for (int i = 1; i < images->n; i++) {
     216        psImage *image = images->data[i]; // Image of interest
     217        PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
     218        PS_ASSERT_IMAGE_SIZE(image, numCols, numRows, NULL);
     219    }
     220    if (errors) {
     221        PS_ASSERT_ARRAYS_SIZE_EQUAL(images, errors, NULL);
     222        for (int i = 0; i < images->n; i++) {
     223            psImage *error = errors->data[i];
     224            PS_ASSERT_IMAGE_SIZE(error, numCols, numRows, NULL);
     225            PS_ASSERT_IMAGE_TYPE(error, PS_TYPE_F32, NULL);
     226        }
     227    }
     228    if (masks) {
     229        PS_ASSERT_ARRAYS_SIZE_EQUAL(images, masks, NULL);
     230        for (int i = 0; i < images->n; i++) {
     231            psImage *mask  = masks->data[i];
     232            PS_ASSERT_IMAGE_SIZE(mask, numCols, numRows, NULL);
     233            PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL);
     234        }
     235    }
     236    PS_ASSERT_INT_POSITIVE(numIter, NULL);
     237    PS_ASSERT_FLOAT_LARGER_THAN(sigmaClip, 0.0, NULL);
    109238
    110239    // Allocate and initialize the combined image, if necessary.
    111     if (combine == NULL) {
     240    if (!combine) {
    112241        combine = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    113         if (pixels != NULL) {
    114             psImageInit(combine, 0.0);
     242        if (pixels) {
     243            // Set everything we're not working on to NAN
     244            psImageInit(combine, NAN);
    115245        }
    116246    }
    … …  
    120250    // struct for each image.
    121251    //
    122     if (*questionablePixels == NULL) {
    123         *questionablePixels = psArrayAlloc(numImages);
    124     } else if ((*questionablePixels)->n != numImages) {
    125         *questionablePixels = psArrayRealloc(*questionablePixels, numImages);
    126     }
    127     for (psS32 im = 0 ; im < numImages ; im++) {
    128         psFree((*questionablePixels)->data[im]);
    129         ((*questionablePixels)->data[im]) =
    130             (psPtr *) psPixelsAlloc(PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH);
    131         ((psPixels *) ((*questionablePixels)->data[im]))->n = 0;
    132     }
    133     //
    134     // qpPtr is used to maintain a count of the questionable pixels for each image.
    135     //
    136     psVector *qpPtr = psVectorAlloc(numImages, PS_TYPE_S32);
    137     psVectorInit(qpPtr, 0);
    138     //
    139     // Allocate the necessary psVectors for the call to psVectorStats().
    140     // These vectors will be used whether we are combining a list of pixels,
    141     // or every pixel in the input images.
    142     //
    143     psVector *pixelData = psVectorAlloc(numImages, PS_TYPE_F32);
    144     psVector *pixelMask = NULL;
    145     if (masks != NULL) {
    146         pixelMask = psVectorAlloc(numImages, PS_TYPE_U8);
    147         psVectorInit(pixelMask, 0);
    148     }
    149 
    150     psVector *pixelErrors = NULL;
    151     if (errors != NULL) {
    152         pixelErrors = psVectorAlloc(numImages, PS_TYPE_F32);
    153         psVectorInit(pixelErrors, 1.0);
    154     }
    155 
    156     if (pixels != NULL) {
     252    if (questionablePixels) {
     253        if (*questionablePixels == NULL) {
     254            *questionablePixels = psArrayAlloc(numImages);
     255        } else if ((*questionablePixels)->n != numImages) {
     256            *questionablePixels = psArrayRealloc(*questionablePixels, numImages);
     257        }
     258        for (int i = 0; i < numImages; i++) {
     259            psFree((*questionablePixels)->data[i]);
     260            (*questionablePixels)->data[i] = psPixelsAlloc(PIXEL_LIST_BUFFER);
     261        }
     262    }
     263
     264    combineBuffer *buffer = combineBufferAlloc(numImages); // Buffer for combination
     265
     266    if (pixels) {
    157267        // Only those specified pixels should be combined.
    158268
    159         psStats *stdevStats = psStatsAlloc(PS_STAT_SAMPLE_STDEV);
    160 
    161         for (psS32 p = 0 ; p < pixels->n ; p++) {
    162             // Must initialize the mask to 0 for every pixel.
    163             psVectorInit(pixelMask, 0);
    164             psS32 col = (pixels->data[p]).x;
    165             psS32 row = (pixels->data[p]).y;
    166 
    167             //
    168             // Loop through each image, extract the pixel/mask/error data
    169             // into psVectors.
    170             //
    171             for (psS32 im = 0 ; im < numImages ; im++) {
    172                 // Set the pixel data
    173                 pixelData->data.F32[im] =       ((psImage *) images->data[im])->data.F32[row][col];
    174                 // Set the pixel mask data, if necessary
    175                 if (masks != NULL) {
    176                     pixelMask->data.U8[im] =   ((psImage *) masks->data[im])->data.F32[row][col];
    177                 }
    178 
    179                 // Set the pixel error data, if necessary
    180                 if (errors != NULL) {
    181                     pixelErrors->data.F32[im] = ((psImage *) errors->data[im])->data.F32[row][col];
    182                 }
     269        for (int p = 0; p < pixels->n; p++) {
     270            int x = pixels->data[p].x; // Column of interest
     271            int y = pixels->data[p].y; // Row of interest
     272
     273            if (!combinePixels(combine, questionablePixels ? *questionablePixels : NULL, x, y,
     274                               images, errors, masks, maskVal, numIter, sigmaClip, buffer)) {
     275                // Bad pixel --- no big deal
     276                psErrorClear();
    183277            }
    184 
    185             //
    186             // Iterate on the pixels, rejecting outliers
    187             //
    188             for (psS32 iter = 0 ; iter < numIter ; iter++) {
    189                 // Combine all the pixels, using the specified stat.
    190 
    191                 psVectorStats(stats, pixelData, pixelErrors, pixelMask, maskVal);
    192                 psF64 combinedPixel = psStatsGetValue(stats, statistic);
    193 
    194                 if (iter == 0) {
    195                     combine->data.F32[row][col] = (psF32) combinedPixel;
    196                 }
    197 
    198                 //
    199                 // Reject all pixels that lie more that sigmaClip standard deviations from
    200                 // the combined pixel value.
    201                 //
    202                 psS32 numRejects = 0;
    203                 for (psS32 im = 0 ; im < numImages ; im++) {
    204                     psVectorStats(stdevStats, pixelData, pixelErrors, pixelMask, maskVal);
    205                     psF64 stdev = stdevStats->sampleStdev;
    206 
    207                     if (!(pixelMask->data.U8[im] & maskVal)) {
    208                         if (fabs(pixelData->data.F32[im]) >
    209                                 (fabs(sigmaClip * stdev) + fabs(combinedPixel))) {
    210                             numRejects++;
    211                             pixelMask->data.U8[im] = maskVal;
    212                             //
    213                             // XXX: These data structures indirections are getting complicated.
    214                             //
    215                             psS32 ptr = qpPtr->data.S32[im];
    216                             psPixels *pixelListPtr = ((psPixels *) ((*questionablePixels)->data[im]));
    217 
    218                             if (ptr >= pixelListPtr->nalloc) {
    219                                 (*questionablePixels)->data[im] =
    220                                     (psPtr *) psPixelsRealloc(((psPixels *) ((*questionablePixels)->data[im])),
    221                                                               ((((psPixels *) ((*questionablePixels)->data[im]))->nalloc) +
    222                                                                PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC));
    223                                 // XXX: Can the realloc() fail?  Must we check for NULL?
    224                             }
    225                             ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].x = col;
    226                             ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].y = row;
    227                             (qpPtr->data.S32[im])++;
    228                             // XXX: this pixel ->n increment is wierd
    229                             ((psPixels *) ((*questionablePixels)->data[im]))->n = qpPtr->data.S32[im];
    230                         }
    231                     }
    232                 }
    233 
    234                 //
    235                 // If the number of rejected pixels is zero, then there's no point
    236                 // continuing the loop.
    237                 //
    238                 if (numRejects == 0) {
    239                     break;
    240                 }
    241                 psS32 totalRejects = 0;
    242                 for (psS32 im = 0 ; im < numImages ; im++) {
    243                     if (pixelMask->data.U8[im] & maskVal) {
    244                         totalRejects++;
    245                     }
    246                 }
    247 
    248                 //
    249                 // XXX: Is it possible to have all pixels rejected?  If so, we should
    250                 // exit the loop.
    251                 //
    252                 if (totalRejects == numImages) {
    253                     break;
    254                 }
    255             }
    256         }
    257         psFree(stdevStats);
     278        }
    258279    } else {
    259280        //
    … …  
    267288        // combine image.
    268289        //
    269         for (psS32 row = 0 ; row < numRows ; row++) {
    270             for (psS32 col = 0 ; col < numCols ; col++) {
    271                 for (psS32 im = 0 ; im < numImages ; im++) {
    272                     // Set the pixel data
    273                     pixelData->data.F32[im] =       ((psImage *) images->data[im])->data.F32[row][col];
    274 
    275                     // Set the pixel mask data, if necessary
    276                     if (masks != NULL) {
    277                         pixelMask->data.U8[im] =   ((psImage *) masks->data[im])->data.F32[row][col];
    278                     }
    279 
    280                     // Set the pixel error data, if necessary
    281                     if (errors != NULL) {
    282                         pixelErrors->data.F32[im] = ((psImage *) errors->data[im])->data.F32[row][col];
    283                     }
    284 
     290        for (int y = 0; y < numRows; y++) {
     291            for (int x = 0; x < numCols; x++) {
     292                if (!combinePixels(combine, questionablePixels ? *questionablePixels : NULL, x, y,
     293                                   images, errors, masks, maskVal, numIter, sigmaClip, buffer)) {
     294                    // Bad pixel --- no big deal
     295                    psErrorClear();
    285296                }
    286                 // Combine all the pixels, using the specified stat.
    287                 psVectorStats(stats, pixelData, pixelErrors, pixelMask, maskVal);
    288                 combine->data.F32[row][col] = psStatsGetValue(stats, statistic);
    289297            }
    290298        }
    291299    }
    292300
    293     psFree(pixelData);
    294     psFree(pixelMask);
    295     psFree(pixelErrors);
    296     psFree(qpPtr);
     301    psFree(buffer);
    297302
    298303    psTrace("psModules.imcombine", 3, "Exiting pmCombineImages(%ld)\n", images->n);
    299     return(combine);
     304    return combine;
    300305}
    301306
Note: See TracChangeset for help on using the changeset viewer.