/** @file  psImageManip.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.35 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-01-18 03:15:03 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */

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

#include "psError.h"
#include "psImage.h"
#include "psStats.h"
#include "psMemory.h"
#include "psImageExtraction.h"
#include "psConstants.h"
#include "psImageErrors.h"

psS32 psImageClip(psImage* input,
                  psF64 min,
                  psF64 vmin,
                  psF64 max,
                  psF64 vmax)
{
    psS32 numClipped = 0;
    psU32 numRows;
    psU32 numCols;

    if (input == NULL) {
        return 0;
    }

    if (max < min) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImageManip_MAXMIN,
                (double)min,(double)max);
        return 0;
    }

    numRows = input->numRows;
    numCols = input->numCols;

    switch (input->type.type) {

        #define psImageClipCase(type) \
    case PS_TYPE_##type: { \
            if (vmin < PS_MIN_##type || vmin > PS_MAX_##type) { \
                psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                        PS_ERRORTEXT_psImage_PARAMETER_OUTOF_TYPERANGE, \
                        "vmin",vmin, PS_TYPE_##type##_NAME, \
                        (psF64)PS_MIN_##type,(psF64)PS_MAX_##type); \
            } \
            if (vmax > PS_MAX_##type || vmax < PS_MIN_##type) { \
                psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                        PS_ERRORTEXT_psImage_PARAMETER_OUTOF_TYPERANGE, \
                        "vmax",vmax, PS_TYPE_##type##_NAME, \
                        (psF64)PS_MIN_##type,(psF64)PS_MAX_##type); \
            } \
            for (psU32 row = 0;row<numRows;row++) { \
                ps##type* inputRow = input->data.type[row]; \
                for (psU32 col = 0; col < numCols; col++) { \
                    if ((psF64)inputRow[col] < min) { \
                        inputRow[col] = (ps##type)vmin; \
                        numClipped++; \
                    } else if ((psF64)inputRow[col] > max) { \
                        inputRow[col] = (ps##type)vmax; \
                        numClipped++; \
                    } \
                } \
            } \
        } \
        break;

        #define psImageClipCaseComplex(type,absfcn)\
    case PS_TYPE_##type: { \
            if (vmin < PS_MIN_##type || vmin > PS_MAX_##type) { \
                psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                        PS_ERRORTEXT_psImage_PARAMETER_OUTOF_TYPERANGE, \
                        "vmin",vmin, PS_TYPE_##type##_NAME, \
                        (psF64)PS_MIN_##type,(psF64)PS_MAX_##type); \
            } \
            if (vmax > PS_MAX_##type || vmax < PS_MIN_##type) { \
                psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                        PS_ERRORTEXT_psImage_PARAMETER_OUTOF_TYPERANGE, \
                        "vmax",vmax, PS_TYPE_##type##_NAME, \
                        (psF64)PS_MIN_##type,(psF64)PS_MAX_##type); \
            } \
            for (psU32 row = 0;row<numRows;row++) { \
                ps##type* inputRow = input->data.type[row]; \
                for (psU32 col = 0; col < numCols; col++) { \
                    if (absfcn(inputRow[col]) < min) { \
                        inputRow[col] = (ps##type)vmin; \
                        numClipped++; \
                    } else if (absfcn(inputRow[col]) > max) { \
                        inputRow[col] = (ps##type)vmax; \
                        numClipped++; \
                    } \
                } \
            } \
        } \
        break;

        psImageClipCase(S8)
        psImageClipCase(S16)
        psImageClipCase(S32)            // Not a requirement
        psImageClipCase(S64)            // Not a requirement
        psImageClipCase(U8)
        psImageClipCase(U16)
        psImageClipCase(U32)            // Not a requirement
        psImageClipCase(U64)            // Not a requirement
        psImageClipCase(F32)
        psImageClipCase(F64)
        psImageClipCaseComplex(C32, cabsf)
        psImageClipCaseComplex(C64, cabs)

    default: {
            char* typeStr;
            PS_TYPE_NAME(typeStr,input->type.type);
            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                    PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
                    typeStr);
        }
    }

    return numClipped;
}

psS32 psImageClipNaN(psImage* input,
                     psF64 value)
{
    psS32 numClipped = 0;
    psU32 numRows;
    psU32 numCols;

    if (input == NULL) {
        return 0;
    }
    numRows = input->numRows;
    numCols = input->numCols;

    switch (input->type.type) {

        #define psImageClipNaNCase(type) \
    case PS_TYPE_##type: \
        for (psU32 row = 0;row<numRows;row++) { \
            ps##type* inputRow = input->data.type[row]; \
            for (psU32 col = 0; col < numCols; col++) { \
                if (! isfinite(inputRow[col])) { \
                    inputRow[col] = (ps##type)value; \
                    numClipped++; \
                } \
            } \
        } \
        break;

        psImageClipNaNCase(F32)
        psImageClipNaNCase(F64)
        psImageClipNaNCase(C32)
        psImageClipNaNCase(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);
        }
    }

    return numClipped;
}

psS32 psImageOverlaySection(psImage* image,
                            const psImage* overlay,
                            psS32 col0,
                            psS32 row0,
                            const char *op)
{
    psU32 imageNumRows;
    psU32 imageNumCols;
    psU32 overlayNumRows;
    psU32 overlayNumCols;
    psU32 imageRowLimit;
    psU32 imageColLimit;
    psElemType type;
    psU32 pixelsOverlaid = 0;

    if (image == NULL || overlay == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        return pixelsOverlaid;
    }

    if (op == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImageManip_OPERATION_NULL);
        return pixelsOverlaid;
    }

    type = image->type.type;

    if (type != overlay->type.type) {
        char* typeStr;
        char* typeStrOverlay;
        PS_TYPE_NAME(typeStr,type);
        PS_TYPE_NAME(typeStrOverlay,overlay->type.type);
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psImageManip_OVERLAY_TYPE_MISMATCH,
                typeStrOverlay, typeStr);
        return pixelsOverlaid;
    }

    imageNumRows = image->numRows;
    imageNumCols = image->numCols;
    overlayNumRows = overlay->numRows;
    overlayNumCols = overlay->numCols;
    imageRowLimit = row0 + overlayNumRows;
    imageColLimit = col0 + overlayNumCols;

    /* check to see if overlay is within the input image */
    if ( row0 < 0 ||
            col0 < 0 ||
            imageRowLimit > imageNumRows ||
            imageColLimit > imageNumCols) {

        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_SUBSET_RANGE_INVALID,
                col0, imageColLimit, row0, imageRowLimit,
                imageNumCols, imageNumRows);
        return pixelsOverlaid;
    }


    #define psImageOverlayLoop(DATATYPE,OP) { \
        for (int row=row0;row<imageRowLimit;row++) { \
            ps##DATATYPE* imageRow = image->data.DATATYPE[row]; \
            ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-row0]; \
            for (int col=col0;col<imageColLimit;col++) { \
                imageRow[col] OP overlayRow[col-col0]; \
            } \
        } \
        pixelsOverlaid += (imageRowLimit - row0) * (imageColLimit - col0); \
    }

    #define psImageOverlayCase(DATATYPE) \
case PS_TYPE_##DATATYPE: \
    switch (*op) { \
    case '+': \
        psImageOverlayLoop(DATATYPE,+=); \
        break; \
    case '-': \
        psImageOverlayLoop(DATATYPE,-=); \
        break; \
    case '*': \
        psImageOverlayLoop(DATATYPE,*=); \
        break; \
    case '/': \
        psImageOverlayLoop(DATATYPE,/=); \
        break; \
    case '=': \
        psImageOverlayLoop(DATATYPE,=); \
        break; \
    default: \
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                PS_ERRORTEXT_psImageManip_OVERLAY_OPERATOR_INVALID, \
                op); \
        return pixelsOverlaid; \
    } \
    break;

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

    default: {
            char* typeStr;
            PS_TYPE_NAME(typeStr,type);
            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                    PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
                    typeStr);
            return pixelsOverlaid;
        }
    }

    return pixelsOverlaid;
}

psS32 psImageClipComplexRegion(psImage* input,
                               psC64 min,
                               psC64 vmin,
                               psC64 max,
                               psC64 vmax)
{
    psS32 numClipped = 0;
    psU32 numRows;
    psU32 numCols;
    psF64 realMin = creal(min);
    psF64 imagMin = cimag(min);
    psF64 realMax = creal(max);
    psF64 imagMax = cimag(max);

    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psImage_IMAGE_NULL);
        return 0;
    }

    if (realMax < realMin) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImageManip_MAXMIN_REAL,
                (double)realMin, (double)realMax);
        return 0;
    }
    if (imagMax < imagMin) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImageManip_MAXMIN_IMAG,
                (double)imagMin, (double)imagMax);
        return 0;
    }

    numRows = input->numRows;
    numCols = input->numCols;

    #define psImageClipComplexRegionCase(type,realfcn,imagfcn) \
case PS_TYPE_##type: { \
        if (realfcn(vmin) < PS_MIN_##type || imagfcn(vmin) < PS_MIN_##type || \
                realfcn(vmin) > PS_MAX_##type || imagfcn(vmin) > PS_MAX_##type ) { \
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                    PS_ERRORTEXT_psImageManip_CLIP_VALUE_INVALID, \
                    "vmin", creal(vmin), cimag(vmin), \
                    PS_TYPE_##type##_NAME, PS_MIN_##type, PS_MAX_##type); \
            break; \
        } \
        if (realfcn(vmax) > PS_MAX_##type || imagfcn(vmax) > PS_MAX_##type || \
                realfcn(vmax) < PS_MIN_##type || imagfcn(vmax) < PS_MIN_##type ) { \
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                    PS_ERRORTEXT_psImageManip_CLIP_VALUE_INVALID, \
                    "vmax", creal(vmax), cimag(vmax), \
                    PS_TYPE_##type##_NAME, PS_MIN_##type, PS_MAX_##type); \
            break; \
        } \
        for (psU32 row = 0;row<numRows;row++) { \
            ps##type* inputRow = input->data.type[row]; \
            for (psU32 col = 0; col < numCols; col++) { \
                if ( (realfcn(inputRow[col]) > realMax) || (imagfcn(inputRow[col]) > imagMax) ) { \
                    inputRow[col] = (ps##type)vmax; \
                    numClipped++; \
                } else if ( (realfcn(inputRow[col]) < realMin) || (imagfcn(inputRow[col]) < imagMin) ){ \
                    inputRow[col] = (ps##type)vmin; \
                    numClipped++; \
                } \
            } \
        } \
    } \
    break;

    switch (input->type.type) {

        psImageClipComplexRegionCase(C32, crealf, cimagf)
        psImageClipComplexRegionCase(C64, creal, cimag)

    default: {
            char* typeStr;
            PS_TYPE_NAME(typeStr,input->type.type);
            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                    PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
                    typeStr);
        }
    }

    return numClipped;
}

psImage* psImageRebin(psImage* out,
                      const psImage* in,
                      const psImage* restrict mask,
                      psMaskType maskVal,
                      psU32 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);
        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,
                         psS32 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_FLAT || mode >= PS_INTERPOLATE_NUM_MODES) {
        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* in,
                     psS32 dx,
                     psS32 dy)
{
    psS32 outRows;
    psS32 outCols;
    psS32 elementSize;

    if (in == 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 = in->numRows;
    outCols = in->numCols;
    elementSize = PSELEMTYPE_SIZEOF(in->type.type);
    out = psImageRecycle(out, outCols, outRows, in->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 = in->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* in,
                       float angle,
                       psC64 unexposedValue,
                       psImageInterpolateMode mode)
{
    if (in == 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*PS_PI) * floor(angle / (2.0*PS_PI)));

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

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

        #define PSIMAGE_ROTATE_LEFT_90(TYPE) \
    case PS_TYPE_##TYPE: { \
            ps##TYPE** inData = in->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 - PS_PI) < FLT_EPSILON) {
        // perform 1/2 rotate
        psS32 numRows = in->numRows;
        psS32 lastRow = numRows - 1;
        psS32 numCols = in->numCols;
        psS32 lastCol = numCols - 1;
        psElemType type = in->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 = in->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 - (PS_PI+PS_PI_2)) < FLT_EPSILON) {
        // perform 1/4 rotate clockwise
        psS32 numRows = in->numCols;
        psS32 lastRow = numRows - 1;
        psS32 numCols = in->numRows;
        psElemType type = in->type.type;

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

        #define PSIMAGE_ROTATE_RIGHT_90(TYPE) \
    case PS_TYPE_##TYPE: { \
            ps##TYPE** inData = in->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, in, in->type.type);
    } else {
        psElemType type = in->type.type;
        psS32 numRows = in->numRows;
        psS32 numCols = in->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(unexposedValue) < PS_MIN_##TYPE || \
                    creal(unexposedValue) > PS_MAX_##TYPE || \
                    cimag(unexposedValue) < PS_MIN_##TYPE || \
                    cimag(unexposedValue) > PS_MAX_##TYPE) { \
                psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                        PS_ERRORTEXT_psImageManip_CLIP_VALUE_INVALID, \
                        "unexposedValue", \
                        creal(unexposedValue),cimag(unexposedValue), \
                        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(in,inX,inY,NULL,0,unexposedValue); \
                    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);
        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* in,
                      float dx,
                      float dy,
                      psC64 unexposedValue,
                      psImageInterpolateMode mode)
{
    psS32 outRows;
    psS32 outCols;
    psS32 elementSize;
    psElemType type;

    if (in == 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 = in->numRows;
    outCols = in->numCols;
    type = in->type.type;
    elementSize = PSELEMTYPE_SIZEOF(type);
    out = psImageRecycle(out, outCols, outRows, type);

    #define PSIMAGE_SHIFT_CASE(MODE,TYPE) \
case PS_TYPE_##TYPE: \
    if (creal(unexposedValue) < PS_MIN_##TYPE || \
            creal(unexposedValue) > PS_MAX_##TYPE || \
            cimag(unexposedValue) < PS_MIN_##TYPE || \
            cimag(unexposedValue) > PS_MAX_##TYPE) { \
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
                PS_ERRORTEXT_psImageManip_CLIP_VALUE_INVALID, \
                "unexposedValue", \
                creal(unexposedValue),cimag(unexposedValue), \
                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( \
                          in,dx+(float)col,y,NULL,0,unexposedValue); \
        } \
    } \
    break;

    #define PSIMAGE_SHIFT_ARBITRARY_CASE(MODE) \
case PS_INTERPOLATE_##MODE: \
    switch (in->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);
    default:
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED,
                mode);
        psFree(out);
        out = NULL;
    }

    return out;
}
