Commit 0facef56 authored by David Frank's avatar David Frank
Browse files

Add view for cartesian index space

Using this CartesianIndices, you can create a n-dimensional recetengular
shape and iterate over the region. This should be quite a powerful but
easy to use abstraction over the volume.

The provided free functions, can be used to visit neighbouring voxels
from a given position.
parent 1232fb05
Pipeline #843818 passed with stages
in 32 minutes and 40 seconds
...@@ -4,6 +4,7 @@ set(MODULE_HEADERS ...@@ -4,6 +4,7 @@ set(MODULE_HEADERS
Backtrace.h Backtrace.h
Cloneable.h Cloneable.h
Utilities/Badge.hpp Utilities/Badge.hpp
Utilities/CartesianIndices.h
Utilities/DataContainerFormatter.hpp Utilities/DataContainerFormatter.hpp
Utilities/FormatConfig.h Utilities/FormatConfig.h
Utilities/Math.hpp Utilities/Math.hpp
...@@ -35,6 +36,7 @@ set(MODULE_HEADERS ...@@ -35,6 +36,7 @@ set(MODULE_HEADERS
set(MODULE_SOURCES set(MODULE_SOURCES
elsaDefines.cpp elsaDefines.cpp
Backtrace.cpp Backtrace.cpp
Utilities/CartesianIndices.cpp
Descriptors/DataDescriptor.cpp Descriptors/DataDescriptor.cpp
Descriptors/VolumeDescriptor.cpp Descriptors/VolumeDescriptor.cpp
Descriptors/PlanarDetectorDescriptor.cpp Descriptors/PlanarDetectorDescriptor.cpp
......
#include "CartesianIndices.h"
#include "TypeCasts.hpp"
#include <iostream>
#include <iterator>
#include <numeric>
namespace elsa
{
template <typename Fn>
auto map(std::vector<IndexRange> v, Fn fn)
{
IndexVector_t result(v.size());
std::transform(v.begin(), v.end(), result.begin(), std::move(fn));
return result;
}
auto CartesianIndices::dims() const -> index_t { return as<index_t>(idxrange_.size()); }
auto CartesianIndices::size() const -> index_t
{
// TODO: Once the coverage CI image is updated, change this to:
// clang-format off
// std::transform_reduce(idxrange_.begin(), idxrange_.end(), 1, std::multiplies<>{},
// [](auto p) { return p.second - p.first; }));
// clang-format on
std::vector<index_t> tmp;
tmp.reserve(idxrange_.size());
std::transform(idxrange_.begin(), idxrange_.end(), std::back_inserter(tmp),
[](auto p) { return p.second - p.first; });
return std::accumulate(tmp.begin(), tmp.end(), 1, std::multiplies<>{});
}
auto CartesianIndices::range(index_t i) const -> IndexRange { return idxrange_[asUnsigned(i)]; }
auto CartesianIndices::first() -> IndexVector_t
{
return map(idxrange_, [](auto range) { return range.first; });
}
auto CartesianIndices::first() const -> IndexVector_t
{
return map(idxrange_, [](auto range) { return range.first; });
}
auto CartesianIndices::last() -> IndexVector_t
{
return map(idxrange_, [](auto range) { return range.second; });
}
auto CartesianIndices::last() const -> IndexVector_t
{
return map(idxrange_, [](auto range) { return range.second; });
}
auto CartesianIndices::begin() -> iterator { return iterator(first(), first(), last()); }
auto CartesianIndices::begin() const -> iterator { return {first(), first(), last()}; }
auto CartesianIndices::end() -> iterator { return {as_sentinel_t, first(), last()}; }
auto CartesianIndices::end() const -> iterator { return {as_sentinel_t, first(), last()}; }
/****** Iterator implementation ******/
CartesianIndices::iterator::iterator(as_sentinel, const IndexVector_t& begin,
const IndexVector_t& end)
: cur_(end), begins_(end), ends_(end)
{
cur_.tail(cur_.size() - 1) = begin.tail(cur_.size() - 1);
}
CartesianIndices::iterator::iterator(IndexVector_t vec) : cur_(vec), begins_(vec), ends_(vec) {}
CartesianIndices::iterator::iterator(IndexVector_t cur, IndexVector_t begin, IndexVector_t end)
: cur_(cur), begins_(begin), ends_(end)
{
}
auto CartesianIndices::iterator::operator*() const -> IndexVector_t { return cur_; }
auto CartesianIndices::iterator::operator++() -> iterator&
{
inc();
return *this;
}
auto CartesianIndices::iterator::operator++(int) -> iterator
{
auto copy = *this;
++(*this);
return copy;
}
auto CartesianIndices::iterator::operator--() -> iterator&
{
prev();
return *this;
}
auto CartesianIndices::iterator::operator--(int) -> iterator
{
auto copy = *this;
--(*this);
return copy;
}
auto CartesianIndices::iterator::operator+=(difference_type n) -> iterator&
{
advance(n);
return *this;
}
auto CartesianIndices::iterator::operator-=(difference_type n) -> iterator&
{
advance(-n);
return *this;
}
auto operator+(const CartesianIndices::iterator& iter,
CartesianIndices::iterator::difference_type n) -> CartesianIndices::iterator
{
auto copy = iter;
copy += n;
return copy;
}
auto operator+(CartesianIndices::iterator::difference_type n,
const CartesianIndices::iterator& iter) -> CartesianIndices::iterator
{
return iter + n;
}
auto operator-(const CartesianIndices::iterator& iter,
CartesianIndices::iterator::difference_type n) -> CartesianIndices::iterator
{
auto copy = iter;
copy -= n;
return copy;
}
auto operator-(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
-> CartesianIndices::iterator::difference_type
{
return rhs.distance_to(lhs);
}
auto CartesianIndices::iterator::operator[](difference_type n) const -> value_type
{
// TODO: Make this more efficient
auto copy = *this;
copy += n;
return *copy;
}
auto operator<(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
-> bool
{
// Not sure why I have to instantiate them this way, else it doesn't work...
using Array = Eigen::Array<index_t, Eigen::Dynamic, 1>;
const Array a1 = lhs.cur_.array();
const Array a2 = rhs.cur_.array();
return (a1 < a2).any();
}
auto operator>(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
-> bool
{
return rhs < lhs;
}
auto operator<=(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
-> bool
{
return !(lhs > rhs);
}
auto operator>=(const CartesianIndices::iterator& lhs, const CartesianIndices::iterator& rhs)
-> bool
{
return !(lhs < rhs);
}
auto CartesianIndices::iterator::at_end() -> bool { return cur_[0] == ends_[0]; }
auto CartesianIndices::iterator::at_end() const -> bool { return cur_[0] == ends_[0]; }
auto CartesianIndices::iterator::inc_recursive(index_t N) -> void
{
auto& iter = cur_[N];
if (++iter == ends_[N] && N > 0) {
iter = begins_[N];
inc_recursive(N - 1);
}
}
auto CartesianIndices::iterator::inc() -> void { inc_recursive(cur_.size() - 1); }
auto CartesianIndices::iterator::prev_recursive(index_t N) -> void
{
auto& iter = cur_[N];
if (iter == begins_[N] && N > 0) {
iter = ends_[N];
inc_recursive(N - 1);
}
--iter;
}
auto CartesianIndices::iterator::prev() -> void { prev_recursive(cur_.size() - 1); }
auto CartesianIndices::iterator::distance_to_recusrive(const iterator& other,
difference_type N) const
-> difference_type
{
if (N == 0) {
return other.cur_[N] - cur_[N];
} else {
const auto d = distance_to_recusrive(other, N - 1);
const auto scale = ends_[N] - begins_[N];
const auto increment = other.cur_[N] - cur_[N];
return difference_type{(d * scale) + increment};
}
}
auto CartesianIndices::iterator::distance_to(const iterator& other) const -> difference_type
{
auto N = cur_.size() - 1;
return distance_to_recusrive(other, N);
}
auto CartesianIndices::iterator::advance_recursive(difference_type n, index_t N) -> void
{
if (n == 0) {
return;
}
auto& iter = cur_[N];
const auto size = ends_[N] - begins_[N];
const auto first = begins_[N];
auto const idx = as<difference_type>(iter - first);
n += idx;
auto div = size ? n / size : 0;
auto mod = size ? n % size : 0;
if (N != 0) {
if (mod < 0) {
mod += size;
div--;
}
advance_recursive(div, N - 1);
} else {
if (div > 0) {
mod = size;
}
}
iter = first + mod;
}
auto CartesianIndices::iterator::advance(difference_type n) -> void
{
advance_recursive(n, cur_.size() - 1);
}
CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist)
{
Eigen::IOFormat format(4, 0, ", ", "", "", "", "[", "]");
IndexVector_t from = pos;
from.tail(pos.size() - 1).array() -= dist.array();
IndexVector_t to = pos;
to.tail(pos.size() - 1).array() += dist.array() + 1;
to[0] += 1; // Be sure this is incremented, so we actually iterate over it
return {from, to};
}
CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist)
{
return neighbours_in_slice(pos, IndexVector_t::Constant(pos.size() - 1, dist));
}
CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist,
const IndexVector_t& lower, const IndexVector_t& upper)
{
IndexVector_t from = pos;
from.tail(pos.size() - 1).array() -= dist.array();
from = from.cwiseMax(lower);
IndexVector_t to = pos;
to.tail(pos.size() - 1).array() += dist.array() + 1;
to[0] += 1;
to = to.cwiseMin(upper);
return {from, to};
}
CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist,
const IndexVector_t& lower, const IndexVector_t& upper)
{
return neighbours_in_slice(pos, IndexVector_t::Constant(pos.size() - 1, dist), lower,
upper);
}
} // namespace elsa
#include "TypeCasts.hpp"
#include "elsaDefines.h"
#include "Error.h"
#include <Eigen/Core>
#include <algorithm>
#include <vector>
#include "spdlog/fmt/fmt.h"
namespace elsa
{
using IndexRange = std::pair<int, int>;
/*
* @brief This is a iterable view over a Cartesian index space. This defines a rectangular
* region of integer indices. This is mostly equivalent to a n-dimensional loop, e.g.:
*
* ```cpp
* for(int i = istart; i < istop; ++i) {
* for(int j = jstart; j < jstop; ++j) {
* // nested arbitrarily deep
* }
* }
* ```
*
* If you want to iterate over all indices form `0` to `upper` in a n-dimensional setting you
* can:
*
* ```cpp
* for(auto pos : CartesianIndex(upper)) {
* // nested arbitrarily deep
* }
* ```
*
* Note, that this is as typically an exclusive range for the end, i.e. [lower, upper).
* At no point any coefficient of `pos` take the value of `upper` (lower[i] <= pos[i] < upper[i]
* for all i in 0, ... pos.size()).
*
* Assuming `x` is a vector of integers. If you don't want to start at `0`. Pass the lower bound
* index as the first argument:
*
* ```cpp
* for(auto pos : CartesianIndex(lower, upper)) {
* // nested arbitrarily deep
* }
* ```
*
* The behaviour is quite similar to `std::views::iota` in C++20, but extended to n dimensions.
* In this sense it is quite similar to Julies, check this
* [blob](https://julialang.org/blog/2016/02/iteration/) for more information. This was also
* the initial inspiration for this implementation.
*
*/
class CartesianIndices
{
private:
std::vector<IndexRange> idxrange_{};
// Tag struct, just to indicate the end iterator as sentinel. TODO: With a switch to C++
// 20, remove this and let `end()` return a proper sentinel type
struct as_sentinel {
};
static constexpr as_sentinel as_sentinel_t{};
public:
/// Vector of pairs, which represent the range for each dimension
CartesianIndices(std::vector<IndexRange> ranges);
/// Create an index space ranging from `0`, to `to` with the dimension of `to`
template <typename Vector>
// requires ForwardRange + Sized
CartesianIndices(Vector to)
{
idxrange_.reserve(to.size());
using std::begin;
using std::end;
std::transform(begin(to), end(to), std::back_inserter(idxrange_), [](auto idx) {
return std::pair{0, idx};
});
}
/// Create an index space ranging from `from`, to `to` with the dimension of `ranges`
template <typename Vector>
// requires ForwardRange + Sized
CartesianIndices(Vector from, Vector to)
{
if (from.size() != to.size()) {
throw InvalidArgumentError("CartesianIndices: vectors must be of same size");
}
idxrange_.reserve(asUnsigned(from.size()));
using std::begin;
using std::end;
std::transform(begin(from), end(from), begin(to), std::back_inserter(idxrange_),
[](auto start, auto end) {
return std::pair{start, end};
});
}
/// Return the dimension of index space
auto dims() const -> index_t;
/// Return the number of coordinates of the space
auto size() const -> index_t;
/// return a range for a given dimension (0 <= i < dims())
auto range(index_t i) const -> IndexRange;
/// Return the lower bound of the index space
auto first() -> IndexVector_t;
/// @overload
auto first() const -> IndexVector_t;
/// Return the upper bound of the index space
auto last() -> IndexVector_t;
/// @overload
auto last() const -> IndexVector_t;
/// Random Access Iterator for index space
struct iterator {
private:
IndexVector_t cur_;
IndexVector_t begins_;
IndexVector_t ends_;
public:
using value_type = IndexVector_t;
using reference = value_type&;
using pointer = value_type*;
using difference_type = std::ptrdiff_t;
using iterator_category = std::random_access_iterator_tag;
iterator(as_sentinel, const IndexVector_t& begin, const IndexVector_t& end);
explicit iterator(IndexVector_t vec);
iterator(IndexVector_t cur, IndexVector_t begin, IndexVector_t end);
// Dereference
IndexVector_t operator*() const;
// Increment
iterator& operator++();
iterator operator++(int);
// decrement
iterator& operator--();
iterator operator--(int);
// advance
auto operator+=(difference_type n) -> CartesianIndices::iterator&;
auto operator-=(difference_type n) -> CartesianIndices::iterator&;
friend auto operator+(const iterator& iter, difference_type n) -> iterator;
friend auto operator+(difference_type n, const iterator& iter) -> iterator;
friend auto operator-(const iterator& iter, difference_type n) -> iterator;
friend auto operator-(const iterator& lhs, const iterator& rhs) -> difference_type;
auto operator[](difference_type n) const -> value_type;
// comparison
friend auto operator==(const iterator& lhs, const iterator& rhs) -> bool
{
return lhs.cur_ == rhs.cur_;
}
friend auto operator!=(const iterator& lhs, const iterator& rhs) -> bool
{
return !(lhs == rhs);
}
// relational operators
friend auto operator<(const iterator& lhs, const iterator& rhs) -> bool;
friend auto operator>(const iterator& lhs, const iterator& rhs) -> bool;
friend auto operator<=(const iterator& lhs, const iterator& rhs) -> bool;
friend auto operator>=(const iterator& lhs, const iterator& rhs) -> bool;
private:
auto at_end() -> bool;
auto at_end() const -> bool;
// TODO: Prefer a for loop, and remove recursive function
auto inc_recursive(index_t N) -> void;
auto inc() -> void;
auto prev_recursive(index_t N) -> void;
auto prev() -> void;
auto distance_to_recusrive(const iterator& iter, difference_type N) const
-> difference_type;
auto distance_to(const iterator& iter) const -> difference_type;
auto advance_recursive(difference_type n, index_t N) -> void;
auto advance(difference_type n) -> void;
};
/// Return the begin iterator to the index space
auto begin() -> iterator;
/// Return the end iterator to the index space
auto end() -> iterator;
/// @overload
auto begin() const -> iterator;
/// @overload
auto end() const -> iterator;
};
/// @brief Visit all neighbours with a certain distance `dist`, which have the same x-axis (i.e.
/// they lie in the same y-z plane)
///
/// The returned position can lie outside of the volume or have negative indices, so be careful
CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist);
/// @overload
CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist);
/// @brief Visit all neighbours with a certain distance `dist`, which have the same x-axis (i.e.
/// they lie in the same y-z plane), plus make sure we stay in bounds of `lower`, and `upper`.
///
/// Given a certain voxel in a volume, you can visit all close by voxels, which are still in the
/// volume
///
/// ```cpp
/// auto volmin = {0, 0, 0};
/// auto volmax = {128, 128, 128};
/// auto current_pos = ...;
/// for(auto neighbours : neighbours_in_slice(current_pos, 1, volmin, volmax)) {
/// // ...
/// }
/// ```
/// Note, that the current_pos is also visited!
CartesianIndices neighbours_in_slice(const IndexVector_t& pos, const IndexVector_t& dist,
const IndexVector_t& lower, const IndexVector_t& upper);
/// @overload
CartesianIndices neighbours_in_slice(const IndexVector_t& pos, index_t dist,
const IndexVector_t& lower, const IndexVector_t& upper);
} // namespace elsa
template <>
struct fmt::formatter<elsa::CartesianIndices> {
template <typename ParseContext>
constexpr auto parse(ParseContext& ctx)
{
return ctx.begin();
}
template <typename FormatContext>
auto format(const elsa::CartesianIndices& idx, FormatContext& ctx)
{
fmt::format_to(ctx.out(), "(");
for (int i = 0; i < idx.dims() - 1; ++i) {
auto p = idx.range(i);
fmt::format_to(ctx.out(), "{}:{}, ", p.first, p.second);
}
auto p = idx.range(idx.dims() - 1);
return fmt::format_to(ctx.out(), "{}:{})", p.first, p.second);
}
};
...@@ -27,3 +27,4 @@ ELSA_DOCTEST(DataHandlers) ...@@ -27,3 +27,4 @@ ELSA_DOCTEST(DataHandlers)
ELSA_DOCTEST(DataHandlerMap) ELSA_DOCTEST(DataHandlerMap)
ELSA_DOCTEST(DataContainer)