PERF: N4BiasFieldCorrectionImageFilter implementation using ImageRange

[N-Dekker]

Adapted N4BiasFieldCorrectionImageFilter to start using ImageRange internally.
This commit yields a significant performance improvement, especially because
it replaces the (rather expensive) itk::Image::GetPixel(index) calls by
more efficient ImageRange[indexValue] calls.

Added checks to GenerateData() that the image sizes of the mask image and
confidence image (if provided) are equal to the size of the input image.

Ensured that the entire output image is computed at once, by overriding
EnlargeOutputRequestedRegion(DataObject*), as suggested by Bradley (@blowekamp).

Using VS2015 (Release), ~10% reduction of the run-time duration was observed on
a 3-D MRI volume + mask (453x340x20 voxels), provided by Denis (@dpshamonin),
at LKEB, Leiden University Medical Center.

[N-Dekker]

Please note: I will remove the "WIP" tag from the commit message when I think the commit is ready to be merged. For now, I still want to have another look at the code, and try it out!

[ntustison]

Hey, this is great. Hope it works out.

[N-Dekker]

Thanks for your encouragements @ntustison and @gdevenyi !

I just did some experiments (using MeVisLab), and observed a performance improvement of about 10%, with this commit. (Going from 77 sec to 70 sec., on a 453x340x20 MRI volume + mask, provided by my LKEB/LUMC colleague Denis, @dpshamonin ) It may not be a huge performance improvement, but it is significant, and reproducible.

[blowekamp]

This is a great addition! I see some change that need to happen and a couple things that could be improved...

[blowekamp]

This check should not be needed. This should be the "PrimaryInput" according to the ProcessObject, and it should by default be marked as required. This should be verified in the base class: https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Core/Common/src/itkProcessObject.cxx#L1419-L1429

[N-Dekker]

@blowekamp Thanks for your feedback. I'll remove itkExceptionMacro("Failed precondition: Input image should not be null") with the next amend.

[blowekamp]

The documentation for VerifyPreconditions says it is called before UpdateOutputInformation is called. The is a step in the ITK pipeline before the filters GenerateData step is done. So, the LargestPossibleRegion is not even updated at this point in the pipeline.

[N-Dekker]

@blowekamp Can you please explain? The maskImage is an input of the filter. So it should be checked before GenerateData. OK?

[blowekamp]

Similar here. Please review section 8.3.1 of the ITK Software guide explaining the step and how the ITK pipeline works by default. By default this filter's request output region is propagated to all the inputs in the substeps of PropagateRequestedRegion(). I'm not sure what the issue was to add these checks, but knowing that we can then work on a proper solution with the ITK pipeline.

[blowekamp]

There seems to be some regions issue with this filter related to the pipeline propagation of the requested region. We could ask the questions can the filter produce a sub-region of the output, or what input regions are needed to produce that output. The answers would indicate output enlarge the output requested regions and the propagate it to the input.

However, it looks like this filter assumes to produce the whole output and consumable the whole input. This can be done with the Pipeline by overloading EnlargeOutputRequestedRegion to set the output to the LargestPossibleRegion. Then the default implementation of GenerateInputRequestedRegion() will propagate that enlarged output to all the inputs, are require that they are present before this filter is executed.

I hope that makes since and helps...

[N-Dekker]

@blowekamp Thanks for your extensive analysis! So if I understand correctly, the question is: Was the N4 filter designed to support sub-regions? Or does it always operate on the entire input image and the entire output image?

If it always operates on the entire image buffer (its buffered region), 10% performance improvement can be achieved by this commit, using ImageRange. In this case EnlargeOutputRequestedRegion should still be called, to be sure that the output region is correct.

However, if ImageRange cannot be used according to this pull request, I would very much like to see an actual test case for this filter that produces a test failure on this pull request.

[blowekamp]

@N-Dekker I believe this filter most certainly needs all if its inputs that this filter can only compute the full bias field. Its possible to not apply the bias field to the whole image in the final expAndDivFilter step, but that would not be significant.

The EnlargeOutputRequestedRegion is already implicitly called in the pipeline. This filter need to overload the method and set it output to the largest possible region as done here:
https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Core/Common/include/itkImportImageFilter.hxx#L125-L142

After that you should not need these regions checks.

---

My questions what does the ImageRange operate on? is it always the entire Buffered image?

[N-Dekker]

My questions what does the ImageRange operate on? is it always the entire Buffered image?

Yes indeed, ImageRange always operates on the entire buffer. (I once considered renaming the class to something like "ImageBufferPixelRange", but so far "ImageRange" was found clear enough.)

I am considering to also propose a range to iterate over the pixels of a region (I think I will call it "ImageRegionRange"). But for this particular filter, it seems preferable to use the existing ImageRange. Because such an "ImageRegionRange" can never be as fast as the existing ImageRange.

[blowekamp]

Yes indeed, ImageRange always operates on the entire buffer. (I once considered renaming the class to something like "ImageBufferPixelRange", but so far "ImageRange" was found clear enough.)

OK... I think I am staring to under stand the requirements of the new ImageRange iterator and how they may need to interact with the ITK pipeline. The pipeline will guarantees that the requested regions is available in the GenerateData method, but it does not guarantee that the requested regions is the same as the buffered region.

So perhaps these checks are needed. The check would need to occur in the GenerateData method.

[N-Dekker]

@blowekamp

So perhaps these checks are needed. The check would need to occur in the GenerateData method.

With the current amend/version of the commit, the image sizes are checked within the added VerifyPreconditions() override. Isn't that sufficient?

Note that I also included your suggestion to override EnlargeOutputRequestedRegion.

[blowekamp]

With the current amend/version of the commit, the image sizes are checked within the added VerifyPreconditions() override. Isn't that sufficient?

The VerifyPreconditions methods is designed to execute early, before any pipeline methods are execute, so that a complex pipeline can fail very quickly if the parameters and input are set incorrectly. It is not expected to have updated meta-data or image data available.

The VerifyPreconditions method is not the right time for this check. The buffered region of an output is only updated/created in the GenerateData method, while VerifyInputInformation is called before even the PropagateRequestedRegion method ( please see section 8.3.1 of the ITK Software guide. Your implementation of VerifyPreconditions should be moved to the GenerateData method.

The N4 test has a bit of "style" to it. After every filter is constructed it is updated, and then the output is disconnected. It does not really exercise the ITK pipeline. You would see this issue if the filters inputs were connected from another filter ( not already updated ).

Note that I also included your suggestion to override EnlargeOutputRequestedRegion.
👍 Nicely done.

[N-Dekker]

@blowekamp OK: With the latest amend, the image size checks are in GenerateData. Hope it's good enough by now, because I'm going home now! Have a good day!

[blowekamp]

The discussion of the ImageRange iterator and what is needed for the iterator to be suitable and when it can safely be used in the ITK pipeline has moved to discourse.

[blowekamp]

Sorry wrong click...

[hjmjohnson]

@blowekamp It appears that your concerns have been addressed in the latest code. Would you mind reviewing your block on this PR?

[blowekamp]

From the discourse discussion, I don't believe we have come to an agreement on how to make the ImageRange more generally useful, how to make the assumption of the iterator safe, or when the iterating on the BufferedRegion is appropriate and safe. Perhaps we can discuss this today at the TCON.

[N-Dekker]

@blowekamp @hjmjohnson @thewtex In this specific case (the N4 bias filter), there's an alternative to this PR, that would yield the very same performance gain: just call itk::Image::GetBufferPointer() and access the pixels directly through the raw pointer. In practice that should be fine, for this particular filter. However, the ImageRange is more safe, as it "knows" its own begin and end, and has debug checks. And it is more generic, as it uses PixelAccessor when necessary (especially for itk::VectorImage).

[thewtex]

This is a great addition! Thank you @N-Dekker !

ImageRange is much better than raw pointer manipulation. While there may be other refactoring to refine ImageRange in the future, we could merge this as-is and make those changes in the future. The fundamental changes here are good.

Good first step, could use refactor in a subsequent PR to make the iteration Region explicit in MakeImageRange to ensure consistency between Ranges.

[N-Dekker]

@blowekamp wrote:

Good first step, could use refactor in a subsequent PR to make the iteration Region explicit in MakeImageRange to ensure consistency between Ranges.

It would be helpful to me if this PR could already be merged, indeed, as a starting point for further improvements.