IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 23, 2007, 12:25:41 PM (19 years ago)
Author:
Paul Price
Message:

Bug fixes and algorithmic improvements to get it to a working state. Seems to work OK now, though hasn't been tested rigourously.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r13457 r13487  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2007-05-22 03:59:32 $
     10 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2007-05-23 22:25:41 $
    1212 *
    1313 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    4848}
    4949
    50 static combineBuffer *combineBufferAlloc(long numImages // Number of images that will be combined
     50static combineBuffer *combineBufferAlloc(long numImages, // Number of images that will be combined
     51                                         psStatsOptions stat // Statistic to use
    5152    )
    5253{
    … …  
    5657    buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32);
    5758    buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK);
    58     buffer->stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     59    buffer->stats = psStatsAlloc(stat);
    5960
    6061    return buffer;
    … …  
    7879                          const psArray *inputs, // Stack data
    7980                          const psVector *weights, // Weights for data, or NULL
     81                          const psVector *reject, // Indices of pixels to reject, or NULL
    8082                          int x, int y, // Coordinates of interest
    8183                          psMaskType maskVal, // Value to mask
    … …  
    9698    int num = inputs->n;          // Number of images to combine
    9799
    98     if (buffer) {
    99         psMemIncrRefCounter(buffer);
    100     } else {
    101         buffer = combineBufferAlloc(num);
    102     }
    103 
    104100    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
    105101    psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest
    106102    psStats *stats = buffer->stats;     // Statistics
     103
    107104
    108105    // Extract the pixel and mask data
    … …  
    114111        pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[y][x];
    115112    }
     113    if (reject) {
     114        for (int i = 0; i < reject->n; i++) {
     115            pixelMasks->data.PS_TYPE_MASK_DATA[reject->data.U16[i]] |= maskVal;
     116        }
     117    }
    116118
    117119    // Record the value derived with no clipping, because pixels rejected using the harsh clipping applied in
    118120    // the first pass might later be accepted.
    119     if (!psVectorStats(stats, pixelData, pixelMasks, weights, maskVal)) {
     121    stats->options |= PS_STAT_SAMPLE_MEAN;
     122    if (!psVectorStats(stats, pixelData, weights, pixelMasks, maskVal)) {
    120123        // Can't do anything about an error except give it a NAN and mask
    121124        psErrorClear();
    … …  
    127130    mask->data.PS_TYPE_MASK_DATA[y][x] = 0;
    128131
     132    stats->options &= ~ PS_STAT_SAMPLE_MEAN;
     133
    129134    long numClipped = LONG_MAX;         // Number of pixels clipped
    130135    psMaskType ignore = maskVal | bad;  // Ignore values with this mask value
    131136    for (int i = 0; i < numIter && numClipped > 0; i++) {
     137#if 1
     138        float mean = stats->sampleMedian;
     139        float stdev = 0.74 * (stats->sampleUQ - stats->sampleLQ); // Rough estimate of the standard deviation
     140#else
     141        float mean = stats->robustMedian;
     142        float stdev = stats->robustStdev;
     143#endif
     144        float limit = rej * stdev; // Rejection limit
    132145        numClipped = 0;
    133         float limit = rej * stats->sampleStdev; // Rejection limit
    134146        for (int j = 0; j < num; j++) {
    135147            if (pixelMasks->data.PS_TYPE_MASK_DATA[j] & ignore) {
    136148                continue;
    137149            }
    138             float diff = fabsf(pixelData->data.F32[j] - stats->sampleMean); // Absolute difference from mean
    139             if (diff > limit) {
    140                 pmStackData *data = inputs->data[i]; // Stack data of interest
     150            float diff = pixelData->data.F32[j] - mean;
     151            if (fabsf(diff) > limit) {
     152                pmStackData *data = inputs->data[j]; // Stack data of interest
    141153                data->pixels = psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y);
    142154                pixelMasks->data.PS_TYPE_MASK_DATA[j] |= bad;
    … …  
    144156            }
    145157        }
    146         if (!psVectorStats(stats, pixelData, pixelMasks, weights, maskVal)) {
     158        if (!psVectorStats(stats, pixelData, weights, pixelMasks, maskVal)) {
    147159            // Can't do anything about an error except give it a NAN and mask
    148160            psErrorClear();
    … …  
    286298    }
    287299
    288     if (!psVectorStats(stats, values, valuesMask, NULL, MASK_BAD | MASK_SUSPECT)) {
     300    if (!psVectorStats(stats, values, NULL, valuesMask, MASK_BAD | MASK_SUSPECT)) {
    289301        psFree(stats);
    290302        psFree(values);
    … …  
    397409    PS_ASSERT_INT_NONNEGATIVE(numIter, false);
    398410    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
    399 
    400 
     411    if (havePixels) {
     412        // This is a subsequent combination, so expect that the image and mask already exist
     413        PS_ASSERT_IMAGE_NON_NULL(combined->image, false);
     414        PS_ASSERT_IMAGE_TYPE(combined->image, PS_TYPE_F32, false);
     415        PS_ASSERT_IMAGE_NON_NULL(combined->mask, false);
     416        PS_ASSERT_IMAGE_TYPE(combined->mask, PS_TYPE_MASK, false);
     417        PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false);
     418    }
    401419    if (!haveDetector && !haveSky) {
    402420        psError(PS_ERR_UNEXPECTED_NULL, true, "Nothing to combine!");
    403421        return false;
    404422    }
     423
    405424
    406425    if (!haveSky) {
    … …  
    415434    }
    416435
    417     // Pull out the weights
     436    // Pull out the weightings
    418437    psVector *weights = psVectorAlloc(num, PS_TYPE_F32);
    419438    for (int i = 0; i < num; i++) {
    … …  
    422441    }
    423442
    424     combineBuffer *buffer = combineBufferAlloc(num); // Buffer for combination
     443    // Buffer for combination
     444    combineBuffer *buffer = combineBufferAlloc(num, numIter == 0 ? PS_STAT_SAMPLE_MEAN :
     445                                               PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN);
    425446
    426447    if (havePixels) {
    427448        // Only combine select pixels
    428 
     449        psImage *combinedImage = combined->image; // Combined image
     450        psImage *combinedMask = combined->mask; // Combined mask
     451
     452        psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source
     453        psPixels *pixels = NULL;            // Total list of pixels, with no duplicates
     454        for (int i = 0; i < num; i++) {
     455            pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels
     456            pixels = psPixelsConcatenate(pixels, data->pixels);
     457            data->pixels = psPixelsRealloc(data->pixels, PIXEL_LIST_BUFFER); // Just in case more rejection
     458        }
     459        for (int i = 0; i < pixels->n; i++) {
     460            int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest
     461            psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely
     462            combinePixels(combinedImage, combinedMask, input, weights, reject, x, y,
     463                          maskVal, bad, numIter, rej, buffer);
     464        }
     465        psFree(pixels);
     466        psFree(pixelMap);
     467        psTrace("psModules.imcombine", 5, "Additional %ld pixels fixed.\n", pixels->n);
    429468    } else {
    430469        // Pull the products out, allocate if necessary
    … …  
    440479        }
    441480
     481        // Generate the pixel lists in which to place the rejected pixels
     482        if (numIter != 0) {
     483            for (int i = 0; i < num; i++) {
     484                pmStackData *data = input->data[i]; // Stack data for this input
     485                data->pixels = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
     486            }
     487        }
     488
    442489        for (int y = 0; y < numRows; y++) {
    443490            for (int x = 0; x < numCols; x++) {
    444                 combinePixels(combinedImage, combinedMask, input, weights, x, y,
     491                combinePixels(combinedImage, combinedMask, input, weights, NULL, x, y,
    445492                              maskVal, bad, numIter, rej, buffer);
    446493            }
    447494        }
    448495
     496        if (psTraceGetLevel("psModules.imcombine") >= 5) {
     497            for (int i = 0; i < num; i++) {
     498                pmStackData *data = input->data[i]; // Stack data for this input
     499                psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->pixels->n);
     500            }
     501        }
    449502    }
    450503
    … …  
    492545    }
    493546    (void)psBinaryOp(seeing, seeing, "*", seeing);
    494     (void)psBinaryOp(seeing, seeing, "-", psScalarAlloc(PS_SQR(seeingMax), PS_TYPE_F32));
     547    (void)psBinaryOp(seeing, psScalarAlloc(PS_SQR(seeingMax), PS_TYPE_F32), "-", seeing);
    495548    (void)psUnaryOp(seeing, seeing, "sqrt");
    496549
    … …  
    504557
    505558        // Can daisy-chain multiple tests here
     559        reject->n = 0;
    506560
    507561        if (!convolveTest(reject, inspect, input, x, y, seeing, maskVal, extent, threshold)) {
    … …  
    519573            psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y);
    520574        }
    521 
     575    }
     576
     577    if (psTraceGetLevel("psModules.imcombine") >= 5) {
     578        for (int i = 0; i < num; i++) {
     579            pmStackData *data = input->data[i]; // Stack data for this input
     580            psTrace("psModules.imcombine", 5, "Image %d: %ld pixels rejected.\n", i, data->pixels->n);
     581        }
    522582    }
    523583
    … …  
    525585}
    526586
    527 
    528 
    529 
    530 
    531 
    532 
    533 
    534 
    535 #if 0
    536 /******************************************************************************
    537 XXX: Directly from Paul Price
    538  *****************************************************************************/
    539 static psF32 CalcGradient(
    540     psImage *image,
    541     psImage *imageMask,
    542     psS32 x,
    543     psS32 y
    544 )
    545 {
    546     psTrace("psModules.imcombine", 4, "Calling CalcGradient(%d, %d)\n", x, y);
    547     int num = 0;
    548     psVector *pixels = psVectorAlloc(8, PS_TYPE_F32); // Array of pixels
    549     psVector *mask = psVectorAlloc(8, PS_TYPE_U8); // Corresponding mask
    550 
    551     // Get limits
    552     int xMin = PS_MAX(x - 1, 0);
    553     int xMax = PS_MIN(x + 1, image->numCols - 1);
    554     int yMin = PS_MAX(y - 1, 0);
    555     int yMax = PS_MIN(y + 1, image->numRows - 1);
    556     if (imageMask != NULL) {
    557         for (int j = yMin; j <= yMax; j++) {
    558             for (int i = xMin; i <= xMax; i++) {
    559                 if ((i != x) && (j != y) && (0 == imageMask->data.U8[j][i])) {
    560                     pixels->data.F32[num] = image->data.F32[j][i];
    561                     mask->data.U8[num] = 0;
    562                     num++;
    563                 } else {
    564                     mask->data.U8[num] = 1;
    565                 }
    566             }
    567         }
    568     } else {
    569         //
    570         // This code is simply the previous loop without the imageMask.
    571         // XXX: Consider restructuring this.
    572         //
    573         for (int j = yMin; j <= yMax; j++) {
    574             for (int i = xMin; i <= xMax; i++) {
    575                 if ((i != x) && (j != y)) {
    576                     pixels->data.F32[num] = image->data.F32[j][i];
    577                     mask->data.U8[num] = 0;
    578                     num++;
    579                 } else {
    580                     mask->data.U8[num] = 1;
    581                 }
    582             }
    583         }
    584     }
    585 
    586     pixels->n = num;
    587     mask->n = num;
    588 
    589     // Get the median
    590     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    591     psVectorStats(stats, pixels, NULL, mask, 1);
    592     float median = stats->sampleMedian;
    593     psFree(stats);
    594     psFree(pixels);
    595     psFree(mask);
    596 
    597     psTrace("psModules.imcombine", 4, "Exiting CalcGradient(%d, %d)\n", x, y);
    598     return(median / image->data.F32[y][x]);
    599 }
    600 
    601 /******************************************************************************
    602 DetermineRegion(image, myOutToIn): for a psImage and a psPlaneTransform to that
    603 image, this routine determines the size of the input image which maps to that
    604 image, and returns the result in a psRegion struct.
    605 
    606 XXX: Basically, this routine is only guaranteed to work if the transform is
    607 linear.
    608 
    609 XXX: Shouldn't this functionality be part of psImageTransform()?
    610  *****************************************************************************/
    611 static psRegion DetermineRegion(psImage *image,
    612                                 psPlaneTransform *myOutToIn)
    613 {
    614     psTrace("psModules.imcombine", 4, "Calling DetermineRegion()\n");
    615     psRegion myRegion;
    616     myRegion.x0 = PS_MAX_F32;
    617     myRegion.x1 = PS_MIN_F32;
    618     myRegion.y0 = PS_MAX_F32;
    619     myRegion.y1 = PS_MIN_F32;
    620     psPlane in;
    621     psPlane out;
    622 
    623     in.x = 0.0;
    624     in.y = 0.0;
    625 
    626     psPlaneTransformApply(&out, myOutToIn, &in);
    627     if (out.x < myRegion.x0) {
    628         myRegion.x0 = out.x;
    629     }
    630     if (out.x > myRegion.x1) {
    631         myRegion.x1 = out.x;
    632     }
    633     if (out.y < myRegion.y0) {
    634         myRegion.y0 = out.y;
    635     }
    636     if (out.y > myRegion.y1) {
    637         myRegion.y1 = out.y;
    638     }
    639 
    640     in.x = (psF32) (image->numCols);
    641     in.y = 0.0;
    642     psPlaneTransformApply(&out, myOutToIn, &in);
    643     if (out.x < myRegion.x0) {
    644         myRegion.x0 = out.x;
    645     }
    646     if (out.x > myRegion.x1) {
    647         myRegion.x1 = out.x;
    648     }
    649     if (out.y < myRegion.y0) {
    650         myRegion.y0 = out.y;
    651     }
    652     if (out.y > myRegion.y1) {
    653         myRegion.y1 = out.y;
    654     }
    655 
    656     in.x = (psF32) (image->numCols);
    657     ;
    658     in.y = 0.0;
    659     psPlaneTransformApply(&out, myOutToIn, &in);
    660     if (out.x < myRegion.x0) {
    661         myRegion.x0 = out.x;
    662     }
    663     if (out.x > myRegion.x1) {
    664         myRegion.x1 = out.x;
    665     }
    666     if (out.y < myRegion.y0) {
    667         myRegion.y0 = out.y;
    668     }
    669     if (out.y > myRegion.y1) {
    670         myRegion.y1 = out.y;
    671     }
    672 
    673     in.x = (psF32) (image->numCols);
    674     in.y = (psF32) (image->numRows);
    675     psPlaneTransformApply(&out, myOutToIn, &in);
    676     if (out.x < myRegion.x0) {
    677         myRegion.x0 = out.x;
    678     }
    679     if (out.x > myRegion.x1) {
    680         myRegion.x1 = out.x;
    681     }
    682     if (out.y < myRegion.y0) {
    683         myRegion.y0 = out.y;
    684     }
    685     if (out.y > myRegion.y1) {
    686         myRegion.y1 = out.y;
    687     }
    688 
    689     psTrace("psModules.imcombine", 4, "Exiting DetermineRegion()\n");
    690     return(myRegion);
    691 }
    692 
    693 /******************************************************************************
    694 XXX: Don't we have a psLib function for this?
    695  *****************************************************************************/
    696 static psImage *ImageConvertF32(psImage *image)
    697 {
    698     psTrace("psModules.imcombine", 4, "Calling ImageConvertF32()\n");
    699     psImage *imgF32 = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
    700 
    701     for (psS32 i = 0 ; i < image->numRows ; i++) {
    702         for (psS32 j = 0 ; j < image->numCols ; j++) {
    703             imgF32->data.F32[i][j] = (psF32) image->data.U8[i][j];
    704         }
    705     }
    706 
    707     psTrace("psModules.imcombine", 4, "Exiting ImageConvertF32()\n");
    708     return(imgF32);
    709 }
    710 
    711 
    712 //
    713 // The following macros define how big the initial pixel list will be, and
    714 // how much it should be incremented when realloc'ed.
    715 //
    716 #define PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH 100
    717 #define PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC 100
    718 /******************************************************************************
    719 pmRejectPixels(images, errors, inToOut, outToIn, rejThreshold,
    720 gradLimit)
    721 
    722 XXX: Optimization: we don't need to transform the entire mask image.
    723 XXX: The inToOut and outToIn transforms are confusing.  Verify that what
    724      I think they mean syncs with PWP.
    725  *****************************************************************************/
    726 psArray *pmRejectPixels(
    727     const psArray *images,              ///< Array of input images
    728     const psArray *masks,               ///< Array of input image masks
    729     const psArray *errors,              ///< The pixels which were rejected in the combination
    730     const psArray *inToOut,             ///< Transformation from input to output system
    731     const psArray *outToIn,             ///< Transformation from output to input system
    732     psF32 rejThreshold,                 ///< Rejection threshold
    733     psF32 gradLimit                     ///< Gradient limit
    734 )
    735 {
    736     psTrace("psModules.imcombine", 3, "Calling pmRejectPixels()\n");
    737     PS_ASSERT_PTR_NON_NULL(images, NULL);
    738     for (psS32 im = 0 ; im < images->n ; im++) {
    739         psImage *tmpImage = (psImage *) images->data[im];
    740         PS_ASSERT_IMAGE_NON_NULL(tmpImage, NULL);
    741         PS_ASSERT_IMAGE_NON_EMPTY(tmpImage, NULL);
    742         PS_ASSERT_IMAGE_TYPE(tmpImage, PS_TYPE_F32, NULL);
    743         if (masks != NULL) {
    744             PS_ASSERT_INT_EQUAL(images->n, masks->n, NULL);
    745             psImage *tmpMask = (psImage *) masks->data[im];
    746             PS_ASSERT_IMAGE_NON_NULL(tmpMask, NULL);
    747             PS_ASSERT_IMAGE_NON_EMPTY(tmpMask, NULL);
    748             PS_ASSERT_IMAGE_TYPE(tmpMask, PS_TYPE_F32, NULL);
    749             PS_ASSERT_IMAGES_SIZE_EQUAL(tmpImage, tmpMask, NULL);
    750         }
    751         PS_ASSERT_IMAGES_SIZE_EQUAL(((psImage *) images->data[0]), tmpImage, NULL);
    752     }
    753     PS_ASSERT_PTR_NON_NULL(errors, NULL);
    754     PS_ASSERT_PTR_NON_NULL(inToOut, NULL);
    755     PS_ASSERT_PTR_NON_NULL(outToIn, NULL);
    756     // Ensure that the psArray parameters have an element for each image.
    757     psS32 numImages = images->n;
    758     PS_ASSERT_INT_EQUAL(numImages, errors->n, NULL);
    759     PS_ASSERT_INT_EQUAL(numImages, inToOut->n, NULL);
    760     PS_ASSERT_INT_EQUAL(numImages, outToIn->n, NULL);
    761 
    762     //
    763     // Create the psArray of psPixelLists, one for each image, for rejected pixels.
    764     //
    765     psArray *rejects = psArrayAlloc(numImages);
    766     for (psS32 im = 0 ; im < numImages ; im++) {
    767         rejects->data[im] = (psPtr *) psPixelsAlloc(PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH);
    768         ((psPixels *)(rejects->data[im]))->n = ((psPixels *)(rejects->data[im]))->nalloc;
    769         psPixels *pixels = (psPixels *) rejects->data[im];
    770         pixels->n = 0;
    771     }
    772     //
    773     // rPtr is used to maintain a count of the questionable pixels for each image.
    774     //
    775     psVector *rPtr = psVectorAlloc(numImages, PS_TYPE_S32);
    776     psVectorInit(rPtr, 0);
    777 
    778     psS32 numCols = ((psImage *) images->data[0])->numCols;
    779     psS32 numRows = ((psImage *) images->data[0])->numRows;
    780     psRegion myRegion = psRegionSet(0, numCols-1, 0, numRows-1);
    781     psU32 maskVal = 1;  // XXX: Is this appropriate?
    782 
    783     psPlane *inCoords = psAlloc(sizeof(psPlane));
    784     psPlane *outCoords = psAlloc(sizeof(psPlane));
    785 
    786     for (psS32 im = 0 ; im < numImages ; im++) {
    787         //
    788         // Extract data from psArrays.
    789         //
    790         psPixels *pixelList = (psPixels *) errors->data[im];
    791 
    792         psImage *currImage = (psImage *) images->data[im];
    793         myRegion.x0 = 0;
    794         myRegion.x1 = currImage->numCols;
    795         myRegion.y0 = 0;
    796         myRegion.y1 = currImage->numRows;
    797         psPlaneTransform *myInToOut = (psPlaneTransform *) inToOut->data[im];
    798         psPlaneTransform *myOutToIn = (psPlaneTransform *) outToIn->data[im];
    799 
    800         //
    801         // Create a psU8 mask image from the list of cosmic pixels.
    802         //
    803         psImage *maskImage = NULL;
    804         maskImage = psPixelsToMask(maskImage, pixelList, myRegion, maskVal);
    805         psImage *maskImageF32 = ImageConvertF32(maskImage);
    806 
    807         //
    808         // Transform that mask image into detector coordinate space
    809         //
    810         psRegion myRegionXForm = DetermineRegion(maskImageF32, myOutToIn);
    811         psImage *transformedImage = psImageTransform(NULL, NULL, maskImageF32, NULL,
    812                                     0, myOutToIn, myRegionXForm, NULL,
    813                                     PS_INTERPOLATE_BILINEAR, 0);
    814 
    815         //
    816         // Loop over all cosmic pixels.  Transform their coords to detector space.
    817         // If the value of the transformed mask is larger than rejThreshold, then
    818         // this might be a cosmic ray pixel.  We then calculate the mean gradient
    819         // in other images.
    820         //
    821 
    822         psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR,
    823                                                                            transformedImage, NULL, NULL,
    824                                                                            0, 0.0, 0.0, 0, 0, 0.0);
    825 
    826         for (psS32 p = 0 ; p < pixelList->n ; p++) {
    827             inCoords->x = 0.5 + (psF32) (pixelList->data[p]).x;
    828             inCoords->y = 0.5 + (psF32) (pixelList->data[p]).y;
    829             psPlaneTransformApply(outCoords, myInToOut, inCoords);
    830             double maskVal;
    831             if (!psImageInterpolate(&maskVal, NULL, NULL, outCoords->x, outCoords->y, interp)) {
    832                 psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
    833                 psFree(interp);
    834                 psFree(maskImage);
    835                 psFree(maskImageF32);
    836                 psFree(transformedImage);
    837                 psFree(inCoords);
    838                 psFree(outCoords);
    839                 psFree(rejects);
    840                 return NULL;
    841             }
    842             if (maskVal > rejThreshold) {
    843 
    844                 // This is a possible cosmic array pixel.  We must calculate the gradient
    845                 // at this location in all input images.
    846                 psF32 meanGrads = 0.0;
    847                 psS32 numGrads = 0;
    848                 //
    849                 // Loop through all other images, calculate their mean gradient.
    850                 //
    851                 for (psS32 otherImg = 0 ; otherImg < numImages ; otherImg++) {
    852                     if (im != otherImg) {
    853                         // Map the outCoords to inCoords that for otherImg space.
    854                         psImage *tmpMask = NULL;
    855                         if (masks != NULL) {
    856                             tmpMask = masks->data[otherImg];
    857                         }
    858                         psPlaneTransformApply(inCoords,
    859                                               (psPlaneTransform * )outToIn->data[otherImg],
    860                                               outCoords);
    861                         psS32 xPix = (int)(inCoords->x + 0.5);
    862                         psS32 yPix = (int)(inCoords->y + 0.5);
    863                         if ((xPix >= 0) && (xPix <= ((psImage*)(images->data[otherImg]))->numCols - 1) &&
    864                                 (yPix >= 0) && (yPix <= ((psImage*)(images->data[otherImg]))->numRows - 1)) {
    865                             meanGrads += CalcGradient(images->data[otherImg], tmpMask, xPix, yPix);
    866                             numGrads++;
    867                         }
    868                     }
    869                 }
    870                 if (numGrads > 0) {
    871                     meanGrads /= (psF32) numGrads;
    872                 } else {
    873                     // XXX: my idea.  Verify with PWP:
    874                     meanGrads = 1.0 + gradLimit;
    875                 }
    876 
    877                 // XXX: The SDRS and the prototype code differ significantly here:
    878                 // if (CalcGradient(inputs->data[i], pixelList->data.x, pixelList->data.y) < (gradLimit * meanGrads)) {
    879                 if (meanGrads < gradLimit) {
    880                     //
    881                     // Add this to the list of questionable pixels.  We must ensure that the
    882                     // pixelList is large enough; if not, we realloc()
    883                     //
    884                     psS32 ptr = rPtr->data.S32[im];
    885                     psPixels *pixelListPtr = (psPixels *) rejects->data[im];
    886                     if (ptr >= pixelListPtr->nalloc) {
    887                         rejects->data[im] = (psPtr *) psPixelsRealloc(((psPixels *) rejects->data[im]),
    888                                             ((((psPixels *) rejects->data[im])->nalloc) + PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC));
    889                         // XXX: Can the realloc() fail?  Must we check for NULL?
    890                     }
    891 
    892                     ((psPixels *) rejects->data[im])->data[ptr].x = (pixelList->data[p]).x;
    893                     ((psPixels *) rejects->data[im])->data[ptr].y = (pixelList->data[p]).y;
    894                     (rPtr->data.S32[im])++;
    895                     // XXX: this pixel ->n increment is wierd
    896                     (((psPixels *) rejects->data[im])->n)++;
    897                 }
    898             }
    899         }
    900 
    901         psFree(interp);
    902         psFree(maskImage);
    903         psFree(maskImageF32);
    904         psFree(transformedImage);
    905     }
    906 
    907     psFree(inCoords);
    908     psFree(outCoords);
    909     psTrace("psModules.imcombine", 3, "Exiting pmRejectPixels()\n");
    910     return(rejects);
    911 }
    912 #endif
Note: See TracChangeset for help on using the changeset viewer.