/*  @file  psImageConvolve.c
 *
 *  @brief Contains FFT transform related functions for psImage.
 *
 *  @author Robert DeSonia, MHPCC
 *
 *  @version $Revision: 1.42 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2007-01-19 04:30:33 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
 */

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#include <string.h>
#include <math.h>
#include "psImageConvolve.h"
#include "psImageFFT.h"
#include "psImageStructManip.h"
#include "psBinaryOp.h"
#include "psMemory.h"
#include "psLogMsg.h"
#include "psError.h"
#include "psAssert.h"



#define FOURIER_PADDING 32 /* padding amount in every side of the image for fourier convolution */

static void freeKernel(psKernel* ptr);

psKernel* psKernelAlloc(int xMin, int xMax, int yMin, int yMax)
{
    psKernel* result;
    psS32 numRows;
    psS32 numCols;

    // following is explicitly spelled out in the SDRS as a requirement
    if (yMin > yMax) {
        psLogMsg(__func__, PS_LOG_WARN,
                 "Specified yMin, %d, was greater than yMax, %d.  Values swapped.",
                 yMin, yMax);

        psS32 temp = yMin;
        yMin = yMax;
        yMax = temp;
    }

    // following is explicitly spelled out in the SDRS as a requirement
    if (xMin > xMax) {
        psLogMsg(__func__, PS_LOG_WARN,
                 "Specified xMin, %d, was greater than xMax, %d.  Values swapped.",
                 xMin, xMax);

        psS32 temp = xMin;
        xMin = xMax;
        xMax = temp;
    }

    numRows = yMax - yMin + 1;
    numCols = xMax - xMin + 1;

    result = psAlloc(sizeof(psKernel));
    result->xMin = xMin;
    result->xMax = xMax;
    result->yMin = yMin;
    result->yMax = yMax;
    result->image = psImageAlloc(numCols,numRows,PS_TYPE_KERNEL);
    memset((result->image->p_rawDataBuffer),0,numCols*numRows*PSELEMTYPE_SIZEOF(PS_TYPE_KERNEL));
    result->p_kernelRows = psAlloc(sizeof(float*)*numRows);

    float** kernelRows = result->p_kernelRows;
    float** imageRows = result->image->data.PS_TYPE_KERNEL_DATA;
    for (psS32 i = 0; i < numRows; i++) {
        kernelRows[i] = imageRows[i] - xMin;
    }
    result->kernel = kernelRows - yMin;

    psMemSetDeallocator(result,(psFreeFunc)freeKernel);

    return result;
}

void freeKernel(psKernel* ptr)
{
    if (ptr != NULL) {
        psFree(ptr->image);
        psFree(ptr->p_kernelRows);
    }
}


bool psMemCheckKernel(psPtr ptr)
{
    PS_ASSERT_PTR(ptr, false);
    return ( psMemGetDeallocator(ptr) == (psFreeFunc)freeKernel );
}


psKernel* psKernelGenerate(const psVector* tShifts,
                           const psVector* xShifts,
                           const psVector* yShifts,
                           bool relative)
{
    psS32 lastX;
    psS32 lastY;
    psS32 lastT;
    psS32 x;
    psS32 y;
    psS32 t;
    psS32 xMin = 0;
    psS32 xMax = 0;
    psS32 yMin = 0;
    psS32 yMax = 0;
    psS32 length = 0;
    float normalizeTime = 1.0;  // fraction of total time for each shift clock
    psKernel* result = NULL;
    float** kernel = NULL;

    // got non-NULL vectors?
    if (tShifts == NULL || xShifts == NULL || yShifts == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                _("Specified shift vectors can not be NULL."));
        return NULL;
    }

    // types match?
    if (xShifts->type.type != yShifts->type.type ||
            tShifts->type.type != xShifts->type.type) {
        char* typeXStr;
        char* typeYStr;
        char* typeTStr;
        PS_TYPE_NAME(typeXStr,xShifts->type.type);
        PS_TYPE_NAME(typeYStr,yShifts->type.type);
        PS_TYPE_NAME(typeTStr,tShifts->type.type);
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                _("Input t-, x-, and y-shift vector types (%s/%s/%s) must match."),
                typeTStr, typeXStr, typeYStr);
        return NULL;
    }

    // sizes match?
    length = xShifts->n;
    if (length != yShifts->n ||
            length != tShifts->n) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                "Shift vectors can not be of different sizes.");
        return NULL;
    }

    // if no shifts, the kernel is just a 1 at 0,0
    if (length < 1) {
        result = psKernelAlloc(0,0,0,0);
        result->kernel[0][0] = 1;
        return result;
    }

    #define KERNEL_GENERATE_CASE(TYPE) \
case PS_TYPE_##TYPE: { \
        ps##TYPE *tShiftData = tShifts->data.TYPE; \
        ps##TYPE *xShiftData = xShifts->data.TYPE; \
        ps##TYPE *yShiftData = yShifts->data.TYPE; \
        lastX = xShiftData[length-1]; \
        lastY = yShiftData[length-1]; \
        lastT = tShiftData[length-1]; \
        \
        for (int lcv = 0; lcv < length; lcv++) { \
            x = lastX - xShiftData[lcv]; \
            y = lastY - yShiftData[lcv]; \
            \
            if (x < xMin) { \
                xMin = x; \
            } else if (x > xMax) { \
                xMax = x; \
            } \
            if (y < yMin) { \
                yMin = y; \
            } else if (y > yMax) { \
                yMax = y; \
            } \
        } \
        \
        normalizeTime = 1.0 / (float)(tShiftData[length-1]); \
        result = psKernelAlloc(xMin,xMax,yMin,yMax); \
        kernel = result->kernel; \
        \
        psS32 prevT = 0; \
        for (int i = 0; i < length; i++) { \
            t = tShiftData[i] - prevT; \
            x = lastX - xShiftData[i]; \
            y = lastY - yShiftData[i]; \
            \
            kernel[y][x] += (float)t / (float)lastT; \
            prevT = tShiftData[i]; \
        } \
        break; \
    }

    #define RELATIVE_KERNEL_GENERATE_CASE(TYPE) \
case PS_TYPE_##TYPE: { \
        ps##TYPE *tShiftData = tShifts->data.TYPE; \
        ps##TYPE *xShiftData = xShifts->data.TYPE; \
        ps##TYPE *yShiftData = yShifts->data.TYPE; \
        \
        x = 0; \
        y = 0; \
        t = 0; \
        \
        for (int lcv = length-1; lcv >= 0; lcv--) { \
            t += tShiftData[lcv]; \
            \
            if (x < xMin) { \
                xMin = x; \
            } else if (x > xMax) { \
                xMax = x; \
            } \
            if (y < yMin) { \
                yMin = y; \
            } else if (y > yMax) { \
                yMax = y; \
            } \
            x -= xShiftData[lcv]; \
            y -= yShiftData[lcv]; \
            \
        } \
        result = psKernelAlloc(xMin,xMax,yMin,yMax); \
        kernel = result->kernel; \
        \
        normalizeTime = 1.0 / (float)t; \
        x = 0; \
        y = 0; \
        for (psS32 i = length-1; i >= 0; i--) { \
            kernel[y][x] += (float)(tShiftData[i]) * normalizeTime; \
            x -= xShiftData[i]; \
            y -= yShiftData[i]; \
            \
        } \
        break; \
    }

    if (relative) {
        switch (xShifts->type.type) {
            RELATIVE_KERNEL_GENERATE_CASE(U8);
            RELATIVE_KERNEL_GENERATE_CASE(U16);
            RELATIVE_KERNEL_GENERATE_CASE(U32);
            RELATIVE_KERNEL_GENERATE_CASE(U64);
            RELATIVE_KERNEL_GENERATE_CASE(S8);
            RELATIVE_KERNEL_GENERATE_CASE(S16);
            RELATIVE_KERNEL_GENERATE_CASE(S32);
            RELATIVE_KERNEL_GENERATE_CASE(S64);
            RELATIVE_KERNEL_GENERATE_CASE(F32);
            RELATIVE_KERNEL_GENERATE_CASE(F64);
            RELATIVE_KERNEL_GENERATE_CASE(C32);
            RELATIVE_KERNEL_GENERATE_CASE(C64);

        default: {
                char* typeStr;
                PS_TYPE_NAME(typeStr,xShifts->type.type);
                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                        _("Specified psImage type, %s, is not supported."),
                        typeStr);
            }
        }
    } else {
        switch (xShifts->type.type) {
            KERNEL_GENERATE_CASE(U8);
            KERNEL_GENERATE_CASE(U16);
            KERNEL_GENERATE_CASE(U32);
            KERNEL_GENERATE_CASE(U64);
            KERNEL_GENERATE_CASE(S8);
            KERNEL_GENERATE_CASE(S16);
            KERNEL_GENERATE_CASE(S32);
            KERNEL_GENERATE_CASE(S64);
            KERNEL_GENERATE_CASE(F32);
            KERNEL_GENERATE_CASE(F64);
            KERNEL_GENERATE_CASE(C32);
            KERNEL_GENERATE_CASE(C64);

        default: {
                char* typeStr;
                PS_TYPE_NAME(typeStr,xShifts->type.type);
                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                        _("Specified psImage type, %s, is not supported."),
                        typeStr);
            }
        }
    }

    return result;
}

psImage* psImageConvolve(psImage* out,
                         const psImage* in,
                         const psKernel* kernel,
                         bool direct)
{
    if (in == NULL) {
        psFree(out);
        return NULL;
    }

    if (kernel == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                _("Specified psKernel can not be NULL."));
        psFree(out);
        return NULL;
    }
    psS32 xMin = kernel->xMin;
    psS32 xMax = kernel->xMax;
    psS32 yMin = kernel->yMin;
    psS32 yMax = kernel->yMax;
    float** kData = kernel->kernel;

    // make the output image to the proper size and type
    psS32 numRows = in->numRows;
    psS32 numCols = in->numCols;



    if (direct) {
        // spatial convolution

        #define SPATIAL_CONVOLVE_CASE(TYPE) \
    case PS_TYPE_##TYPE: { \
            ps##TYPE** inData = in->data.TYPE; \
            out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \
            for (psS32 row=0;row<numRows;row++) { \
                ps##TYPE* outRow = out->data.TYPE[row]; \
                for (psS32 col=0;col<numCols;col++) { \
                    ps##TYPE pixel = 0.0; \
                    for (psS32 kRow = yMin; kRow < yMax; kRow++) { \
                        if (row-kRow >= 0 && row-kRow < numRows) { \
                            for (psS32 kCol = xMin; kCol < xMax; kCol++) { \
                                if (col-kCol >= 0 && col-kCol < numCols) { \
                                    pixel += kData[kRow][kCol] * inData[row-kRow][col-kCol]; \
                                } \
                            } \
                        } \
                    } \
                    outRow[col] = pixel; \
                } \
            } \
        } \
        break;

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

        default: {
                char* typeStr;
                PS_TYPE_NAME(typeStr,in->type.type);
                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                        _("Specified psImage type, %s, is not supported."),
                        typeStr);
                psFree(out);
                return NULL;

            }
        }


    } else {
        // fourier convolution
        psS32 paddedCols = numCols+2*FOURIER_PADDING;
        psS32 paddedRows = numRows+2*FOURIER_PADDING;

        // check to see if kernel is smaller, otherwise padding it up will fail.
        psS32 kRows = kernel->image->numRows;
        psS32 kCols = kernel->image->numCols;
        if (kRows >= numRows || kCols >= numCols) {
            psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                    _("Specified psKernel size, %dx%d, can not be larger than input psImage size, %dx%d."),
                    kCols,kRows,
                    numCols, numRows);
            psFree(out);
            return NULL;
        }

        // pad the image
        psImage* paddedImage = psImageAlloc(paddedCols,paddedRows,in->type.type);
        psS32 elementSize = PSELEMTYPE_SIZEOF(in->type.type);

        // zero out padded area on top and bottom
        memset(paddedImage->data.U8[0],0,FOURIER_PADDING*paddedCols*elementSize);
        memset(paddedImage->data.U8[FOURIER_PADDING+numRows-1],0,FOURIER_PADDING*paddedCols*elementSize);

        // fill in the image-containing rows.
        psS32 sidePaddingSize = FOURIER_PADDING*elementSize;
        psS32 imageRowSize = numCols*elementSize;
        psU8* paddedData = paddedImage->data.U8[FOURIER_PADDING];
        for (psS32 row=0;row<numRows;row++) {
            // zero out padded area on left edge.
            memset(paddedData,0,sidePaddingSize);
            paddedData += sidePaddingSize;
            memcpy(paddedData,in->data.U8[row],imageRowSize);
            paddedData += imageRowSize;
            // zero out padded area on right edge.
            memset(paddedData,0,sidePaddingSize);
            paddedData += sidePaddingSize;
        }

        // pad the kernel to the same size of paddedImage
        psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL);
        memset(paddedKernel->data.U8[0],0,sizeof(float)*numCols*numRows); // zero-out image
        psS32 yMax = kernel->yMax;
        psS32 xMax = kernel->xMax;
        for (psS32 row = kernel->yMin; row <= yMax;row++) {
            psS32 padRow = row;
            if (padRow < 0) {
                padRow += paddedRows;
            }
            float* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow];
            float* kernelRow = kernel->kernel[row];
            for (psS32 col = kernel->xMin; col <= xMax; col++) {
                if (col < 0) {
                    padData[col+paddedCols] = kernelRow[col];
                } else {
                    padData[col] = kernelRow[col];
                }
            }
        }

        psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD);
        if (kernelFourier == NULL) {
            psError(PS_ERR_UNKNOWN, false,
                    _("Failed to perform a fourier transform of kernel."));
            psFree(out);
            return NULL;
        }

        psImage* inFourier = psImageFFT(NULL, paddedImage, PS_FFT_FORWARD);
        if (inFourier == NULL) {
            psError(PS_ERR_UNKNOWN, false,
                    _("Failed to perform a fourier transform of input image."));
            psFree(out);
            return NULL;
        }

        // convolution in fourier domain is just a pixel-wise multiplication
        for (int row = 0; row < paddedRows; row++) {
            psC32* inRow = inFourier->data.C32[row];
            psC32* kRow = kernelFourier->data.C32[row];
            for (int col = 0; col < paddedCols; col++) {
                inRow[col] *= kRow[col];
            }
        }

        psImage* complexOut = psImageFFT(NULL, inFourier,
                                         PS_FFT_REVERSE);
        if (complexOut == NULL) {
            psError(PS_ERR_UNKNOWN, false,
                    _("Failed to perform a fourier transform of input image."));
            psFree(out);
            return NULL;
        }

        // subset out the padded area now.
        psImage* complexOutSansPad = psImageSubset(complexOut,
                                     psRegionSet(FOURIER_PADDING, FOURIER_PADDING+numCols,FOURIER_PADDING,FOURIER_PADDING+numRows));

        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
        float factor = 1.0f/(float)paddedCols/(float)paddedRows;
        for (psS32 row = 0; row < numRows; row++) {
            psF32* outRow = out->data.F32[row];
            psC32* resultRow = complexOutSansPad->data.C32[row];
            for (psS32 col = 0; col < numCols; col++) {
                outRow[col] = crealf(resultRow[col])*factor;
            }
        }

        psFree(complexOut); // frees complexOutSansPad, as it is a child.
        psFree(kernelFourier);
        psFree(inFourier);
        psFree(paddedImage);
        psFree(paddedKernel);

    }

    return out;
}

bool psImageSmooth (psImage *image,
                    double  sigma,
                    double  Nsigma)
{
    PS_ASSERT_IMAGE_NON_NULL(image, NULL);

    // relevant terms
    int Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
    int Npixel = 2*Nrange + 1;          // Total number of pixels in convolution kernel
    int Nx = image->numCols;            // Number of columns
    int Ny = image->numRows;            // Number of rows

    #define IMAGESMOOTH_CASE(TYPE) \
case PS_TYPE_##TYPE: { \
        /* generate normalized gaussian */ \
        psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_##TYPE); \
        { \
            double sum = 0.0; \
            double factor = -0.5/(sigma*sigma); \
            for (int i = -Nrange; i < Nrange + 1; i++) { \
                gaussnorm->data.TYPE[i+Nrange] = exp(factor*i*i); \
                sum += gaussnorm->data.TYPE[i+Nrange]; \
            } \
            for (int i = -Nrange; i < Nrange + 1; i++) { \
                gaussnorm->data.TYPE[i+Nrange] /= sum; \
            } \
        } \
        ps##TYPE *gauss = &gaussnorm->data.TYPE[Nrange]; \
        \
        /* Smooth in X direction */ \
        { \
            psVector *calculation = psVectorAlloc(Nx, PS_TYPE_##TYPE); \
            for (int j = 0; j < Ny; j++) { \
                ps##TYPE *vi = image->data.TYPE[j]; \
                ps##TYPE *vo = calculation->data.TYPE; \
                int xMax = PS_MIN(Nrange, Nx); \
                int convRange = PS_MIN(Nrange + 1, Nx); \
                /* Smooth first Nrange pixels, with renorm */ \
                for (int i = 0; i < xMax; i++, vi++, vo++) { \
                    ps##TYPE *vr = vi - i; \
                    ps##TYPE *vg = gauss - i; \
                    double g = 0.0; \
                    double s = 0.0; \
                    for (int n = -i; n < convRange; n++, vr++, vg++) { \
                        s += *vg * *vr; \
                        g += *vg; \
                    } \
                    *vo = s / g; \
                } \
                /* If that's all the pixels we have, then we're done already */ \
                if (Nx > Nrange) { \
                    /* Smooth middle pixels */ \
                    for (int i = Nrange; i < Nx - Nrange; i++, vi++, vo++) { \
                        ps##TYPE *vr = vi - Nrange; \
                        ps##TYPE *vg = gauss - Nrange; \
                        double s = 0; \
                        for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) { \
                            s += *vg * *vr; \
                        } \
                        *vo = s; \
                    } \
                    /* Smooth last Nrange pixels, with renorm */ \
                    /* XXX does this miss the last column? */ \
                    for (int i = Nx - Nrange; i < Nx; i++, vi++, vo++) { \
                        ps##TYPE *vr = vi - Nrange; \
                        ps##TYPE *vg = gauss - Nrange; \
                        double g = 0.0; \
                        double s = 0.0; \
                        for (int n = -Nrange; n < Nx - i; n++, vr++, vg++) { \
                            s += *vg * *vr; \
                            g += *vg; \
                        } \
                        *vo = s / g; \
                    } \
                } \
                memcpy(image->data.TYPE[j], calculation->data.TYPE, Nx*sizeof(ps##TYPE)); \
            } \
            psFree(calculation); \
        } \
        \
        /* Smooth in Y direction */ \
        psArray *rows = psArrayAlloc(Nrange); \
        /* Smooth the first Nrange pixels, with renorm */ \
        int yMax = PS_MIN(Nrange, Ny); \
        int convRange = PS_MIN(Nrange + 1, Ny); \
        for (int j = 0; j < yMax; j++) { \
            psVector *calculation = psVectorAlloc(Nx, PS_TYPE_##TYPE); \
            /* Zero the output row */ \
            memset(calculation->data.TYPE, 0, Nx*sizeof(ps##TYPE)); \
            double sum = 0.0; \
            for (int n = -j; n < convRange; n++) { sum += gauss[n]; } \
            for (int n = -j; n < convRange; n++) { \
                ps##TYPE *vi = image->data.TYPE[j+n]; \
                ps##TYPE *vo = calculation->data.TYPE; \
                double g = gauss[n] / sum; \
                for (int i = 0; i < Nx; i++, vi++, vo++) { \
                    *vo += *vi * g; \
                } \
            } \
            /* Save output rows on temp array of rows */ \
            rows->data[j] = calculation; \
        } \
        if (Ny < Nrange) { \
            /* Need to save the first bit, then we're done */ \
            for (int j = 0; j < Ny; j++) {  \
                psVector *save = rows->data[j]; \
                memcpy(image->data.TYPE[j], save->data.TYPE, Nx*sizeof(ps##TYPE)); \
            } \
        } else { \
            /* Smooth middle pixels */ \
            psVector *calculation = psVectorAlloc(Nx, PS_TYPE_##TYPE); \
            for (int j = Nrange; j < Ny - Nrange; j++) { \
                memset(calculation->data.TYPE, 0, Nx*sizeof(ps##TYPE)); \
                for (int n = -Nrange; n < Nrange + 1; n++) { \
                    ps##TYPE *vi = image->data.TYPE[j+n]; \
                    ps##TYPE *vo = calculation->data.TYPE; \
                    double g = gauss[n]; \
                    for (int i = 0; i < Nx; i++, vi++, vo++) { \
                        *vo += *vi * g; \
                    } \
                } \
                /* Write the output row */ \
                int Nr = j % Nrange; \
                psVector *save = rows->data[Nr]; \
                memcpy(image->data.TYPE[j-Nrange], save->data.TYPE, Nx*sizeof(ps##TYPE)); \
                /* Juggle the pointers, so that next run we use the one we just wrote */ \
                rows->data[Nr] = calculation; \
                calculation = save; \
            } \
            /* Smooth last Nrange pixels, with renorm */ \
            for (int j = Ny - Nrange; j < Ny; j++) { \
                /* save the Nrange-offset output row, then zero */ \
                memset(calculation->data.TYPE, 0, Nx*sizeof(ps##TYPE)); \
                double sum = 0.0; \
                for (int n = -Nrange; n < Ny - j; n++) { sum += gauss[n]; } \
                for (int n = -Nrange; n < Ny - j; n++) { \
                    ps##TYPE *vi = image->data.TYPE[j+n]; \
                    ps##TYPE *vo = calculation->data.TYPE; \
                    double g = gauss[n] / sum; \
                    for (int i = 0; i < Nx; i++, vi++, vo++) { \
                        *vo += *vi * g; \
                    } \
                } \
                /* Write the output row */ \
                int Nr = j % Nrange; \
                psVector *save = rows->data[Nr]; \
                memcpy(image->data.TYPE[j-Nrange], save->data.TYPE, Nx*sizeof(ps##TYPE)); \
                /* Juggle the pointers, so that next run we use the one we just wrote */ \
                rows->data[Nr] = calculation; \
                calculation = save; \
            } \
            psFree(calculation); \
            /* Write the remaining rows */ \
            for (int j = Ny; j < Ny + Nrange; j++) {  \
                int Nr = j % Nrange; \
                psVector *save = rows->data[Nr]; \
                memcpy(image->data.TYPE[j-Nrange], save->data.TYPE, Nx*sizeof(ps##TYPE)); \
            } \
        } \
        psFree(rows); \
        psFree(gaussnorm); \
        break; \
    }

    switch (image->type.type) {
        IMAGESMOOTH_CASE(F32);
        IMAGESMOOTH_CASE(F64);
    default: {
            char* typeStr;
            PS_TYPE_NAME(typeStr,image->type.type);
            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                    _("Specified psImage type, %s, is not supported."),
                    typeStr);
            return false;
        }
    }
    return true;
}

bool psImageSmoothMaskF32 (psImage *image,
                           psImage *mask,
                           psMaskType maskVal,
                           double  sigma,
                           double  Nsigma)
{
    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);

    // relevant terms
    int Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
    int Npixel = 2*Nrange + 1;          // Total number of pixels in convolution kernel
    int Nx = image->numCols;            // Number of columns
    int Ny = image->numRows;            // Number of rows

    /* generate normalized gaussian */
    psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_F32);
    {
        double sum = 0.0;
        double factor = -0.5/(sigma*sigma);
        for (int i = -Nrange; i < Nrange + 1; i++) {
            gaussnorm->data.F32[i+Nrange] = exp(factor*i*i);
            sum += gaussnorm->data.F32[i+Nrange];
        }
        for (int i = -Nrange; i < Nrange + 1; i++) {
            gaussnorm->data.F32[i+Nrange] /= sum;
        }
    }
    psF32 *gauss = &gaussnorm->data.F32[Nrange];

    /** Smooth in X direction **/
    psVector *calculation = psVectorAlloc(Nx, PS_TYPE_F32);
    for (int j = 0; j < Ny; j++) {
        psU8  *vm = mask->data.U8[j];
        psF32 *vi = image->data.F32[j];
        psF32 *vo = calculation->data.F32;
        // loop over all pixels in the row
        for (int i = 0; i < Nx; i++, vi++, vo++, vm++) {
            int offset = PS_MIN (i, Nrange);
            psU8  *sm = vm - offset;
            psF32 *si = vi - offset;
            psF32 *sg = gauss - offset;
            double g = 0.0;
            double s = 0.0;
            // loop over all valid pixels in the smoothing kernel
            int xMin = PS_MAX (i - Nrange, 0);
            int xMax = PS_MIN (i + Nrange + 1, Nx);
            for (int n = xMin; n < xMax; n++, sm++, si++, sg++) {
                if (*sm & maskVal)
                    continue;
                s += *sg * *si;
                g += *sg;
            }
            *vo = (g > 0) ? s / g : 0.0;
        }
        memcpy(image->data.F32[j], calculation->data.F32, Nx*sizeof(psF32));
    }
    psFree(calculation);

    /** Smooth in Y direction  **/
    // allocate and save Nrange extra row vectors for storage. these will be
    // written over the output rows only after we are Nrange rows beyond them
    int Nsave = Nrange + 1;
    psArray *rows = psArrayAlloc(Nsave);
    for (int i = 0; i < Nsave; i++) {
        rows->data[i] = psVectorAlloc(Nx, PS_TYPE_F32);
    }

    psVector *outsum = psVectorAlloc(Nx, PS_TYPE_F32);
    for (int j = 0; j < Ny; j++) {
        psVector *output = rows->data[j % Nsave];
        memset (output->data.F32, 0, Nx*sizeof(psF32));
        memset (outsum->data.F32, 0, Nx*sizeof(psF32));
        int yMin = PS_MAX (j - Nrange, 0);
        int yMax = PS_MIN (j + Nrange + 1, Ny);
        for (int n = yMin; n < yMax; n++) {
            psU8  *vm = mask->data.U8[n];
            psF32 *vi = image->data.F32[n];
            psF32 *vo = output->data.F32;
            psF32 *vs = outsum->data.F32;
            double g = gauss[n - j];
            for (int i = 0; i < Nx; i++, vi++, vo++, vm++, vs++) {
                if (*vm & maskVal)
                    continue;
                *vo += *vi * g;
                *vs += g;
            }
        }
        // renormalize the row
        psF32 *vo = output->data.F32;
        psF32 *vs = outsum->data.F32;
        for (int i = 0; i < Nx; i++, vo++, vs++) {
            *vo = (*vs > 0) ? *vo / *vs : 0.0;
        }

        // Write the output row
        if (j - Nrange >= 0) {
            int Nout = (j - Nrange) % Nsave;
            psVector *save = rows->data[Nout];
            memcpy(image->data.F32[j-Nrange], save->data.F32, Nx*sizeof(psF32));
        }
    }

    // Write the remaining output rows
    for (int j = PS_MAX(0, Ny - Nrange); j < Ny; j++) {
        psVector *save = rows->data[j % Nsave];
        memcpy(image->data.F32[j], save->data.F32, Nx*sizeof(psF32));
    }
    psFree(rows);
    psFree(outsum);
    psFree(gaussnorm);
    return true;
}
