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