Line data    Source code

       1             : #include "TypeCasts.hpp"
       2             : #include "elsaDefines.h"
       3             : #include "Error.h"
       4             : 
       5             : #include <Eigen/Core>
       6             : #include <algorithm>
       7             : #include <vector>
       8             : 
       9             : #include "spdlog/fmt/fmt.h"
      10             : 
      11             : namespace elsa
      12             : {
      13             :     using IndexRange = std::pair<int, int>;
      14             : 
      15             :     /*
      16             :      * @brief This is a iterable view over a Cartesian index space. This defines a rectangular
      17             :      * region of integer indices. This is mostly equivalent to a n-dimensional loop, e.g.:
      18             :      *
      19             :      * ```cpp
      20             :      * for(int i = istart; i < istop; ++i) {
      21             :      *     for(int j = jstart; j < jstop; ++j) {
      22             :      *         // nested arbitrarily deep
      23             :      *     }
      24             :      * }
      25             :      * ```
      26             :      *
      27             :      * If you want to iterate over all indices form `0` to `upper` in a n-dimensional setting you
      28             :      * can:
      29             :      *
      30             :      * ```cpp
      31             :      * for(auto pos : CartesianIndex(upper)) {
      32             :      *     // nested arbitrarily deep
      33             :      * }
      34             :      * ```
      35             :      *
      36             :      * Note, that this is as typically an exclusive range for the end, i.e. [lower, upper).
      37             :      * At no point any coefficient of `pos` take the value of `upper` (lower[i] <= pos[i] < upper[i]
      38             :      * for all i in 0, ... pos.size()).
      39             :      *
      40             :      * Assuming `x` is a vector of integers. If you don't want to start at `0`. Pass the lower bound
      41             :      * index as the first argument:
      42             :      *
      43             :      * ```cpp
      44             :      * for(auto pos : CartesianIndex(lower, upper)) {
      45             :      *     // nested arbitrarily deep
      46             :      * }
      47             :      * ```
      48             :      *
      49             :      * The behaviour is quite similar to `std::views::iota` in C++20, but extended to n dimensions.
      50             :      * In this sense it is quite similar to Julies, check this
      51             :      * [blob](https://julialang.org/blog/2016/02/iteration/) for more information. This was also
      52             :      * the initial inspiration for this implementation.
      53             :      *
      54             :      */
      55             :     class CartesianIndices
      56             :     {
      57             :     private:
      58             :         index_t size_;
      59             :         IndexVector_t first_;
      60             :         IndexVector_t last_;
      61             : 
      62             :         // Tag struct, just to indicate the end iterator as sentinel. TODO: With a switch to C++
      63             :         // 20, remove this and let `end()` return a proper sentinel type
      64             :         struct as_sentinel {
      65             :         };
      66             : 
      67             :         static constexpr as_sentinel as_sentinel_t{};
      68             : 
      69             :     public:
      70             :         /// Vector of pairs, which represent the range for each dimension
      71             :         CartesianIndices(std::vector<IndexRange> ranges);
      72             : 
      73             :         /// Create an index space ranging from `0`, to `to` with the dimension of `to`
      74             :         CartesianIndices(const IndexVector_t& to);
      75             : 
      76             :         template <typename Vector>
      77           0 :         CartesianIndices(const Vector& to)
      78           0 :             : size_(to.size()), first_(IndexVector_t::Zero(size_)), last_(to.size())
      79             :         {
      80             :             using std::begin;
      81           0 :             std::copy_n(begin(to), size_, last_.begin());
      82           0 :         }
      83             : 
      84             :         /// Create an index space ranging from `from`, to `to` with the dimension of `ranges`
      85             :         CartesianIndices(const IndexVector_t& from, const IndexVector_t& to);
      86             : 
      87             :         template <typename Vector>
      88           0 :         CartesianIndices(const Vector& from, const Vector& to)
      89           0 :             : size_(to.size()), first_(from.size()), last_(to.size())
      90             :         {
      91             :             using std::begin;
      92           0 :             std::copy_n(begin(from), size_, first_.begin());
      93           0 :             std::copy_n(begin(to), size_, last_.begin());
      94           0 :         }
      95             : 
      96             :         /// Return the dimension of index space
      97             :         auto dims() const -> index_t;
      98             : 
      99             :         /// Return the number of coordinates of the space
     100             :         auto size() const -> index_t;
     101             : 
     102             :         /// return a range for a given dimension (0 <= i < dims())
     103             :         auto range(index_t i) const -> IndexRange;
     104             : 
     105             :         /// Return the lower bound of the index space
     106             :         auto first() -> IndexVector_t&;
     107             :         /// @overload
     108             :         auto first() const -> const IndexVector_t&;
     109             : 
     110             :         /// Return the upper bound of the index space
     111             :         auto last() -> IndexVector_t&;
     112             :         /// @overload
     113             :         auto last() const -> const IndexVector_t&;
     114             : 
     115             :         /// Random Access Iterator for index space
     116             :         struct iterator {
     117             :         private:
     118             :             IndexVector_t cur_;
     119             :             const IndexVector_t& begins_;
     120             :             const IndexVector_t& ends_;
     121             : 
     122             :         public:
     123             :             using value_type = IndexVector_t;
     124             :             using reference = value_type&;
     125             :             using pointer = value_type*;
     126             :             using difference_type = std::ptrdiff_t;
     127             :             using iterator_category = std::random_access_iterator_tag;
     128             : 
     129             :             iterator(as_sentinel, const IndexVector_t& begin, const IndexVector_t& end);
     130             : 
     131             :             explicit iterator(const IndexVector_t& vec);
     132             : 
     133             :             iterator(const IndexVector_t& cur, const IndexVector_t& begin,
     134             :                      const IndexVector_t& end);
     135             : 
     136             :             // Dereference
     137             :             const IndexVector_t& operator*() const;
     138             : 
     139             :             // Increment
     140             :             iterator& operator++();
     141             :             iterator operator++(int);
     142             : 
     143             :             // decrement
     144             :             iterator& operator--();
     145             :             iterator operator--(int);
     146             : 
     147             :             // advance
     148             :             auto operator+=(difference_type n) -> CartesianIndices::iterator&;
     149             :             auto operator-=(difference_type n) -> CartesianIndices::iterator&;
     150             : 
     151             :             friend auto operator+(const iterator& iter, difference_type n) -> iterator;
     152             :             friend auto operator+(difference_type n, const iterator& iter) -> iterator;
     153             : 
     154             :             friend auto operator-(const iterator& iter, difference_type n) -> iterator;
     155             :             friend auto operator-(const iterator& lhs, const iterator& rhs) -> difference_type;
     156             : 
     157             :             auto operator[](difference_type n) const -> value_type;
     158             : 
     159             :             // comparison
     160           0 :             friend auto operator==(const iterator& lhs, const iterator& rhs) -> bool
     161             :             {
     162           0 :                 return lhs.cur_ == rhs.cur_;
     163             :             }
     164             : 
     165           0 :             friend auto operator!=(const iterator& lhs, const iterator& rhs) -> bool
     166             :             {
     167           0 :                 return !(lhs == rhs);
     168             :             }
     169             : 
     170             :             // relational operators
     171             :             friend auto operator<(const iterator& lhs, const iterator& rhs) -> bool;
     172             :             friend auto operator>(const iterator& lhs, const iterator& rhs) -> bool;
     173             :             friend auto operator<=(const iterator& lhs, const iterator& rhs) -> bool;
     174             :             friend auto operator>=(const iterator& lhs, const iterator& rhs) -> bool;
     175             : 
     176             :         private:
     177             :             auto at_end() -> bool;
     178             :             auto at_end() const -> bool;
     179             : 
     180             :             // TODO: Prefer a for loop, and remove recursive function
     181             :             auto inc_recursive(index_t N) -> void;
     182             :             auto inc() -> void;
     183             : 
     184             :             auto prev_recursive(index_t N) -> void;
     185             :             auto prev() -> void;
     186             : 
     187             :             auto distance_to_recusrive(const iterator& iter, difference_type N) const
     188             :                 -> difference_type;
     189             :             auto distance_to(const iterator& iter) const -> difference_type;
     190             : 
     191             :             auto advance_recursive(difference_type n, index_t N) -> void;
     192             :             auto advance(difference_type n) -> void;
     193             :         };
     194             : 
     195             :         /// Return the begin iterator to the index space
     196             :         auto begin() -> iterator;
     197             :         /// Return the end iterator to the index space
     198             :         auto end() -> iterator;
     199             : 
     200             :         /// @overload
     201             :         auto begin() const -> iterator;
     202             :         /// @overload
     203             :         auto end() const -> iterator;
     204             :     };
     205             : 
     206             :     /// @brief Visit all neighbours with a certain distance `dist`, which have the same x-axis (i.e.
     207             :     /// they lie in the same y-z plane)
     208             :     ///
     209             :     /// The returned position can lie outside of the volume or have negative indices, so be careful
     210             :     CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist);
     211             : 
     212             :     /// @overload
     213             :     CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist);
     214             : 
     215             :     /// @brief Visit all neighbours with a certain distance `dist`, which have the same x-axis (i.e.
     216             :     /// they lie in the same y-z plane), plus make sure we stay in bounds of `lower`, and `upper`.
     217             :     ///
     218             :     /// Given a certain voxel in a volume, you can visit all close by voxels, which are still in the
     219             :     /// volume
     220             :     ///
     221             :     /// ```cpp
     222             :     /// auto volmin = {0, 0, 0};
     223             :     /// auto volmax = {128, 128, 128};
     224             :     /// auto current_pos = ...;
     225             :     /// for(auto neighbours : neighbours_in_slice(current_pos, 1, volmin, volmax)) {
     226             :     ///     // ...
     227             :     /// }
     228             :     /// ```
     229             :     /// Note, that the current_pos is also visited!
     230             :     CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist,
     231             :                                          const IndexVector_t& lower, const IndexVector_t& upper);
     232             : 
     233             :     /// @overload
     234             :     CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist, index_t leadingDim,
     235             :                                          const IndexVector_t& lower, const IndexVector_t& upper);
     236             : } // namespace elsa
     237             : 
     238             : template <>
     239             : struct fmt::formatter<elsa::CartesianIndices> {
     240             :     template <typename ParseContext>
     241           0 :     constexpr auto parse(ParseContext& ctx)
     242             :     {
     243           0 :         return ctx.begin();
     244             :     }
     245             : 
     246             :     template <typename FormatContext>
     247           0 :     auto format(const elsa::CartesianIndices& idx, FormatContext& ctx)
     248             :     {
     249           0 :         fmt::format_to(ctx.out(), "(");
     250             : 
     251           0 :         for (int i = 0; i < idx.dims() - 1; ++i) {
     252           0 :             auto p = idx.range(i);
     253           0 :             fmt::format_to(ctx.out(), "{}:{}, ", p.first, p.second);
     254             :         }
     255           0 :         auto p = idx.range(idx.dims() - 1);
     256           0 :         return fmt::format_to(ctx.out(), "{}:{})", p.first, p.second);
     257             :     }
     258             : };