/** @file  psImageGeomManip.c
 *
 *  @brief Contains basic image pixel and geometry manipulation operations, as
 *         specified in the PSLIB SDRS sections "Image Pixel Manipulations" and
 *         "Image Geometry Manipulations".
 *
 *  @ingroup Image
 *
 *  @author Robert DeSonia, MHPCC
 *  @author Ross Harman, MHPCC
 *
 *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-04-06 22:55:18 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
 */

#include <complex.h>
#include <math.h>                          // for isfinite(), etc.
#include <stdlib.h>
#include <string.h>                        // for memcpy, etc.

#include "psImageGeomManip.h"

#include "psError.h"
#include "psImage.h"
#include "psImageStructManip.h"
#include "psStats.h"
#include "psMemory.h"
#include "psConstants.h"
#include "psErrorText.h"
#include "psCoord.h"

psImage* psImageRebin(psImage* out,
                      const psImage* in,
                      const psImage* mask,
                      psMaskType maskVal,
                      int scale,
                      const psStats* stats)
{
    psS32 inRows;
    psS32 inCols;
    psS32 outRows;
    psS32 outCols;
    psVector* vec;                     // vector to hold the values of a single bin.
    psVector* maskVec = NULL;          // vector to hold the mask of a single bin.
    psMaskType* maskData = NULL;
    psStats* myStats;
    double statVal;

    if (in == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        psFree(out);
        return NULL;
    }

    if (scale < 1) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImageManip_SCALE_NOT_POSITIVE,
                scale);
        psFree(out);
        return NULL;
    }

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

    if (p_psGetStatValue(stats, &statVal) == false) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_BAD_STAT,
                stats->options);
        psFree(out);
        return NULL;
    }

    vec = psVectorAlloc(scale * scale, in->type.type);

    if (mask != 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);
            psFree(vec);
            return NULL;
        }
        maskVec = psVectorAlloc(scale * scale, PS_TYPE_MASK);
        maskVec->n = maskVec->nalloc;
        maskData = maskVec->data.PS_TYPE_MASK_DATA;
    }

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

    // create output image.
    inRows = in->numRows;
    inCols = in->numCols;
    outRows = (inRows + scale - 1) / scale;     // round-up for remainders
    outCols = (inCols + scale - 1) / scale;     // round-up for remainders
    out = psImageRecycle(out, outCols, outRows, in->type.type);

    #define PS_IMAGE_REBIN_CASE(type) \
case PS_TYPE_##type: { \
        ps##type* outRowData; \
        ps##type* vecData = vec->data.type; \
        psMaskType* inRowMask = NULL; \
        for (psS32 row = 0; row < outRows; row++) { \
            outRowData = out->data.type[row]; \
            psS32 inCurrentRow = row*scale; \
            psS32 inNextRow = (row+1)*scale; \
            for (psS32 col = 0; col < outCols; col++) { \
                psS32 inCurrentCol = col*scale; \
                psS32 inNextCol = (col+1)*scale; \
                psS32 n = 0; \
                for (psS32 inRow = inCurrentRow; inRow < inNextRow && inRow < inRows; inRow++) { \
                    ps##type* inRowData = in->data.type[inRow]; \
                    if (mask != NULL) { \
                        inRowMask = mask->data.PS_TYPE_MASK_DATA[inRow]; \
                    } \
                    for (psS32 inCol = inCurrentCol; inCol < inNextCol && inCol < inCols; inCol++) { \
                        if (maskData != NULL) { \
                            maskData[n] = inRowMask[inCol]; \
                        } \
                        vecData[n++] = inRowData[inCol]; \
                    } \
                } \
                vec->n = n; \
                myStats = psVectorStats(myStats,vec,NULL,maskVec,maskVal); \
                p_psGetStatValue(myStats,&statVal); \
                outRowData[col] = (ps##type)statVal; \
            } \
        } \
    } \
    break;

    switch (in->type.type) {
        //        PS_IMAGE_REBIN_CASE(U8);       Not valid since psVectorStats doesn't allow
        PS_IMAGE_REBIN_CASE(U16);
        PS_IMAGE_REBIN_CASE(U32);      // Not a requirement
        PS_IMAGE_REBIN_CASE(U64);      // Not a requirement
        PS_IMAGE_REBIN_CASE(S8);
        //        PS_IMAGE_REBIN_CASE(S16);      Not valid since psVectorStats doesn't allow
        PS_IMAGE_REBIN_CASE(S32);      // Not a requirement
        PS_IMAGE_REBIN_CASE(S64);      // Not a requirement
        PS_IMAGE_REBIN_CASE(F32);
        PS_IMAGE_REBIN_CASE(F64);
        //        PS_IMAGE_REBIN_CASE(C32);      Not valid since psVectorStats doesn't allow
        //        PS_IMAGE_REBIN_CASE(C64);      Not valid since psVectorStats doesn't allow

    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;
        }
    }

    psFree(vec);
    psFree(maskVec);
    psFree(myStats);

    return out;
}

psImage* psImageResample(psImage* out,
                         const psImage* in,
                         int scale,
                         psImageInterpolateMode mode)
{
    psS32 outRows;
    psS32 outCols;
    float invScale;

    if (in == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        psFree(out);
        return NULL;
    }

    if (scale < 1) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImageManip_SCALE_NOT_POSITIVE,
                scale);
        psFree(out);
        return NULL;
    }

    if (mode > PS_INTERPOLATE_LANCZOS4_VARIANCE ) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED,
                mode);
        psFree(out);
        return NULL;
    }

    // create an output image of the same size
    // and type
    outRows = in->numRows * scale;
    outCols = in->numCols * scale;
    invScale = 1.0f / (float)scale;

    #define PSIMAGE_RESAMPLE_CASE(TYPE) \
case PS_TYPE_##TYPE: { \
        out = psImageRecycle(out,outCols, outRows, PS_TYPE_##TYPE); \
        for (psS32 row=0;row<outRows;row++) { \
            ps##TYPE* rowData = out->data.TYPE[row]; \
            float inRow = (float)row * invScale; \
            for (psS32 col=0;col<outCols;col++) { \
                rowData[col] = psImagePixelInterpolate(in,(float)col*invScale,inRow,NULL,0,0,mode); \
            } \
        }  \
        break; \
    }

    switch (in->type.type) {
        PSIMAGE_RESAMPLE_CASE(U8)
        PSIMAGE_RESAMPLE_CASE(U16)
        PSIMAGE_RESAMPLE_CASE(U32)
        PSIMAGE_RESAMPLE_CASE(U64)
        PSIMAGE_RESAMPLE_CASE(S8)
        PSIMAGE_RESAMPLE_CASE(S16)
        PSIMAGE_RESAMPLE_CASE(S32)
        PSIMAGE_RESAMPLE_CASE(S64)
        PSIMAGE_RESAMPLE_CASE(F32)
        PSIMAGE_RESAMPLE_CASE(F64)
        PSIMAGE_RESAMPLE_CASE(C32)
        PSIMAGE_RESAMPLE_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;
}

psImage* psImageRoll(psImage* out,
                     const psImage* input,
                     int dx,
                     int dy)
{
    psS32 outRows;
    psS32 outCols;
    psS32 elementSize;

    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        psFree(out);
        return NULL;
    }
    // create an output image of the same size
    // and type
    outRows = input->numRows;
    outCols = input->numCols;
    elementSize = PSELEMTYPE_SIZEOF(input->type.type);
    out = psImageRecycle(out, outCols, outRows, input->type.type);

    // make dx and dy between 0 and outCols or
    // outRows, respectively
    dx = dx % outCols;
    dy = dy % outRows;
    if (dx < 0) {
        dx += outCols;
    }
    if (dy < 0) {
        dy += outRows;
    }

    psS32 segment1Size = elementSize * (outCols - dx);
    psS32 segment2Size = elementSize * dx;

    for (psS32 row = 0; row < outRows; row++) {
        psS32 inRowNumber = row + dy;

        if (inRowNumber >= outRows) {
            inRowNumber -= outRows;
        }
        psU8* inRow = input->data.U8[inRowNumber]; // use byte arithmetic for all types
        psU8* outRow = out->data.U8[row];

        memcpy(outRow, inRow + segment2Size, segment1Size);
        memcpy(outRow + segment1Size, inRow, segment2Size);
    }

    return out;
}

psImage* psImageRotate(psImage* out,
                       const psImage* input,
                       float angle,
                       double complex exposed,
                       psImageInterpolateMode mode)
{
    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        psFree(out);
        return NULL;
    }
    // put the angle in the range of 0...2PI.
    angle = (float)((double)angle - (2.0*M_PI) * floor(angle / (2.0*M_PI)));

    if (fabsf(angle - M_PI_2) < FLT_EPSILON) {
        // perform 1/4 rotate counter-clockwise
        psS32 numRows = input->numCols;
        psS32 numCols = input->numRows;
        psS32 lastCol = numCols - 1;
        psElemType type = input->type.type;

        out = psImageRecycle(out, numCols, numRows, type);

        #define PSIMAGE_ROTATE_LEFT_90(TYPE) \
    case PS_TYPE_##TYPE: { \
            ps##TYPE** inData = input->data.TYPE; \
            for (psS32 row=0;row<numRows;row++) { \
                ps##TYPE* outRow = out->data.TYPE[row]; \
                for (psS32 col=0;col<numCols;col++) { \
                    outRow[col] = inData[lastCol-col][row]; \
                } \
            } \
        } \
        break;

        switch (type) {
            PSIMAGE_ROTATE_LEFT_90(U8);
            PSIMAGE_ROTATE_LEFT_90(U16);
            PSIMAGE_ROTATE_LEFT_90(U32);    //  Not a requirement
            PSIMAGE_ROTATE_LEFT_90(U64);    //  Not a requirement
            PSIMAGE_ROTATE_LEFT_90(S8);
            PSIMAGE_ROTATE_LEFT_90(S16);
            PSIMAGE_ROTATE_LEFT_90(S32);    //  Not a requirement
            PSIMAGE_ROTATE_LEFT_90(S64);    //  Not a requirement
            PSIMAGE_ROTATE_LEFT_90(F32);
            PSIMAGE_ROTATE_LEFT_90(F64);
            PSIMAGE_ROTATE_LEFT_90(C32);
            PSIMAGE_ROTATE_LEFT_90(C64);

        default: {
                char* typeStr;
                PS_TYPE_NAME(typeStr,type);
                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                        PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
                        typeStr);
                psFree(out);
                return NULL;
            }
        }
    } else if (fabsf(angle - M_PI) < FLT_EPSILON) {
        // perform 1/2 rotate
        psS32 numRows = input->numRows;
        psS32 lastRow = numRows - 1;
        psS32 numCols = input->numCols;
        psS32 lastCol = numCols - 1;
        psElemType type = input->type.type;

        out = psImageRecycle(out, numCols, numRows, type);

        #define PSIMAGE_ROTATE_180_CASE(TYPE) \
    case PS_TYPE_##TYPE: { \
            for (psS32 row=0;row<numRows;row++) { \
                ps##TYPE* outRow = out->data.TYPE[row]; \
                ps##TYPE* inRow = input->data.TYPE[lastRow-row]; \
                for (psS32 col=0;col<numCols;col++) { \
                    outRow[col] = inRow[lastCol - col]; \
                } \
            } \
        } \
        break;

        switch (type) {
            PSIMAGE_ROTATE_180_CASE(U8);
            PSIMAGE_ROTATE_180_CASE(U16);
            PSIMAGE_ROTATE_180_CASE(U32);    // Not a requirement
            PSIMAGE_ROTATE_180_CASE(U64);    // Not a requirement
            PSIMAGE_ROTATE_180_CASE(S8);
            PSIMAGE_ROTATE_180_CASE(S16);
            PSIMAGE_ROTATE_180_CASE(S32);    // Not a requirement
            PSIMAGE_ROTATE_180_CASE(S64);    // Not a requirement
            PSIMAGE_ROTATE_180_CASE(F32);
            PSIMAGE_ROTATE_180_CASE(F64);
            PSIMAGE_ROTATE_180_CASE(C32);
            PSIMAGE_ROTATE_180_CASE(C64);

        default: {
                char* typeStr;
                PS_TYPE_NAME(typeStr,type);
                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                        PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
                        typeStr);
                psFree(out);
                return NULL;
            }
        }
    } else if (fabsf(angle - (M_PI+M_PI_2)) < FLT_EPSILON) {
        // perform 1/4 rotate clockwise
        psS32 numRows = input->numCols;
        psS32 lastRow = numRows - 1;
        psS32 numCols = input->numRows;
        psElemType type = input->type.type;

        out = psImageRecycle(out, numCols, numRows, type);

        #define PSIMAGE_ROTATE_RIGHT_90(TYPE) \
    case PS_TYPE_##TYPE: { \
            ps##TYPE** inData = input->data.TYPE; \
            for (psS32 row=0;row<numRows;row++) { \
                ps##TYPE* outRow = out->data.TYPE[row]; \
                for (psS32 col=0;col<numCols;col++) { \
                    outRow[col] = inData[col][lastRow-row]; \
                } \
            } \
        } \
        break;

        switch (type) {
            PSIMAGE_ROTATE_RIGHT_90(U8);
            PSIMAGE_ROTATE_RIGHT_90(U16);
            PSIMAGE_ROTATE_RIGHT_90(U32);     // Not a requirement
            PSIMAGE_ROTATE_RIGHT_90(U64);     // Not a requirement
            PSIMAGE_ROTATE_RIGHT_90(S8);
            PSIMAGE_ROTATE_RIGHT_90(S16);
            PSIMAGE_ROTATE_RIGHT_90(S32);     // Not a requirement
            PSIMAGE_ROTATE_RIGHT_90(S64);     // Not a requirement
            PSIMAGE_ROTATE_RIGHT_90(F32);
            PSIMAGE_ROTATE_RIGHT_90(F64);
            PSIMAGE_ROTATE_RIGHT_90(C32);
            PSIMAGE_ROTATE_RIGHT_90(C64);

        default: {
                char* typeStr;
                PS_TYPE_NAME(typeStr,type);
                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                        PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
                        typeStr);
                psFree(out);
                return NULL;
            }
        }
    } else if (fabsf(angle) < FLT_EPSILON) {
        out = psImageCopy(out, input, input->type.type);
    } else {
        psElemType type = input->type.type;
        psS32 numRows = input->numRows;
        psS32 numCols = input->numCols;
        float centerX = (float)(numCols) / 2.0f;
        float centerY = (float)(numRows) / 2.0f;
        double cosT = cosf(angle);
        double sinT = sinf(angle);

        // calculate the corners of the rotated image so we know the proper output image size.
        // x' = x cos(t) + y sin(t); i.e, x' = (x-centerX)*cosT + (y-centerY)*sinT;
        // y' = y cos(t) - x sin(t); i.e. y' = (y-centerY)*cosT - (x-centerX)*sinT;

        psS32 outCols = ceil(abs(numCols * cosT) + abs(numRows * sinT)) + 1;
        psS32 outRows = ceil(abs(numCols * sinT) + abs(numRows * cosT)) + 1;
        float minX = (float)outCols / -2.0f;
        psS32 intMinY = outRows / -2;

        out = psImageRecycle(out, outCols, outRows, type);

        /* optimized public domain rotation routine by Karl Lager
         *
         * float cosT,sinT;
         * cosT = cos(t);
         * sinT = sin(t);
         * for (y = min_y; y <= max_y; y++) {
         *     x' = min_x * cosT + y * sinT + x1';
         *     y' = y * cosT - min_x * sinT + y1';
         *     for (x = min_x; x <= max_x; x++) {
         *         if (x', y') is in the bounds of the bitmap, get pixel
         *            (x', y') and plot the pixel to (x, y) on screen.
         *         x' += cosT;
         *         y' -= sinT;
         *     }
         * }
         */

        // precalculate some figures that are used within loop
        float minXTimesCosTPlusCenterX = minX * cosT + centerX;
        float CenterYMinusminXTimesSinT = centerY - minX * sinT;

        #define PSIMAGE_ROTATE_ARBITRARY_LOOP(TYPE,MODE) { \
            if (creal(exposed) < PS_MIN_##TYPE || \
                    creal(exposed) > PS_MAX_##TYPE || \
                    cimag(exposed) < PS_MIN_##TYPE || \
                    cimag(exposed) > PS_MAX_##TYPE) { \
                psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                        PS_ERRORTEXT_psImageManip_CLIP_VALUE_INVALID, \
                        "exposed", \
                        creal(exposed),cimag(exposed), \
                        PS_TYPE_##TYPE##_NAME,  \
                        (double)PS_MIN_##TYPE,(double)PS_MAX_##TYPE); \
                psFree(out); \
                out = NULL; \
                break; \
            } \
            float inX; \
            float inY; \
            ps##TYPE* outRow; \
            for (psS32 y = 0; y < outRows; y++) { \
                inX = minXTimesCosTPlusCenterX + (y+intMinY) * sinT; \
                inY = CenterYMinusminXTimesSinT + (y+intMinY) * cosT; \
                outRow = out->data.TYPE[y]; \
                for (psS32 x = 0; x < outCols; x++) { \
                    outRow[x] = p_psImagePixelInterpolate##MODE##_##TYPE(input,inX,inY,NULL,0,exposed); \
                    inX += cosT; \
                    inY -= sinT; \
                } \
            } \
        }

        #define PSIMAGE_ROTATE_ARBITRARY_CASE(MODE) \
    case PS_INTERPOLATE_##MODE: \
        switch (type) { \
        case PS_TYPE_U8: \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(U8,MODE); \
            break; \
        case PS_TYPE_U16: \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(U16,MODE); \
            break; \
        case PS_TYPE_U32:   /* Not a requirement */ \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(U32,MODE); \
            break; \
        case PS_TYPE_U64:   /* Not a requirement */ \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(U64,MODE); \
            break;  \
        case PS_TYPE_S8: \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(S8,MODE); \
            break; \
        case PS_TYPE_S16: \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(S16,MODE); \
            break; \
        case PS_TYPE_S32:   /* Not a requirement */ \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(S32,MODE); \
            break; \
        case PS_TYPE_S64:   /* Not a requirement */ \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(S64,MODE); \
            break; \
        case PS_TYPE_F32: \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(F32,MODE); \
            break; \
        case PS_TYPE_F64: \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(F64,MODE); \
            break; \
        case PS_TYPE_C32: \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(C32,MODE); \
            break; \
        case PS_TYPE_C64: \
            PSIMAGE_ROTATE_ARBITRARY_LOOP(C64,MODE); \
            break; \
        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; \
            } \
        } \
        break;

        switch (mode) {
            PSIMAGE_ROTATE_ARBITRARY_CASE(FLAT);
            PSIMAGE_ROTATE_ARBITRARY_CASE(BILINEAR);
            PSIMAGE_ROTATE_ARBITRARY_CASE(BILINEAR_VARIANCE);
        default:
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED,
                    mode);
            psFree(out);
            out = NULL;
        }
    }

    return out;
}

psImage* psImageShift(psImage* out,
                      const psImage* input,
                      float dx,
                      float dy,
                      double complex exposed,
                      psImageInterpolateMode mode)
{
    psS32 outRows;
    psS32 outCols;
    psS32 elementSize;
    psElemType type;

    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        psFree(out);
        return NULL;
    }
    // create an output image of the same size
    // and type
    outRows = input->numRows;
    outCols = input->numCols;
    type = input->type.type;
    elementSize = PSELEMTYPE_SIZEOF(type);
    out = psImageRecycle(out, outCols, outRows, type);

    #define PSIMAGE_SHIFT_CASE(MODE,TYPE) \
case PS_TYPE_##TYPE: \
    if (creal(exposed) < PS_MIN_##TYPE || \
            creal(exposed) > PS_MAX_##TYPE || \
            cimag(exposed) < PS_MIN_##TYPE || \
            cimag(exposed) > PS_MAX_##TYPE) { \
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                PS_ERRORTEXT_psImageManip_CLIP_VALUE_INVALID, \
                "exposed", \
                creal(exposed),cimag(exposed), \
                PS_TYPE_##TYPE##_NAME,  \
                (double)PS_MIN_##TYPE,(double)PS_MAX_##TYPE); \
        psFree(out); \
        out = NULL; \
        break; \
    } \
    for (psS32 row=0;row<outRows;row++) { \
        ps##TYPE* outRow = out->data.TYPE[row]; \
        float y = dy+(float)row; \
        for (psS32 col=0;col<outCols;col++) { \
            outRow[col] = p_psImagePixelInterpolate##MODE##_##TYPE( \
                          input,dx+(float)col,y,NULL,0,exposed); \
        } \
    } \
    break;

    #define PSIMAGE_SHIFT_ARBITRARY_CASE(MODE) \
case PS_INTERPOLATE_##MODE: \
    switch (input->type.type) { \
        PSIMAGE_SHIFT_CASE(MODE,U8);  \
        PSIMAGE_SHIFT_CASE(MODE,U16); \
        PSIMAGE_SHIFT_CASE(MODE,U32);     /* Not a requirement */ \
        PSIMAGE_SHIFT_CASE(MODE,U64);     /* Not a requirement */ \
        PSIMAGE_SHIFT_CASE(MODE,S8);  \
        PSIMAGE_SHIFT_CASE(MODE,S16); \
        PSIMAGE_SHIFT_CASE(MODE,S32);    /* Not a requirement */ \
        PSIMAGE_SHIFT_CASE(MODE,S64);    /* Not a requirement */ \
        PSIMAGE_SHIFT_CASE(MODE,F32); \
        PSIMAGE_SHIFT_CASE(MODE,F64); \
        PSIMAGE_SHIFT_CASE(MODE,C32); \
        PSIMAGE_SHIFT_CASE(MODE,C64); \
        \
    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; \
        } \
    } \
    break;

    switch (mode) {
        PSIMAGE_SHIFT_ARBITRARY_CASE(FLAT);
        PSIMAGE_SHIFT_ARBITRARY_CASE(BILINEAR);
        PSIMAGE_SHIFT_ARBITRARY_CASE(BILINEAR_VARIANCE);
    default:
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED,
                mode);
        psFree(out);
        out = NULL;
    }

    return out;
}

psImage* psImageTransform(psImage *output,
                          psPixels** blankPixels,
                          const psImage *input,
                          const psImage *inputMask,
                          psMaskType inputMaskVal,
                          const psPlaneTransform *outToIn,
                          psRegion region,
                          const psPixels* pixels,
                          psImageInterpolateMode mode,
                          double exposedValue)
{
    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        psFree(output);
        return NULL;
    }
    psElemType type = input->type.type;

    if (blankPixels != NULL && *blankPixels == NULL) {
        *blankPixels = psPixelsAlloc(0);
    }

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

    if (outToIn == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImageManip_TRANSFORM_NULL);
        return NULL;
    }

    int row0;
    int row1;
    int col0;
    int col1;
    int numRows;
    int numCols;
    if (output == NULL) { // output image size is determined by psRegion
        row0 = region.y0;
        row1 = region.y1;
        col0 = region.x0;
        col1 = region.x1;
        //If [0,0,0,0] specified, the whole image is to be included
        if (row0 == 0 && col0 == 0 && row1 == 0 && col1 == 0) {
            row0 = input->row0;
            col0 = input->col0;
            row1 = input->row0 + input->numRows - 1;
            col1 = input->col0 + input->numCols - 1;
        }
        if (col1 < 1) {
            col1 += input->col0 + input->numCols;
        }

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

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

        if (numRows < 1 || numCols < 1) {
            psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                    "The specified region is invalid.");
            psFree(output);
            return NULL;
        }
        // create the output image.
        output = psImageRecycle(output, numCols, numRows, input->type.type);
        *(psS32*)&output->col0 = region.x0;
        *(psS32*)&output->row0 = region.y0;
    } else { // size of output is determined by output parameter
        numRows = output->numRows;
        numCols = output->numCols;
        row0 = output->row0;
        col0 = output->col0;
        row1 = row0+numRows -1;
        col1 = col0+numCols -1;
    }

    // loop through the output image using the domain above and transform
    // each output pixel to input coordinates and use psImagePixelInterpolate
    // to determine the pixel value.
    psPlane outPosition;
    psPlane* inPosition = NULL;

    #define PSIMAGE_TRANSFORM_DOTRANSFORM(TYPE,MODE) \
    /* apply the transform to get the position in the input image */ \
    inPosition = psPlaneTransformApply(inPosition, outToIn, &outPosition); \
    \
    if (inPosition == NULL) { \
        psError(PS_ERR_UNKNOWN, false, \
                "Failed to apply the transform"); \
        psFree(output); \
        return NULL; \
    } \
    /* interpolate the cooresponding input pixel to get the output pixel value. */ \
    ps##TYPE value = p_psImagePixelInterpolate##MODE##_##TYPE(input, \
                     inPosition->x, inPosition->y, \
                     inputMask, inputMaskVal, NAN); \
    /*    psFree(inPosition); */\
    if (isnan(value)) { \
        if (blankPixels != NULL) { \
            p_psPixelsAppend(*blankPixels, (*blankPixels)->nalloc, outPosition.x, outPosition.y); \
        } \
        value = exposedValue; \
    } \

    #define PSIMAGE_TRANSFORM_LOOP(TYPE, MODE) { \
        for (int row = 0; row < numRows; row++) { \
            outPosition.y = row+row0; \
            ps##TYPE* outputData=output->data.TYPE[row]; \
            for (int col = 0; col < numCols; col++) { \
                outPosition.x = col+col0; \
                PSIMAGE_TRANSFORM_DOTRANSFORM(TYPE,MODE) \
                outputData[col] = value; \
            } \
        } \
    }

    #define PSIMAGE_TRANSFORM_FROMLIST(TYPE, MODE) { \
        int n = pixels->n; \
        for (int i= 0; i < n; i++) { \
            int x = pixels->data[i].x; \
            int y = pixels->data[i].y; \
            if (x >= col0 && x < col1 && y >= row0 && y < row1) { \
                outPosition.x = x; \
                outPosition.y = y; \
                PSIMAGE_TRANSFORM_DOTRANSFORM(TYPE,MODE) \
                output->data.TYPE[y][x] = value; \
            } \
        } \
    }

    #define PSIMAGE_TRANSFORM_CASE(MODE) \
case PS_INTERPOLATE_##MODE: \
    switch (type) { \
    case PS_TYPE_F32: \
        PSIMAGE_TRANSFORM_LOOP(F32,MODE); \
        break; \
    case PS_TYPE_F64: \
        PSIMAGE_TRANSFORM_LOOP(F64,MODE); \
        break; \
    default: { \
            char* typeStr; \
            PS_TYPE_NAME(typeStr,type); \
            psError(PS_ERR_BAD_PARAMETER_TYPE, true, \
                    PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, \
                    typeStr); \
            psFree(output); \
            return NULL; \
        } \
    } \
    break;

    switch (mode) {
        PSIMAGE_TRANSFORM_CASE(FLAT);
        PSIMAGE_TRANSFORM_CASE(BILINEAR);
        PSIMAGE_TRANSFORM_CASE(BILINEAR_VARIANCE);
    default:
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED,
                mode);
        psFree(output);
        return NULL;
    }

    psFree(inPosition);

    return output;
}

