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