/** @file  psImagePixelExtract.c
 *
 *  @brief Contains basic image extraction operations, as specified in the
 *         PSLIB SDRS sections "Image Pixel Extractions".
 *
 *  @ingroup Image
 *
 *  @author Robert DeSonia, MHPCC
 *
 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-09-30 23:09:13 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
 *
 */

#include <string.h>

#include "psMemory.h"
#include "psImagePixelExtract.h"
#include "psError.h"

#include "psErrorText.h"


#define FUNC_MACRO_VECTOR_STORE_ROW(TYPE) \
static psVector *vectorStoreRow##TYPE(psVector *vec, const psImage *in, int row) \
{ \
    \
    for (int i = 0; i < in->numCols; i++) \
    { \
        vec->data.TYPE[i] = in->data.TYPE[i][row]; \
    } \
    return vec; \
} \

FUNC_MACRO_VECTOR_STORE_ROW(S8)
FUNC_MACRO_VECTOR_STORE_ROW(S16)
FUNC_MACRO_VECTOR_STORE_ROW(S32)
FUNC_MACRO_VECTOR_STORE_ROW(S64)
FUNC_MACRO_VECTOR_STORE_ROW(U8)
FUNC_MACRO_VECTOR_STORE_ROW(U16)
FUNC_MACRO_VECTOR_STORE_ROW(U32)
FUNC_MACRO_VECTOR_STORE_ROW(U64)
FUNC_MACRO_VECTOR_STORE_ROW(F32)
FUNC_MACRO_VECTOR_STORE_ROW(F64)
FUNC_MACRO_VECTOR_STORE_ROW(C32)
FUNC_MACRO_VECTOR_STORE_ROW(C64)

psVector *psImageRow(psVector *out,
                     const psImage *input,
                     int row)
{
    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psImage_IMAGE_NULL);
        return NULL;
    }
    if (row >= input->numRows) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                "Specified row number is out of range for specified image.\n");
        return NULL;
    }
    if (row < 0) {
        row += input->numRows;
        if (row < 0) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    "Specified row number is out of range for specified image.\n");
            return NULL;
        }
    }
    psVectorRecycle(out, input->numCols, input->type.type);

    switch (input->type.type) {
    case PS_TYPE_S8:
        out = vectorStoreRowS8(out, input, row);
        break;
    case PS_TYPE_S16:
        out = vectorStoreRowS16(out, input, row);
        break;
    case PS_TYPE_S32:
        out = vectorStoreRowS32(out, input, row);
        break;
    case PS_TYPE_S64:
        out = vectorStoreRowS64(out, input, row);
        break;
    case PS_TYPE_U8:
        out = vectorStoreRowU8(out, input, row);
        break;
    case PS_TYPE_U16:
        out = vectorStoreRowU16(out, input, row);
        break;
    case PS_TYPE_U32:
        out = vectorStoreRowU32(out, input, row);
        break;
    case PS_TYPE_U64:
        out = vectorStoreRowU64(out, input, row);
        break;
    case PS_TYPE_F32:
        out = vectorStoreRowF32(out, input, row);
        break;
    case PS_TYPE_F64:
        out = vectorStoreRowF64(out, input, row);
        break;
    case PS_TYPE_C32:
        out = vectorStoreRowC32(out, input, row);
        break;
    case PS_TYPE_C64:
        out = vectorStoreRowC64(out, input, row);
        break;
    default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                "Specified psImage has invalid type for this function.\n");
        return NULL;
    }

    return out;
}


#define FUNC_MACRO_VECTOR_STORE_COL(TYPE) \
static psVector *vectorStoreCol##TYPE(psVector *vec, const psImage *in, int col) \
{ \
    \
    for (int i = 0; i < in->numRows; i++) \
    { \
        vec->data.TYPE[i] = in->data.TYPE[col][i]; \
    } \
    return vec; \
} \

FUNC_MACRO_VECTOR_STORE_COL(S8)
FUNC_MACRO_VECTOR_STORE_COL(S16)
FUNC_MACRO_VECTOR_STORE_COL(S32)
FUNC_MACRO_VECTOR_STORE_COL(S64)
FUNC_MACRO_VECTOR_STORE_COL(U8)
FUNC_MACRO_VECTOR_STORE_COL(U16)
FUNC_MACRO_VECTOR_STORE_COL(U32)
FUNC_MACRO_VECTOR_STORE_COL(U64)
FUNC_MACRO_VECTOR_STORE_COL(F32)
FUNC_MACRO_VECTOR_STORE_COL(F64)
FUNC_MACRO_VECTOR_STORE_COL(C32)
FUNC_MACRO_VECTOR_STORE_COL(C64)


psVector *psImageCol(psVector *out,
                     const psImage *input,
                     int column)
{
    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psImage_IMAGE_NULL);
        return NULL;
    }
    if (column >= input->numCols) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                "Specified column number is out of range for specified image.\n");
        return NULL;
    }
    if (column < 0) {
        column += input->numCols;
        if (column < 0) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    "Specified column number is out of range for specified image.\n");
            return NULL;
        }
    }


    psVectorRecycle(out, input->numRows, input->type.type);

    switch (input->type.type) {
    case PS_TYPE_S8:
        out = vectorStoreColS8(out, input, column);
        break;
    case PS_TYPE_S16:
        out = vectorStoreColS16(out, input, column);
        break;
    case PS_TYPE_S32:
        out = vectorStoreColS32(out, input, column);
        break;
    case PS_TYPE_S64:
        out = vectorStoreColS64(out, input, column);
        break;
    case PS_TYPE_U8:
        out = vectorStoreColU8(out, input, column);
        break;
    case PS_TYPE_U16:
        out = vectorStoreColU16(out, input, column);
        break;
    case PS_TYPE_U32:
        out = vectorStoreColU32(out, input, column);
        break;
    case PS_TYPE_U64:
        out = vectorStoreColU64(out, input, column);
        break;
    case PS_TYPE_F32:
        out = vectorStoreColF32(out, input, column);
        break;
    case PS_TYPE_F64:
        out = vectorStoreColF64(out, input, column);
        break;
    case PS_TYPE_C32:
        out = vectorStoreColC32(out, input, column);
        break;
    case PS_TYPE_C64:
        out = vectorStoreColC64(out, input, column);
        break;
    default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                "Specified psImage has invalid type for this function.\n");
        return NULL;
    }

    return out;

}


psVector* psImageSlice(psVector* out,
                       psVector* coords,
                       const psImage* restrict input,
                       const psImage* restrict mask,
                       psMaskType maskVal,
                       psRegion region,
                       psImageCutDirection direction,
                       const psStats* stats)
{
    double statVal;
    psStats* myStats;
    psElemType type;
    psS32 inRows;
    psS32 inCols;
    psS32 delta = 1;
    psF64* outData;
    psS32 row0 = region.y0;
    psS32 row1 = region.y1;
    psS32 col0 = region.x0;
    psS32 col1 = region.x1;

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

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

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

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

    type = input->type.type;
    inRows = input->numRows;
    inCols = input->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 (coords != NULL) {
            coords = psVectorRecycle(coords, numCols, PS_TYPE_U32);
            outPosition = coords->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 = input->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 = input->type;
        imgVec->n = *(int*)&imgVec->nalloc = numCols;
        if (mask != NULL) {
            maskVec = psAlloc(sizeof(psVector));
            maskVec->type = mask->type;
            maskVec->n = *(int*)&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 (coords != NULL) {
            coords = psVectorRecycle(coords, numRows, PS_TYPE_U32);
            outPosition = coords->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 )(input->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* input,
                     const psImage* mask,
                     psMaskType maskVal,
                     psRegion region,
                     unsigned int nSamples,
                     psImageInterpolateMode mode)
{

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

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

    float startCol = region.x0;
    float startRow = region.y0;
    float endCol = region.x1;
    float endRow = region.y1;
    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 ) {
        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, input->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(input,x,y,mask,maskVal,0,mode); \
        } \
    } \
    break;


    switch (input->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,input->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* input,
                           const psImage* restrict mask,
                           psMaskType maskVal,
                           float x,
                           float y,
                           const psVector* radii,
                           const psStats* stats)
{
    double statVal;

    /* check the parameters */

    if (input == NULL || input->data.V == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        psFree(out);
        return NULL;
    }
    psS32 numCols = input->numCols;
    psS32 numRows = input->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 (x < 0 || x >= numCols ||
            y < 0 || y >= numRows) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_CENTER_NOT_IN_IMAGE,
                x, y,
                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 = y - rSq[numOut];
    psS32 endRow = y + rSq[numOut];
    psS32 startCol = x - rSq[numOut];
    psS32 endCol = x + 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]),
                                    input->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 = input->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 = x - (float)col - 0.5f;
            dY = y - (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;
}
