/** @file  psImagePixelManip.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.11 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-09-01 01:05:41 $
 *
 *  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 "psImagePixelManip.h"

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

int psImageClip(psImage* input,
                double min,
                double vmin,
                double max,
                double 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;
}

int psImageClipNaN(psImage* input,
                   float 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;
}

int psImageOverlaySection(psImage* image,
                          const psImage* overlay,
                          int x0,
                          int y0,
                          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 = y0 + overlayNumRows;
    imageColLimit = x0 + overlayNumCols;

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

        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psImage_SUBSET_RANGE_INVALID,
                x0, imageColLimit, y0, imageRowLimit,
                imageNumCols, imageNumRows);
        return pixelsOverlaid;
    }


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

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

int psImageClipComplexRegion(psImage* input,
                             complex min,
                             complex vmin,
                             complex max,
                             complex 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;
}

// mask the area contained by the region
// the region is defined wrt the parent image
void psImageMaskRegion(psImage *image,
                       psRegion *region,
                       const char *op,
                       psMaskType maskValue)
{

    for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
            if (ix + image->col0 <  region->x0)
                continue;
            if (ix + image->col0 >= region->x1)
                continue;
            if (iy + image->row0 <  region->y0)
                continue;
            if (iy + image->row0 >= region->y1)
                continue;
            if ( !strncmp(op, "&", 2) || !strncmp(op, "AND", 5) ) {
                image->data.U8[iy][ix] &= maskValue;
            } else if ( !strncmp(op, "|", 2) || !strncmp(op, "OR", 5) ) {
                image->data.U8[iy][ix] |= maskValue;
            } else if ( !strncmp(op, "=", 2) || !strncmp(op, "EQUAL", 5) ) {
                image->data.U8[iy][ix] = maskValue;
            } else if ( !strncmp(op, "^", 2) || !strncmp(op, "XOR", 5) ) {
                image->data.U8[iy][ix] ^= maskValue;
            }
        }
    }
}

// mask the area not contained by the region
// the region is defined wrt the parent image
void psImageKeepRegion(psImage *image,
                       psRegion *region,
                       const char *op,
                       psMaskType maskValue)
{

    for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
            if (ix + image->col0 <  region->x0)
                goto maskit;
            if (ix + image->col0 >= region->x1)
                goto maskit;
            if (iy + image->row0 <  region->y0)
                goto maskit;
            if (iy + image->row0 >= region->y1)
                goto maskit;
            continue;
maskit:
            if ( !strncmp(op, "&", 2) || !strncmp(op, "AND", 5) ) {
                image->data.U8[iy][ix] &= maskValue;
            } else if ( !strncmp(op, "|", 2) || !strncmp(op, "OR", 5) ) {
                image->data.U8[iy][ix] |= maskValue;
            } else if ( !strncmp(op, "=", 2) || !strncmp(op, "EQUAL", 5) ) {
                image->data.U8[iy][ix] = maskValue;
            } else if ( !strncmp(op, "^", 2) || !strncmp(op, "XOR", 5) ) {
                image->data.U8[iy][ix] ^= maskValue;
            }
        }
    }
}

// mask the area contained by the region
// the region is defined wrt the parent image
void psImageMaskCircle(psImage *image,
                       double x,
                       double y,
                       double radius,
                       const char *op,
                       psMaskType maskValue)
{

    double dx, dy, r2, R2;

    R2 = PS_SQR(radius);

    for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
            dx = ix + image->col0 - x;
            dy = iy + image->row0 - y;
            r2 = PS_SQR(dx) + PS_SQR(dy);
            if (r2 > R2)
                continue;
            if ( !strncmp(op, "&", 2) || !strncmp(op, "AND", 5) ) {
                image->data.U8[iy][ix] &= maskValue;
            } else if ( !strncmp(op, "|", 2) || !strncmp(op, "OR", 5) ) {
                image->data.U8[iy][ix] |= maskValue;
            } else if ( !strncmp(op, "=", 2) || !strncmp(op, "EQUAL", 5) ) {
                image->data.U8[iy][ix] = maskValue;
            } else if ( !strncmp(op, "^", 2) || !strncmp(op, "XOR", 5) ) {
                image->data.U8[iy][ix] ^= maskValue;
            }
        }
    }
}

// mask the area contained by the region
// the region is defined wrt the parent image
void psImageKeepCircle(psImage *image,
                       double x,
                       double y,
                       double radius,
                       const char *op,
                       psMaskType maskValue)
{

    double dx, dy, r2, R2;

    R2 = PS_SQR(radius);

    for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
            dx = ix + image->col0 - x;
            dy = iy + image->row0 - y;
            r2 = PS_SQR(dx) + PS_SQR(dy);
            if (r2 < R2)
                continue;
            if ( !strncmp(op, "&", 2) || !strncmp(op, "AND", 5) ) {
                image->data.U8[iy][ix] &= maskValue;
            } else if ( !strncmp(op, "|", 2) || !strncmp(op, "OR", 5) ) {
                image->data.U8[iy][ix] |= maskValue;
            } else if ( !strncmp(op, "=", 2) || !strncmp(op, "EQUAL", 5) ) {
                image->data.U8[iy][ix] = maskValue;
            } else if ( !strncmp(op, "^", 2) || !strncmp(op, "XOR", 5) ) {
                image->data.U8[iy][ix] ^= maskValue;
            }
        }
    }
}

