PeriodicBoundaryCondition


(Tobias Wood) #1

Hello,

Is there an example of how to use a PeriodicBoundaryCondition anywhere?

I had hoped it was as simple as changing the definition of my NeighborhoodIterator to ConstNeighborhoodIterator<TImage, PeriodicBoundaryCondition<TImage, TImage>>;, however, this is causing a segmentation fault at line 57 of itkPeriodicBoundaryCondition.h, and the call to GetImagePointer is hitting a null pointer (which makes no sense to me). If I remove the template parameter so the default constant flux boundary condition is used, the code runs fine.

Thanks in advance.


(Dženan Zukić) #2

A short runnable example would kick-start debugging.


(Tobias Wood) #3

Hi,

I’ve attached an example that gives the segmentation fault (compiled with ITK 4.13).

Thanks!

example_periodic.cpp (1.1 KB)


(Tobias Wood) #4

Hi,

I’m wondering if anyone had a chance to look at the code I attached above?

Thanks


(Matt McCormick) #5

Hi Toby,

Thanks for sharing the nice example.

There were two issues to make it work.

  1. Call OverrideBoundaryCondition in the example.
  2. The use of auto has revealed a bug in itk::PeriodicBoundaryCondition.

Please review this patch:

http://review.source.kitware.com/#/c/23222/

And here is the working example code:

#include <iostream>

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

int main(int argc, char **argv) {
    itk::MultiThreaderBase::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 BoundaryConditionType = itk::PeriodicBoundaryCondition< ImageType, ImageType >;
    using IterType = itk::ConstNeighborhoodIterator<ImageType, BoundaryConditionType>;
    IterType::RadiusType radius; radius.Fill(1);
    IterType test_iter(radius, test_image, test_region);
    BoundaryConditionType boundaryCondition;
    test_iter.OverrideBoundaryCondition(&boundaryCondition);
    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;
}

CMakeLists.txt (224 Bytes)


(Tobias Wood) #6

Hi Matt,

Thanks for finding the bug. How detailed should code reviews like this be? The change looks good to me, but would you prefer that I check out that gerrit branch and confirm my own program works?

The only comment I really have is that it would be good to add an example to the docs for PeriodicBoundaryCondition, as I don’t think it’s intuitive to have to use OverrideBoundaryCondition. I’m happy for the above code to be used if necessary. What’s the best way to get that added to the docs?


(Dženan Zukić) #7

Make a pull request to SoftwareGuide and/or edit the doxygen comments in the header and then make a Gerrit patch.


(Matt McCormick) #8

@spinicist Yes, if you could please check out the patch and just verify that it works locally for you, that would be helpful.

Yes, the need to call OverrideBoundaryCondition is not intuitive. As @dzenanz mentioned, it would be nice to have this noted in the Iterators chapter of the Software Guide. We recently made it easier to contribute to the Software Guide on GitHub.


(Pablo Hernandez-Cerdan) #9

Why is OverrideBoundaryCondition needed if we are using the “default” templated BoundaryCondition type?


(Matt McCormick) #10

Yes, if we could make the OverrideBoundaryCondition call optional, that would be nice.


(Bradley Lowekamp) #11

The second template parameter of ConstNeighborhoodIterator is miss leading.

template<typename TImage, typename TBoundaryCondition = ZeroFluxNeumannBoundaryCondition< TImage >>
class itk::ConstNeighborhoodIterator< TImage, TBoundaryCondition >

It represents the default boundary condition that the neighborhood iterator is constructed ( and can be reset to ), and not the only boundary condition that can be used. The boundary condition can set at runtime via the OverrideBoundaryCondition method to another boundary type. This makes the class have a rather ugly mix of runtime, and compile time specification of the boundary condition. We MUST consider both of these use cases!


(Matt McCormick) #12

After investigation, it looks like calling OverrideBoundaryCondition should be / will be optional.


(Tobias Wood) #13

Hi,

I was very busy this week so have completely lost track of the multiple Gerrit threads. Is there something I can do to help move things along?


(Matt McCormick) #14

Hi Toby,

The issue has been resolved on ITK Git master.

Matt


(Tobias Wood) #15

Great, I will hopefully get round to testing it soon.