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

psArray *poisKernelBasisFunctions(const poisConfig *config // Configuration
    )
{
    assert(config);

    psTrace("pois.kernelBasisFunctions", 1, "Setting kernel basis functions....\n");

    int xKernelHalfSize = config->xKernel; // Half size of kernel in x
    int yKernelHalfSize = config->yKernel; // Half size of kernel in y
    int spatialOrder = config->spatialOrder; // Order of spatial variations

    // Number of basis functions
    int nBF = (2 * xKernelHalfSize + 1) * (2 * yKernelHalfSize + 1) * (spatialOrder + 1) *
        (spatialOrder + 2) / 2;

    psArray *array = psArrayAlloc(nBF); // Array to hold the basis functions

    int num = 0;                        // Kernel parameter number

    // Put the (u,v) = (0,0) component right at the start of the list for convenience
    for (int order = 0; order <= spatialOrder; order++) {
        for (int i = 0; i <= order; i++) {
            int j = order - i;
            poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis));
            kernel->u = 0;
            kernel->v = 0;
            kernel->xOrder = i;
            kernel->yOrder = j;

            // Stuff into the array
            array->data[num] = kernel;

            psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", i, kernel->u,
                    kernel->v, kernel->xOrder, kernel->yOrder);
            num++;
        }
    }

    // Iterate over (u,v)
    for (int u = -xKernelHalfSize; u <= xKernelHalfSize; u++) {
        for (int v = -yKernelHalfSize; v <= yKernelHalfSize; v++) {

            // Already did (u,v) = (0,0): it's at the start, so skip it now.
            if ((u != 0) || (v != 0)) {

                // Iterate over spatial order
                for (int order = 0; order <= spatialOrder; order++) {

                    for (int i = 0; i <= order; i++) {
                        int j = order - i;
                        poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis));

                        kernel->u = u;
                        kernel->v = v;
                        kernel->xOrder = i;
                        kernel->yOrder = j;

                        // Stuff into the array
                        array->data[num] = kernel;

                        psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", num,
                                kernel->u, kernel->v, kernel->xOrder, kernel->yOrder);
                        num++;
                    }
                } // Iterating over spatial order
            } // Except (u,v) = (0,0)
        }
    } // Iterating over (u,v)

    return array;
}

