#ifndef INC_PYRAMIDS
#define INC_PYRAMIDS

#include <algorithm>
#include <vector>

#include "FImage.hpp"

// Free-space 2D Green function convolution kernel approximations F_{5,3}
// here h1 = h2 =: h12
static Matrix<float, 1, 5> h12;
static Matrix<float, 5, 1> h12t;
static const Matrix<float, 1, 3> g(0.175, 0.547, 0.175);
static const Matrix<float, 3, 1> gt(0.175, 0.547, 0.175);


// Same-size backtracking memory
static std::vector<FImage> a, b, c;

void poissonPrepare(int w, int h)
{
    // Requires GCC, not portable, oh well
    int nbLevels = 32 - __builtin_clz(std::max(w, h));
    
    a.clear();
    b.clear();
    c.clear();
    a.reserve(nbLevels + 1);
    b.reserve(nbLevels + 1);
    c.reserve(nbLevels + 1);
    for(int level = 0; level < nbLevels + 1; level++)
    {
        w += 10; h += 10;
        a.push_back(FImage(h, w));
        b.push_back(a.back());
        c.push_back(a.back());
        w = (w + 1) / 2; h = (h + 1) / 2;
    }
    
    h12 << 0.15, 0.5, 0.7, 0.5, 0.15;
    h12t << 0.15, 0.5, 0.7, 0.5, 0.15;
}

/**
 * Convolution pyramids for Poisson equation solving on images.
 * Based on Convolution Pyramids, Farbman et. al, 2011.
 */
template <typename Derived>
FImage poissonIntegrate(const DenseBase<Derived> &src)
{
    Array<float, 5, 5> h12f = h12t * h12;
    Array<float, 3, 3> gf = gt * g;
    
    a[0].setZero();
    a[0].block(5, 5, src.rows(), src.cols()) = src;
    
    int nbLevels = a.size() - 1;
    
    // Forward transform
    for(int level = 0; level < nbLevels; level++)
    {
        convolveSepInto<5, 5>(a[level], b[level], h12, h12t);
        // convolveInto<5, 5>(a[level], b[level], h12f);
        a[level + 1].setZero();
        int br = (b[level].rows() + 1) / 2, bc = (b[level].cols() + 1) / 2;
        a[level + 1].block(5, 5, br, bc) = Map<FImage, 0, Stride<Dynamic, 2>>(b[level].data(), br, bc,
            Stride<Dynamic, 2>(b[level].outerStride() * 2, 2));
    }
    // At this point b is free to use and a holds the downsampled convolutions with h1
    
    // Backward transform
    for(int level = 0; level <= nbLevels; level++)
        convolveSepInto<3, 3>(a[level], b[level], g, gt);
        // convolveInto<3, 3>(a[level], b[level], gf);
    // At this point a is free to use and b holds the convolutions with g
    
    for(int level = nbLevels - 1; level >= 0; level--)
    {
        a[level].setZero();
        
        int br = b[level + 1].rows() - 10, bc = b[level + 1].cols() - 10;
        {
            Map<FImage, 0, Stride<Dynamic, 2>> m(a[level].data(), br, bc, Stride<Dynamic, 2>(a[level].outerStride() * 2, 2));
            m = b[level + 1].block(5, 5, br, bc);
        }
        convolveSepInto<5, 5>(a[level], c[level], h12, h12t);
        // convolveInto<5, 5>(a[level], c[level], h12f);
        b[level] += c[level];
    }
    
    return b[0].block(5, 5, b[0].rows() - 10, b[0].cols() - 10);
}

#endif
