Line data Source code
1 : #pragma once
2 :
3 : #include "Intersection.h"
4 : #include "elsaDefines.h"
5 : #include "BoundingBox.h"
6 : #include "Logger.h"
7 : #include "Assertions.h"
8 :
9 : #include <Eigen/Geometry>
10 : #include <cassert>
11 : #include <optional>
12 : #include <cmath>
13 : #include <utility>
14 : #include <variant>
15 :
16 : // TODO: remove
17 : #include <iostream>
18 :
19 : namespace elsa
20 : {
21 : /// Strong type to distinguish transformation of points vs transformation of vectors
22 : template <class data_t>
23 : struct Point {
24 0 : explicit Point(Vector_t<data_t> point) : point_(std::move(point)) {}
25 :
26 : Vector_t<data_t> point_;
27 : };
28 :
29 : /// Strong type to distinguish transformation of points vs transformation of vectors
30 : template <class data_t>
31 : struct Vec {
32 0 : explicit Vec(Vector_t<data_t> vec) : vec_(std::move(vec)) {}
33 :
34 : Vector_t<data_t> vec_;
35 : };
36 :
37 : /**
38 : * @brief Represent a transformation, which will transform any point into the
39 : * traversal coordinate system.
40 : *
41 : * In the traversal coordinate system, the reference direction - form which it is constructed -
42 : * is transformed, such that the leading coefficient (i.e. the coefficient of largest absolute
43 : * value), is in the positive x-direction (i.e. the first component)
44 : *
45 : * The algorithm to determine the rotation in a 2D case is quite simple: Given the vector
46 : * \f$\begin{bmatrix} x & y \end{bmatrix}\f$,
47 : * first determine the leading axis, by finding the index of the component with maximum
48 : * absolute value, i.e. \f$\text{axis} = \text{maxAxis} (\lvert x \rvert, \lvert y \rvert) \f$
49 : * (in this case, the max returns "x", or "y" for our use cases instead of the actual value).
50 : * Assume `axis` is either `"x"` or `"y"` and `leadingCoeff` stores corresponding value
51 : * to the leading axis (but it's signed!). Then in pseudocode the algorithm determines the
52 : * rotation the following way:
53 : *
54 : * ```python
55 : * if axis == "x" and maxCoeff >= 0:
56 : * return rotate_by(0)
57 : * elif axis == "x" and maxCoeff <= 0:
58 : * return rotate_by(180)
59 : * elif axis == "y" and maxCoeff >= 0:
60 : * return rotate_by(90)
61 : * elif axis == "y" and maxCoeff <= 0:
62 : * return rotate_by(270)
63 : * ```
64 : *
65 : * To extent it to 3D one further decision has to be made, do we only rotate, such that the
66 : * leading direction is in the x direction, or do we rotate, such that y is the second largest
67 : * direction and z is the smallest.
68 : *
69 : * TODO: It might be nice to move this to a separate file and module, to wrap some of Eigens
70 : * transformation stuff, such that it's usable for dynamic cases, as we have here. This
71 : * might be nice to have for the Geometry class as well, but not for now.
72 : */
73 : class TransformToTraversal
74 : {
75 : private:
76 : static RealMatrix_t create2DRotation(const real_t leadingCoeff,
77 : const index_t leadingAxisIndex);
78 :
79 : static RealMatrix_t create3DRotation(const real_t leadingCoeff,
80 : const index_t leadingAxisIndex);
81 :
82 : static RealMatrix_t createRotation(RealRay_t ray);
83 :
84 : static RealMatrix_t createTransformation(const RealRay_t& ray,
85 : const RealVector_t& centerOfRotation);
86 :
87 : RealMatrix_t transformation_;
88 : RealVector_t translation_;
89 :
90 : public:
91 : TransformToTraversal(const RealRay_t& ray, const RealVector_t& centerOfRotation);
92 :
93 : RealRay_t toTraversalCoordinates(RealRay_t ray) const;
94 :
95 : BoundingBox toTraversalCoordinates(BoundingBox aabb) const;
96 :
97 : RealMatrix_t transformation() const;
98 :
99 : RealMatrix_t invTransformation() const;
100 :
101 : RealMatrix_t rotation() const;
102 :
103 : RealMatrix_t invRotation() const;
104 :
105 : RealVector_t translation() const;
106 :
107 : RealMatrix_t linear() const;
108 :
109 : RealVector_t operator*(const Point<real_t>& point) const;
110 :
111 : RealVector_t operator*(const Vec<real_t>& vec) const;
112 : };
113 :
114 : /**
115 : * @brief Traverse a volume along the direction of a given ray. Each step it will advance one
116 : * slice further in the direction of the ray. The traversal visits voxels at the center planes,
117 : * i.e. it always evaluates at the center of voxel in the plane of leading direction.
118 : *
119 : * This is a slightly modified version of the algorithm of Amanatides & Woo's "A Fast Voxel
120 : * Traversal Algorithm". The reference visits each and every single voxel along the ray.
121 : * However, for this case we don't need that.
122 : *
123 : * Given a bounding box of the volume and a ray to traverse the volume, one can simply call:
124 : *
125 : * ```cpp
126 : * for(auto [pos, voxel, t] = SliceTraversal(aabb, ray)) {
127 : * // pos: exact position on center plane of voxel
128 : * // voxel: current voxel
129 : * // t: position on the ray (i.e. pos = ray.origin() + t * ray.direction())
130 : * }
131 : * ```
132 : *
133 : * The number of visited voxels is known at construction, therefore moving along the ray is
134 : * basically decrementing a count. However, to return useful information, more bookkeeping is
135 : * needed.
136 : *
137 : * Dereferencing the iterator will return a small struct, which has information about the
138 : * current position in the volume, the current voxel and the t parameter of the ray.
139 : *
140 : * Note: When using voxel returned from iterating over the volume, it can happen that voxel
141 : * on the boundary of the volume are returned. This can trigger wrong behaviour for certain
142 : * applications and therefore should be handled with care.
143 : *
144 : * A note about the implementation: To ease the construction and handling, the incoming ray and
145 : * bounding box are transformed to a traversal internal coordinate space, in which the leading
146 : * direction lies in the positive x-axis. The 'world' space is partitioned in 4 quadrants based
147 : * on the direction. The first ranges for all direction from with a polar angle in the range of
148 : * [-45°, 45°], the second (45°, 135°), the third [135°, 225°], and the fourth (225°, 315°).
149 : * Note that -45° = 315°.
150 : *
151 : * TODO:
152 : * - Make the iterator random access
153 : *
154 : * @author
155 : * - David Frank - initial code
156 : */
157 : class SliceTraversal
158 : {
159 : private:
160 : TransformToTraversal transformation_;
161 :
162 : /// TODO: Do we want to consider points that spawn inside the volume?
163 : index_t startIndex_{0};
164 :
165 : index_t endIndex_{0};
166 :
167 : /// t value for the first intersection of the ray and the bounding box
168 : real_t t_{0.0};
169 :
170 : /// t value to increment each iteration
171 : real_t tDelta_{0.0};
172 :
173 : RealRay_t ray_;
174 :
175 : public:
176 : /// Traversal iterator, models forward iterator, maybe this should actually be an input
177 : /// iterator due to the non reference type of the dereference type. IDK, as we use it, this
178 : /// works, but in the future this might should be different.
179 : struct Iter {
180 : public:
181 : using iterator_category = std::forward_iterator_tag;
182 : using difference_type = std::ptrdiff_t;
183 : using value_type = IndexVector_t;
184 : using pointer = value_type*;
185 : using reference = value_type&;
186 :
187 : private:
188 : index_t pos_{};
189 : RealVector_t cur_{};
190 : RealVector_t dir_{};
191 :
192 : public:
193 : /// Construct iterator
194 : ///
195 : /// @param pos index position of the traversal, used mainly for comparing two iterators
196 : /// @param ray traversed ray used to compute exact position on dereference
197 : /// @param deltat increment of t each increment
198 : /// @param t position along the ray
199 0 : Iter(index_t pos, const RealVector_t& entry, const RealVector_t& dir)
200 0 : : pos_(pos), cur_(entry), dir_(dir)
201 : {
202 0 : }
203 :
204 0 : Iter(index_t pos) : pos_(pos) {}
205 :
206 : /// Dereference iterator
207 : value_type operator*() const;
208 :
209 : /// Pre increment
210 : Iter& operator++();
211 :
212 : /// Post increment
213 : Iter operator++(int);
214 :
215 : /// Equal to operator with iterator
216 : friend bool operator==(const Iter& lhs, const Iter& rhs);
217 :
218 : /// Not equal to operator with iterator
219 : friend bool operator!=(const Iter& lhs, const Iter& rhs);
220 : };
221 :
222 : /// Delete default construction
223 : SliceTraversal() = delete;
224 :
225 : /// Construct traversal from bounding box and ray
226 : SliceTraversal(BoundingBox aabb, RealRay_t ray);
227 :
228 : /// Get the first visited voxel
229 : Iter begin() const;
230 :
231 : /// Get on past the end
232 : Iter end() const;
233 :
234 : /// Get the index of the first visited voxel
235 : index_t startIndex() const;
236 :
237 : /// Get the index of the one past the last visited voxel
238 : index_t endIndex() const;
239 :
240 : /// Get the leading direction of the ray
241 : index_t leadingDirection() const;
242 :
243 : /// Get the initial computed t value for the entry point
244 : real_t t() const;
245 :
246 : /// Get the computed delta t value
247 : real_t tDelta() const;
248 : };
249 : } // namespace elsa
|
