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

#include <stdio.h>
#include <pslib.h>

#include "pmSubtractionKernels.h"

#define RINGS_BUFFER 10                 // Buffer size for RINGS data


// Free function for pmSubtractionKernels
static void subtractionKernelsFree(pmSubtractionKernels *kernels)
{
    psFree(kernels->u);
    psFree(kernels->v);
    psFree(kernels->widths);
    psFree(kernels->uStop);
    psFree(kernels->vStop);
    psFree(kernels->xOrder);
    psFree(kernels->yOrder);
    psFree(kernels->preCalc);
}

// Raise a number to an integer power
static inline long power(int value,     // Value
                         int exp        // Exponent
    )
{
    if (exp == 0) {
        return 1.0;
    }
    long result = value;               // Result to return
    for (int i = 2; i <= exp; i++) {
        result *= value;
    }
    return result;
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Public functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
                                                int size, int spatialOrder)
{
    pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return
    psMemSetDeallocator(kernels, (psFreeFunc)subtractionKernelsFree);

    kernels->type = type;
    kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
    kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
    kernels->widths = NULL;
    kernels->uStop = NULL;
    kernels->vStop = NULL;
    kernels->xOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
    kernels->yOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
    kernels->subIndex = 0;
    kernels->preCalc = NULL;
    kernels->size = size;
    kernels->inner = 0;
    kernels->spatialOrder = spatialOrder;

    return kernels;
}

pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder)
{
    PS_ASSERT_INT_POSITIVE(size, NULL);
    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);

    int num = PS_SQR(2 * size + 1) * (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions

    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS,
                                                              size, spatialOrder); // The kernels

    psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
             size, spatialOrder, num);

    // Generate a set of kernels for each (u,v)
    for (int v = - size, index = 0; v <= size; v++) {
        for (int u = - size; u <= size; u++) {
            // Iterate over spatial order.  This loop creates the terms for
            // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
            for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
                for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
                    kernels->u->data.S32[index] = u;
                    kernels->v->data.S32[index] = v;
                    kernels->xOrder->data.S32[index] = xOrder;
                    kernels->yOrder->data.S32[index] = yOrder;

                    psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
                            u, v, xOrder, yOrder);
                }
            }
        }
    }

    kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2;
    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
           kernels->v->data.S32[kernels->subIndex] == 0 &&
           kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
           kernels->yOrder->data.S32[kernels->subIndex] == 0);

    return kernels;
}

pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
                                               const psVector *sigmas, const psVector *orders)
{
    PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL);
    PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL);
    PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
    PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
    PS_ASSERT_VECTORS_SIZE_EQUAL(sigmas, orders, NULL);
    PS_ASSERT_INT_POSITIVE(size, NULL);
    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);

    int numGaussians = sigmas->n;       // Number of Gaussians

    int num = 0;                        // Number of basis functions
    psString params = NULL;             // List of parameters
    for (int i = 0; i < numGaussians; i++) {
        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
        psStringAppend(&params, "(%.2f,%d)", sigmas->data.F32[i], orders->data.S32[i]);
        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
    }
    num *= (spatialOrder + 1) * (spatialOrder + 2) / 2;

    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS,
                                                              size, spatialOrder); // The kernels

    psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS kernel: %s,%d --> %d elements",
             params, spatialOrder, num);
    psFree(params);

    kernels->widths = psVectorAlloc(num, PS_TYPE_F32);
    kernels->preCalc = psArrayAlloc(num);

    // Set the kernel parameters
    for (int i = 0, index = 0; i < numGaussians; i++) {
        // Iterate over (u,v) order
        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++) {

                // Set the pre-calculated kernel
                psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
                for (int v = -size; v <= size; v++) {
                    for (int u = -size; u <= size; u++) {
                        preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *
                            expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
                    }
                }

                // Iterate over spatial order.  This loop creates the terms for
                // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
                for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
                    for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
                        kernels->widths->data.F32[index] = sigmas->data.F32[i];
                        kernels->u->data.S32[index] = uOrder;
                        kernels->v->data.S32[index] = vOrder;
                        kernels->xOrder->data.S32[index] = xOrder;
                        kernels->yOrder->data.S32[index] = yOrder;
                        kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc);

                        psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %d %d\n", index,
                                sigmas->data.F32[i], uOrder, vOrder, xOrder, yOrder);
                    }
                }

                psFree(preCalc);        // Drop reference
            }
        }
    }

    kernels->subIndex = 0;
    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
           kernels->v->data.S32[kernels->subIndex] == 0 &&
           kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
           kernels->yOrder->data.S32[kernels->subIndex] == 0);

    if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
        for (int i = 0; i < num; i++) {
            psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
            psString kernelName = NULL;
            psStringAppend(&kernelName, "kernel%03d.fits", i);
            psFits *kernelFile = psFitsOpen(kernelName, "w");
            psFree(kernelName);
            psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL);
            psFitsClose(kernelFile);
        }
    }

    return kernels;
}

/// Generate SPAM kernels
pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel
                                               int spatialOrder, ///< Order of spatial variations
                                               int inner, ///< Inner radius to preserve unbinned
                                               int binning ///< Kernel binning factor
    )
{
    PS_ASSERT_INT_POSITIVE(size, NULL);
    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
    PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
    PS_ASSERT_INT_LARGER_THAN(size, inner, NULL);
    PS_ASSERT_INT_POSITIVE(binning, NULL);

    // The outer region should be divisible by the "binning"; otherwise allocate remainder to the inner region
    int numOuter = (size - inner) / binning; // Number of summed pixels in the outer region
    int numInner = inner + (size - inner) % binning; // Number of pixels in the inner region
    assert(numOuter * binning + numInner == size);
    int numTotal = numOuter + numInner; // Total number of summed pixels

    psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter);

    int num = PS_SQR(2 * numTotal + 1) *
        (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions

    psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);

    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM,
                                                              size, spatialOrder); // The kernels

    psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements",
             size, inner, binning, spatialOrder, num);

    kernels->uStop = psVectorAlloc(num, PS_TYPE_F32);
    kernels->vStop = psVectorAlloc(num, PS_TYPE_F32);

    psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element
    psVector *widths = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Widths for each kernel element
    locations->data.S32[numTotal] = 0;
    widths->data.S32[numTotal] = 0;
    for (int i = 1; i <= numInner; i++) {
        locations->data.S32[numTotal + i] = i;
        widths->data.S32[numTotal + i] = 0;
        locations->data.S32[numTotal - i] = - i;
        widths->data.S32[numTotal - i] = 0;
    }
    for (int i = numInner + 1; i <= numTotal; i++) {
        locations->data.S32[numTotal + i] = locations->data.S32[numTotal + i - 1] +
            widths->data.S32[numTotal + i - 1] + 1;
        widths->data.S32[numTotal + i] = binning - 1;
        locations->data.S32[numTotal - i] = locations->data.S32[numTotal - i + 1] - binning;
        widths->data.S32[numTotal - i] = binning - 1;
    }

    if (psTraceGetLevel("psModules.imcombine") >= 10) {
        for (int i = 0; i < 2 * numTotal + 1; i++) {
            psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, locations->data.S32[i],
                    locations->data.S32[i] + widths->data.S32[i]);
        }
    }

    // Set the kernel parameters
    for (int i = - numTotal, index = 0; i <= numTotal; i++) {
        int u = locations->data.S32[numTotal + i]; // Location of pixel
        int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel

        for (int j = - numTotal; j <= numTotal; j++) {
            int v = locations->data.S32[numTotal + j]; // Location of pixel
            int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel

            // Iterate over spatial order.  This loop creates the terms for
            // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
            for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
                for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
                    kernels->u->data.S32[index] = u;
                    kernels->v->data.S32[index] = v;
                    kernels->uStop->data.S32[index] = uStop;
                    kernels->vStop->data.S32[index] = vStop;
                    kernels->xOrder->data.S32[index] = xOrder;
                    kernels->yOrder->data.S32[index] = yOrder;

                    psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index,
                            u, uStop, v, vStop, xOrder, yOrder);
                }
            }
        }
    }

    kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2;
    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
           kernels->v->data.S32[kernels->subIndex] == 0 &&
           kernels->uStop->data.S32[kernels->subIndex] == 0 &&
           kernels->vStop->data.S32[kernels->subIndex] == 0 &&
           kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
           kernels->yOrder->data.S32[kernels->subIndex] == 0);

    psFree(locations);
    psFree(widths);

    return kernels;
}


/// Generate FRIES kernels
pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel
                                                int spatialOrder, ///< Order of spatial variations
                                                int inner ///< Inner radius to preserve unbinned
    )
{
    PS_ASSERT_INT_POSITIVE(size, NULL);
    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
    PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
    PS_ASSERT_INT_LARGER_THAN(size, inner, NULL);

    int fibNum = 0;                     // Number of Fibonacci values
    int fibLast = 1, fibTotal = 2;      // Fibonacci sequence
    while (fibTotal < size - inner) {
        int temp = fibTotal;
        fibTotal += fibLast;
        fibLast = temp;
        fibNum++;
    }

    int numInner = inner;               // Number of pixels in the inner region
    int numOuter = fibNum;              // Number of summed pixels in the outer region
    int numTotal = numOuter + numInner; // Total number of summed pixels

    psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter);

    int num = PS_SQR(2 * numTotal + 1) *
        (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions

    psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);

    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES,
                                                              size, spatialOrder); // The kernels

    psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements",
             size, inner, spatialOrder, num);

    kernels->uStop = psVectorAlloc(num, PS_TYPE_F32);
    kernels->vStop = psVectorAlloc(num, PS_TYPE_F32);

    psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32);
    psVector *stop = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32);
    start->data.S32[numTotal] = 0;
    stop->data.S32[numTotal] = 0;
    for (int i = 1; i <= numInner; i++) {
        start->data.S32[numTotal + i] = i;
        stop->data.S32[numTotal + i] = i;
        start->data.S32[numTotal - i] = -i;
        stop->data.S32[numTotal - i] = -i;
    }
    for (int i = numInner + 1, fibLast = 1, fib = 2, temp; i <= numTotal;
         i++, fib = (temp = fib) + fibLast, fibLast = temp) {
        start->data.S32[numTotal + i] = stop->data.S32[numTotal + i - 1] + 1;
        stop->data.S32[numTotal + i] = PS_MIN(start->data.S32[numTotal + i] + fib - 1, size);
        start->data.S32[numTotal - i] = - stop->data.S32[numTotal + i];
        stop->data.S32[numTotal - i] = - start->data.S32[numTotal + i];
    }

    if (psTraceGetLevel("psModules.imcombine") >= 10) {
        for (int i = 0; i < 2 * numTotal + 1; i++) {
            psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, start->data.S32[i], stop->data.S32[i]);
        }
    }

    // Set the kernel parameters
    for (int i = - numTotal, index = 0; i <= numTotal; i++) {
        int u = start->data.S32[numTotal + i]; // Location of pixel
        int uStop = stop->data.S32[numTotal + i]; // Width of pixel
        for (int j = - numTotal; j <= numTotal; j++) {
            int v = start->data.S32[numTotal + j]; // Location of pixel
            int vStop = stop->data.S32[numTotal + j]; // Width of pixel

            // Iterate over spatial order.  This loop creates the terms for
            // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
            for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
                for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
                    kernels->u->data.S32[index] = u;
                    kernels->v->data.S32[index] = v;
                    kernels->uStop->data.S32[index] = uStop;
                    kernels->vStop->data.S32[index] = vStop;
                    kernels->xOrder->data.S32[index] = xOrder;
                    kernels->yOrder->data.S32[index] = yOrder;

                    psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index,
                            u, uStop, v, vStop, xOrder, yOrder);
                }
            }
        }
    }

    kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2;
    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
           kernels->v->data.S32[kernels->subIndex] == 0 &&
           kernels->uStop->data.S32[kernels->subIndex] == 0 &&
           kernels->vStop->data.S32[kernels->subIndex] == 0 &&
           kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
           kernels->yOrder->data.S32[kernels->subIndex] == 0);

    psFree(start);
    psFree(stop);

    return kernels;
}

// Grid United with Normal Kernel
pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *sigmas,
                                               const psVector *orders, int inner)
{
    PS_ASSERT_INT_POSITIVE(size, NULL);
    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
    PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL);
    PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL);
    PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
    PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
    PS_ASSERT_VECTORS_SIZE_EQUAL(sigmas, orders, NULL);
    PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
    PS_ASSERT_INT_LESS_THAN(inner, size, NULL);

    int numGaussians = sigmas->n;       // Number of Gaussians
    int numGaussianVars = 0;            // Number of Gaussian variant functions in the kernel
    psString params = NULL;             // List of params
    for (int i = 0; i < numGaussians; i++) {
        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
        numGaussianVars += (gaussOrder + 1) * (gaussOrder + 2) / 2;
        psStringAppend(&params, "(%.2f,%d)", sigmas->data.F32[i], orders->data.S32[i]);
    }

    int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements
    int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel

    int num = (numGaussianVars + numInner) * numSpatial; // Total number of basis functions

    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK,
                                                              size, spatialOrder); // The kernels

    psLogMsg("psModules.imcombine", PS_LOG_INFO, "GUNK kernel: %d,%s,%d --> %d elements",
             inner, params, spatialOrder, num);
    psFree(params);

    kernels->widths = psVectorAlloc(numGaussianVars * numSpatial, PS_TYPE_F32);
    kernels->preCalc = psArrayAlloc(numGaussianVars * numSpatial);
    kernels->inner = numGaussianVars * numSpatial;

    // Set the Gaussian kernel parameters
    for (int i = 0, index = 0; i < numGaussians; i++) {
        // Iterate over (u,v) order
        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++) {


                // Set the pre-calculated kernel
                psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
                for (int v = -size; v <= size; v++) {
                    for (int u = -size; u <= size; u++) {
                        preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *
                            expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
                    }
                }

                // Iterate over spatial order.  This loop creates the terms for
                // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
                for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
                    for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
                        kernels->widths->data.F32[index] = sigmas->data.F32[i];
                        kernels->u->data.S32[index] = uOrder;
                        kernels->v->data.S32[index] = vOrder;
                        kernels->xOrder->data.S32[index] = xOrder;
                        kernels->yOrder->data.S32[index] = yOrder;
                        kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc);

                        psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %d %d\n", index,
                                sigmas->data.F32[i], uOrder, vOrder, xOrder, yOrder);
                    }
                }

                psFree(preCalc);        // Drop reference
            }
        }
    }


    // Generate a grid set of kernels for each (u,v)
    for (int v = - inner, index = kernels->inner; v <= inner; v++) {
        for (int u = - inner; u <= inner; u++) {
            // Iterate over spatial order.  This loop creates the terms for
            // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
            for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
                for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
                    kernels->u->data.S32[index] = u;
                    kernels->v->data.S32[index] = v;
                    kernels->xOrder->data.S32[index] = xOrder;
                    kernels->yOrder->data.S32[index] = yOrder;

                    psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
                            u, v, xOrder, yOrder);
                }
            }
        }
    }

    kernels->subIndex = kernels->inner + (numInner - 1) * numSpatial / 2;
    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
           kernels->v->data.S32[kernels->subIndex] == 0 &&
           kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
           kernels->yOrder->data.S32[kernels->subIndex] == 0);

    return kernels;
}

// RINGS --- just what it says
pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder)
{
    PS_ASSERT_INT_POSITIVE(size, NULL);
    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
    PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
    PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
    PS_ASSERT_INT_NONNEGATIVE(ringsOrder, NULL);

    int fibNum = 0;                     // Number of Fibonacci values
    {
        int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters
        int radius = inner;
        while (radius < size) {
            radius++;
            int fibNew = fibIndex + fibIndexMinus1;
            fibIndexMinus1 = fibIndex;
            fibIndex = fibNew;
            radius += fibIndex;
            fibNum++;
            printf("%d ", radius);
        }
        printf("\n");
    }

    int numInner = inner;               // Number of pixels in the inner region
    int numOuter = fibNum;              // Number of summed pixels in the outer region

    int numRings = numOuter + numInner; // Number of rings (not including the central pixel)
    int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring
    int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel

    int num = (numRings * numPoly + 1) * numSpatial; // Total number of basis functions

    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS,
                                                              size, spatialOrder); // The kernels

    psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements",
             size, inner, ringsOrder, spatialOrder, num);

    kernels->preCalc = psArrayAlloc(num);

    // Set the Gaussian kernel parameters
    int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters
    int radiusLast = 0;                 // Last radius
    for (int i = 0, index = 0; i < numRings + 1; i++) {

        float lower2;           // Lower limit of radius^2
        float upper2;           // Upper limit of radius^2
        if (i == 0) {
            lower2 = 0;
            upper2 = PS_SQR(0.5);
        } else if (i <= numInner) {
            // A ring every pixel width
            float radius = i;
            lower2 = PS_SQR(radius - 0.5);
            upper2 = PS_SQR(radius + 0.5);
            radiusLast = i;
        } else {
            // Rings Fibonacci distributed (2, 3, 5...)
            int fibNew = fibIndex + fibIndexMinus1;
            fibIndexMinus1 = fibIndex;
            fibIndex = fibNew;

            float radiusLower = radiusLast + 1;
            radiusLast = radiusLower + fibIndex;
            float radiusUpper = radiusLast;

            lower2 = PS_SQR(radiusLower - 0.5);
            upper2 = PS_SQR(radiusUpper + 0.5);
        }

        psTrace("psModules.imcombine", 8, "Radius limits: %f --> %f\n", sqrtf(lower2), sqrtf(upper2));

        // Iterate over (u,v) order
        for (int uOrder = 0; uOrder <= (i == 0 ? 0 : ringsOrder); uOrder++) {
            for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder) - uOrder; vOrder++) {

                psArray *data = psArrayAlloc(3); // Container for data
                psVector *uCoords = data->data[0] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // u coords
                psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords
                psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial

                if (i == 0) {
                    // Central pixel is easy
                    uCoords->data.S32[0] = vCoords->data.S32[0] = 0;
                    poly->data.F32[0] = 1.0;
                    uCoords->n = vCoords->n = poly->n = 1;
                    radiusLast = 0;
                } else {
                    int j = 0;          // Index for data
                    for (int v = -size; v <= size; v++) {
                        int v2 = PS_SQR(v);   // Square of v
                        float vPoly = power(v, vOrder); // Value of v^vOrder

                        for (int u = -size; u <= size; u++) {
                            int u2 = PS_SQR(u); // Square of u
                            int distance2 = u2 + v2; // Distance from the centre
                            if (distance2 > lower2 && distance2 < upper2) {
                                float uPoly = power(u, uOrder); // Value of u^uOrder

                                uCoords->data.S32[j] = u;
                                vCoords->data.S32[j] = v;
                                poly->data.F32[j] = uPoly * vPoly;

                                psVectorExtend(uCoords, RINGS_BUFFER, 1);
                                psVectorExtend(vCoords, RINGS_BUFFER, 1);
                                psVectorExtend(poly, RINGS_BUFFER, 1);

                                psTrace("psModules.imcombine", 9, "u = %d, v = %d, poly = %f\n",
                                        u, v, poly->data.F32[j]);

                                j++;
                            }
                        }
                    }
                }

                psTrace("psModules.imcombine", 8, "%ld pixels in ring\n", uCoords->n);

                // Iterate over spatial order.  This loop creates the terms for
                // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
                for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
                    for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
                        kernels->preCalc->data[index] = psMemIncrRefCounter(data);
                        kernels->u->data.S32[index] = uOrder;
                        kernels->v->data.S32[index] = vOrder;
                        kernels->xOrder->data.S32[index] = xOrder;
                        kernels->yOrder->data.S32[index] = yOrder;

                        psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d\n", index,
                                i, uOrder, vOrder, xOrder, yOrder);
                    }
                }
                psFree(data);
            }
        }
    }

    kernels->subIndex = 0;
    assert(kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
           kernels->yOrder->data.S32[kernels->subIndex] == 0);

    return kernels;
}

pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
                                                   const psVector *sigmas, const psVector *orders, int inner,
                                                   int binning, int ringsOrder)
{
    switch (type) {
      case PM_SUBTRACTION_KERNEL_POIS:
        return pmSubtractionKernelsPOIS(size, spatialOrder);
      case PM_SUBTRACTION_KERNEL_ISIS:
        return pmSubtractionKernelsISIS(size, spatialOrder, sigmas, orders);
      case PM_SUBTRACTION_KERNEL_SPAM:
        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning);
      case PM_SUBTRACTION_KERNEL_FRIES:
        return pmSubtractionKernelsFRIES(size, spatialOrder, inner);
      case PM_SUBTRACTION_KERNEL_GUNK:
        return pmSubtractionKernelsGUNK(size, spatialOrder, sigmas, orders, inner);
      case PM_SUBTRACTION_KERNEL_RINGS:
        return pmSubtractionKernelsRINGS(size, ringsOrder, inner, spatialOrder);
      default:
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
        return NULL;
    }
}
