#include <iostream>
#include <boost/array.hpp>
#include <boost/multi_array.hpp>

using namespace boost;
using namespace std;

const int ndim = 3;
typedef multi_array_ref<float, ndim> MArray3D;

int main() {
    const int nx = 512;
    const int ny = nx;
    const int nz = 51;
    const int nkx = 11;
    const int nky = nkx;
    const int nkz = nkx;
    float* resultdat = new float[nx*ny*nz];
    for (int i = 0; i < nx*ny*nz; i++) resultdat[i] = 0.f;
    float* imagedat = new float[nx*ny*nz];
    float* kerneldat = new float[nkx*nky*nkz];
    for (int i = 0; i < nkx*nky*nkz; i++) kerneldat[i] = 0.5f;
    array<size_t,ndim> dims = {{nx, ny, nz}};
    MArray3D image(imagedat, dims, fortran_storage_order());
    double sum = 0.;
    for (int iz = 0; iz < nz; iz++)
        for (int iy = 0; iy < ny; iy++)
            for (int ix = 0; ix < nx; ix++) {
                image[ix][iy][iz] = iz + 1 - nz/2;
                sum += image[ix][iy][iz];
            }
    cout << sum << endl;

    MArray3D result(resultdat, dims, fortran_storage_order());
    array<size_t,ndim> kdims = {{nkx, nky, nkz}};
    MArray3D kernel(kerneldat, kdims, fortran_storage_order());
    array<size_t,ndim> bases = {{-nkx/2, -nky/2, -nkz/2}};
    kernel.reindex(bases);
    // kernel corners, need to check for degenerate case
    int kxmin = -nkx/2; int kymin = -nky/2; int kzmin = -nkz/2;
    int kxmax = -kxmin, kymax = -kymin, kzmax = -kzmin;
    // interior boundaries, need to check for degenerate cases
    int izmin = 0, izmax = 0, iymin = 0, iymax = 0, ixmin = 0, ixmax = 0;
    if (1 != nz) {
        izmin = -kzmin;
        izmax = nz - 1 - kzmax;
    }
    if (1 != ny) {
        iymin = -kymin;
        iymax = ny - 1 - kymax;
    }
    if (1 != nx) {
        ixmin = -kxmin;
        ixmax = nx - 1 - kxmax;
    }
    // interior (no boundary condition issues here)
    float totsum = 0.f;
    for (int iz = izmin; iz <= izmax; iz++) {
        for (int iy = iymin; iy <= iymax; iy++) {
            for (int ix = ixmin; ix <= ixmax; ix++) {
                double sum = 0.f;
                for (int kz = kzmin; kz <= kzmax; kz++) {
                    for (int ky = kymin; ky <= kymax; ky++) {
                        for (int kx = kxmin; kx <= kxmax; kx++) {
                            float Kp = kernel[kx][ky][kz];
                            float fp = image[ix-kx][iy-ky][iz-kz];
                            sum += double(Kp)*fp;
                        }
                    }
                }
                result[ix][iy][iz] = float(sum);
                totsum += result[ix][iy][iz];
            }
        }
    }
    cout << totsum << endl;
    delete [] imagedat;
    delete [] kerneldat;
    delete [] resultdat;
}