/** @file  psImageExtraction.c
 *
 *  @brief Contains basic image extraction operations, as specified in the 
 *         PSLIB SDRS sections "Image Pixel Extractions" and "Image Structure
 *         Manipulation".
 *
 *  @ingroup Image
 *
 *  @author Robert DeSonia, MHPCC
 *
 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-03-22 21:52:49 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
 *
 */

#include <string.h>

#include "psMemory.h"
#include "psImageExtraction.h"
#include "psError.h"

#include "psImageErrors.h"

psImage* imageSubset(psImage* out,
                     psImage* image,
                     psS32 col0,
                     psS32 row0,
                     psS32 col1,
                     psS32 row1)
{
    psU32 elementSize;          // size of image element in bytes
    psU32 inputColOffset;       // offset in bytes to first subset pixel in input row

    if (image == NULL || image->data.V == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        return NULL;
    }

    if (image->type.dimen != PS_DIMEN_IMAGE) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psImage_NOT_AN_IMAGE);
        return NULL;
    }

    if (col1 < 1) {
        col1 = image->numCols + col1;
    }
    if (row1 < 1) {
        row1 = image->numRows + row1;
    }

    if (    col1 <= col0 ||
            row1 <= row0 ||
            col0 >= image->numCols ||
            row0 >= image->numRows ||
            col1 > image->numCols ||
            row1 > image->numRows ) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_SUBSET_RANGE_INVALID,
                col0, col1-1, row0, row1-1,
                image->numCols-1, image->numRows-1);
        return NULL;
    }
    psS32 numRows = row1-row0;
    psS32 numCols = col1-col0;

    elementSize = PSELEMTYPE_SIZEOF(image->type.type);

    if (image->parent != NULL) { // if this is a child, we need to start working with parent.
        col0 += image->col0;
        col1 += image->col0;
        row0 += image->row0;
        row1 += image->row0;
        image = (psImage*)image->parent;
    }

    // increment the raw data buffer before freeing anything in the 'out'
    psPtr rawData = psMemIncrRefCounter(image->rawDataBuffer);

    if (out != NULL) {
        // if a child, need to orphan (disassociate from parent) first
        if (out->parent != NULL) {
            psArrayRemove(out->parent->children,out); // remove from parent's knowledge
            out->parent = NULL; // break link to parent
        }

        psFree(out->rawDataBuffer); // free the previous data reference
    } else {
        out = psAlloc(sizeof(psImage));
        out->data.V = NULL;
    }

    out->data.V = psRealloc(out->data.V,sizeof(psPtr)*numRows); // resize row pointer array
    *(psType*)&out->type = image->type;
    *(psU32*)&out->numCols = numCols;
    *(psU32*)&out->numRows = numRows;
    *(psS32*)&out->row0 = row0;
    *(psS32*)&out->col0 = col0;
    out->parent = image;
    out->children = NULL;
    out->rawDataBuffer = rawData;

    // set the new psImage's deallocator to the same as the input image
    p_psMemSetDeallocator(out,p_psMemGetDeallocator(image));

    inputColOffset = elementSize * col0;
    for (psS32 row = 0; row < numRows; row++) {
        out->data.V[row] = image->data.U8[row0 + row] + inputColOffset;
    }

    // add output image as a child of the input image.
    psS32 n = 0;
    psArray* children = image->children;
    if (children == NULL) {
        children = psArrayAlloc(16); // start with a reasonable size for growth
    } else if (children->nalloc == children->n) { // full?
        n = children->n;
        children = psArrayRealloc(children,n*2); // double the array size
    } else {
        n = children->n;
    }
    children->data[n] = out;
    children->n = n+1;
    image->children = children; // push back any change (esp. if children==NULL before)

    return (out);
}

psImage* psImageSubset(psImage* image,
                       psS32 col0,
                       psS32 row0,
                       psS32 col1,
                       psS32 row1)
{
    return imageSubset(NULL,image,col0,row0,col1,row1);
}

psImage* psImageSubsection(psImage* image,
                           const char* section)
{
    psS32 col0;
    psS32 col1;
    psS32 row0;
    psS32 row1;

    // section should be of the form '[col0:col1,row0:row1]'
    if (section == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_SUBSECTION_NULL);
        return NULL;
    }

    if (sscanf(section,"[%d:%d,%d:%d]",&col0,&col1,&row0,&row1) < 4) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_SUBSECTION_INVALID,
                section);
        return NULL;
    }

    if (col0 > col1 || row0 > row1) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_SUBSET_RANGE_MALFORMED,
                col0,col1,row0,row1);
        return NULL;
    }

    return imageSubset(NULL,image,col0,row0,col1+1,row1+1);
}

psImage* psImageTrim(psImage* image, psS32 col0, psS32 row0, psS32 col1, psS32 row1)
{
    if (image == NULL || image->data.V == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        return NULL;
    }

    if (image->parent != NULL) {
        return imageSubset(image,
                           (psImage*)image->parent,
                           col0+image->col0,
                           row0+image->row0,
                           col1+image->col0,
                           row1+image->row0);
    }

    if (col1 < 1) {
        col1 += image->numCols;
    }

    if (row1 < 1) {
        row1 += image->numRows;
    }

    if (    col0 < 0 ||
            row0 < 0 ||
            col1 > image->numCols ||
            row1 > image->numRows ||
            col0 >= col1 ||
            row0 >= row1 ) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_SUBSET_RANGE_INVALID,
                col0, col1-1, row0, row1-1,
                image->numCols-1, image->numRows-1);
        psFree(image);
        return NULL;
    }

    psImageFreeChildren(image);

    psU32 elementSize = PSELEMTYPE_SIZEOF(image->type.type);
    psU32 numCols = col1-col0;
    psU32 numRows = row1-row0;
    psU32 rowSize = elementSize*numCols;
    psU32 colOffset = elementSize * col0;
    psU8* imageData = image->rawDataBuffer;
    for (psS32 row = row0; row < row1; row++) {
        memmove(imageData,image->data.U8[row] + colOffset,rowSize);
        imageData += rowSize;
    }

    *(psU32*)&image->numRows = numRows;
    *(psU32*)&image->numCols = numCols;

    // XXX: should I really resize the buffers?
    image->data.V = psRealloc(image->data.V,sizeof(psPtr)*numRows);
    image->rawDataBuffer = psRealloc(image->rawDataBuffer,rowSize*numRows);

    image->data.V[0] = image->rawDataBuffer;
    for (psS32 r = 1; r < numRows; r++) {
        image->data.U8[r] = image->data.U8[r-1] + rowSize;
    }

    return (image);
}

psVector* psImageSlice(psVector* out,
                       psVector* slicePositions,
                       const psImage* restrict in,
                       const psImage* restrict mask,
                       psU32 maskVal,
                       psS32 col0,
                       psS32 row0,
                       psS32 col1,
                       psS32 row1,
                       psImageCutDirection direction,
                       const psStats* stats)
{
    double statVal;
    psStats* myStats;
    psElemType type;
    psS32 inRows;
    psS32 inCols;
    psS32 delta = 1;
    psF64* outData;

    if (in == NULL || in->data.V == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        psFree(out);
        return NULL;
    }

    if (col1 < 1) {
        col1 += in->numCols;
    }

    if (row1 < 1) {
        row1 += in->numRows;
    }

    if (    col0 < 0 ||
            row0 < 0 ||
            col1 > in->numCols ||
            row1 > in->numRows ||
            col0 >= col1 ||
            row0 >= row1) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_SUBSET_RANGE_INVALID,
                col0, col1, row0, row1,
                in->numCols, in->numRows);
        psFree(out);
        return NULL;
    }

    type = in->type.type;
    inRows = in->numRows;
    inCols = in->numCols;

    if (direction == PS_CUT_X_NEG || direction == PS_CUT_Y_NEG) {
        delta = -1;
    }

    if (mask != NULL) {
        if (inRows != mask->numRows || inCols != mask->numCols) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    PS_ERRORTEXT_psImage_IMAGE_MASK_SIZE,
                    mask->numCols,mask->numRows,
                    inCols, inRows);
            psFree(out);
            return NULL;
        }
        if (mask->type.type != PS_TYPE_MASK) {
            char* typeStr;
            PS_TYPE_NAME(typeStr,mask->type.type);
            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                    PS_ERRORTEXT_psImage_IMAGE_MASK_TYPE,
                    typeStr, PS_TYPE_MASK_NAME);
            psFree(out);
            return NULL;
        }
    }

    if (stats == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_STAT_NULL);
        psFree(out);
        return NULL;
    }

    // verify that the stats struct specifies a
    // single stats operation
    if (p_psGetStatValue(stats, &statVal) == false) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                PS_ERRORTEXT_psImage_BAD_STAT,stats->options);
        psFree(out);
        return NULL;
    }
    // since stats input is const, I need to
    // create a 'scratch' stats struct
    myStats = psAlloc(sizeof(psStats));
    *myStats = *stats;

    psS32 numCols = col1-col0;
    psS32 numRows = row1-row0;

    if (direction == PS_CUT_X_POS || direction == PS_CUT_X_NEG) {
        psVector* imgVec = psVectorAlloc(numRows, type);
        psVector* maskVec = NULL;
        psMaskType* maskData = NULL;
        psU32* outPosition = NULL;

        // recycle output to make a proper sized/type output structure
        // n.b. type is double as that is the type given for all stats is
        // psStats.
        out = psVectorRecycle(out, numCols, PS_TYPE_F64);
        if (slicePositions != NULL) {
            slicePositions = psVectorRecycle(slicePositions, numCols, PS_TYPE_U32);
            outPosition = slicePositions->data.U32;
        }
        outData = out->data.F64;
        if (delta < 0) {
            outData += numCols - 1;
            if (outPosition != NULL) {
                outPosition += numCols - 1;
            }
        }

        if (mask != NULL) {
            maskVec = psVectorAlloc(numRows, mask->type.type);
        }
        #define PSIMAGE_CUT_VERTICAL(TYPE) \
    case PS_TYPE_##TYPE: { \
            psMaskType* maskVecData = NULL; \
            for (psS32 c=col0;c<col1;c++) { \
                ps##TYPE *imgData = in->data.TYPE[row0] + c; \
                ps##TYPE *imgVecData = imgVec->data.TYPE; \
                if (maskVec != NULL) { \
                    maskVecData = maskVec->data.U8; \
                    maskData = (psMaskType* )(mask->data.U8[row0]) + c; \
                } \
                for (psS32 r=row0;r<row1;r++) { \
                    *(imgVecData++) = *imgData; \
                    imgData += inCols; \
                    if (maskVecData != NULL) { \
                        *(maskVecData++) = *maskData; \
                        maskData += inCols; \
                    } \
                } \
                myStats = psVectorStats(myStats,imgVec,NULL,maskVec,maskVal); \
                (void)p_psGetStatValue(myStats,&statVal); \
                *outData = statVal; \
                if (outPosition != NULL) { \
                    *outPosition = c; \
                    outPosition += delta; \
                } \
                outData += delta; \
            } \
            break; \
        }

        switch (type) {
            PSIMAGE_CUT_VERTICAL(U8);  // Not a requirement
            PSIMAGE_CUT_VERTICAL(U16);
            PSIMAGE_CUT_VERTICAL(U32); // Not a requirement
            PSIMAGE_CUT_VERTICAL(U64); // Not a requirement
            PSIMAGE_CUT_VERTICAL(S8);
            PSIMAGE_CUT_VERTICAL(S16); // Not a requirement
            PSIMAGE_CUT_VERTICAL(S32); // Not a requirement
            PSIMAGE_CUT_VERTICAL(S64); // Not a requirement
            PSIMAGE_CUT_VERTICAL(F32);
            PSIMAGE_CUT_VERTICAL(F64);
            PSIMAGE_CUT_VERTICAL(C32); // Not a requirement
            PSIMAGE_CUT_VERTICAL(C64); // Not a requirement
        default: {
                char* typeStr;
                PS_TYPE_NAME(typeStr,type);
                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                        PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
                        typeStr);
                psFree(out);
                out = NULL;
            }
        }
        psFree(imgVec);
        psFree(maskVec);
    } else if (direction == PS_CUT_Y_POS || direction == PS_CUT_Y_NEG) {
        // Cut in Y direction
        psVector* imgVec = NULL;
        psVector* maskVec = NULL;
        psS32 elementSize = PSELEMTYPE_SIZEOF(type);
        psU32* outPosition = NULL;

        // fill in psVector to fake out the statistics functions.
        imgVec = psAlloc(sizeof(psVector));
        imgVec->type = in->type;
        imgVec->n = imgVec->nalloc = numCols;
        if (mask != NULL) {
            maskVec = psAlloc(sizeof(psVector));
            maskVec->type = mask->type;
            maskVec->n = maskVec->nalloc = numCols;
        }
        // recycle output to make a proper sized/type output structure
        // n.b. type is double as that is the type given for all stats in
        // psStats.
        out = psVectorRecycle(out, numRows, PS_TYPE_F64);
        if (slicePositions != NULL) {
            slicePositions = psVectorRecycle(slicePositions, numRows, PS_TYPE_U32);
            outPosition = slicePositions->data.U32;
        }
        outData = out->data.F64;
        if (delta < 0) {
            outData += numRows-1;
            if (outPosition != NULL) {
                outPosition += numRows-1;
            }
        }

        for (psS32 r = row0; r < row1; r++) {
            // point the vector struct to the
            // data to calculate the stats
            imgVec->data.U8 = (psPtr )(in->data.U8[r] + col0 * elementSize);
            if (maskVec != NULL) {
                maskVec->data.U8 = (psPtr )(mask->data.U8[r] + col0 * sizeof(psMaskType));
            }
            myStats = psVectorStats(myStats, imgVec, NULL, maskVec, maskVal);
            (void)p_psGetStatValue(myStats, &statVal);  // we know it works cause we tested it above
            *outData = statVal;
            if (outPosition != NULL) {
                *outPosition = r;
                outPosition += delta;

            }
            outData += delta;
        }
        psFree(imgVec);
        psFree(maskVec);
    } else { // don't know what the direction flag is
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_SLICE_DIRECTION_INVALID,
                direction);
        psFree(out);
        out = NULL;
    }

    psFree(myStats);

    return out;
}

psVector* psImageCut(psVector* out,
                     psVector* cutCols,
                     psVector* cutRows,
                     const psImage* in,
                     const psImage* restrict mask,
                     psU32 maskVal,
                     float startCol,
                     float startRow,
                     float endCol,
                     float endRow,
                     psU32 nSamples,
                     psImageInterpolateMode mode)
{
    if (in == NULL || in->data.V == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        psFree(out);
        return NULL;
    }
    psS32 numCols = in->numCols;
    psS32 numRows = in->numRows;

    if (nSamples < 2) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_nSamples_TOOSMALL,
                nSamples);
        psFree(out);
        return NULL;
    }

    if (startCol < 0 || startCol >= numCols ||
            startRow < 0 || startRow >= numRows ||
            endCol < 0 || endCol >= numCols ||
            endRow < 0 || endRow >= numRows) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_LINE_NOT_IN_IMAGE,
                startCol,startRow,endCol,endRow,
                numCols-1,numRows-1);
        psFree(out);
        return NULL;
    }

    if (mode < PS_INTERPOLATE_FLAT || mode >= PS_INTERPOLATE_NUM_MODES) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED,
                mode);
        psFree(out);
        return NULL;
    }

    if (mask != NULL) {
        if (numRows != mask->numRows || numCols != mask->numCols) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    PS_ERRORTEXT_psImage_IMAGE_MASK_SIZE,
                    mask->numCols,mask->numRows,
                    numCols-1, numRows);
            psFree(out);
            return NULL;
        }
        if (mask->type.type != PS_TYPE_MASK) {
            char* typeStr;
            PS_TYPE_NAME(typeStr,mask->type.type);
            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                    PS_ERRORTEXT_psImage_IMAGE_MASK_TYPE,
                    typeStr, PS_TYPE_MASK_NAME);
            psFree(out);
            return NULL;
        }
    }

    //resize the vectors for the coordinate output
    psF32* cutColsData = NULL;
    psF32* cutRowsData = NULL;
    if (cutCols != NULL) {
        cutCols = psVectorRecycle(cutCols, nSamples, PS_TYPE_F32);
        cutColsData = cutCols->data.F32;
    }
    if (cutRows != NULL) {
        cutRows = psVectorRecycle(cutRows, nSamples, PS_TYPE_F32);
        cutRowsData = cutRows->data.F32;
    }

    out = psVectorRecycle(out, nSamples, in->type.type);

    float dX = (endCol - startCol) / (float)(nSamples-1);
    float dY = (endRow - startRow) / (float)(nSamples-1);

    #define LINEAR_CUT_CASE(TYPE) \
case PS_TYPE_##TYPE: { \
        ps##TYPE* outData = out->data.TYPE; \
        for (psS32 i = 0; i < nSamples; i++) { \
            float x = startCol + (float)i*dX; \
            float y = startRow + (float)i*dY; \
            /* store off the location of the sample. */ \
            if (cutColsData != NULL) { \
                cutColsData[i] = x; \
            } \
            if (cutRowsData != NULL) { \
                cutRowsData[i] = y; \
            } \
            outData[i] = psImagePixelInterpolate(in,x,y,mask,maskVal,0,mode); \
        } \
    } \
    break;


    switch (in->type.type) {
        LINEAR_CUT_CASE(U8);
        LINEAR_CUT_CASE(U16);
        LINEAR_CUT_CASE(U32);
        LINEAR_CUT_CASE(U64);
        LINEAR_CUT_CASE(S8);
        LINEAR_CUT_CASE(S16);
        LINEAR_CUT_CASE(S32);
        LINEAR_CUT_CASE(S64);
        LINEAR_CUT_CASE(F32);
        LINEAR_CUT_CASE(F64);
        LINEAR_CUT_CASE(C32);
        LINEAR_CUT_CASE(C64);

    default: {
            char* typeStr;
            PS_TYPE_NAME(typeStr,in->type.type);
            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                    PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
                    typeStr);
            psFree(out);
            out = NULL;
        }
    }

    return out;
}

psVector* psImageRadialCut(psVector* out,
                           const psImage* in,
                           const psImage* restrict mask,
                           psU32 maskVal,
                           float centerCol,
                           float centerRow,
                           const psVector* radii,
                           const psStats* stats)
{
    double statVal;

    /* check the parameters */

    if (in == NULL || in->data.V == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        psFree(out);
        return NULL;
    }
    psS32 numCols = in->numCols;
    psS32 numRows = in->numRows;

    if (mask != NULL) {
        if (numRows != mask->numRows || numCols != mask->numCols) {
            psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                    PS_ERRORTEXT_psImage_IMAGE_MASK_SIZE,
                    mask->numCols,mask->numRows,
                    numCols, numRows);
            psFree(out);
            return NULL;
        }
        if (mask->type.type != PS_TYPE_MASK) {
            char* typeStr;
            PS_TYPE_NAME(typeStr,mask->type.type);
            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                    PS_ERRORTEXT_psImage_IMAGE_MASK_TYPE,
                    typeStr, PS_TYPE_MASK_NAME);
            psFree(out);
            return NULL;
        }
    }

    if (centerCol < 0 || centerCol >= numCols ||
            centerRow < 0 || centerRow >= numRows) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_CENTER_NOT_IN_IMAGE,
                centerCol, centerRow,
                numCols-1, numRows-1);
        psFree(out);
        return NULL;
    }

    if (radii == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_RADII_VECTOR_NULL);
        psFree(out);
        return NULL;
    }

    if (radii->n < 2) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                PS_ERRORTEXT_psImage_RADII_VECTOR_TOOSMALL,
                radii->n);
        psFree(out);
        return NULL;
    }

    if (stats == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_STAT_NULL);
        psFree(out);
        return NULL;
    }

    // verify that the stats struct specifies a
    // single stats operation
    if (p_psGetStatValue(stats, &statVal) == false) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                PS_ERRORTEXT_psImage_BAD_STAT,
                stats->options);
        psFree(out);
        return NULL;
    }

    /* completed checking the parameters */

    // size the output vector to proper size.
    psS32 numOut = radii->n - 1;
    out = psVectorRecycle(out, numOut, PS_TYPE_F64);
    psF64* outData = out->data.F64;

    psVector* rSqVec = psVectorCopy(NULL, radii, PS_TYPE_F32);
    psF32* rSq = rSqVec->data.F32;

    psS32 startRow = centerRow - rSq[numOut];
    psS32 endRow = centerRow + rSq[numOut];
    psS32 startCol = centerCol - rSq[numOut];
    psS32 endCol = centerCol + rSq[numOut];

    if (startRow < 0) {
        startRow = 0;
    }

    if (startCol < 0) {
        startCol = 0;
    }

    if (endRow >= numRows) {
        endRow = numRows - 1;
    }

    if (endCol >= numCols) {
        endCol = numCols - 1;
    }

    // Square the radii data
    for (psS32 d = 0; d <= numOut; d++) {
        rSq[d] *= rSq[d];
    }

    // create temporary vectors for the data binning step
    psVector** buffer = psAlloc(sizeof(psVector*)*numOut);
    psVector** bufferMask = psAlloc(sizeof(psVector*)*numOut);
    for (psS32 lcv = 0; lcv < numOut; lcv++) {
        // n.b. alloc enough for the data by making the vectors slightly larger
        // than the area of the region of interest.
        buffer[lcv] = psVectorAlloc(1+4*(rSq[lcv+1]-rSq[lcv]),
                                    in->type.type);
        buffer[lcv]->n = 0;

        bufferMask[lcv] = NULL;
        if (mask != NULL) {
            bufferMask[lcv] = psVectorAlloc(1+4*(rSq[lcv+1]-rSq[lcv]),
                                            PS_TYPE_MASK);
            bufferMask[lcv]->n = 0;
        }
    }

    float dX;
    float dY;
    float dist;
    for (psS32 row=startRow; row <= endRow; row++) {
        psF32* inRow = in->data.F32[row];
        psMaskType* maskRow = NULL;
        if (mask != NULL) {
            maskRow = mask->data.PS_TYPE_MASK_DATA[row];
        }
        for (psS32 col=startCol; col <= endCol; col++) {
            dX = centerCol - (float)col - 0.5f;
            dY = centerRow - (float)row - 0.5f;
            dist = dX*dX+dY*dY;
            for (psS32 r = 0; r < numOut; r++) {
                if (rSq[r] < dist && dist < rSq[r+1]) {
                    psS32 n = buffer[r]->n;
                    if (n == buffer[r]->nalloc) { // in case buffers already full, expand
                        buffer[r] = psVectorRealloc(buffer[r], n*2);
                        if (bufferMask[r] != NULL) {
                            bufferMask[r] = psVectorRealloc(bufferMask[r], n*2);
                        }
                    }

                    buffer[r]->data.F32[n] = inRow[col];
                    buffer[r]->n = n+1;

                    if (maskRow != NULL) {
                        bufferMask[r]->data.PS_TYPE_MASK_DATA[n] = maskRow[col];
                        bufferMask[r]->n = n+1;
                    }

                    break;
                }
            }
        }
    }

    psStats* myStats = psAlloc(sizeof(psStats));
    *myStats = *stats;

    for (psS32 r = 0; r < numOut; r++) {
        myStats = psVectorStats(myStats,buffer[r], NULL, bufferMask[r],maskVal);
        (void)p_psGetStatValue(myStats,&statVal);
        outData[r] = statVal;
    }

    psFree(myStats);

    for (psS32 lcv = 0; lcv < numOut; lcv++) {
        psFree(buffer[lcv]);
        psFree(bufferMask[lcv]);
    }
    psFree(buffer);
    psFree(bufferMask);
    psFree(rSqVec);
    return out;
}
