/// @file  psImageConvolve.c
///
/// @brief Contains FFT transform related functions for psImage.
///
/// @author Robert DeSonia, MHPCC
/// @author Paul Price, IfA
/// @author Eugene Magnier, IfA
///
/// @version $Revision: 1.66 $ $Name: not supported by cvs2svn $
/// @date $Date: 2008-04-07 20:30:38 $
///
/// Copyright 2004-2007 Institute for Astronomy, University of Hawaii
///

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

#include <string.h>
#include <math.h>
#include "psAbort.h"
#include "psMemory.h"
#include "psLogMsg.h"
#include "psError.h"
#include "psAssert.h"
#include "psScalar.h"
#include "psBinaryOp.h"
#include "psImageFFT.h"
#include "psImageStructManip.h"
#include "psImagePixelManip.h"
#include "psTrace.h"

#include "psImageConvolve.h"

static void kernelFree(psKernel *kernel)
{
    if (kernel) {
        psFree(kernel->image);
        psFree(kernel->p_kernelRows);
    }
    return;
}

// Set up indirections, so we can refer to kernel->kernel[-1][-3] for the (-1,-1) element, instead of
// kernel->image[kernel->yMax - kernel->yMin + 1][kernel->yMax - kernel->yMin + 3] (yuk!).
static void kernelRedirects(psKernel *kernel, // Kernel for which to set up the redirects
                            int numRows // Number of rows
    )
{
    int xMin = kernel->xMin, yMin = kernel->yMin; // Minimum values

    kernel->p_kernelRows = psAlloc(sizeof(float*)*numRows);
    for (int i = 0; i < numRows; i++) {
        kernel->p_kernelRows[i] = kernel->image->data.PS_TYPE_KERNEL_DATA[i] - xMin;
    }
    kernel->kernel = kernel->p_kernelRows - yMin;
    return;
}

psKernel *psKernelAlloc(int xMin, int xMax, int yMin, int yMax)
{
    // Check the inputs to make sure max > min; if not, switch.
    // following is explicitly spelled out in the SDRS as a requirement
    if (yMin > yMax) {
        psWarning("Specified yMin, %d, was greater than yMax, %d.  Values swapped.",
                 yMin, yMax);

        int temp = yMin;
        yMin = yMax;
        yMax = temp;
    }
    if (xMin > xMax) {
        psWarning("Specified xMin, %d, was greater than xMax, %d.  Values swapped.",
                 xMin, xMax);

        int temp = xMin;
        xMin = xMax;
        xMax = temp;
    }

    int numRows = yMax - yMin + 1;      // Number of rows for kernel image
    int numCols = xMax - xMin + 1;      // Number of columns for kernel image

    psKernel *kernel = psAlloc(sizeof(psKernel)); // The kernel, to be returned
    psMemSetDeallocator(kernel,(psFreeFunc)kernelFree);

    kernel->xMin = xMin;
    kernel->xMax = xMax;
    kernel->yMin = yMin;
    kernel->yMax = yMax;
    kernel->image = psImageAlloc(numCols, numRows, PS_TYPE_KERNEL);
    psImageInit(kernel->image, 0.0);

    kernelRedirects(kernel, numRows);

    return kernel;
}

psKernel *psKernelAllocFromImage(psImage *image, int x0, int y0)
{
    psKernel *kernel = psAlloc(sizeof(psKernel)); // The kernel, to be returned
    psMemSetDeallocator(kernel,(psFreeFunc)kernelFree);

    int numCols = image->numCols, numRows = image->numRows; // Size of image

    kernel->xMin = - x0;
    kernel->xMax = numCols - 1 - x0;
    kernel->yMin = - y0;
    kernel->yMax = numRows - 1 - y0;
    kernel->image = psMemIncrRefCounter(image);

    kernelRedirects(kernel, numRows);

    return kernel;
}

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


psKernel *psKernelGenerate(const psVector *tShifts,
                           const psVector *xShifts,
                           const psVector *yShifts,
                           float totalTime,
                           bool xyRelative)
{
    PS_ASSERT_VECTOR_NON_NULL(tShifts, NULL);
    PS_ASSERT_VECTOR_NON_NULL(xShifts, NULL);
    PS_ASSERT_VECTOR_NON_NULL(yShifts, NULL);
    PS_ASSERT_VECTORS_SIZE_EQUAL(tShifts, xShifts, NULL);
    PS_ASSERT_VECTORS_SIZE_EQUAL(tShifts, yShifts, NULL);
    PS_ASSERT_VECTOR_TYPE(tShifts, PS_TYPE_F32, NULL);
    PS_ASSERT_VECTOR_TYPE(xShifts, PS_TYPE_S32, NULL);
    PS_ASSERT_VECTOR_TYPE(yShifts, PS_TYPE_S32, NULL);

    if (isnan(totalTime)) {
        // It's more expensive to check for NAN than 0.0
        totalTime = 0.0;
    }

    // If there are no shifts, the kernel is just a 1 at 0,0
    long num = tShifts->n;              // Number of shifts
    if (num == 0) {
        psKernel *kernel = psKernelAlloc(0,0,0,0);
        kernel->kernel[0][0] = 1;
        return kernel;
    }

    // Get dimensions and scaling
    int xMin, xMax, yMin, yMax;         // Range of values for kernel
    int xLast, yLast;                   // Last location, for relative shifts
    float tSum = tShifts->data.F32[0];   // Sum of the times
    xLast = xMin = xMax = xShifts->data.S32[0];
    yLast = yMin = yMax = yShifts->data.S32[0];
    int x0, y0;                         // Final location; everything is relative to this
    x0 = xShifts->data.S32[num - 1];
    y0 = yShifts->data.S32[num - 1];
    for (long i = 1; i < num; i++) {
        int x = xShifts->data.S32[i] - x0; // x position in kernel
        int y = yShifts->data.S32[i] - y0; // y position in kernel
        if (xyRelative) {
            x += xLast;
            y += yLast;
            xLast = x;
            yLast = y;
        }
        if (x < xMin) {
            xMin = x;
        }
        if (x > xMax) {
            xMax = x;
        }
        if (y < yMin) {
            yMin = y;
        }
        if (y > yMax) {
            yMax = y;
        }

        if (totalTime <= 0) {
            tSum += tShifts->data.F32[i];
        }
    }

    psTrace("psLib.imageops", 5, "Kernel range: %d:%d,%d:%d\n", xMin, xMax, yMin, yMax);

    if (totalTime > 0) {
        // Then the total time is simply the final value
        // NB: We assume the counter starts at zero!
        tSum = totalTime;
    }

    // One more pass through to set the kernel
    psKernel *kernel = psKernelAlloc(xMin, xMax, yMin, yMax); // The kernel
    xLast = xShifts->data.S32[0];
    yLast = yShifts->data.S32[0];
    float tLast = 0.0;                  // Last value for t
    for (int i = 0; i < num; i++) {
        int x = xShifts->data.S32[i] - x0; // x position in kernel
        int y = yShifts->data.S32[i] - y0; // y position in kernel
        if (xyRelative) {
            x += xLast;
            y += yLast;
            xLast = x;
            yLast = y;
        }
        float t = tShifts->data.F32[i];
        if (totalTime > 0) {
            t -= tLast;
            tLast = tShifts->data.F32[i];
        }

        kernel->kernel[y][x] += t;
    }

    // Normalise the kernel by the total time (kernel sum should be unity)
    psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(1.0 / tSum, PS_TYPE_F32));

    return kernel;
}

psImage *psImageConvolveDirect(psImage *out,
                               const psImage *in,
                               const psKernel *kernel)
{
    PS_ASSERT_IMAGE_NON_NULL(in, NULL);
    if (in->type.type != PS_TYPE_S32 && in->type.type != PS_TYPE_F32 && in->type.type != PS_TYPE_F64) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                "Unallowable operation: psImage %s is not of type S32 or F{32,64}.", "in");
        return(NULL);
    }
    PS_ASSERT_PTR_NON_NULL(kernel, NULL);
    PS_ASSERT_PTR_NON_NULL(kernel->kernel, NULL);

    // Pull out kernel parameters, for convenience
    int xMin = kernel->xMin;
    int xMax = kernel->xMax;
    int yMin = kernel->yMin;
    int yMax = kernel->yMax;
    float **kernelData = kernel->kernel;

    int numRows = in->numRows;          // Number of rows
    int numCols = in->numCols;          // Number of columns

#if 1

    // This is the usual way of doing the convolution, but it wastes time if there's a large number of zero
    // pixels in the kernel.
#define SPATIAL_CONVOLVE_CASE(TYPE) \
    case PS_TYPE_##TYPE: { \
        ps##TYPE **inData = in->data.TYPE; /* Dereference input data */ \
        out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \
        for (int row = 0; row < numRows; row++) { \
            ps##TYPE *outRow = out->data.TYPE[row]; \
            int kRowMin = PS_MAX(yMin, row + 1 - numRows); \
            int kRowMax = PS_MIN(yMax, row); \
            for (int col = 0; col < numCols; col++) { \
                int kColMin = PS_MAX(xMin, col + 1 - numCols); \
                int kColMax = PS_MIN(xMax, col); \
                ps##TYPE pixel = 0.0; \
                for (int kRow = kRowMin; kRow <= kRowMax; kRow++) { \
                    for (int kCol = kColMin; kCol <= kColMax; kCol++) { \
                        pixel += kernelData[kRow][kCol] * inData[row - kRow][col - kCol]; \
                    } \
                } \
                outRow[col] = pixel; \
            } \
        } \
    } \
    break;

#else

    // This turns the convolution inside-out, allowing us to skip kernel pixels that have no contribution.
#define SPATIAL_CONVOLVE_CASE(TYPE) \
    case PS_TYPE_##TYPE: { \
        ps##TYPE **inData = in->data.TYPE; /* Dereference input data */ \
        out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \
        psImageInit(out, 0.0); \
        for (int ky = yMin; ky <= yMax; ky++) { \
            for (int kx = xMin; kx <= xMax; kx++) { \
                float kValue = kernelData[ky][kx]; /* Kernel value */ \
                if (kValue == 0.0) { \
                    continue; \
                } \
                for (int y = PS_MAX(ky, 0); y < PS_MIN(numRows, numRows + ky); y++) { \
                    for (int x = PS_MAX(kx, 0); x < PS_MIN(numCols, numCols + kx); x++) { \
                        out->data.TYPE[y][x] += kValue * inData[y - ky][x - kx]; \
                    } \
                } \
            } \
        } \
    } \
    break;

#endif

    switch (in->type.type) {
        SPATIAL_CONVOLVE_CASE(S32);
        SPATIAL_CONVOLVE_CASE(F32);
        SPATIAL_CONVOLVE_CASE(F64);
      default:
        psAbort("Should never get here: bad type that was asserted on previously.");
    }

    return out;
}

psImage *psImageConvolveMaskDirect(psImage *out, const psImage *mask, psMaskType maskVal,
                                   psMaskType setVal, int xMin, int xMax, int yMin, int yMax)
{
    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL);
    if (out == mask && ((maskVal & setVal) || !setVal)) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                "Can't convolve mask in-place if values to set contains values to convolve.");
        return NULL;
    }

    if (yMin > yMax) {
        psWarning("Specified yMin, %d, was greater than yMax, %d.  Values swapped.",
                 yMin, yMax);

        int temp = yMin;
        yMin = yMax;
        yMax = temp;
    }
    if (xMin > xMax) {
        psWarning("Specified xMin, %d, was greater than xMax, %d.  Values swapped.",
                 xMin, xMax);

        int temp = xMin;
        xMin = xMax;
        xMax = temp;
    }

    int numRows = mask->numRows;        // Number of rows
    int numCols = mask->numCols;        // Number of columns

    if (!out) {
        // Propagate the non-masked values
        out = (psImage*)psBinaryOp(NULL, (const psPtr)mask, "&", psScalarAlloc(~maskVal, PS_TYPE_MASK));
    }

    // Dereference mask images
    psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA;
    psMaskType **outData = out->data.PS_TYPE_MASK_DATA;

    if (setVal) {
        // Grow any pixels matching maskVal, setting setVal
        for (int row = 0; row < numRows; row++) {
            for (int col = 0; col < numCols; col++) {
                if (maskData[row][col] & maskVal) {
                    for (int kRow = PS_MAX(yMin, -row); kRow <= PS_MIN(yMax, numRows - row - 1); kRow++) {
                        for (int kCol = PS_MAX(xMin, -col); kCol <= PS_MIN(xMax, numCols - col - 1); kCol++) {
                            outData[row + kRow][col + kCol] |= setVal;
                        }
                    }
                }
            }
        }
    } else {
        // Each pixel receives any maskVal bits set within the convolution window
        for (int row = 0; row < numRows; row++) {
            for (int col = 0; col < numCols; col++) {
                psMaskType pixel = outData[row][col]; // Pixel value to set
                if (pixel & maskVal) {
                    // Already done this one
                    continue;
                }
                for (int kRow = PS_MAX(yMin, -row); kRow <= PS_MIN(yMax, numRows - row - 1); kRow++) {
                    for (int kCol = PS_MAX(xMin, -col); kCol <= PS_MIN(xMax, numCols - col - 1); kCol++) {
                        pixel |= maskData[row][col] & maskVal;
                    }
                }
                outData[row][col] = pixel;
            }
        }
    }

    return out;
}


psImage *psImageConvolveMaskFFT(psImage *out, const psImage *mask, psMaskType maskVal,
                                psMaskType setVal, int xMin, int xMax, int yMin, int yMax, float thresh)
{
    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL);
    PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, NULL);
    PS_ASSERT_FLOAT_LESS_THAN(thresh, 1.0, NULL);
    if (out == mask && (maskVal & setVal)) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                "Can't convolve mask in-place if values to set contains values to convolve.");
        return NULL;
    }

    if (yMin > yMax) {
        psWarning("Specified yMin, %d, was greater than yMax, %d.  Values swapped.",
                 yMin, yMax);

        int temp = yMin;
        yMin = yMax;
        yMax = temp;
    }
    if (xMin > xMax) {
        psWarning("Specified xMin, %d, was greater than xMax, %d.  Values swapped.",
                 xMin, xMax);

        int temp = xMin;
        xMin = xMax;
        xMax = temp;
    }

    int numRows = mask->numRows, numCols = mask->numCols; // Size of image

    psImage *onoff = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Pixels on or off
    psImageInit(onoff, 0);
    for (int y = 0; y < numRows; y++) {
        for (int x = 0; x < numCols; x++) {
            if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) {
                onoff->data.F32[y][x] = 1.0;
            }
        }
    }

    psKernel *kernel = psKernelAlloc(xMin, xMax, yMin, yMax);
    psImageInit(kernel->image, 1.0);
    psImage *convolved = psImageConvolveFFT(NULL, onoff, NULL, 0, kernel);
    psFree(onoff);
    psFree(kernel);
    if (!convolved) {
        psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask.");
        return NULL;
    }

    if (!setVal) {
        setVal = maskVal;
    }

    if (!out) {
        out = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
    }
    for (int y = 0; y < numRows; y++) {
        for (int x = 0; x < numCols; x++) {
            out->data.PS_TYPE_MASK_DATA[y][x] = (convolved->data.F32[y][x] >= thresh) ?
                (mask->data.PS_TYPE_MASK_DATA[y][x] | setVal) : mask->data.PS_TYPE_MASK_DATA[y][x];
        }
    }

    psFree(convolved);

    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.25) ? s / g : NAN;
        }
        memcpy(image->data.F32[j], calculation->data.F32, Nx*sizeof(psF32));
    }
    psFree(calculation);

    // XXX test
    if (0) {
        psFree(gaussnorm);
        return true;
    }

    /** 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;
                if (!isfinite(*vi)) 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.25) ? *vo / *vs : NAN;
        }

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