OpenVDB  0.104.0
MapsUtil.h
Go to the documentation of this file.
1 
2 //
3 // Copyright (c) 2012 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 //
32 
33 #ifndef OPENVDB_UTIL_MAPSUTIL_HAS_BEEN_INCLUDED
34 #define OPENVDB_UTIL_MAPSUTIL_HAS_BEEN_INCLUDED
35 
36 #include <openvdb/math/Maps.h>
37 
38 
39 namespace openvdb {
41 namespace OPENVDB_VERSION_NAME {
42 namespace util {
43 
44 // Utility methods for calculating bounding boxes
45 
49 template<typename MapType>
50 inline void
51 calculateBounds(const MapType& map, const BBoxd& in, BBoxd& out)
52 {
53  const Vec3d& min = in.min();
54  const Vec3d& max = in.max();
55 
56  // the pre-image of the 8 corners of the box
57  Vec3d corners[8];
58  corners[0] = in.min();;
59  corners[1] = Vec3d(min(0), min(1), min(2));
60  corners[2] = Vec3d(max(0), max(1), min(2));
61  corners[3] = Vec3d(min(0), max(1), min(2));
62  corners[4] = Vec3d(min(0), min(1), max(2));
63  corners[5] = Vec3d(max(0), min(1), max(2));
64  corners[6] = max;
65  corners[7] = Vec3d(min(0), max(1), max(2));
66 
67  Vec3d pre_image;
68  Vec3d& out_min = out.min();
69  Vec3d& out_max = out.max();
70  out_min = map.applyInverseMap(corners[0]);
71  out_max = min;
72  for (int i = 1; i < 8; ++i) {
73  pre_image = map.applyInverseMap(corners[i]);
74  for (int j = 0; j < 3; ++j) {
75  out_min(j) = std::min( out_min(j), pre_image(j));
76  out_max(j) = std::max( out_max(j), pre_image(j));
77  }
78  }
79 }
80 
81 
84 template<typename MapType>
85 inline void
86 calculateBounds(const MapType& map, const Vec3d& center, const Real radius, BBoxd& out)
87 {
88  // On return, out gives a bounding box in continuous index space
89  // that encloses the sphere.
90  //
91  // the image of a sphere under the inverse of the linearMap will be an ellipsoid.
92 
94  // I want to find extrema for three functions f(x', y', z') = x', or = y', or = z'
95  // with the constraint that g = (x-xo)^2 + (y-yo)^2 + (z-zo)^2 = r^2.
96  // Where the point x,y,z is the image of x',y',z'
97  // Solve: \lambda Grad(g) = Grad(f) and g = r^2.
98  // Note: here (x,y,z) is the image of (x',y',z'), and the gradient
99  // is w.r.t the (') space.
100  //
101  // This can be solved exactly: e_a^T (x' -xo') =\pm r\sqrt(e_a^T J^(-1)J^(-T)e_a)
102  // where e_a is one of the three unit vectors. - djh.
103 
105  Vec3d center_pre_image = map.applyInverseMap(center);
106 
107  std::vector<Vec3d> coordinate_units;
108  coordinate_units.push_back(Vec3d(1,0,0));
109  coordinate_units.push_back(Vec3d(0,1,0));
110  coordinate_units.push_back(Vec3d(0,0,1));
111 
112  Vec3d& out_min = out.min();
113  Vec3d& out_max = out.max();
114  for (int direction = 0; direction < 3; ++direction) {
115  Vec3d temp = map.applyIJT(coordinate_units[direction]);
116  double offset =
117  radius * sqrt(temp.x()*temp.x() + temp.y()*temp.y() + temp.z()*temp.z());
118  out_min(direction) = center_pre_image(direction) - offset;
119  out_max(direction) = center_pre_image(direction) + offset;
120  }
121 
122  } else {
123  // This is some unknown map type. In this case, we form an axis-aligned
124  // bounding box for the sphere in world space and find the pre-images of
125  // the corners in index space. From these corners we compute an axis-aligned
126  // bounding box in index space.
127  BBoxd bounding_box(center - radius*Vec3d(1,1,1), center + radius*Vec3d(1,1,1));
128  calculateBounds<MapType>(map, bounding_box, out);
129  }
130 }
131 
132 
133 namespace { // anonymous namespace for this helper function
134 
140 inline int
141 findTangentPoints(const double g, const double xo, const double zo,
142  const double r, double& xp, double& zp, double& xm, double& zm)
143 {
144  double x2 = xo * xo;
145  double r2 = r * r;
146  double xd = g * xo;
147  double xd2 = xd*xd;
148  double zd = g * zo + 1.;
149  double zd2 = zd*zd;
150  double rd2 = r2*g*g;
151 
152  double distA = xd2 + zd2;
153  double distB = distA - rd2;
154 
155  if (distB > 0) {
156  double discriminate = sqrt(distB);
157 
158  xp = xo - xo*rd2/distA + r * zd *discriminate / distA;
159  xm = xo - xo*rd2/distA - r * zd *discriminate / distA;
160 
161  zp = (zo*zd2 + zd*g*(x2 - r2) - xo*xo*g - r*xd*discriminate) / distA;
162  zm = (zo*zd2 + zd*g*(x2 - r2) - xo*xo*g + r*xd*discriminate) / distA;
163 
164  return 2;
165 
166  } if (0 >= distB && distB >= -1e-9) {
167  // the circle touches the focal point (x=0, z = -1/g)
168  xp = 0; xm = 0;
169  zp = -1/g; zm = -1/g;
170 
171  return 1;
172  }
173 
174  return 0;
175 }
176 
177 } // end anonymous namespace
178 
179 
183 template<>
184 inline void
185 calculateBounds<math::NonlinearFrustumMap>(const math::NonlinearFrustumMap& map,
186  const Vec3d& center, const Real radius, BBoxd& out)
187 {
188  // The frustum is a nonlinear map followed by a uniform scale, rotation, translation.
189  // First we invert the translation, rotation and scale to find the spherical pre-image
190  // of the sphere in "local" coordinates where the frustum is aligned with the near plane
191  // on the z=0 plane and the "camera" is located at (x=0, y=0, z=-1/g).
192 
193  // check that the internal map has no shear.
194  const math::AffineMap& secondMap = map.secondMap();
195  // test if the linear part has shear or non-uniform scaling
196  if (!map.hasSimpleAffine()) {
197  // In this case, we form an axis-aligned bounding box for sphere in world space
198  // and find the pre_images of the corners in voxel space. From these corners we
199  // compute an axis-algined bounding box in voxel-spae
200  BBoxd bounding_box(center - radius*Vec3d(1,1,1), center + radius*Vec3d(1,1,1));
201  calculateBounds<math::NonlinearFrustumMap>(map, bounding_box, out);
202  return;
203  }
204 
205  // for convenience
206  Vec3d& out_min = out.min();
207  Vec3d& out_max = out.max();
208 
209  Vec3d centerLS = secondMap.applyInverseMap(center);
210  Vec3d voxelSize = secondMap.voxelSize();
211  // all the voxels have the same size since we know this is a simple affine map
212  double radiusLS = radius/ voxelSize(0);
213 
214  double gamma = map.getGamma();
215  double xp;
216  double zp;
217  double xm;
218  double zm;
219  int soln_number;
220 
221  // the bounding box in index space for the points in the frustum
222  const BBoxd& bbox = map.getBBox();
223  // initialize min and max
224  const double x_min = bbox.min().x();
225  const double y_min = bbox.min().y();
226  const double z_min = bbox.min().z();
227 
228  const double x_max = bbox.max().x();
229  const double y_max = bbox.max().y();
230  const double z_max = bbox.max().z();
231 
232  out_min.x() = x_min;
233  out_max.x() = x_max;
234  out_min.y() = y_min;
235  out_max.y() = y_max;
236 
237  Vec3d extreme;
238  Vec3d extreme2;
239  Vec3d pre_image;
240  // find the x-range
241  soln_number = findTangentPoints(gamma, centerLS.x(), centerLS.z(), radiusLS, xp, zp, xm, zm);
242  if (soln_number == 2) {
243  extreme.x() = xp;
244  extreme.y() = centerLS.y();
245  extreme.z() = zp;
246 
247  extreme2 = secondMap.applyMap(extreme);
248  // back into voxel space
249  pre_image = map.applyInverseMap(extreme2);
250  out_max.x() = std::max(x_min, std::min(x_max, pre_image.x()));
251 
252  extreme.x() = xm;
253  extreme.y() = centerLS.y();
254  extreme.z() = zm;
255  extreme2 = secondMap.applyMap(extreme);
256  pre_image = map.applyInverseMap(extreme2);
257  out_min.x() = std::max(x_min, std::min(x_max, pre_image.x()));
258  } else if (soln_number == 1) {
259  // the circle was tangent at the focal point
260  } else if (soln_number == 0) {
261  // the focal point was inside the circle
262  }
263 
264  // find the y-range
265  soln_number = findTangentPoints(gamma, centerLS.y(), centerLS.z(), radiusLS, xp, zp, xm, zm);
266  if (soln_number == 2) {
267  extreme.x() = centerLS.x();
268  extreme.y() = xp;
269  extreme.z() = zp;
270 
271  extreme2 = secondMap.applyMap(extreme);
272 
273  pre_image = map.applyInverseMap(extreme2);
274  out_max.y() = std::max(y_min, std::min(y_max, pre_image.y()));
275 
276  extreme.x() = centerLS.x();
277  extreme.y() = xm;
278  extreme.z() = zm;
279  extreme2 = secondMap.applyMap(extreme);
280  pre_image = map.applyInverseMap(extreme);
281  out_min.y() = std::max(y_min, std::min(y_max, pre_image.y()));
282 
283  } else if (soln_number == 1) {
284  // the circle was tangent at the focal point
285  } else if (soln_number == 0) {
286  // the focal point was inside the circle
287  }
288 
289  // the near and far
290  // the closest point. The front of the frustum is at 0
291  double near_dist = std::max(centerLS.z() - radiusLS, 0.);
292  // the farthest point. The back of the frustum is at mDepth
293  double far_dist = std::min(centerLS.z() + radiusLS, map.getDepth() );
294 
295  Vec3d near_point(0, 0, near_dist);
296  Vec3d far_point(0, 0, far_dist);
297 
298  out_min.z() = std::max(z_min, map.applyInverseMap(secondMap.applyMap(near_point)).z());
299  out_max.z() = std::min(z_max, map.applyInverseMap(secondMap.applyMap(far_point)).z());
300 }
301 
302 } // namespace util
303 } // namespace OPENVDB_VERSION_NAME
304 } // namespace openvdb
305 
306 #endif // OPENVDB_UTIL_MAPSUTIL_HAS_BEEN_INCLUDED
307 
308 // Copyright (c) 2012 DreamWorks Animation LLC
309 // All rights reserved. This software is distributed under the
310 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )