Line data Source code
1 : #include "SliceTraversal.h"
2 : #include "Error.h"
3 : #include "elsaDefines.h"
4 : #include "spdlog/fmt/fmt.h"
5 : #include "spdlog/fmt/ostr.h"
6 :
7 : namespace elsa
8 : {
9 0 : RealMatrix_t TransformToTraversal::create2DRotation(const real_t leadingCoeff,
10 : const index_t leadingAxisIndex)
11 : {
12 0 : auto createRotationMatrix = [](real_t radian) {
13 0 : real_t c = std::cos(radian);
14 0 : real_t s = std::sin(radian);
15 0 : return RealMatrix_t({{c, -s}, {s, c}});
16 : };
17 :
18 : // TODO: Does a closed form solution exist?
19 : // Use conversion to polar coordinates
20 0 : if (leadingAxisIndex == 0 && leadingCoeff >= 0) {
21 : // Already identity, so do nothing
22 0 : return createRotationMatrix(0);
23 0 : } else if (leadingAxisIndex == 0 && leadingCoeff < 0) {
24 : // create a 2D 180° rotation matrix
25 0 : return createRotationMatrix(pi_t);
26 0 : } else if (leadingAxisIndex == 1 && leadingCoeff >= 0) {
27 : // create a 2D 270° rotation matrix counter clockwise
28 0 : return createRotationMatrix(3 * pi_t / 2);
29 0 : } else if (leadingAxisIndex == 1 && leadingCoeff <= 0) {
30 : // create a 2D 90° rotation matrix counter clockwise
31 0 : return createRotationMatrix(pi_t / 2);
32 : }
33 0 : return createRotationMatrix(0);
34 : }
35 :
36 0 : RealMatrix_t TransformToTraversal::create3DRotation(const real_t leadingCoeff,
37 : const index_t leadingAxisIndex)
38 : {
39 0 : auto identity = []() { return RealMatrix_t({{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}); };
40 :
41 : auto rotationX = [](real_t radian) {
42 : real_t c = std::cos(radian);
43 : real_t s = std::sin(radian);
44 : return RealMatrix_t({{1, 0, 0}, {0, c, -s}, {0, s, c}});
45 : };
46 :
47 0 : auto rotationY = [](real_t radian) {
48 0 : real_t c = std::cos(radian);
49 0 : real_t s = std::sin(radian);
50 0 : return RealMatrix_t({{c, 0, s}, {0, 1, 0}, {-s, 0, c}});
51 : };
52 :
53 0 : auto rotationZ = [](real_t radian) {
54 0 : real_t c = std::cos(radian);
55 0 : real_t s = std::sin(radian);
56 0 : return RealMatrix_t({{c, -s, 0}, {s, c, 0}, {0, 0, 1}});
57 : };
58 :
59 : // TODO: Does a closed form solution exist?
60 : // Use conversion to polar coordinates
61 0 : if (leadingAxisIndex == 0 && leadingCoeff >= 0) {
62 : // Already identity, so do nothing
63 0 : return identity();
64 0 : } else if (leadingAxisIndex == 0 && leadingCoeff < 0) {
65 0 : return rotationY(pi_t);
66 0 : } else if (leadingAxisIndex == 1 && leadingCoeff >= 0) {
67 0 : return rotationZ(-0.5 * pi_t);
68 0 : } else if (leadingAxisIndex == 1 && leadingCoeff <= 0) {
69 0 : return rotationZ(0.5 * pi_t);
70 0 : } else if (leadingAxisIndex == 2 && leadingCoeff >= 0) {
71 0 : return rotationY(0.5 * pi_t);
72 0 : } else if (leadingAxisIndex == 2 && leadingCoeff <= 0) {
73 0 : return rotationY(-0.5 * pi_t);
74 : }
75 :
76 0 : return identity();
77 : }
78 :
79 0 : RealMatrix_t TransformToTraversal::createRotation(RealRay_t ray)
80 : {
81 : // Get the leading axis, by absolute value
82 0 : index_t leadingAxisIndex = 0;
83 0 : ray.direction().array().abs().maxCoeff(&leadingAxisIndex);
84 :
85 : // Get the signed leading coefficient
86 0 : const auto leadingCoeff = ray.direction()(leadingAxisIndex);
87 :
88 0 : if (ray.dim() == 2) {
89 0 : return create2DRotation(leadingCoeff, leadingAxisIndex);
90 0 : } else if (ray.dim() == 3) {
91 0 : return create3DRotation(leadingCoeff, leadingAxisIndex);
92 : }
93 0 : throw Error("Can not create a {}-dimensional transformation for the slice traversal",
94 0 : ray.dim());
95 : }
96 :
97 0 : RealMatrix_t TransformToTraversal::createTransformation(const RealRay_t& ray,
98 : const RealVector_t& centerOfRotation)
99 : {
100 0 : const auto dim = ray.dim();
101 :
102 : // Setup translation
103 0 : RealMatrix_t translation(dim + 1, dim + 1);
104 0 : translation.setIdentity();
105 0 : translation.block(0, dim, dim, 1) = -centerOfRotation;
106 :
107 : // Setup rotation
108 0 : RealMatrix_t rotation(dim + 1, dim + 1);
109 0 : rotation.setIdentity();
110 0 : rotation.block(0, 0, dim, dim) = createRotation(ray);
111 :
112 0 : return rotation * translation;
113 0 : }
114 :
115 0 : TransformToTraversal::TransformToTraversal(const RealRay_t& ray,
116 0 : const RealVector_t& centerOfRotation)
117 0 : : transformation_(createTransformation(ray, centerOfRotation)),
118 0 : translation_(-centerOfRotation)
119 : {
120 0 : }
121 :
122 0 : RealRay_t TransformToTraversal::toTraversalCoordinates(RealRay_t ray) const
123 : {
124 0 : ray.origin() = *this * Point(ray.origin());
125 0 : ray.direction() = *this * Vec(ray.direction());
126 0 : return ray;
127 : }
128 :
129 0 : BoundingBox TransformToTraversal::toTraversalCoordinates(BoundingBox aabb) const
130 : {
131 : // Only translate, as the rotations are always 90, 180, 270 degrees,
132 : // TODO: this might be wrong, idk, what about non square things
133 0 : aabb._min = *this * Point(aabb._min);
134 0 : aabb._max = *this * Point(aabb._max);
135 0 : aabb.recomputeBounds();
136 0 : return aabb;
137 : }
138 :
139 0 : RealMatrix_t TransformToTraversal::transformation() const { return transformation_; }
140 :
141 0 : RealMatrix_t TransformToTraversal::invTransformation() const
142 : {
143 0 : return transformation().inverse();
144 : }
145 :
146 0 : RealMatrix_t TransformToTraversal::rotation() const
147 : {
148 0 : const auto dim = transformation_.rows();
149 0 : return transformation_.block(0, 0, dim - 1, dim - 1);
150 : }
151 :
152 0 : RealMatrix_t TransformToTraversal::invRotation() const { return rotation().transpose(); }
153 :
154 0 : RealVector_t TransformToTraversal::translation() const { return translation_; }
155 :
156 0 : RealMatrix_t TransformToTraversal::linear() const { return rotation(); }
157 :
158 0 : RealVector_t TransformToTraversal::operator*(const Point<real_t>& point) const
159 : {
160 0 : return (transformation() * point.point_.homogeneous()).hnormalized();
161 : }
162 :
163 0 : RealVector_t TransformToTraversal::operator*(const Vec<real_t>& vec) const
164 : {
165 0 : return linear() * vec.vec_;
166 : }
167 :
168 0 : SliceTraversal::SliceTraversal(BoundingBox aabb, RealRay_t ray)
169 0 : : transformation_(ray, aabb._max.array() / 2), ray_(ray)
170 : {
171 : // Transform ray and aabb to traversal coordinate space
172 0 : ray = transformation_.toTraversalCoordinates(ray);
173 0 : aabb = transformation_.toTraversalCoordinates(aabb);
174 :
175 : // TODO: We only want to shift in x direction by 0.5, not in y
176 0 : auto hit = Intersection::xPlanesWithRay(aabb, ray);
177 :
178 : // Default init is sufficient if we don't hit, TODO: test that!
179 0 : if (hit) {
180 0 : const auto firstVoxel = ray.pointAt(hit->_tmin);
181 0 : const auto lastVoxel = ray.pointAt(hit->_tmax);
182 :
183 0 : endIndex_ = std::max<index_t>(
184 0 : 1, static_cast<index_t>(std::ceil(lastVoxel[0] - firstVoxel[0] + 0.5f)));
185 :
186 0 : tDelta_ = 1 / ray.direction()[0];
187 :
188 0 : t_ = hit->_tmin;
189 0 : }
190 0 : }
191 :
192 0 : index_t SliceTraversal::leadingDirection() const
193 : {
194 0 : index_t idx = 0;
195 0 : ray_.direction().array().cwiseAbs().maxCoeff(&idx);
196 0 : return idx;
197 : }
198 :
199 0 : real_t SliceTraversal::t() const { return t_; }
200 :
201 0 : real_t SliceTraversal::tDelta() const { return tDelta_; }
202 :
203 0 : SliceTraversal::Iter SliceTraversal::begin() const
204 : {
205 0 : return {startIndex_, ray_.pointAt(t_), ray_.direction() * tDelta_};
206 : }
207 :
208 0 : SliceTraversal::Iter SliceTraversal::end() const { return {endIndex_}; }
209 :
210 0 : index_t SliceTraversal::startIndex() const { return startIndex_; }
211 :
212 0 : index_t SliceTraversal::endIndex() const { return endIndex_; }
213 :
214 0 : SliceTraversal::Iter::value_type SliceTraversal::Iter::operator*() const
215 : {
216 0 : return cur_.template cast<index_t>();
217 : }
218 :
219 0 : SliceTraversal::Iter& SliceTraversal::Iter::operator++()
220 : {
221 0 : ++pos_;
222 0 : cur_ += dir_;
223 0 : return *this;
224 : }
225 :
226 0 : SliceTraversal::Iter SliceTraversal::Iter::operator++(int)
227 : {
228 0 : auto copy = *this;
229 0 : ++(*this);
230 0 : return copy;
231 0 : }
232 :
233 0 : bool operator==(const SliceTraversal::Iter& lhs, const SliceTraversal::Iter& rhs)
234 : {
235 0 : return lhs.pos_ == rhs.pos_;
236 : }
237 :
238 0 : bool operator!=(const SliceTraversal::Iter& lhs, const SliceTraversal::Iter& rhs)
239 : {
240 0 : return !(lhs == rhs);
241 : }
242 :
243 : } // namespace elsa
|