#include <iostream>

using namespace std;



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* result = new float[nx*ny*nz];
    for (int i = 0; i < nx*ny*nz; i++) result[i] = 0.f;
    float* image = new float[nx*ny*nz];
    float* kernel = new float[nkx*nky*nkz];
    for (int i = 0; i < nkx*nky*nkz; i++) kernel[i] = 0.5f;
    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*ny)*nx] = iz + 1 - nz/2;
                sum += image[ix+(iy+iz*ny)*nx];
            }
    cout << sum << endl;

    // 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;
                int resindx = ix+(iy+iz*ny)*nx;
                for (int kz = kzmin; kz <= kzmax; kz++) {
                    for (int ky = kymin; ky <= kymax; ky++) {
                        for (int kx = kxmin; kx <= kxmax; kx++) {
                            int kindx = (kx-kxmin)+((ky-kymin)+(kz-kzmin)*nky)*nkx;
                            float Kp = kernel[kindx];
                            int newx = ix - kx, newy = iy - ky, newz = iz-kz;
                            int newindx = newx+(newy+newz*ny)*nx;
                            float fp = image[newindx];
                            sum += double(Kp)*fp;
                        }
                    }
                }
                result[resindx] = float(sum);
                totsum += result[resindx];
            }
        }
    }
    cout << totsum << endl;
    delete [] image;
    delete [] kernel;
    delete [] result;
}