// Copyright (C) 2006  Paul A. Price (price@ifa.hawaii.edu)
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <assert.h>
#include "pslib.h"
#include "pois.h"


/* Convolve an image by the kernel */
psImage *poisConvolveImage(const psImage *input, // Input image, to be convolved
                           const psImage *mask, // Mask image
                           const psVector *solution, // The solution vector, containing the coefficients
                           const psArray *kernelParams, // The kernel parameters
                           const poisConfig *config // The configuration
    )
{
    assert(solution->n == kernelParams->n + 1); // Extra parameter for background
    assert(solution->type.type == PS_TYPE_F64);
    assert(input->numCols == mask->numCols && input->numRows == mask->numRows);
    assert(mask->type.type == PS_TYPE_U8);
    assert(input->type.type == PS_TYPE_F32);

    int nx = input->numCols;            // The size of the image in x
    int ny = input->numRows;            // The size of the image in y
    float nxHalf = 0.5 * (float)nx;     // Half the size of the image
    float nyHalf = 0.5 * (float)ny;     // Half the size of the image
    int xKernel = config->xKernel;      // The half-size of the kernel in x
    int yKernel = config->yKernel;      // The half-size of the kernel in y
    float background = solution->data.F64[solution->n-1]; // The difference in background
    int numParams = kernelParams->n;    // Number of kernel parameters

    psImage *convolved = psImageAlloc(nx, ny, PS_TYPE_F32); // The convolved image, to be returned

    psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, config->spatialOrder + 1,
                                               config->spatialOrder + 1); // Polynomial
    for (int i = 0; i < config->spatialOrder + 1; i++) {
        for (int j = 0; j < config->spatialOrder + 1; j++) {
            poly->coeff[i][j] = 1.0;
            poly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
        }
    }

    // Convolve only the middle part
    for (int y = yKernel; y < ny - yKernel; y++) {
        for (int x = xKernel; x < nx - xKernel; x++) {
            if (! mask->data.U8[y][x] & (POIS_MASK_BAD | POIS_MASK_NEAR_BAD)) {
                double conv = background; // Convolved value

                float imageX = (x - nxHalf) / nxHalf; // Normalised position in x
                float imageY = (y - nyHalf) / nyHalf; // Normalised position in y

                // Iterate over the kernel basis functions
                for (int k = 0; k < numParams; k++) {
                    poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function
                    int u = kernel->u;
                    int v = kernel->v;
                    int xOrder = kernel->xOrder;
                    int yOrder = kernel->yOrder;

                    double polyVal = 1.0;
                    if (xOrder != 0 || yOrder != 0) {
                        // Evaluate the polynomial
                        poly->mask[xOrder][yOrder] = 0;
                        polyVal = psPolynomial2DEval(poly, imageX, imageY);
                        poly->mask[xOrder][yOrder] = 1;
                    }

                    if (k == 0) {
                        conv += solution->data.F64[k] * input->data.F32[y - v][x - u] * polyVal;
                    } else {
                        conv += solution->data.F64[k] *
                            (input->data.F32[y - v][x - u] * polyVal - input->data.F32[y][x]);
                    }
                }
                convolved->data.F32[y][x] = (float)conv;
            }
        }

        /* Pad the rest with zeros */
        for (int x = 0; x < xKernel; x++) {
            convolved->data.F32[y][x] = 0.0;
        }
        for (int x = nx - xKernel; x < nx; x++) {
            convolved->data.F32[y][x] = 0.0;
        }
    }

    for (int y = 0; y < yKernel; y++) {
        for (int x = 0; x < nx; x++) {
            convolved->data.F32[y][x] = 0.0;
        }
    }
    for (int y = ny - yKernel; y < ny; y++) {
        for (int x = 0; x < nx; x++) {
            convolved->data.F32[y][x] = 0.0;
        }
    }

    psFree(poly);

    return convolved;
}


