#include <iostream>

#include "itkImage.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkPeriodicBoundaryCondition.h"

int main(int argc, char **argv) {
    itk::MultiThreader::SetGlobalMaximumNumberOfThreads(4);

    typedef itk::Image<float, 3> ImageType;

    auto test_image = ImageType::New();
    ImageType::RegionType test_region;
    test_region.GetModifiableSize() = {{64, 64, 64}};
    test_image->SetRegions(test_region);
    test_image->Allocate();

    using IterType = itk::ConstNeighborhoodIterator<ImageType, itk::PeriodicBoundaryCondition<ImageType, ImageType>>;
    IterType::RadiusType radius; radius.Fill(1);
    IterType test_iter(radius, test_image, test_region);
    std::vector<IterType::OffsetType> back = { {{-1, 0, 0}}, {{ 0,-1, 0}}, {{ 0, 0,-1}} };
    test_iter.SetNeedToUseBoundaryCondition(true);
    test_iter.GoToBegin();
    while (!test_iter.IsAtEnd()) {
        float sum = 0;
        for (auto j = 0; j < back.size(); j++) {
            const float d = test_iter.GetPixel(back[j]);
            sum += d*d;
        }
        ++test_iter;
    }
    std::cout << "FINISHED" << std::endl;
    return EXIT_SUCCESS;
}