Line data Source code
1 : #include "Intersection.h"
2 :
3 : namespace elsa
4 : {
5 : namespace detail
6 : {
7 : template <class data_t, int Dim>
8 : std::tuple<data_t, data_t, data_t, data_t>
9 : intersect(const BoundingBox& aabb, const Eigen::ParametrizedLine<data_t, Dim>& r)
10 243366 : {
11 243366 : data_t invDir = 1 / r.direction()(0);
12 :
13 243366 : data_t t1 = (aabb.min()(0) - r.origin()(0)) * invDir;
14 243366 : data_t t2 = (aabb.max()(0) - r.origin()(0)) * invDir;
15 :
16 243366 : data_t tmin = invDir >= 0 ? t1 : t2;
17 243366 : data_t tmax = invDir >= 0 ? t2 : t1;
18 243366 : const auto txmin = tmin;
19 243366 : const auto txmax = tmax;
20 :
21 500415 : for (int i = 1; i < aabb.min().rows(); ++i) {
22 257049 : invDir = 1 / r.direction()(i);
23 :
24 257049 : t1 = (aabb.min()(i) - r.origin()(i)) * invDir;
25 257049 : t2 = (aabb.max()(i) - r.origin()(i)) * invDir;
26 :
27 257049 : tmin = maxNum(tmin, invDir >= 0 ? t1 : t2);
28 257049 : tmax = minNum(tmax, invDir >= 0 ? t2 : t1);
29 257049 : }
30 :
31 243366 : return {tmin, tmax, txmin, txmax};
32 243366 : }
33 : } // namespace detail
34 :
35 : template <class data_t, int Dim>
36 : std::optional<IntersectionResult<data_t>>
37 : intersectRay(const BoundingBox& aabb, const Eigen::ParametrizedLine<data_t, Dim>& r)
38 236798 : {
39 236798 : const auto [tmin, tmax, txmin, txmax] = detail::intersect(aabb, r);
40 :
41 236798 : if (tmax == 0.0 && tmin == 0.0)
42 0 : return std::nullopt;
43 236798 : if (tmax >= maxNum(tmin, 0.0f)) // hit
44 218039 : return std::make_optional<IntersectionResult<data_t>>(tmin, tmax);
45 18759 : return std::nullopt;
46 18759 : }
47 :
48 : template <class data_t, int Dim>
49 : std::optional<IntersectionResult<data_t>>
50 : intersectXPlanes(BoundingBox aabb, const Eigen::ParametrizedLine<data_t, Dim>& r)
51 6856 : {
52 : // Slightly increase the size of the bounding box, such that center of voxels are actually
53 : // inside the bounding box, parallel rays which directly go through the center of the
54 : // voxel, will be recognized with this cheapo hack
55 6856 : aabb.min()[0] += 0.5f;
56 6856 : aabb.max()[0] -= 0.5f;
57 :
58 6856 : auto [tmin, tmax, txmin, txmax] = detail::intersect(aabb, r);
59 :
60 6856 : auto dist_to_integer = [](auto real) {
61 384 : const auto aabbmin = static_cast<data_t>(static_cast<int>(std::round(real)));
62 384 : return std::abs(real - aabbmin);
63 384 : };
64 :
65 6856 : auto fractional = [](auto real) {
66 384 : data_t trash = 0;
67 384 : return std::modf(real, &trash);
68 384 : };
69 :
70 6856 : auto advance_to_next_voxel = [&](auto aabb, auto& tmin) -> data_t {
71 : // the x-axis coord for voxel centers
72 304 : const auto xvoxelcenter = dist_to_integer(aabb.min()[0]);
73 :
74 304 : Eigen::Vector<data_t, Dim> entry = r.pointAt(tmin);
75 :
76 : // Calculate the distance from entry.x to the x coord of the next voxel
77 304 : const auto tmp = std::abs(fractional(entry[0] - xvoxelcenter));
78 304 : const auto frac = entry[0] < xvoxelcenter ? tmp : 1 - tmp;
79 :
80 : // Calculate distance in t, but don't advance if it's too small
81 304 : return tmp > 0.00001 && frac > 0.00001 ? frac / r.direction()[0] : 0;
82 304 : };
83 :
84 6856 : auto retreat_to_last_voxel = [&](auto aabb, auto& tmax) -> data_t {
85 80 : const auto xvoxelcenter = dist_to_integer(aabb.min()[0]);
86 :
87 80 : Eigen::Vector<data_t, Dim> exit = r.pointAt(tmax);
88 :
89 : // calculate the distance from entry.x to the x coord of the previous voxel
90 80 : const auto tmp = std::abs(fractional(exit[0] - xvoxelcenter));
91 80 : const auto frac = exit[0] < xvoxelcenter ? 1 - tmp : tmp;
92 :
93 : // Calculate distance in t, but don't advance if it's too small
94 80 : return tmp > 0.00001 && frac > 0.00001 ? frac / r.direction()[0] : 0;
95 80 : };
96 :
97 : // Sometimes we hit the y axis first, not the x axis, but we don't want
98 : // that, we always want the intersection with the x-plane
99 6856 : if (txmin < tmin) {
100 304 : tmin += advance_to_next_voxel(aabb, tmin);
101 304 : }
102 :
103 : // This can also happen if we leave at the top of bottom
104 6856 : if (txmax > tmax) {
105 80 : tmax -= retreat_to_last_voxel(aabb, tmax);
106 80 : }
107 :
108 6856 : if (tmax == 0.0 && tmin == 0.0)
109 0 : return std::nullopt;
110 :
111 6856 : if (tmax - maxNum(tmin, 0.0f) > -0.000001) // hit
112 6852 : return std::make_optional<IntersectionResult<data_t>>(tmin, tmax);
113 :
114 4 : return std::nullopt;
115 4 : }
116 :
117 : // ------------------------------------------
118 : // explicit template instantiation
119 : template struct IntersectionResult<float>;
120 : template struct IntersectionResult<double>;
121 :
122 : namespace detail
123 : {
124 : #define ELSA_INSTANTIATE_INTERSECT(type, dim) \
125 : template std::tuple<type, type, type, type> intersect<type, dim>( \
126 : const BoundingBox& aabb, const Eigen::ParametrizedLine<type, dim>& r);
127 :
128 : ELSA_INSTANTIATE_INTERSECT(float, 2)
129 : ELSA_INSTANTIATE_INTERSECT(double, 2)
130 : ELSA_INSTANTIATE_INTERSECT(float, 3)
131 : ELSA_INSTANTIATE_INTERSECT(double, 3)
132 : ELSA_INSTANTIATE_INTERSECT(float, Eigen::Dynamic)
133 : ELSA_INSTANTIATE_INTERSECT(double, Eigen::Dynamic)
134 : } // namespace detail
135 :
136 : #define ELSA_INSTANTIATE_INTERSECTRAY(type, dim) \
137 : template std::optional<IntersectionResult<type>> intersectRay<type, dim>( \
138 : const BoundingBox& aabb, const Eigen::ParametrizedLine<type, dim>& r);
139 :
140 : ELSA_INSTANTIATE_INTERSECTRAY(float, Eigen::Dynamic)
141 : ELSA_INSTANTIATE_INTERSECTRAY(double, Eigen::Dynamic)
142 :
143 : #define ELSA_INSTANTIATE_INTERSECTXPLANES(type, dim) \
144 : template std::optional<IntersectionResult<type>> intersectXPlanes<type, dim>( \
145 : BoundingBox aabb, const Eigen::ParametrizedLine<type, dim>& r);
146 :
147 : ELSA_INSTANTIATE_INTERSECTXPLANES(float, Eigen::Dynamic)
148 : ELSA_INSTANTIATE_INTERSECTXPLANES(double, Eigen::Dynamic)
149 :
150 : #undef ELSA_INSTANTIATE_INTERSECT
151 : #undef ELSA_INSTANTIATE_INTERSECTRAY
152 : #undef ELSA_INSTANTIATE_INTERSECTXPLANES
153 : } // namespace elsa
|
