GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/util/MapsUtil.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 90 96 93.8%
Functions: 4 5 80.0%
Branches: 20 36 55.6%
Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @file MapsUtil.h
5
6 #ifndef OPENVDB_UTIL_MAPSUTIL_HAS_BEEN_INCLUDED
7 #define OPENVDB_UTIL_MAPSUTIL_HAS_BEEN_INCLUDED
8
9 #include <openvdb/math/Maps.h>
10 #include <algorithm> // for std::min(), std::max()
11 #include <cmath>
12 #include <vector>
13
14
15 namespace openvdb {
16 OPENVDB_USE_VERSION_NAMESPACE
17 namespace OPENVDB_VERSION_NAME {
18 namespace util {
19
20 // Utility methods for calculating bounding boxes
21
22 /// @brief Calculate an axis-aligned bounding box in the given map's domain
23 /// (e.g., index space) from an axis-aligned bounding box in its range
24 /// (e.g., world space)
25 template<typename MapType>
26 inline void
27 4 calculateBounds(const MapType& map, const BBoxd& in, BBoxd& out)
28 {
29 const Vec3d& min = in.min();
30 const Vec3d& max = in.max();
31
32 // the pre-image of the 8 corners of the box
33
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
 36 Vec3d corners[8];
34 4 corners[0] = in.min();;
35 4 corners[1] = Vec3d(min(0), min(1), min(2));
36 4 corners[2] = Vec3d(max(0), max(1), min(2));
37 4 corners[3] = Vec3d(min(0), max(1), min(2));
38 4 corners[4] = Vec3d(min(0), min(1), max(2));
39 4 corners[5] = Vec3d(max(0), min(1), max(2));
40 4 corners[6] = max;
41 4 corners[7] = Vec3d(min(0), max(1), max(2));
42
43 Vec3d pre_image;
44 Vec3d& out_min = out.min();
45 Vec3d& out_max = out.max();
46 4 out_min = map.applyInverseMap(corners[0]);
47 4 out_max = min;
48
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 2 times.
 32 for (int i = 1; i < 8; ++i) {
49 28 pre_image = map.applyInverseMap(corners[i]);
50
2/2
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 14 times.
 112 for (int j = 0; j < 3; ++j) {
51
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 39 times.
 84 out_min(j) = std::min( out_min(j), pre_image(j));
52 84 out_max(j) = std::max( out_max(j), pre_image(j));
53 }
54 }
55 }
56
57
58 /// @brief Calculate an axis-aligned bounding box in the given map's domain
59 /// from a spherical bounding box in its range.
60 template<typename MapType>
61 inline void
62 4 calculateBounds(const MapType& map, const Vec3d& center, const Real radius, BBoxd& out)
63 {
64 // On return, out gives a bounding box in continuous index space
65 // that encloses the sphere.
66 //
67 // the image of a sphere under the inverse of the linearMap will be an ellipsoid.
68
69 if (math::is_linear<MapType>::value) {
70 // I want to find extrema for three functions f(x', y', z') = x', or = y', or = z'
71 // with the constraint that g = (x-xo)^2 + (y-yo)^2 + (z-zo)^2 = r^2.
72 // Where the point x,y,z is the image of x',y',z'
73 // Solve: \lambda Grad(g) = Grad(f) and g = r^2.
74 // Note: here (x,y,z) is the image of (x',y',z'), and the gradient
75 // is w.r.t the (') space.
76 //
77 // This can be solved exactly: e_a^T (x' -xo') =\pm r\sqrt(e_a^T J^(-1)J^(-T)e_a)
78 // where e_a is one of the three unit vectors. - djh.
79
80 /// find the image of the center of the sphere
81 Vec3d center_pre_image = map.applyInverseMap(center);
82
83 std::vector<Vec3d> coordinate_units;
84
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 coordinate_units.push_back(Vec3d(1,0,0));
85
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 coordinate_units.push_back(Vec3d(0,1,0));
86
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 4 coordinate_units.push_back(Vec3d(0,0,1));
87
88 Vec3d& out_min = out.min();
89 Vec3d& out_max = out.max();
90
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 4 times.
 16 for (int direction = 0; direction < 3; ++direction) {
91 12 Vec3d temp = map.applyIJT(coordinate_units[direction]);
92 12 double offset =
93 12 radius * sqrt(temp.x()*temp.x() + temp.y()*temp.y() + temp.z()*temp.z());
94 12 out_min(direction) = center_pre_image(direction) - offset;
95 12 out_max(direction) = center_pre_image(direction) + offset;
96 }
97
98 } else {
99 // This is some unknown map type. In this case, we form an axis-aligned
100 // bounding box for the sphere in world space and find the pre-images of
101 // the corners in index space. From these corners we compute an axis-aligned
102 // bounding box in index space.
103 BBoxd bounding_box(center - radius*Vec3d(1,1,1), center + radius*Vec3d(1,1,1));
104 calculateBounds<MapType>(map, bounding_box, out);
105 }
106 4 }
107
108
109 namespace { // anonymous namespace for this helper function
110
111 /// @brief Find the intersection of a line passing through the point
112 /// (<I>x</I>=0,&nbsp;<I>z</I>=&minus;1/<I>g</I>) with the circle
113 /// (<I>x</I> &minus; <I>xo</I>)&sup2; + (<I>z</I> &minus; <I>zo</I>)&sup2; = <I>r</I>&sup2;
114 /// at a point tangent to the circle.
115 /// @return 0 if the focal point (0, -1/<I>g</I>) is inside the circle,
116 /// 1 if the focal point touches the circle, or 2 when both points are found.
117 inline int
118 2 findTangentPoints(const double g, const double xo, const double zo,
119 const double r, double& xp, double& zp, double& xm, double& zm)
120 {
121 2 double x2 = xo * xo;
122 2 double r2 = r * r;
123 2 double xd = g * xo;
124 2 double xd2 = xd*xd;
125 2 double zd = g * zo + 1.;
126 2 double zd2 = zd*zd;
127 2 double rd2 = r2*g*g;
128
129 2 double distA = xd2 + zd2;
130 2 double distB = distA - rd2;
131
132
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (distB > 0) {
133 2 double discriminate = sqrt(distB);
134
135 2 xp = xo - xo*rd2/distA + r * zd *discriminate / distA;
136 2 xm = xo - xo*rd2/distA - r * zd *discriminate / distA;
137
138 2 zp = (zo*zd2 + zd*g*(x2 - r2) - xo*xo*g - r*xd*discriminate) / distA;
139 2 zm = (zo*zd2 + zd*g*(x2 - r2) - xo*xo*g + r*xd*discriminate) / distA;
140
141 2 return 2;
142
143 ✗ } if (0 >= distB && distB >= -1e-9) {
144 // the circle touches the focal point (x=0, z = -1/g)
145 ✗ xp = 0; xm = 0;
146 ✗ zp = -1/g; zm = -1/g;
147
148 ✗ return 1;
149 }
150
151 return 0;
152 }
153
154 } // end anonymous namespace
155
156
157 /// @brief Calculate an axis-aligned bounding box in index space
158 /// from a spherical bounding box in world space.
159 /// @note This specialization is optimized for a frustum map
160 template<>
161 inline void
162
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 calculateBounds<math::NonlinearFrustumMap>(const math::NonlinearFrustumMap& frustum,
163 const Vec3d& center, const Real radius, BBoxd& out)
164 {
165 // The frustum is a nonlinear map followed by a uniform scale, rotation, translation.
166 // First we invert the translation, rotation and scale to find the spherical pre-image
167 // of the sphere in "local" coordinates where the frustum is aligned with the near plane
168 // on the z=0 plane and the "camera" is located at (x=0, y=0, z=-1/g).
169
170 // check that the internal map has no shear.
171 const math::AffineMap& secondMap = frustum.secondMap();
172 // test if the linear part has shear or non-uniform scaling
173
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (!frustum.hasSimpleAffine()) {
174
175 // In this case, we form an axis-aligned bounding box for sphere in world space
176 // and find the pre_images of the corners in voxel space. From these corners we
177 // compute an axis-algined bounding box in voxel-spae
178 ✗ BBoxd bounding_box(center - radius*Vec3d(1,1,1), center + radius*Vec3d(1,1,1));
179 ✗ calculateBounds<math::NonlinearFrustumMap>(frustum, bounding_box, out);
180 return;
181 }
182
183 // for convenience
184 Vec3d& out_min = out.min();
185 Vec3d& out_max = out.max();
186
187 Vec3d centerLS = secondMap.applyInverseMap(center);
188 Vec3d voxelSize = secondMap.voxelSize();
189
190 // all the voxels have the same size since we know this is a simple affine map
191 1 double radiusLS = radius / voxelSize(0);
192
193 double gamma = frustum.getGamma();
194 double xp;
195 double zp;
196 double xm;
197 double zm;
198 int soln_number;
199
200 // the bounding box in index space for the points in the frustum
201 const BBoxd& bbox = frustum.getBBox();
202 // initialize min and max
203 1 const double x_min = bbox.min().x();
204 1 const double y_min = bbox.min().y();
205 1 const double z_min = bbox.min().z();
206
207 1 const double x_max = bbox.max().x();
208 1 const double y_max = bbox.max().y();
209 1 const double z_max = bbox.max().z();
210
211 1 out_min.x() = x_min;
212 1 out_max.x() = x_max;
213 1 out_min.y() = y_min;
214 1 out_max.y() = y_max;
215
216 Vec3d extreme;
217 Vec3d extreme2;
218 Vec3d pre_image;
219 // find the x-range
220 1 soln_number = findTangentPoints(gamma, centerLS.x(), centerLS.z(), radiusLS, xp, zp, xm, zm);
221
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (soln_number == 2) {
222 1 extreme.x() = xp;
223 1 extreme.y() = centerLS.y();
224 1 extreme.z() = zp;
225
226 // location in world space of the tangent point
227 1 extreme2 = secondMap.applyMap(extreme);
228 // convert back to voxel space
229 1 pre_image = frustum.applyInverseMap(extreme2);
230 1 out_max.x() = std::max(x_min, std::min(x_max, pre_image.x()));
231
232 1 extreme.x() = xm;
233 1 extreme.y() = centerLS.y();
234 1 extreme.z() = zm;
235 // location in world space of the tangent point
236 1 extreme2 = secondMap.applyMap(extreme);
237
238 // convert back to voxel space
239 1 pre_image = frustum.applyInverseMap(extreme2);
240 1 out_min.x() = std::max(x_min, std::min(x_max, pre_image.x()));
241
242 } else if (soln_number == 1) {
243 // the circle was tangent at the focal point
244 } else if (soln_number == 0) {
245 // the focal point was inside the circle
246 }
247
248 // find the y-range
249 1 soln_number = findTangentPoints(gamma, centerLS.y(), centerLS.z(), radiusLS, xp, zp, xm, zm);
250
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (soln_number == 2) {
251 1 extreme.x() = centerLS.x();
252 1 extreme.y() = xp;
253 1 extreme.z() = zp;
254
255 // location in world space of the tangent point
256 1 extreme2 = secondMap.applyMap(extreme);
257 // convert back to voxel space
258 1 pre_image = frustum.applyInverseMap(extreme2);
259 1 out_max.y() = std::max(y_min, std::min(y_max, pre_image.y()));
260
261 1 extreme.x() = centerLS.x();
262 1 extreme.y() = xm;
263 1 extreme.z() = zm;
264 1 extreme2 = secondMap.applyMap(extreme);
265
266 // convert back to voxel space
267 1 pre_image = frustum.applyInverseMap(extreme2);
268 1 out_min.y() = std::max(y_min, std::min(y_max, pre_image.y()));
269
270 } else if (soln_number == 1) {
271 // the circle was tangent at the focal point
272 } else if (soln_number == 0) {
273 // the focal point was inside the circle
274 }
275
276 // the near and far
277 // the closest point. The front of the frustum is at 0 in index space
278
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 double near_dist = std::max(centerLS.z() - radiusLS, 0.);
279 // the farthest point. The back of the frustum is at mDepth in index space
280
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 double far_dist = std::min(centerLS.z() + radiusLS, frustum.getDepth() );
281
282 Vec3d near_point(0.f, 0.f, near_dist);
283 Vec3d far_point(0.f, 0.f, far_dist);
284
285 2 out_min.z() = std::max(z_min, frustum.applyInverseMap(secondMap.applyMap(near_point)).z());
286 2 out_max.z() = std::min(z_max, frustum.applyInverseMap(secondMap.applyMap(far_point)).z());
287
288 }
289
290 } // namespace util
291 } // namespace OPENVDB_VERSION_NAME
292 } // namespace openvdb
293
294 #endif // OPENVDB_UTIL_MAPSUTIL_HAS_BEEN_INCLUDED
295