A modern C range for neighborhood pixel iteration

A modern C++ range for neighborhood pixel iteration with Insight Toolkit Niels Dekker LKEB Division of Image Processing Department of Radiology Leiden University Medical Center

Background • Research project at LUMC, Radiology + Ophthalmology • MRI (Magnetic Resonance Imaging) • Image segmentation of the eyeball • Using Insight Toolkit (ITK): Open source C++ library for image analysis

Edge enhancement by ITK’s Gaussian Derivative Input from MRI scanner Result of Gaussian Derivative

Gaussian derivative image function: a pixel neighborhood based operation x x y input image Array of kernel values output image

ITK implementation Inside itk: : Gaussian. Derivative. Image. Function: : Evaluate. At. Index(index) // “iterator” is an itk: : Const. Neighborhood. Iterator iterator. Set. Location(index); // index has N coordinates Neighborhood. Inner. Product (iterator, kernel. Values); Inside itk: : Neighborhood. Inner. Product: sum = 0. 0; for (unsigned i = 0; i < number. Of. Neighbors; ++i) sum += kernel. Values[i] * iterator. Get. Pixel(i);

ITK neighborhood iterator class hierarchy Neighborhood Const. Neighborhood. Iterator. With. Only. Index Neighborhood. Iterator private inheritance Const. Shaped. Neighborhood. Iterator

itk: : Const. Neighborhood. Iterator Image location iteration: • Set. Location(position) • Go. To. Begin(), Go. To. End() • Is. At. Begin(), Is. At. End() • operator++(), operator--() Neighbor pixel access: • Get. Pixel(i) // i < number. Of. Neighbors . . Neighborhood pixel pointers

Performance bottlenecks • Iterator: : Get. Pixel() was virtual (and not final) • Performance penalty 15% • Memory reallocated with each Iterator: : Set. Location(index) • Performance penalty 10% • Calculation neighborhood pixel pointers with each Set. Location • mutable data (caching “in bounds” property) makes multi-threading more difficult • No noexcept Please, none of these!!!

Neighborhood. Range: modern C++ style iteration! Neighborhood. Range: : begin() Neighborhood. Range: : end() Neighborhood. Range: : iterator it = range. begin(); pixel. Value x = *it; ++it;

Modern ways to get the inner product Standard C++ Library (<numeric>): sum = std: : inner_product(kernel. Values. begin(), kernel. Values. end(), neighborhood. Range. begin(), 0. 0); Range library (Eric Niebler): sum = ranges: : inner_product(kernel. Values, neighborhood. Range, 0. 0); Range-based for-loop: sum = 0. 0; for(auto neighbor. Pixel : neighborhood. Range) sum += kernel. Values[i++] * neighbor. Pixel; Manual iteration: it = neighborhood. Range. begin(); sum = 0. 0; for(const auto& kernel. Value: kernel. Values) { sum += kernel. Value * (*it); ++it; } Fastest!

Cheat sheet Iterator Examples ++it ; *it Input Iterator it 1 != it 2; v = *it; it++; it->mem; Output Iterator Forward Iterator Bidirectional Iterator Random Access Iterator Contiguous Iterator (C++17) *it = v; it++; istream_iterator ostream_iterator Default. Constructible Multi-pass guarantee forward_list --it; it--; list, set, map it + n; it[i]; it 1 < it 2; and eight more deque *(it+n) equals *(addressof(*it) + n) vector, array

Which iterator category or concept? • For std: : inner_product, Input. Iterator is OK • For range-based for, Input. Iterator + Output. Iterator • Forward. Iterator almost for free • Multi-pass: it 1 == it 2 implies ++it 1 == ++it 2 • To replace the old Iterator: : Get. Pixel(i): Random. Access. Iterator • Contiguous. Iterator is not feasible • neighborhood pixels may not be contiguously located in memory • Conclusion: Random. Access. Iterator

Neighborhood. Range<Image. Type> - initial version // Constructor center location Neighborhood. Range( Image. Type&, // (Possibly const) reference to N-dimensional image const Index. Type& center. Location, // Index of N integer coordinates const Radius. Size. Type& radius); // Radius value for each of the N image dimensions const_iterator cbegin() const; const_iterator cend() const; iterator begin() const; iterator end() const; const here means: the range does not “own” the pixels // Const or non-const qualified iterator template <bool Is. Const> class Qualified. Iterator; using const_iterator = Qualified. Iterator<true>; using iterator = Qualified. Iterator<false>; radius[0] radius[1]

Qualified. Iterator& Qualified. Iterator: : operator++() { for (unsigned i = 0; i < Num. Dimensions-1; ++i) { ++m_Current. Index[i]; if (m_Current. Index[i] <= m_Max. Index[i]) return *this; Iterator object has three N-D index data members: 1. m_Current. Index 2. m_Min. Index 3. m_Max. Index Min. Index m_Current. Index[i] = m_Min. Index[i]; } assert(m_Current. Index[Image. Dimension - 1] <= m_Max. Index[Image. Dimension - 1]); ++m_Current. Index[Image. Dimension - 1]; return *this; } Max. Index

Pixel access clamped within image borders reference operator*() const { Pixel. Type* ptr = m_Image. Buffer. Pointer; for (unsigned i = 0; i < Image. Dimension; ++i) // Sorry @Martin. Moene, cannot yet use C++17 std: : clamp! ptr += Clamp. Index. Value(m_Current. Index[i], m_Image. Size[i]) * m_Offset. Table[i]; return *ptr; } Me and Mr. Clamp (Martin Moene) @ C++Meetup, March 2016

Neigborhoods of any custom shape • Like the old (2003) itk: : Shaped. Neighborhood. Iterator • Triggered by ITK developer Pablo Hernandez-Cerdan (hyper-)rectangular shape + shape x shape arbitrary shape, arbitrary number of neighbors • Represented by offsets relative to neighborhood location + shape: {{0, -1}, {-1, 0}, {0, 1}} x shape: {{-1, -1}, {-1, 1}, {1, 1}}

A shaped neighborhood, specified by offsets Neighborhood. Range( Constructor Image. Type& image, const Index. Type& location, const Offset. Type* shape. Offsets, // Not owned or copied by the range! (C++20 std: : span<Offset. Type>? ) const size_t number. Of. Neigborhood. Pixels); // Equal to the number of offsets. template <typename TContainer. Offsets> Neighborhood. Range( Convenience overload Image. Type& image, const Index. Type& location, TContainer. Offsets&& shape. Offsets) : overprotective? Neighborhood. Range{ image, location, shape. Offsets. data(), shape. Offsets. size() } { static_assert( ! std: : is_rvalue_reference<decltype(shape. Offsets)>: : value, “No temporary please!”); }

Iterator increment became extremely simple! Qualified. Iterator& operator++() { // Just go to the next offset! ++m_Current. Offset; return *this; } Remember slide 14: Qualified. Iterator& Qualified. Iterator: : operator++() { for (unsigned i = 0; i < Num. Dimensions-1; ++i) { ++m_Current. Index[i]; if (m_Current. Index[i] <= m_Max. Index[i]) return *this; m_Current. Index[i] = m_Min. Index[i]; } assert(m_Current. Index[Image. Dimension - 1] <= m_Max. Index[Image. Dimension - 1]); ++m_Current. Index[Image. Dimension - 1]; return *this; }

ITK’s Image Accessor Bradley Lowekamp (ITK): “You are not using the Image Accessor”! “ITK image classes support various data layouts, internal and external pixel types [. . . ] Utilizing these to access the image is part of the fundamental design of ITK's image data. ” Example: itk: : Vector. Image image. Get. Neighborhood. Accessor() returns a Neighborhood. Accessor. Functor<Image. Type> object: template< typename TImage > class Neighborhood. Accessor. Functor { Pixel. Type Get(const Internal. Pixel. Type *pixel. Pointer) const; void Set(Internal. Pixel. Type * const pixel. Pointer, const Pixel. Type & p) const; };

class Pixel. Proxy // aims to act like a reference to a pixel { Internal. Pixel. Type* m_Internal. Pixel; Neighborhood. Accessor. Functor. Type m_Neighborhood. Accessor; public: …. // Converts proxy to pixel value (pixel. Value = proxy) operator Pixel. Type() const { return m_Neighborhood. Accessor. Get( m_Internal. Pixel ); } // Assigns pixel value to proxy (proxy = pixel. Value) Pixel. Proxy& operator=(const Pixel. Type& pixel. Value) // non-const, OK? { m_Neighborhood. Accessor. Set(m_Internal. Pixel, pixel. Value); return *this; } };

Using Pixel. Proxy as iterator: : reference type • • using reference = Pixel. Proxy; reference iterator: : operator*() const; But what about operator->() ? Use case (RGB color encoded pixels): it->r = red. Value; it->g = green. Value; it->b = blue. Value • However… auto operator->() const noexcept { return &*(*this); } warning C 4172: returning address of local variable or temporary

ACCU mailing list, April 19: A C++ iterator without operator->(), OK? Jonathan Wakely and Anthony Williams replied: • Support for it->mem is required, but • std: : istreambuf_iterator: : operator->() removed from Standard (LWG issue 2790) • Standard Library implementations does not really need to do it->mem • Return type of iterator: : operator*() must be T& • Not for std: : vector<bool>: : iterator, but “std: : vector<bool> is evil”.

Proxy Iterators – Eric Niebler • Proposal: Proxy Iterators for the Ranges Extensions (2015 - 2016) • The One Ranges Proposal (w/ Casey Carter and Christopher Di Bella) • Iterators++ blog, 2015 (Part 1, 2, 3) at http: //ericniebler. com Me: “Do you think I’m just lucky that [my proxy iterator] happens to work well? ” Eric: “I’ll say that you’re just getting lucky. Now, it’s more than just luck; STL implementors try to support proxy iterators, I’m told, because of bastard types like vector<bool>. But without the proper iter_swap and iter_move customization points I describe in my blog posts, the support can be half-baked at best. ”

Proxy iterators and range-based for-loops for( Pixel. Type x: range ) // OK (copies each pixel value) for( Pixel. Type& x: range ) // ERROR! for( const Pixel. Type& x: range ) // OK (read-only access) for( auto&& x: range )// OK: x is a proxy object. But still a warning: Clang warning: loop variable ‘x' is always a copy because the range does not return a reference [-Wrange-loop-analysis] for( auto x: range ) // Note: mutable pixel access (x is a proxy object) for( Neighborhood. Range: : reference x: range ) // More clear, maybe!

Accepted in ITK • Renamed the class as discussed with Matt Mc. Cormick (ITK): • itk: : Experimental: : Shaped. Image. Neighborhood. Range • April 25, 2018, after 20 revisions (15 days): accepted! • To be included with ITK 5. 0 (out soon)! • Possible future extensions: • Operator overloadings (+, -, +=, …) to ease use of pixel proxy • More utilities to generate shapes • Custom handling of pixel access outside image borders

Handling pixel access outside image borders ? • Currently clamping pixel access between image borders: ? • 1 1 1 |1 2 3 4 5 6 7|7 7 7 ? • Other possible approaches • • • Constant/zero-padding: Periodic (circular wrapping): Mirror-reflecting/symmetric: Linear extrapolation: Undefined behavior: 0 0 3 4 3 2 -2 -1 X X 0|1 5|1 1|1 0|1 X|1 2 2 2 3 3 3 4 4 4 5|0 5|1 5|5 5|6 5|X • Proposed: allow customizing pixel access by policy class: 0 2 4 7 X 0 3 3 8 X • Neighborhood. Range<class Image. Type, class Pixel. Access. Policy> • How to name the current “clamping” approach? Two proposals, “Border. Replicating. Pixel. Access. Policy” and “void”, not accepted, so far.

What’s in the name? (Mind the default!) ITK Image Boundary Conditions Open. CV Border. Types Mat. Lab imfilter Boundary Options Zero. Flux. Neumann. Boundary. Condition (default) BORDER_REPLICATE 'replicate' Constant. Boundary. Condition BORDER_CONSTANT X (default: X = 0) Periodic. Boundary. Condition BORDER_WRAP 'circular' Apply. Mirror. Boundary. Conditions (private) BORDER_REFLECT 'symmetric' BORDER_REFLECT 101 (default) BORDER_TRANSPARENT BORDER_ISOLATED

Conclusion • Shaped. Image. Neighborhood. Range to be released with ITK 5. 0 • Significant performance improvement • Modern C++ style iteration • Proxy iterators seem to work with std algorithms, but no guarantee! • Custom pixel access at borders still under discussion • Thanks to: • my colleagues at LKEB, LUMC • the guys from ITK (Kitware, Inc. ) • Ray Burgemeestre, Copernica, and all of you! talk. end()