/** @file  pmImageCombine.c
 *
 *  This file will perform image combination of several images of the
 *  same field, produce a list of questionable pixels, then tag some
 *  of those pixels as cosmic rays.
 *
 *  @author Paul Price, IfA (original prototype)
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-12-08 00:29:53 $
 *
 *  XXX: pmRejectPixels() has a known bug with the pmImageTransform() call.
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#include <config.h>
#include <stdio.h>
#include <math.h>
#include "pslib.h"

//
// The following macros define how big the initial pixel list will be, and
// how much it should be incremented when realloc'ed.
//
#define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH 100
#define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC 100
#define PS_COMBINE_IMAGE_MAX_QUESTIONABLE_PIXELS 1000
/******************************************************************************
pmCombineImages(combine, questionablePixels, images, errors, masks, maskVal,
                pixels, numIter, sigmaClip, stats)
 
XXX: Allocate a dummy psStats structure so that we don't destroy away its data.
 *****************************************************************************/
psImage *pmCombineImages(
    psImage *combine,                   ///< Combined image (output)
    psArray **questionablePixels,       ///< Array of rejection masks
    const psArray *images,              ///< Array of input images
    const psArray *errors,              ///< Array of input error images
    const psArray *masks,               ///< Array of input masks
    psU32 maskVal,                      ///< Mask value
    const psPixels *pixels,             ///< Pixels to combine
    psS32 numIter,                      ///< Number of rejection iterations
    psF32 sigmaClip,                    ///< Number of standard deviations at which to reject
    const psStats *stats)               ///< Statistics to use in the combination
{

    PS_ASSERT_PTR_NON_NULL(images, combine);
    psU32 numImages = images->n;
    psTrace("ImageCombine.pmCombineImages", 3, "Calling pmCombineImages(%d)\n", images->n);

    if (errors != NULL) {
        if (images->n != errors->n) {
            psError(PS_ERR_UNKNOWN, true, "images and errors args must have same length (%d != %d)\n",
                    images->n, errors->n);
            return(combine);
        }
    }
    if (masks != NULL) {
        if (images->n != masks->n) {
            psError(PS_ERR_UNKNOWN, true, "images and masks args must have same length (%d != %d)\n",
                    images->n, masks->n);
            return(combine);
        }
    }

    psImage *tmpImg = (psImage *) images->data[0];
    psU32 numRows = tmpImg->numRows;
    psU32 numCols = tmpImg->numCols;

    //
    // Check that all images have the appropriate size and type.
    //
    for (psS32 i = 0 ; i < numImages ; i++) {
        psImage *tmpDataImg = (psImage *) images->data[i];
        PS_ASSERT_IMAGE_NON_NULL(tmpDataImg, combine);
        PS_ASSERT_IMAGE_TYPE(tmpDataImg, PS_TYPE_F32, combine);
        if ((tmpDataImg->numRows != numRows) || (tmpDataImg->numCols != numCols)) {
            psError(PS_ERR_UNKNOWN, true, "image %d has size (%d, %d); should be (%d, %d).\n",
                    i, tmpDataImg->numRows, tmpDataImg->numCols, numRows, numCols);
        }

        if (errors != NULL) {
            psImage *tmpErrorImg = (psImage *) errors->data[i];
            PS_ASSERT_IMAGE_NON_NULL(tmpErrorImg, combine);
            PS_ASSERT_IMAGE_TYPE(tmpErrorImg, PS_TYPE_F32, combine);
            PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpErrorImg, NULL);
        }

        if (masks != NULL) {
            psImage *tmpMaskImg = (psImage *) masks->data[i];
            PS_ASSERT_IMAGE_NON_NULL(tmpMaskImg, combine);
            PS_ASSERT_IMAGE_TYPE(tmpMaskImg, PS_TYPE_U8, combine);
            PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpMaskImg, NULL);
        }
    }
    PS_ASSERT_PTR_NON_NULL(stats, combine);

    // Allocate and initialize the combined image, if necessary.
    if (combine == NULL) {
        combine = psImageAlloc(numCols, numRows, PS_TYPE_F32);
        if (pixels != NULL) {
            PS_IMAGE_SET_F32(combine, 0.0);
        }
    }

    //
    // Allocate the questionablePixels psArray, if necesssary, then create a psPixels
    // struct for each image.
    //
    if (*questionablePixels == NULL) {
        *questionablePixels = psArrayAlloc(numImages);
    } else if ((*questionablePixels)->n != numImages) {
        *questionablePixels = psArrayRealloc(*questionablePixels, numImages);
    }
    for (psS32 im = 0 ; im < numImages ; im++) {
        psFree((*questionablePixels)->data[im]);
        ((*questionablePixels)->data[im]) = (psPtr *) psPixelsAlloc(PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH);
        ((psPixels *) ((*questionablePixels)->data[im]))->n = 0;
    }
    //
    // qpPtr is used to maintain a count of the questionable pixels for each image.
    //
    psVector *qpPtr = psVectorAlloc(numImages, PS_TYPE_S32);
    PS_VECTOR_SET_S32(qpPtr, 0);

    //
    // Allocate the necessary psVectors for the call to psVectorStats().
    // These vectors will be used whether we are combining a list of pixels,
    // or every pixel in the input images.
    //
    psVector *pixelData = psVectorAlloc(numImages, PS_TYPE_F32);

    psVector *pixelMask = NULL;
    if (masks != NULL) {
        pixelMask = psVectorAlloc(numImages, PS_TYPE_U8);
        PS_VECTOR_SET_U8(pixelMask, 0);
    }

    psVector *pixelErrors = NULL;
    if (errors != NULL) {
        pixelErrors = psVectorAlloc(numImages, PS_TYPE_F32);
        PS_VECTOR_SET_F32(pixelErrors, 1.0);
    }

    if (pixels != NULL) {
        // Only those specified pixels should be combined.

        psStats *stdevStats = psStatsAlloc(PS_STAT_SAMPLE_STDEV);

        for (psS32 p = 0 ; p < pixels->n ; p++) {
            // Must initialize the mask to 0 for every pixel.
            PS_VECTOR_SET_U8(pixelMask, 0);
            psS32 col = (pixels->data[p]).x;
            psS32 row = (pixels->data[p]).y;

            //
            // Loop through each image, extract the pixel/mask/error data
            // into psVectors.
            //
            for (psS32 im = 0 ; im < numImages ; im++) {
                // Set the pixel data
                pixelData->data.F32[im] =       ((psImage *) images->data[im])->data.F32[row][col];
                // Set the pixel mask data, if necessary
                if (masks != NULL) {
                    pixelMask->data.U8[im] =   ((psImage *) masks->data[im])->data.F32[row][col];
                }

                // Set the pixel error data, if necessary
                if (errors != NULL) {
                    pixelErrors->data.F32[im] = ((psImage *) errors->data[im])->data.F32[row][col];
                }
            }

            //
            // Iterate on the pixels, rejecting outliers
            //
            for (psS32 iter = 0 ; iter < numIter ; iter++) {
                // Combine all the pixels, using the specified stat.

                stats = psVectorStats((psStats *) stats, pixelData, pixelErrors, pixelMask, maskVal);
                psF64 combinedPixel;
                psBool rc = p_psGetStatValue(stats, &combinedPixel);
                if (rc != true) {
                    psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not combine pixels (%d, %d) with the specified stat.\n", row, col);
                }
                if (iter == 0) {
                    combine->data.F32[row][col] = (psF32) combinedPixel;
                }

                //
                // Reject all pixels that lie more that sigmaClip standard deviations from
                // the combined pixel value.
                //
                psS32 numRejects = 0;
                for (psS32 im = 0 ; im < numImages ; im++) {
                    stdevStats = psVectorStats(stdevStats, pixelData, pixelErrors, pixelMask, maskVal);
                    psF64 stdev;
                    psBool rc = p_psGetStatValue(stdevStats, &stdev);
                    if (rc != true) {
                        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not compute the standard deviation of pixel (%d, %d).\n", row, col);
                        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not combine pixels (%d, %d) with the specified stat.\n", row, col);
                    }

                    if (!(pixelMask->data.U8[im] & maskVal)) {
                        if (fabs(pixelData->data.F32[im]) >
                                (fabs(sigmaClip * stdev) + fabs(combinedPixel))) {
                            numRejects++;
                            pixelMask->data.U8[im] = maskVal;
                            //
                            // XXX: These data structures indirections are getting complicated.
                            //
                            psS32 ptr = qpPtr->data.S32[im];
                            psPixels *pixelListPtr = ((psPixels *) ((*questionablePixels)->data[im]));

                            if (ptr >= pixelListPtr->nalloc) {
                                (*questionablePixels)->data[im] =
                                    (psPtr *) psPixelsRealloc(((psPixels *) ((*questionablePixels)->data[im])),
                                                              ((((psPixels *) ((*questionablePixels)->data[im]))->nalloc) +
                                                               PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC));
                                // XXX: Can the realloc() fail?  Must we check for NULL?
                            }
                            ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].x = col;
                            ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].y = row;
                            (qpPtr->data.S32[im])++;
                            // XXX: this pixel ->n increment is wierd
                            ((psPixels *) ((*questionablePixels)->data[im]))->n = qpPtr->data.S32[im];
                        }
                    }
                }

                //
                // If the number of rejected pixels is zero, then there's no point
                // continuing the loop.
                //
                if (numRejects == 0) {
                    break;
                }
                psS32 totalRejects = 0;
                for (psS32 im = 0 ; im < numImages ; im++) {
                    if (pixelMask->data.U8[im] & maskVal) {
                        totalRejects++;
                    }
                }

                //
                // XXX: Is it possible to have all pixels rejected?  If so, we should
                // exit the loop.
                //
                if (totalRejects == numImages) {
                    break;
                }
            }
        }
        psFree(stdevStats);
    } else {
        //
        // We get here if there is a NULL list of pixels to combine.
        // Therefore, we combine all pixels in all images.
        //

        //
        // Loop over all pixels in all images, set the appropriate data, mask,
        // error vectors, call psVectorStats(), and set the result in the
        // combine image.
        //
        for (psS32 row = 0 ; row < numRows ; row++) {
            for (psS32 col = 0 ; col < numCols ; col++) {
                for (psS32 im = 0 ; im < numImages ; im++) {
                    // Set the pixel data
                    pixelData->data.F32[im] =       ((psImage *) images->data[im])->data.F32[row][col];

                    // Set the pixel mask data, if necessary
                    if (masks != NULL) {
                        pixelMask->data.U8[im] =   ((psImage *) masks->data[im])->data.F32[row][col];
                    }

                    // Set the pixel error data, if necessary
                    if (errors != NULL) {
                        pixelErrors->data.F32[im] = ((psImage *) errors->data[im])->data.F32[row][col];
                    }

                }
                // Combine all the pixels, using the specified stat.
                stats = psVectorStats((psStats *) stats, pixelData, pixelErrors, pixelMask, maskVal);
                psF64 tmpF64;
                psBool rc = p_psGetStatValue(stats, &tmpF64);
                combine->data.F32[row][col] = (psF32) tmpF64;
                if (rc != true) {
                    psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not combine pixels (%d, %d) with the specified stat.\n", row, col);
                }
            }
        }
    }

    psFree(pixelData);
    psFree(pixelMask);
    psFree(pixelErrors);
    psFree(qpPtr);

    psTrace("ImageCombine.pmCombineImages", 3, "Exiting pmCombineImages(%d)\n", images->n);
    return(combine);
}


/******************************************************************************
XXX: Directly from Paul Price
 *****************************************************************************/
static psF32 CalcGradient(
    psImage *image,
    psImage *imageMask,
    psS32 x,
    psS32 y
)
{
    psTrace("ImageCombine.CalcGradient", 4, "Calling CalcGradient(%d, %d)\n", x, y);
    int num = 0;
    psVector *pixels = psVectorAlloc(8, PS_TYPE_F32); // Array of pixels
    psVector *mask = psVectorAlloc(8, PS_TYPE_U8); // Corresponding mask

    // Get limits
    int xMin = PS_MAX(x - 1, 0);
    int xMax = PS_MIN(x + 1, image->numCols - 1);
    int yMin = PS_MAX(y - 1, 0);
    int yMax = PS_MIN(y + 1, image->numRows - 1);
    if (imageMask != NULL) {
        for (int j = yMin; j <= yMax; j++) {
            for (int i = xMin; i <= xMax; i++) {
                if ((i != x) && (j != y) && (0 == imageMask->data.U8[j][i])) {
                    pixels->data.F32[num] = image->data.F32[j][i];
                    mask->data.U8[num] = 0;
                    num++;
                } else {
                    mask->data.U8[num] = 1;
                }
            }
        }
    } else {
        //
        // This code is simply the previous loop without the imageMask.
        // XXX: Consider restructuring this.
        //
        for (int j = yMin; j <= yMax; j++) {
            for (int i = xMin; i <= xMax; i++) {
                if ((i != x) && (j != y)) {
                    pixels->data.F32[num] = image->data.F32[j][i];
                    mask->data.U8[num] = 0;
                    num++;
                } else {
                    mask->data.U8[num] = 1;
                }
            }
        }
    }

    pixels->n = num;
    mask->n = num;

    // Get the median
    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    (void)psVectorStats(stats, pixels, NULL, mask, 1);
    float median = stats->sampleMedian;
    psFree(stats);
    psFree(pixels);
    psFree(mask);

    psTrace("ImageCombine.CalcGradient", 4, "Exiting CalcGradient(%d, %d)\n", x, y);
    return(median / image->data.F32[y][x]);
}

/******************************************************************************
DetermineRegion(image, myOutToIn): for a psImage and a psPlaneTransform to that
image, this routine determines the size of the input image which maps to that
image, and returns the result in a psRegion struct.
 
XXX: Basically, this routine is only guaranteed to work if the transform is
linear.
 
XXX: Shouldn't this functionality be part of psImageTransform()?
 *****************************************************************************/
static psRegion DetermineRegion(psImage *image,
                                psPlaneTransform *myOutToIn)
{
    psTrace("ImageCombine.DetermineRegion", 4, "Calling DetermineRegion()\n");
    psRegion myRegion;
    myRegion.x0 = PS_MAX_F32;
    myRegion.x1 = PS_MIN_F32;
    myRegion.y0 = PS_MAX_F32;
    myRegion.y1 = PS_MIN_F32;
    psPlane in;
    psPlane out;

    in.x = 0.0;
    in.y = 0.0;

    psPlaneTransformApply(&out, myOutToIn, &in);
    if (out.x < myRegion.x0) {
        myRegion.x0 = out.x;
    }
    if (out.x > myRegion.x1) {
        myRegion.x1 = out.x;
    }
    if (out.y < myRegion.y0) {
        myRegion.y0 = out.y;
    }
    if (out.y > myRegion.y1) {
        myRegion.y1 = out.y;
    }

    in.x = (psF32) (image->numCols);
    in.y = 0.0;
    psPlaneTransformApply(&out, myOutToIn, &in);
    if (out.x < myRegion.x0) {
        myRegion.x0 = out.x;
    }
    if (out.x > myRegion.x1) {
        myRegion.x1 = out.x;
    }
    if (out.y < myRegion.y0) {
        myRegion.y0 = out.y;
    }
    if (out.y > myRegion.y1) {
        myRegion.y1 = out.y;
    }

    in.x = (psF32) (image->numCols);
    ;
    in.y = 0.0;
    psPlaneTransformApply(&out, myOutToIn, &in);
    if (out.x < myRegion.x0) {
        myRegion.x0 = out.x;
    }
    if (out.x > myRegion.x1) {
        myRegion.x1 = out.x;
    }
    if (out.y < myRegion.y0) {
        myRegion.y0 = out.y;
    }
    if (out.y > myRegion.y1) {
        myRegion.y1 = out.y;
    }

    in.x = (psF32) (image->numCols);
    in.y = (psF32) (image->numRows);
    psPlaneTransformApply(&out, myOutToIn, &in);
    if (out.x < myRegion.x0) {
        myRegion.x0 = out.x;
    }
    if (out.x > myRegion.x1) {
        myRegion.x1 = out.x;
    }
    if (out.y < myRegion.y0) {
        myRegion.y0 = out.y;
    }
    if (out.y > myRegion.y1) {
        myRegion.y1 = out.y;
    }

    psTrace("ImageCombine.DetermineRegion", 4, "Exiting DetermineRegion()\n");
    return(myRegion);
}

/******************************************************************************
XXX: Don't we have a psLib function for this?
 *****************************************************************************/
static psImage *ImageConvertF32(psImage *image)
{
    psTrace("ImageCombine.ImageConvertF32", 4, "Calling ImageConvertF32()\n");
    psImage *imgF32 = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);

    for (psS32 i = 0 ; i < image->numRows ; i++) {
        for (psS32 j = 0 ; j < image->numCols ; j++) {
            imgF32->data.F32[i][j] = (psF32) image->data.U8[i][j];
        }
    }

    psTrace("ImageCombine.ImageConvertF32", 4, "Exiting ImageConvertF32()\n");
    return(imgF32);
}


//
// The following macros define how big the initial pixel list will be, and
// how much it should be incremented when realloc'ed.
//
#define PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH 100
#define PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC 100
/******************************************************************************
pmRejectPixels(images, errors, inToOut, outToIn, rejThreshold,
gradLimit)
 
XXX: Optimization: we don't need to transform the entire mask image.
XXX: The inToOut and outToIn transforms are confusing.  Verify that what
     I think they mean syncs with PWP.
 *****************************************************************************/
psArray *pmRejectPixels(
    const psArray *images,              ///< Array of input images
    const psArray *masks,               ///< Array of input image masks
    const psArray *errors,              ///< The pixels which were rejected in the combination
    const psArray *inToOut,             ///< Transformation from input to output system
    const psArray *outToIn,             ///< Transformation from output to input system
    psF32 rejThreshold,                 ///< Rejection threshold
    psF32 gradLimit                     ///< Gradient limit
)
{
    psTrace("ImageCombine.pmRejectPixels", 3, "Calling pmRejectPixels()\n");
    PS_ASSERT_PTR_NON_NULL(images, NULL);
    for (psS32 im = 0 ; im < images->n ; im++) {
        psImage *tmpImage = (psImage *) images->data[im];
        PS_ASSERT_IMAGE_NON_NULL(tmpImage, NULL);
        PS_ASSERT_IMAGE_NON_EMPTY(tmpImage, NULL);
        PS_ASSERT_IMAGE_TYPE(tmpImage, PS_TYPE_F32, NULL);
        if (masks != NULL) {
            PS_ASSERT_INT_EQUAL(images->n, masks->n, NULL);
            psImage *tmpMask = (psImage *) masks->data[im];
            PS_ASSERT_IMAGE_NON_NULL(tmpMask, NULL);
            PS_ASSERT_IMAGE_NON_EMPTY(tmpMask, NULL);
            PS_ASSERT_IMAGE_TYPE(tmpMask, PS_TYPE_F32, NULL);
            PS_ASSERT_IMAGES_SIZE_EQUAL(tmpImage, tmpMask, NULL);
        }
        PS_ASSERT_IMAGES_SIZE_EQUAL(((psImage *) images->data[0]), tmpImage, NULL);
    }
    PS_ASSERT_PTR_NON_NULL(errors, NULL);
    PS_ASSERT_PTR_NON_NULL(inToOut, NULL);
    PS_ASSERT_PTR_NON_NULL(outToIn, NULL);
    // Ensure that the psArray parameters have an element for each image.
    psS32 numImages = images->n;
    PS_ASSERT_INT_EQUAL(numImages, errors->n, NULL);
    PS_ASSERT_INT_EQUAL(numImages, inToOut->n, NULL);
    PS_ASSERT_INT_EQUAL(numImages, outToIn->n, NULL);

    //
    // Create the psArray of psPixelLists, one for each image, for rejected pixels.
    //
    psArray *rejects = psArrayAlloc(numImages);
    for (psS32 im = 0 ; im < numImages ; im++) {
        rejects->data[im] = (psPtr *) psPixelsAlloc(PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH);
        psPixels *pixels = (psPixels *) rejects->data[im];
        pixels->n = 0;
    }
    //
    // rPtr is used to maintain a count of the questionable pixels for each image.
    //
    psVector *rPtr = psVectorAlloc(numImages, PS_TYPE_S32);
    PS_VECTOR_SET_S32(rPtr, 0);

    psS32 numCols = ((psImage *) images->data[0])->numCols;
    psS32 numRows = ((psImage *) images->data[0])->numRows;
    psRegion myRegion = psRegionSet(0, numCols-1, 0, numRows-1);
    psU32 maskVal = 1;  // XXX: Is this appropriate?

    psPlane *inCoords = psAlloc(sizeof(psPlane));
    psPlane *outCoords = psAlloc(sizeof(psPlane));

    for (psS32 im = 0 ; im < numImages ; im++) {
        //
        // Extract data from psArrays.
        //
        psPixels *pixelList = (psPixels *) errors->data[im];

        psImage *currImage = (psImage *) images->data[im];
        myRegion.x0 = 0;
        myRegion.x1 = currImage->numCols;
        myRegion.y0 = 0;
        myRegion.y1 = currImage->numRows;
        psPlaneTransform *myInToOut = (psPlaneTransform *) inToOut->data[im];
        psPlaneTransform *myOutToIn = (psPlaneTransform *) outToIn->data[im];

        //
        // Create a psU8 mask image from the list of cosmic pixels.
        //
        psImage *maskImage = NULL;
        maskImage = psPixelsToMask(maskImage, pixelList, myRegion, maskVal);
        psImage *maskImageF32 = ImageConvertF32(maskImage);

        //
        // Transform that mask image into detector coordinate space
        //
        psRegion myRegionXForm = DetermineRegion(maskImageF32, myOutToIn);
        psImage *transformedImage = psImageTransform(NULL, NULL, maskImageF32, NULL,
                                    0, myOutToIn, myRegionXForm, NULL,
                                    PS_INTERPOLATE_BILINEAR, 0);

        //
        // Loop over all cosmic pixels.  Transform their coords to detector space.
        // If the value of the transformed mask is larger than rejThreshold, then
        // this might be a cosmic ray pixel.  We then calculate the mean gradient
        // in other images.
        //
        for (psS32 p = 0 ; p < pixelList->n ; p++) {
            inCoords->x = 0.5 + (psF32) (pixelList->data[p]).x;
            inCoords->y = 0.5 + (psF32) (pixelList->data[p]).y;
            psPlaneTransformApply(outCoords, myInToOut, inCoords);
            psF32 maskVal = (psF32) psImagePixelInterpolate(transformedImage, outCoords->x, outCoords->y,
                            NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
            if (maskVal > rejThreshold) {

                // This is a possible cosmic array pixel.  We must calculate the gradient
                // at this location in all input images.
                psF32 meanGrads = 0.0;
                psS32 numGrads = 0;
                //
                // Loop through all other images, calculate their mean gradient.
                //
                for (psS32 otherImg = 0 ; otherImg < numImages ; otherImg++) {
                    if (im != otherImg) {
                        // Map the outCoords to inCoords that for otherImg space.
                        psImage *tmpMask = NULL;
                        if (masks != NULL) {
                            tmpMask = masks->data[otherImg];
                        }
                        psPlaneTransformApply(inCoords,
                                              (psPlaneTransform * )outToIn->data[otherImg],
                                              outCoords);
                        psS32 xPix = (int)(inCoords->x + 0.5);
                        psS32 yPix = (int)(inCoords->y + 0.5);
                        if ((xPix >= 0) && (xPix <= ((psImage*)(images->data[otherImg]))->numCols - 1) &&
                                (yPix >= 0) && (yPix <= ((psImage*)(images->data[otherImg]))->numRows - 1)) {
                            meanGrads += CalcGradient(images->data[otherImg], tmpMask, xPix, yPix);
                            numGrads++;
                        }
                    }
                }
                if (numGrads > 0) {
                    meanGrads /= (psF32) numGrads;
                } else {
                    // XXX: my idea.  Verify with PWP:
                    meanGrads = 1.0 + gradLimit;
                }

                // XXX: The SDRS and the prototype code differ significantly here:
                // if (CalcGradient(inputs->data[i], pixelList->data.x, pixelList->data.y) < (gradLimit * meanGrads)) {
                if (meanGrads < gradLimit) {
                    //
                    // Add this to the list of questionable pixels.  We must ensure that the
                    // pixelList is large enough; if not, we realloc()
                    //
                    psS32 ptr = rPtr->data.S32[im];
                    psPixels *pixelListPtr = (psPixels *) rejects->data[im];
                    if (ptr >= pixelListPtr->nalloc) {
                        rejects->data[im] = (psPtr *) psPixelsRealloc(((psPixels *) rejects->data[im]),
                                            ((((psPixels *) rejects->data[im])->nalloc) + PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC));
                        // XXX: Can the realloc() fail?  Must we check for NULL?
                    }

                    ((psPixels *) rejects->data[im])->data[ptr].x = (pixelList->data[p]).x;
                    ((psPixels *) rejects->data[im])->data[ptr].y = (pixelList->data[p]).y;
                    (rPtr->data.S32[im])++;
                    // XXX: this pixel ->n increment is wierd
                    (((psPixels *) rejects->data[im])->n)++;
                }
            }
        }

        psFree(maskImage);
        psFree(maskImageF32);
        psFree(transformedImage);
    }

    psFree(inCoords);
    psFree(outCoords);
    psTrace("ImageCombine.pmRejectPixels", 3, "Exiting pmRejectPixels()\n");
    return(rejects);
}
