OpenVDB  7.0.0
MapsUtil.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
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 {
17 namespace OPENVDB_VERSION_NAME {
18 namespace util {
19 
20 // Utility methods for calculating bounding boxes
21 
25 template<typename MapType>
26 inline void
27 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  Vec3d corners[8];
34  corners[0] = in.min();;
35  corners[1] = Vec3d(min(0), min(1), min(2));
36  corners[2] = Vec3d(max(0), max(1), min(2));
37  corners[3] = Vec3d(min(0), max(1), min(2));
38  corners[4] = Vec3d(min(0), min(1), max(2));
39  corners[5] = Vec3d(max(0), min(1), max(2));
40  corners[6] = max;
41  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  out_min = map.applyInverseMap(corners[0]);
47  out_max = min;
48  for (int i = 1; i < 8; ++i) {
49  pre_image = map.applyInverseMap(corners[i]);
50  for (int j = 0; j < 3; ++j) {
51  out_min(j) = std::min( out_min(j), pre_image(j));
52  out_max(j) = std::max( out_max(j), pre_image(j));
53  }
54  }
55 }
56 
57 
60 template<typename MapType>
61 inline void
62 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 
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 
81  Vec3d center_pre_image = map.applyInverseMap(center);
82 
83  std::vector<Vec3d> coordinate_units;
84  coordinate_units.push_back(Vec3d(1,0,0));
85  coordinate_units.push_back(Vec3d(0,1,0));
86  coordinate_units.push_back(Vec3d(0,0,1));
87 
88  Vec3d& out_min = out.min();
89  Vec3d& out_max = out.max();
90  for (int direction = 0; direction < 3; ++direction) {
91  Vec3d temp = map.applyIJT(coordinate_units[direction]);
92  double offset =
93  radius * sqrt(temp.x()*temp.x() + temp.y()*temp.y() + temp.z()*temp.z());
94  out_min(direction) = center_pre_image(direction) - offset;
95  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 }
107 
108 
109 namespace { // anonymous namespace for this helper function
110 
117 inline int
118 findTangentPoints(const double g, const double xo, const double zo,
119  const double r, double& xp, double& zp, double& xm, double& zm)
120 {
121  double x2 = xo * xo;
122  double r2 = r * r;
123  double xd = g * xo;
124  double xd2 = xd*xd;
125  double zd = g * zo + 1.;
126  double zd2 = zd*zd;
127  double rd2 = r2*g*g;
128 
129  double distA = xd2 + zd2;
130  double distB = distA - rd2;
131 
132  if (distB > 0) {
133  double discriminate = sqrt(distB);
134 
135  xp = xo - xo*rd2/distA + r * zd *discriminate / distA;
136  xm = xo - xo*rd2/distA - r * zd *discriminate / distA;
137 
138  zp = (zo*zd2 + zd*g*(x2 - r2) - xo*xo*g - r*xd*discriminate) / distA;
139  zm = (zo*zd2 + zd*g*(x2 - r2) - xo*xo*g + r*xd*discriminate) / distA;
140 
141  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 
160 template<>
161 inline void
162 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  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  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  const double x_min = bbox.min().x();
204  const double y_min = bbox.min().y();
205  const double z_min = bbox.min().z();
206 
207  const double x_max = bbox.max().x();
208  const double y_max = bbox.max().y();
209  const double z_max = bbox.max().z();
210 
211  out_min.x() = x_min;
212  out_max.x() = x_max;
213  out_min.y() = y_min;
214  out_max.y() = y_max;
215 
216  Vec3d extreme;
217  Vec3d extreme2;
218  Vec3d pre_image;
219  // find the x-range
220  soln_number = findTangentPoints(gamma, centerLS.x(), centerLS.z(), radiusLS, xp, zp, xm, zm);
221  if (soln_number == 2) {
222  extreme.x() = xp;
223  extreme.y() = centerLS.y();
224  extreme.z() = zp;
225 
226  // location in world space of the tangent point
227  extreme2 = secondMap.applyMap(extreme);
228  // convert back to voxel space
229  pre_image = frustum.applyInverseMap(extreme2);
230  out_max.x() = std::max(x_min, std::min(x_max, pre_image.x()));
231 
232  extreme.x() = xm;
233  extreme.y() = centerLS.y();
234  extreme.z() = zm;
235  // location in world space of the tangent point
236  extreme2 = secondMap.applyMap(extreme);
237 
238  // convert back to voxel space
239  pre_image = frustum.applyInverseMap(extreme2);
240  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  soln_number = findTangentPoints(gamma, centerLS.y(), centerLS.z(), radiusLS, xp, zp, xm, zm);
250  if (soln_number == 2) {
251  extreme.x() = centerLS.x();
252  extreme.y() = xp;
253  extreme.z() = zp;
254 
255  // location in world space of the tangent point
256  extreme2 = secondMap.applyMap(extreme);
257  // convert back to voxel space
258  pre_image = frustum.applyInverseMap(extreme2);
259  out_max.y() = std::max(y_min, std::min(y_max, pre_image.y()));
260 
261  extreme.x() = centerLS.x();
262  extreme.y() = xm;
263  extreme.z() = zm;
264  extreme2 = secondMap.applyMap(extreme);
265 
266  // convert back to voxel space
267  pre_image = frustum.applyInverseMap(extreme2);
268  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  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  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  out_min.z() = std::max(z_min, frustum.applyInverseMap(secondMap.applyMap(near_point)).z());
286  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
Vec3d applyInverseMap(const Vec3d &in) const override
Return the pre-image of in under the map.
Definition: Maps.h:416
T & z()
Definition: Vec3.h:85
Map traits.
Definition: Maps.h:55
Vec3< double > Vec3d
Definition: Vec3.h:662
This map is composed of three steps. First it will take a box of size (Lx X Ly X Lz) defined by a mem...
Definition: Maps.h:1876
A general linear transform using homogeneous coordinates to perform rotation, scaling, shear and translation.
Definition: Maps.h:298
Vec3d voxelSize() const override
Return the lengths of the images of the segments (0,0,0)-(1,0,0), (0,0,0)-(0,1,0) and (0...
Definition: Maps.h:464
void calculateBounds(const MapType &map, const Vec3d &center, const Real radius, BBoxd &out)
Calculate an axis-aligned bounding box in the given map&#39;s domain from a spherical bounding box in its...
Definition: MapsUtil.h:62
const Vec3T & min() const
Return a const reference to the minimum point of this bounding box.
Definition: BBox.h:62
T & y()
Definition: Vec3.h:84
double Real
Definition: Types.h:37
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
Vec3d applyMap(const Vec3d &in) const override
Return the image of in under the map.
Definition: Maps.h:414
Definition: Exceptions.h:13
const Vec3T & max() const
Return a const reference to the maximum point of this bounding box.
Definition: BBox.h:64
T & x()
Reference to the component, e.g. v.x() = 4.5f;.
Definition: Vec3.h:83
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:102
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:106