// 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"

/* Create a kernel image from the solution vector */
psImage *poisExtractKernel(const psVector *solution, // Vector containing the kernel elements
                           const psArray *kernelParams, // The kernel parameters
                           float x,     // The relative (-1 to 1) x position for which to get the kernel
                           float y,     // The relative (-1 to 1) y position for which to get the kernel
                           const poisConfig *config // Configuration
    )
{
    assert(solution->n == kernelParams->n + 1);
    assert(solution->type.type == PS_TYPE_F64);
    assert(config);
    assert(x >= -1.0 && x <= 1.0);
    assert(y >= -1.0 && y <= 1.0);

    int xKernel = config->xKernel;      // Half-size of kernel in x
    int yKernel = config->yKernel;      // Half-size of kernel in y
    float kernelSum = 0.0;              // Sum of the kernel

    /* Allocate the image */
    psImage *image = psImageAlloc(2*xKernel + 1, 2*yKernel + 1, PS_TYPE_F32); // The kernel image
    for (int j = 0; j < (2*yKernel + 1); j++) {
        for (int i = 0; i < (2*xKernel + 1); i++) {
            image->data.F32[j][i] = 0.0;
        }
    }

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

    psTrace("pois.extractKernel", 4, "Evaluating kernel at %f,%f\n", x, y);
    /* Iterate over the kernel parameters */
    for (int k = 0; k < kernelParams->n; 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;

        // Evaluate polynomial
        poly->mask[xOrder][yOrder] = 0;
        double evaluation = solution->data.F64[k] * psPolynomial2DEval(poly, x, y);
        poly->mask[xOrder][yOrder] = 1;

        image->data.F32[v + yKernel][u + xKernel] += (float)evaluation;
        kernelSum += (float)evaluation;
    }

    // Subtract from the centre pixel to remove the effect of the flux-conservation treatment
    // This assumes that the "reference" is the centre pixel: u = 0, v = 0, xOrder = 0, yOrder = 0
    image->data.F32[yKernel][xKernel] -= kernelSum - (float)solution->data.F64[0];

    psTrace("pois.extractKernel", 4, "Kernel sum: %f\n", kernelSum);

    psFree(poly);

    return image;
}

