HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MapsUtil.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2018 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 //
29 ///////////////////////////////////////////////////////////////////////////
30 //
31 /// @file MapsUtil.h
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 #include <algorithm> // for std::min(), std::max()
38 #include <cmath>
39 #include <vector>
40 
41 
42 namespace openvdb {
44 namespace OPENVDB_VERSION_NAME {
45 namespace util {
46 
47 // Utility methods for calculating bounding boxes
48 
49 /// @brief Calculate an axis-aligned bounding box in the given map's domain
50 /// (e.g., index space) from an axis-aligned bounding box in its range
51 /// (e.g., world space)
52 template<typename MapType>
53 inline void
54 calculateBounds(const MapType& map, const BBoxd& in, BBoxd& out)
55 {
56  const Vec3d& min = in.min();
57  const Vec3d& max = in.max();
58 
59  // the pre-image of the 8 corners of the box
60  Vec3d corners[8];
61  corners[0] = in.min();;
62  corners[1] = Vec3d(min(0), min(1), min(2));
63  corners[2] = Vec3d(max(0), max(1), min(2));
64  corners[3] = Vec3d(min(0), max(1), min(2));
65  corners[4] = Vec3d(min(0), min(1), max(2));
66  corners[5] = Vec3d(max(0), min(1), max(2));
67  corners[6] = max;
68  corners[7] = Vec3d(min(0), max(1), max(2));
69 
70  Vec3d pre_image;
71  Vec3d& out_min = out.min();
72  Vec3d& out_max = out.max();
73  out_min = map.applyInverseMap(corners[0]);
74  out_max = min;
75  for (int i = 1; i < 8; ++i) {
76  pre_image = map.applyInverseMap(corners[i]);
77  for (int j = 0; j < 3; ++j) {
78  out_min(j) = std::min( out_min(j), pre_image(j));
79  out_max(j) = std::max( out_max(j), pre_image(j));
80  }
81  }
82 }
83 
84 
85 /// @brief Calculate an axis-aligned bounding box in the given map's domain
86 /// from a spherical bounding box in its range.
87 template<typename MapType>
88 inline void
89 calculateBounds(const MapType& map, const Vec3d& center, const Real radius, BBoxd& out)
90 {
91  // On return, out gives a bounding box in continuous index space
92  // that encloses the sphere.
93  //
94  // the image of a sphere under the inverse of the linearMap will be an ellipsoid.
95 
97  // I want to find extrema for three functions f(x', y', z') = x', or = y', or = z'
98  // with the constraint that g = (x-xo)^2 + (y-yo)^2 + (z-zo)^2 = r^2.
99  // Where the point x,y,z is the image of x',y',z'
100  // Solve: \lambda Grad(g) = Grad(f) and g = r^2.
101  // Note: here (x,y,z) is the image of (x',y',z'), and the gradient
102  // is w.r.t the (') space.
103  //
104  // This can be solved exactly: e_a^T (x' -xo') =\pm r\sqrt(e_a^T J^(-1)J^(-T)e_a)
105  // where e_a is one of the three unit vectors. - djh.
106 
107  /// find the image of the center of the sphere
108  Vec3d center_pre_image = map.applyInverseMap(center);
109 
110  std::vector<Vec3d> coordinate_units;
111  coordinate_units.push_back(Vec3d(1,0,0));
112  coordinate_units.push_back(Vec3d(0,1,0));
113  coordinate_units.push_back(Vec3d(0,0,1));
114 
115  Vec3d& out_min = out.min();
116  Vec3d& out_max = out.max();
117  for (int direction = 0; direction < 3; ++direction) {
118  Vec3d temp = map.applyIJT(coordinate_units[direction]);
119  double offset =
120  radius * sqrt(temp.x()*temp.x() + temp.y()*temp.y() + temp.z()*temp.z());
121  out_min(direction) = center_pre_image(direction) - offset;
122  out_max(direction) = center_pre_image(direction) + offset;
123  }
124 
125  } else {
126  // This is some unknown map type. In this case, we form an axis-aligned
127  // bounding box for the sphere in world space and find the pre-images of
128  // the corners in index space. From these corners we compute an axis-aligned
129  // bounding box in index space.
130  BBoxd bounding_box(center - radius*Vec3d(1,1,1), center + radius*Vec3d(1,1,1));
131  calculateBounds<MapType>(map, bounding_box, out);
132  }
133 }
134 
135 
136 namespace { // anonymous namespace for this helper function
137 
138 /// @brief Find the intersection of a line passing through the point
139 /// (<I>x</I>=0,&nbsp;<I>z</I>=&minus;1/<I>g</I>) with the circle
140 /// (<I>x</I> &minus; <I>xo</I>)&sup2; + (<I>z</I> &minus; <I>zo</I>)&sup2; = <I>r</I>&sup2;
141 /// at a point tangent to the circle.
142 /// @return 0 if the focal point (0, -1/<I>g</I>) is inside the circle,
143 /// 1 if the focal point touches the circle, or 2 when both points are found.
144 inline int
145 findTangentPoints(const double g, const double xo, const double zo,
146  const double r, double& xp, double& zp, double& xm, double& zm)
147 {
148  double x2 = xo * xo;
149  double r2 = r * r;
150  double xd = g * xo;
151  double xd2 = xd*xd;
152  double zd = g * zo + 1.;
153  double zd2 = zd*zd;
154  double rd2 = r2*g*g;
155 
156  double distA = xd2 + zd2;
157  double distB = distA - rd2;
158 
159  if (distB > 0) {
160  double discriminate = sqrt(distB);
161 
162  xp = xo - xo*rd2/distA + r * zd *discriminate / distA;
163  xm = xo - xo*rd2/distA - r * zd *discriminate / distA;
164 
165  zp = (zo*zd2 + zd*g*(x2 - r2) - xo*xo*g - r*xd*discriminate) / distA;
166  zm = (zo*zd2 + zd*g*(x2 - r2) - xo*xo*g + r*xd*discriminate) / distA;
167 
168  return 2;
169 
170  } if (0 >= distB && distB >= -1e-9) {
171  // the circle touches the focal point (x=0, z = -1/g)
172  xp = 0; xm = 0;
173  zp = -1/g; zm = -1/g;
174 
175  return 1;
176  }
177 
178  return 0;
179 }
180 
181 } // end anonymous namespace
182 
183 
184 /// @brief Calculate an axis-aligned bounding box in index space
185 /// from a spherical bounding box in world space.
186 /// @note This specialization is optimized for a frustum map
187 template<>
188 inline void
189 calculateBounds<math::NonlinearFrustumMap>(const math::NonlinearFrustumMap& frustum,
190  const Vec3d& center, const Real radius, BBoxd& out)
191 {
192  // The frustum is a nonlinear map followed by a uniform scale, rotation, translation.
193  // First we invert the translation, rotation and scale to find the spherical pre-image
194  // of the sphere in "local" coordinates where the frustum is aligned with the near plane
195  // on the z=0 plane and the "camera" is located at (x=0, y=0, z=-1/g).
196 
197  // check that the internal map has no shear.
198  const math::AffineMap& secondMap = frustum.secondMap();
199  // test if the linear part has shear or non-uniform scaling
200  if (!frustum.hasSimpleAffine()) {
201 
202  // In this case, we form an axis-aligned bounding box for sphere in world space
203  // and find the pre_images of the corners in voxel space. From these corners we
204  // compute an axis-algined bounding box in voxel-spae
205  BBoxd bounding_box(center - radius*Vec3d(1,1,1), center + radius*Vec3d(1,1,1));
206  calculateBounds<math::NonlinearFrustumMap>(frustum, bounding_box, out);
207  return;
208  }
209 
210  // for convenience
211  Vec3d& out_min = out.min();
212  Vec3d& out_max = out.max();
213 
214  Vec3d centerLS = secondMap.applyInverseMap(center);
215  Vec3d voxelSize = secondMap.voxelSize();
216 
217  // all the voxels have the same size since we know this is a simple affine map
218  double radiusLS = radius / voxelSize(0);
219 
220  double gamma = frustum.getGamma();
221  double xp;
222  double zp;
223  double xm;
224  double zm;
225  int soln_number;
226 
227  // the bounding box in index space for the points in the frustum
228  const BBoxd& bbox = frustum.getBBox();
229  // initialize min and max
230  const double x_min = bbox.min().x();
231  const double y_min = bbox.min().y();
232  const double z_min = bbox.min().z();
233 
234  const double x_max = bbox.max().x();
235  const double y_max = bbox.max().y();
236  const double z_max = bbox.max().z();
237 
238  out_min.x() = x_min;
239  out_max.x() = x_max;
240  out_min.y() = y_min;
241  out_max.y() = y_max;
242 
243  Vec3d extreme;
244  Vec3d extreme2;
245  Vec3d pre_image;
246  // find the x-range
247  soln_number = findTangentPoints(gamma, centerLS.x(), centerLS.z(), radiusLS, xp, zp, xm, zm);
248  if (soln_number == 2) {
249  extreme.x() = xp;
250  extreme.y() = centerLS.y();
251  extreme.z() = zp;
252 
253  // location in world space of the tangent point
254  extreme2 = secondMap.applyMap(extreme);
255  // convert back to voxel space
256  pre_image = frustum.applyInverseMap(extreme2);
257  out_max.x() = std::max(x_min, std::min(x_max, pre_image.x()));
258 
259  extreme.x() = xm;
260  extreme.y() = centerLS.y();
261  extreme.z() = zm;
262  // location in world space of the tangent point
263  extreme2 = secondMap.applyMap(extreme);
264 
265  // convert back to voxel space
266  pre_image = frustum.applyInverseMap(extreme2);
267  out_min.x() = std::max(x_min, std::min(x_max, pre_image.x()));
268 
269  } else if (soln_number == 1) {
270  // the circle was tangent at the focal point
271  } else if (soln_number == 0) {
272  // the focal point was inside the circle
273  }
274 
275  // find the y-range
276  soln_number = findTangentPoints(gamma, centerLS.y(), centerLS.z(), radiusLS, xp, zp, xm, zm);
277  if (soln_number == 2) {
278  extreme.x() = centerLS.x();
279  extreme.y() = xp;
280  extreme.z() = zp;
281 
282  // location in world space of the tangent point
283  extreme2 = secondMap.applyMap(extreme);
284  // convert back to voxel space
285  pre_image = frustum.applyInverseMap(extreme2);
286  out_max.y() = std::max(y_min, std::min(y_max, pre_image.y()));
287 
288  extreme.x() = centerLS.x();
289  extreme.y() = xm;
290  extreme.z() = zm;
291  extreme2 = secondMap.applyMap(extreme);
292 
293  // convert back to voxel space
294  pre_image = frustum.applyInverseMap(extreme2);
295  out_min.y() = std::max(y_min, std::min(y_max, pre_image.y()));
296 
297  } else if (soln_number == 1) {
298  // the circle was tangent at the focal point
299  } else if (soln_number == 0) {
300  // the focal point was inside the circle
301  }
302 
303  // the near and far
304  // the closest point. The front of the frustum is at 0 in index space
305  double near_dist = std::max(centerLS.z() - radiusLS, 0.);
306  // the farthest point. The back of the frustum is at mDepth in index space
307  double far_dist = std::min(centerLS.z() + radiusLS, frustum.getDepth() );
308 
309  Vec3d near_point(0.f, 0.f, near_dist);
310  Vec3d far_point(0.f, 0.f, far_dist);
311 
312  out_min.z() = std::max(z_min, frustum.applyInverseMap(secondMap.applyMap(near_point)).z());
313  out_max.z() = std::min(z_max, frustum.applyInverseMap(secondMap.applyMap(far_point)).z());
314 
315 }
316 
317 } // namespace util
318 } // namespace OPENVDB_VERSION_NAME
319 } // namespace openvdb
320 
321 #endif // OPENVDB_UTIL_MAPSUTIL_HAS_BEEN_INCLUDED
322 
323 // Copyright (c) 2012-2018 DreamWorks Animation LLC
324 // All rights reserved. This software is distributed under the
325 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
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:493
IMF_EXPORT IMATH_NAMESPACE::V3f direction(const IMATH_NAMESPACE::Box2i &dataWindow, const IMATH_NAMESPACE::V2f &pixelPosition)
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:847
GLboolean GLboolean g
Definition: glcorearb.h:1221
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:1905
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:189
png_uint_32 i
Definition: png.h:2877
GLfloat f
Definition: glcorearb.h:1925
const Vec3T & min() const
Return a const reference to the minimum point of this bounding box.
Definition: BBox.h:89
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:133
A general linear transform using homogeneous coordinates to perform rotation, scaling, shear and translation.
Definition: Maps.h:327
GLintptr offset
Definition: glcorearb.h:664
Vec3d applyInverseMap(const Vec3d &in) const override
Return the pre-image of in under the map.
Definition: Maps.h:445
const Vec3T & max() const
Return a const reference to the maximum point of this bounding box.
Definition: BBox.h:91
Vec3d applyMap(const Vec3d &in) const override
Return the image of in under the map.
Definition: Maps.h:443
void calculateBounds(const MapType &map, const BBoxd &in, BBoxd &out)
Calculate an axis-aligned bounding box in the given map's domain (e.g., index space) from an axis-ali...
Definition: MapsUtil.h:54
T & x()
Reference to the component, e.g. v.x() = 4.5f;.
Definition: Vec3.h:110
GLboolean r
Definition: glcorearb.h:1221
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:129
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:135