OpenVDB  0.104.0
Maps.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_MATH_MAPS_HAS_BEEN_INCLUDED
34 #define OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED
35 
36 #include "Mat4.h"
37 #include "Vec3.h"
38 #include "BBox.h"
39 #include "Coord.h"
40 #include <openvdb/util/Name.h>
41 #include <openvdb/Types.h>
42 #include <boost/shared_ptr.hpp>
43 #include <tbb/mutex.h>
44 #include <map>
45 
46 namespace openvdb {
48 namespace OPENVDB_VERSION_NAME {
49 namespace math {
50 
51 
53 
55 
56 class MapBase;
57 class ScaleMap;
58 class TranslationMap;
59 class ScaleTranslateMap;
60 class UniformScaleMap;
61 class UniformScaleTranslateMap;
62 class AffineMap;
63 class UnitaryMap;
64 class NonlinearFrustumMap;
65 
66 template<typename T1, typename T2> class CompoundMap;
67 
73 
74 
76 
78 
79 template<typename T> struct is_linear { static const bool value = false; };
80 template<> struct is_linear<AffineMap> { static const bool value = true; };
81 template<> struct is_linear<ScaleMap> { static const bool value = true; };
82 template<> struct is_linear<UniformScaleMap> { static const bool value = true; };
83 template<> struct is_linear<UnitaryMap> { static const bool value = true; };
84 template<> struct is_linear<TranslationMap> { static const bool value = true; };
85 template<> struct is_linear<ScaleTranslateMap> { static const bool value = true; };
86 template<> struct is_linear<UniformScaleTranslateMap> { static const bool value = true; };
87 
88 template<typename T1, typename T2> struct is_linear<CompoundMap<T1, T2> > {
89  static const bool value = is_linear<T1>::value && is_linear<T2>::value;
90 };
91 
92 
93 template<typename T> struct is_uniform_scale { static const bool value = false; };
94 template<> struct is_uniform_scale<UniformScaleMap> { static const bool value = true; };
95 
96 template<typename T> struct is_uniform_scale_translate { static const bool value = false; };
97 template<> struct is_uniform_scale_translate<TranslationMap> { static const bool value = true; };
99  static const bool value = true;
100 };
101 
102 
103 template<typename T> struct is_scale { static const bool value = false; };
104 template<> struct is_scale<ScaleMap> { static const bool value = true; };
105 
106 template<typename T> struct is_scale_translate { static const bool value = false; };
107 template<> struct is_scale_translate<ScaleTranslateMap> { static const bool value = true; };
108 
109 
110 template<typename T> struct is_uniform_diagonal_jacobian {
112 };
113 
114 template<typename T> struct is_diagonal_jacobian {
115  static const bool value = is_scale<T>::value || is_scale_translate<T>::value;
116 };
117 
118 
120 
122 
125 OPENVDB_API boost::shared_ptr<SymmetricMap> createSymmetricMap(const Mat3d& m);
126 
127 
130 OPENVDB_API boost::shared_ptr<FullyDecomposedMap> createFullyDecomposedMap(const Mat4d& m);
131 
132 
143 OPENVDB_API boost::shared_ptr<PolarDecomposedMap> createPolarDecomposedMap(const Mat3d& m);
144 
145 
147 OPENVDB_API boost::shared_ptr<MapBase> simplify(boost::shared_ptr<AffineMap> affine);
148 
149 
151 
152 
155 {
156 public:
157  typedef boost::shared_ptr<MapBase> Ptr;
158  typedef boost::shared_ptr<const MapBase> ConstPtr;
159  typedef Ptr (*MapFactory)();
160 
161  virtual ~MapBase(){}
162 
163  virtual boost::shared_ptr<AffineMap> getAffineMap() const = 0;
164 
166  virtual Name type() const = 0;
167 
169  template<typename MapT> bool isType() const { return this->type() == MapT::mapType(); }
170 
172  virtual bool isEqual(const MapBase& other) const = 0;
173 
175  virtual bool isLinear() const = 0;
177  virtual bool hasUniformScale() const = 0;
178 
179  virtual Vec3d applyMap(const Vec3d& in) const = 0;
180  virtual Vec3d applyInverseMap(const Vec3d& in) const = 0;
181 
182  virtual Vec3d applyIJT(const Vec3d& in) const = 0;
183  virtual Vec3d applyIJT(const Vec3d& in, const Vec3d& pos) const = 0;
184  virtual Mat3d applyIJC(const Mat3d& m) const = 0;
185  virtual Mat3d applyIJC(const Mat3d& m, const Vec3d& v, const Vec3d& pos) const = 0;
186 
187 
188  virtual double determinant() const = 0;
189  virtual double determinant(const Vec3d&) const = 0;
190 
191 
193 
194  virtual Vec3d voxelSize() const = 0;
195  virtual Vec3d voxelSize(const Vec3d&) const = 0;
197 
199 
200 
201  virtual Vec3d voxelDimensions() const = 0;
202  virtual Vec3d voxelDimensions(const Vec3d&) const = 0;
204 
205  virtual void read(std::istream&) = 0;
206  virtual void write(std::ostream&) const = 0;
207 
208  virtual std::string str() const = 0;
209 
210  virtual MapBase::Ptr copy() const = 0;
211 
213 
214  virtual MapBase::Ptr preRotate(double radians, Axis axis = X_AXIS) const = 0;
215  virtual MapBase::Ptr preTranslate(const Vec3d&) const = 0;
216  virtual MapBase::Ptr preScale(const Vec3d&) const = 0;
217  virtual MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const = 0;
218 
219  virtual MapBase::Ptr postRotate(double radians, Axis axis = X_AXIS) const = 0;
220  virtual MapBase::Ptr postTranslate(const Vec3d&) const = 0;
221  virtual MapBase::Ptr postScale(const Vec3d&) const = 0;
222  virtual MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const = 0;
224 
225 
227  virtual MapBase::Ptr rotate(double radians, Axis axis = X_AXIS) const = 0;
229  virtual MapBase::Ptr translate(const Vec3d&) const = 0;
231  virtual MapBase::Ptr scale(const Vec3d&) const = 0;
233  virtual MapBase::Ptr shear(double shear, Axis axis0, Axis axis1) const = 0;
234 
235 protected:
236  MapBase() {}
237 
238  template<typename MapT>
239  static bool isEqualBase(const MapT& self, const MapBase& other)
240  {
241  return other.isType<MapT>() && (self == *static_cast<const MapT*>(&other));
242  }
243 };
244 
245 
247 
248 
252 {
253 public:
254  typedef std::map<Name, MapBase::MapFactory> MapDictionary;
255  typedef tbb::mutex Mutex;
256  typedef Mutex::scoped_lock Lock;
257 
258  static MapRegistry* instance();
259 
261  static MapBase::Ptr createMap(const Name&);
262 
264  static bool isRegistered(const Name&);
265 
267  static void registerMap(const Name&, MapBase::MapFactory);
268 
270  static void unregisterMap(const Name&);
271 
273  static void clear();
274 
275 private:
276  MapRegistry() {}
277  static MapRegistry* mInstance;
278 
279  mutable Mutex mMutex;
280  MapDictionary mMap;
281 };
282 
283 
285 
286 
290 {
291 public:
292  typedef boost::shared_ptr<AffineMap> Ptr;
293  typedef boost::shared_ptr<const AffineMap> ConstPtr;
294 
296  mMatrix(Mat4d::identity()),
297  mMatrixInv(Mat4d::identity()),
298  mJacobianInv(Mat3d::identity()),
299  mDeterminant(1),
300  mVoxelSize(Vec3d(1,1,1)),
301  mIsDiagonal(true),
302  mIsIdentity(true)
303  // the default constructor for translation is zero
304  {
305  }
306 
307  AffineMap(const Mat3d& m)
308  {
309  Mat4d mat4(Mat4d::identity());
310  mat4.setMat3(m);
311  mMatrix = mat4;
312  updateAcceleration();
313  }
314 
315  AffineMap(const Mat4d& m): mMatrix(m)
316  {
317  if (!isAffine(m)) {
319  "Tried to initialize an affine transform from a non-affine 4x4 matrix");
320  }
321  updateAcceleration();
322  }
323 
324  AffineMap(const AffineMap& other):
325  MapBase(other),
326  mMatrix(other.mMatrix),
327  mMatrixInv(other.mMatrixInv),
328  mJacobianInv(other.mJacobianInv),
329  mDeterminant(other.mDeterminant),
330  mVoxelSize(other.mVoxelSize),
331  mIsDiagonal(other.mIsDiagonal),
332  mIsIdentity(other.mIsIdentity)
333  {
334  }
335 
337  AffineMap(const AffineMap& first, const AffineMap& second):
338  mMatrix(first.mMatrix * second.mMatrix)
339  {
340  updateAcceleration();
341  }
342 
344 
346  static MapBase::Ptr create() { return MapBase::Ptr(new AffineMap()); }
348  MapBase::Ptr copy() const { return MapBase::Ptr(new AffineMap(*this)); }
349 
350  static bool isRegistered() { return MapRegistry::isRegistered(AffineMap::mapType()); }
351 
352  static void registerMap()
353  {
354  MapRegistry::registerMap(
355  AffineMap::mapType(),
356  AffineMap::create);
357  }
358 
359  Name type() const { return mapType(); }
360  static Name mapType() { return Name("AffineMap"); }
361 
363  bool isLinear() const { return true; }
364 
366  bool hasUniformScale() const { return isUnitary(mMatrix.getMat3());}
367 
368  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
369 
370  bool operator==(const AffineMap& other) const
371  {
372  // the Mat.eq() is approximate
373  if (!mMatrix.eq(other.mMatrix)) { return false; }
374  if (!mMatrixInv.eq(other.mMatrixInv)) { return false; }
375  return true;
376  }
377 
378  bool operator!=(const AffineMap& other) const { return !(*this == other); }
379 
380  AffineMap& operator=(const AffineMap& other)
381  {
382  mMatrix = other.mMatrix;
383  mMatrixInv = other.mMatrixInv;
384 
385  mJacobianInv = other.mJacobianInv;
386  mDeterminant = other.mDeterminant;
387  mVoxelSize = other.mVoxelSize;
388  mIsDiagonal = other.mIsDiagonal;
389  mIsIdentity = other.mIsIdentity;
390  return *this;
391  }
393  Vec3d applyMap(const Vec3d& in) const { return in * mMatrix; }
395  Vec3d applyInverseMap(const Vec3d& in) const {return in * mMatrixInv; }
396 
398  Vec3d applyIJT(const Vec3d& in, const Vec3d&) const { return applyIJT(in); }
400  Vec3d applyIJT(const Vec3d& in) const { return in * mJacobianInv; }
402  Mat3d applyIJC(const Mat3d& m) const {
403  return mJacobianInv.transpose()* m * mJacobianInv;
404  }
405  Mat3d applyIJC(const Mat3d& in, const Vec3d& , const Vec3d& ) const {
406  return applyIJC(in);
407  }
409  double determinant(const Vec3d& ) const { return determinant(); }
411  double determinant() const { return mDeterminant; }
412 
414 OPENVDB_DEPRECATED Vec3d voxelDimensions() const { return voxelSize(); }
416 OPENVDB_DEPRECATED Vec3d voxelDimensions(const Vec3d&) const { return voxelSize(); }
417 
419 
420 
421  Vec3d voxelSize() const { return mVoxelSize; }
422  Vec3d voxelSize(const Vec3d&) const { return voxelSize(); }
424 
426  bool isIdentity() const { return mIsIdentity; }
428  bool isDiagonal() const { return mIsDiagonal; }
430  bool isScale() const { return isDiagonal(); }
432  bool isScaleTranslate() const { return math::isDiagonal(mMatrix.getMat3()); }
433 
434 
435  // Methods that modify the existing affine map
436 
438 
439 
440  void accumPreRotation(Axis axis, double radians)
441  {
442  mMatrix.preRotate(axis, radians);
443  updateAcceleration();
444  }
445  void accumPreScale(const Vec3d& v)
446  {
447  mMatrix.preScale(v);
448  updateAcceleration();
449  }
450  void accumPreTranslation(const Vec3d& v)
451  {
452  mMatrix.preTranslate(v);
453  updateAcceleration();
454  }
455  void accumPreShear(Axis axis0, Axis axis1, double shear)
456  {
457  mMatrix.preShear(axis0, axis1, shear);
458  updateAcceleration();
459  }
461 
462 
464 
465 
466  void accumPostRotation(Axis axis, double radians)
467  {
468  mMatrix.postRotate(axis, radians);
469  updateAcceleration();
470  }
471  void accumPostScale(const Vec3d& v)
472  {
473  mMatrix.postScale(v);
474  updateAcceleration();
475  }
476  void accumPostTranslation(const Vec3d& v)
477  {
478  mMatrix.postTranslate(v);
479  updateAcceleration();
480  }
481  void accumPostShear(Axis axis0, Axis axis1, double shear)
482  {
483  mMatrix.postShear(axis0, axis1, shear);
484  updateAcceleration();
485  }
487 
488 
490  void read(std::istream& is)
491  {
492  mMatrix.read(is);
493  updateAcceleration();
494  }
495 
497  void write(std::ostream& os) const
498  {
499  mMatrix.write(os);
500  }
501 
503  std::string str() const
504  {
505  std::ostringstream buffer;
506  buffer << " - mat4:\n" << mMatrix.str() << std::endl;
507  buffer << " - voxel dimensions: " << mVoxelSize << std::endl;
508  return buffer.str();
509  }
510 
512  boost::shared_ptr<FullyDecomposedMap> createDecomposedMap()
513  {
514  return createFullyDecomposedMap(mMatrix);
515  }
516 
518  AffineMap::Ptr getAffineMap() const { return AffineMap::Ptr(new AffineMap(*this)); }
519 
521  AffineMap::Ptr inverse() const { return AffineMap::Ptr(new AffineMap(mMatrixInv)); }
522 
523 
525 
526 
527  MapBase::Ptr preRotate(double radians, Axis axis = X_AXIS) const
528  {
529  AffineMap::Ptr affineMap = getAffineMap();
530  affineMap->accumPreRotation(axis, radians);
531  return simplify(affineMap);
532  }
533  MapBase::Ptr preTranslate(const Vec3d& t) const
534  {
535  AffineMap::Ptr affineMap = getAffineMap();
536  affineMap->accumPreTranslation(t);
537  return boost::static_pointer_cast<MapBase, AffineMap>(affineMap);
538  }
539  MapBase::Ptr preScale(const Vec3d& s) const
540  {
541  AffineMap::Ptr affineMap = getAffineMap();
542  affineMap->accumPreScale(s);
543  return boost::static_pointer_cast<MapBase, AffineMap>(affineMap);
544  }
545  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
546  {
547  AffineMap::Ptr affineMap = getAffineMap();
548  affineMap->accumPreShear(axis0, axis1, shear);
549  return simplify(affineMap);
550  }
552 
553 
555 
556 
557  MapBase::Ptr postRotate(double radians, Axis axis = X_AXIS) const
558  {
559  AffineMap::Ptr affineMap = getAffineMap();
560  affineMap->accumPostRotation(axis, radians);
561  return simplify(affineMap);
562  }
563  MapBase::Ptr postTranslate(const Vec3d& t) const
564  {
565  AffineMap::Ptr affineMap = getAffineMap();
566  affineMap->accumPostTranslation(t);
567  return boost::static_pointer_cast<MapBase, AffineMap>(affineMap);
568  }
569  MapBase::Ptr postScale(const Vec3d& s) const
570  {
571  AffineMap::Ptr affineMap = getAffineMap();
572  affineMap->accumPostScale(s);
573  return boost::static_pointer_cast<MapBase, AffineMap>(affineMap);
574  }
575  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
576  {
577  AffineMap::Ptr affineMap = getAffineMap();
578  affineMap->accumPostShear(axis0, axis1, shear);
579  return simplify(affineMap);
580  }
582 
583 
585 OPENVDB_DEPRECATED void accumRotation(Axis axis, double radians)
586  { accumPreRotation(axis, radians); }
588 OPENVDB_DEPRECATED void accumScale(const Vec3d& v)
589  { accumPreScale(v); }
591 OPENVDB_DEPRECATED void accumTranslation(const Vec3d& v)
592  { accumPostTranslation(v); }
594 OPENVDB_DEPRECATED void accumShear(Axis axis0, Axis axis1, double shear)
595  { accumPreShear(axis0, axis1, shear); }
596 
598 OPENVDB_DEPRECATED MapBase::Ptr rotate(double radians, Axis axis = X_AXIS) const
599  { return preRotate(radians, axis); }
601 OPENVDB_DEPRECATED MapBase::Ptr translate(const Vec3d& t) const
602  { return postTranslate(t); }
605  { return preScale(s); }
607 OPENVDB_DEPRECATED MapBase::Ptr shear(double shear, Axis axis0, Axis axis1) const
608  { return preShear(shear, axis0, axis1); }
609 
611  Mat4d getMat4() const { return mMatrix;}
612  const Mat4d& getConstMat4() const {return mMatrix;}
613  const Mat3d& getConstJacobianInv() const {return mJacobianInv;}
614 private:
615  void updateAcceleration() {
616  mDeterminant = mMatrix.getMat3().det();
617 
618  if (std::abs(mDeterminant) < (3.0 * tolerance<double>::value())) {
620  "Tried to initialize an affine transform from a nearly signular matrix");
621  }
622  mMatrixInv = mMatrix.inverse();
623  mJacobianInv = mMatrixInv.getMat3().transpose();
624  mIsDiagonal = math::isDiagonal(mMatrix);
625  mIsIdentity = math::isIdentity(mMatrix);
626  Vec3d pos = applyMap(Vec3d(0,0,0));
627  mVoxelSize(0) = (applyMap(Vec3d(1,0,0)) - pos).length();
628  mVoxelSize(1) = (applyMap(Vec3d(0,1,0)) - pos).length();
629  mVoxelSize(2) = (applyMap(Vec3d(0,0,1)) - pos).length();
630  }
631 
632  // the underlying matrix
633  Mat4d mMatrix;
634 
635  // stored for acceleration
636  Mat4d mMatrixInv;
637  Mat3d mJacobianInv;
638  double mDeterminant;
639  Vec3d mVoxelSize;
640  bool mIsDiagonal, mIsIdentity;
641 }; // class AffineMap
642 
643 
645 
646 
650 {
651 public:
652  typedef boost::shared_ptr<ScaleMap> Ptr;
653  typedef boost::shared_ptr<const ScaleMap> ConstPtr;
654 
655  ScaleMap(): MapBase(), mScaleValues(Vec3d(1,1,1)), mVoxelSize(Vec3d(1,1,1)),
656  mInvScaleSqr(1,1,1), mInvTwiceScale(0.5,0.5,0.5){}
657 
659  MapBase(),
660  mScaleValues(scale),
661  mVoxelSize(Vec3d(std::abs(scale(0)),std::abs(scale(1)), std::abs(scale(2))))
662  {
663  if (scale.eq(Vec3d(0.0), 1.0e-10)) {
664  OPENVDB_THROW(ArithmeticError, "Non-zero scale values required");
665  }
666  mScaleValuesInverse = 1.0 / mScaleValues;
667  mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
668  mInvTwiceScale = mScaleValuesInverse / 2;
669  }
670 
671  ScaleMap(const ScaleMap& other):
672  MapBase(),
673  mScaleValues(other.mScaleValues),
674  mVoxelSize(other.mVoxelSize),
675  mScaleValuesInverse(other.mScaleValuesInverse),
676  mInvScaleSqr(other.mInvScaleSqr),
677  mInvTwiceScale(other.mInvTwiceScale)
678  {
679  }
680 
682 
684  static MapBase::Ptr create() { return MapBase::Ptr(new ScaleMap()); }
686  MapBase::Ptr copy() const { return MapBase::Ptr(new ScaleMap(*this)); }
687 
688  static bool isRegistered() { return MapRegistry::isRegistered(ScaleMap::mapType()); }
689 
690  static void registerMap()
691  {
692  MapRegistry::registerMap(
693  ScaleMap::mapType(),
694  ScaleMap::create);
695  }
696 
697  Name type() const { return mapType(); }
698  static Name mapType() { return Name("ScaleMap"); }
699 
701  bool isLinear() const { return true; }
702 
704  bool hasUniformScale() const {
705  bool value;
706  value = isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.y()));
707  value = value && isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.z()));
708  return value; }
709 
711  Vec3d applyMap(const Vec3d& in) const
712  {
713  return Vec3d(
714  in.x() * mScaleValues.x(),
715  in.y() * mScaleValues.y(),
716  in.z() * mScaleValues.z());
717  }
719  Vec3d applyInverseMap(const Vec3d& in) const
720  {
721  return Vec3d(
722  in.x() * mScaleValuesInverse.x(),
723  in.y() * mScaleValuesInverse.y(),
724  in.z() * mScaleValuesInverse.z());
725  }
728  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
730  Vec3d applyIJT(const Vec3d& in) const { return applyInverseMap(in); }
732  Mat3d applyIJC(const Mat3d& in) const {
733  Mat3d tmp;
734  for (int i=0; i<3; i++){
735  tmp.setRow(i, in.row(i)*mScaleValuesInverse(i));
736  }
737  for (int i=0; i<3; i++){
738  tmp.setCol(i, tmp.col(i)*mScaleValuesInverse(i));
739  }
740  return tmp;
741  }
742  Mat3d applyIJC(const Mat3d& in, const Vec3d&, const Vec3d& ) const {
743  return applyIJC(in);
744  }
746  double determinant(const Vec3d& ) const { return determinant(); }
748  double determinant() const { return mScaleValues.x()*mScaleValues.y()*mScaleValues.z(); }
749 
751  const Vec3d& getScale() const {return mScaleValues;}
752 
754  const Vec3d& getInvScaleSqr() const {return mInvScaleSqr;}
756  const Vec3d& getInvTwiceScale() const {return mInvTwiceScale;}
758  const Vec3d& getInvScale() const {return mScaleValuesInverse; }
759 
760 
762  OPENVDB_DEPRECATED Vec3d voxelDimensions() const { return mVoxelSize; }
764  OPENVDB_DEPRECATED Vec3d voxelDimensions(const Vec3d&) const { return voxelSize();}
765 
767 
768 
769 
770 
771  Vec3d voxelSize() const { return mVoxelSize; }
772  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
774 
776  void read(std::istream& is)
777  {
778  mScaleValues.read(is);
779  mVoxelSize.read(is);
780  mScaleValuesInverse.read(is);
781  mInvScaleSqr.read(is);
782  mInvTwiceScale.read(is);
783  }
785  void write(std::ostream& os) const
786  {
787  mScaleValues.write(os);
788  mVoxelSize.write(os);
789  mScaleValuesInverse.write(os);
790  mInvScaleSqr.write(os);
791  mInvTwiceScale.write(os);
792  }
794  std::string str() const
795  {
796  std::ostringstream buffer;
797  buffer << " - scale: " << mScaleValues << std::endl;
798  buffer << " - voxel dimensions: " << mVoxelSize << std::endl;
799  return buffer.str();
800  }
801 
802  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
803 
804  bool operator==(const ScaleMap& other) const
805  {
806  // ::eq() uses a tolerance
807  if (!mScaleValues.eq(other.mScaleValues)) { return false; }
808  return true;
809  }
810 
811  bool operator!=(const ScaleMap& other) const { return !(*this == other); }
812 
814  AffineMap::Ptr getAffineMap() const
815  {
816  return AffineMap::Ptr(new AffineMap(math::scale<Mat4d>(mScaleValues)));
817  }
818 
819 
820 
822 
823 
824  MapBase::Ptr preRotate(double radians, Axis axis) const
825  {
826  AffineMap::Ptr affineMap = getAffineMap();
827  affineMap->accumPreRotation(axis, radians);
828  return simplify(affineMap);
829  }
830 
831  MapBase::Ptr preTranslate(const Vec3d& tr) const;
832 
833  MapBase::Ptr preScale(const Vec3d& v) const;
834 
835  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
836  {
837  AffineMap::Ptr affineMap = getAffineMap();
838  affineMap->accumPreShear(axis0, axis1, shear);
839  return simplify(affineMap);
840  }
842 
843 
845 
846 
847  MapBase::Ptr postRotate(double radians, Axis axis) const
848  {
849  AffineMap::Ptr affineMap = getAffineMap();
850  affineMap->accumPostRotation(axis, radians);
851  return simplify(affineMap);
852  }
853 
854  MapBase::Ptr postTranslate(const Vec3d& tr) const;
855 
856  MapBase::Ptr postScale(const Vec3d& v) const;
857 
858  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
859  {
860  AffineMap::Ptr affineMap = getAffineMap();
861  affineMap->accumPostShear(axis0, axis1, shear);
862  return simplify(affineMap);
863  }
865 
867 OPENVDB_DEPRECATED MapBase::Ptr rotate(double radians, Axis axis) const
868  { return preRotate(radians, axis); }
870 OPENVDB_DEPRECATED MapBase::Ptr translate(const Vec3d& tr) const;
874 OPENVDB_DEPRECATED MapBase::Ptr shear(double shear, Axis axis0, Axis axis1) const
875  { return preShear(shear, axis0, axis1); }
876 
877 private:
878  Vec3d mScaleValues, mVoxelSize, mScaleValuesInverse, mInvScaleSqr, mInvTwiceScale;
879 }; // class ScaleMap
880 
881 
885 {
886 public:
887  typedef boost::shared_ptr<UniformScaleMap> Ptr;
888  typedef boost::shared_ptr<const UniformScaleMap> ConstPtr;
889 
891  UniformScaleMap(double scale): ScaleMap(Vec3d(scale, scale, scale)) {}
892  UniformScaleMap(const UniformScaleMap& other): ScaleMap(other) {}
894 
896  static MapBase::Ptr create() { return MapBase::Ptr(new UniformScaleMap()); }
898  MapBase::Ptr copy() const { return MapBase::Ptr(new UniformScaleMap(*this)); }
899 
900  static bool isRegistered() { return MapRegistry::isRegistered(UniformScaleMap::mapType()); }
901  static void registerMap()
902  {
903  MapRegistry::registerMap(
904  UniformScaleMap::mapType(),
905  UniformScaleMap::create);
906  }
907 
908  Name type() const { return mapType(); }
909  static Name mapType() { return Name("UniformScaleMap"); }
910 
911  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
912 
913  bool operator==(const UniformScaleMap& other) const { return ScaleMap::operator==(other); }
914  bool operator!=(const UniformScaleMap& other) const { return !(*this == other); }
915 
918  MapBase::Ptr preTranslate(const Vec3d& tr) const;
919 
922  MapBase::Ptr postTranslate(const Vec3d& tr) const;
923 
925 OPENVDB_DEPRECATED MapBase::Ptr translate(const Vec3d& tr) const;
926 
927 }; // class UniformScaleMap
928 
930 
931 
932 inline MapBase::Ptr
933 ScaleMap::preScale(const Vec3d& v) const
934 {
935  Vec3d new_scale(v * mScaleValues);
936  if (isApproxEqual(new_scale[0],new_scale[1]) && isApproxEqual(new_scale[0],new_scale[2])) {
937  return MapBase::Ptr(new UniformScaleMap(new_scale[0]));
938  } else {
939  return MapBase::Ptr(new ScaleMap(new_scale));
940  }
941 }
942 
943 
944 inline MapBase::Ptr
945 ScaleMap::postScale(const Vec3d& v) const
946 { // pre-post Scale are the same for a scale map
947  return preScale(v);
948 }
949 
950 
952 ScaleMap::scale(const Vec3d& v) const
953 { return preScale(v); }
954 
957 {
958 public:
959  typedef boost::shared_ptr<TranslationMap> Ptr;
960  typedef boost::shared_ptr<const TranslationMap> ConstPtr;
961 
962  // default constructor is a translation by zero.
963  TranslationMap(): MapBase(),mTranslation(Vec3d(0,0,0)) {}
964 
965  TranslationMap(const Vec3d& t): MapBase(),mTranslation(t)
966  {
967  }
968 
969  TranslationMap(const TranslationMap& other): MapBase(),mTranslation(other.mTranslation)
970  {
971  }
972 
974 
976  static MapBase::Ptr create() { return MapBase::Ptr(new TranslationMap()); }
978  MapBase::Ptr copy() const { return MapBase::Ptr(new TranslationMap(*this)); }
979 
980  static bool isRegistered() { return MapRegistry::isRegistered(TranslationMap::mapType()); }
981 
982  static void registerMap()
983  {
984  MapRegistry::registerMap(
985  TranslationMap::mapType(),
986  TranslationMap::create);
987  }
988 
989  Name type() const { return mapType(); }
990  static Name mapType() { return Name("TranslationMap"); }
991 
993  bool isLinear() const { return true; }
994 
996  bool hasUniformScale() const { return true; }
997 
999  Vec3d applyMap(const Vec3d& in) const { return in + mTranslation; }
1001  Vec3d applyInverseMap(const Vec3d& in) const { return in - mTranslation; }
1002 
1005  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
1008  Vec3d applyIJT(const Vec3d& in) const {return in;}
1010  Mat3d applyIJC(const Mat3d& mat) const {return mat;}
1011  Mat3d applyIJC(const Mat3d& mat, const Vec3d&, const Vec3d&) const {
1012  return applyIJC(mat);
1013  }
1014 
1016  double determinant(const Vec3d& ) const { return determinant(); }
1018  double determinant() const { return 1.0; }
1019 
1021  Vec3d voxelSize() const { return Vec3d(1,1,1);}
1025  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
1027 OPENVDB_DEPRECATED Vec3d voxelDimensions(const Vec3d&) const { return voxelSize();}
1028 
1030  const Vec3d& getTranslation() const { return mTranslation; }
1032  void read(std::istream& is) { mTranslation.read(is); }
1034  void write(std::ostream& os) const { mTranslation.write(os); }
1035 
1037  std::string str() const
1038  {
1039  std::ostringstream buffer;
1040  buffer << " - translation: " << mTranslation << std::endl;
1041  return buffer.str();
1042  }
1043 
1044  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1045 
1046  bool operator==(const TranslationMap& other) const
1047  {
1048  // ::eq() uses a tolerance
1049  if (!mTranslation.eq(other.mTranslation)) { return false; }
1050  return true;
1051  }
1052 
1053  bool operator!=(const TranslationMap& other) const { return !(*this == other); }
1054 
1056  AffineMap::Ptr getAffineMap() const
1057  {
1058  Mat4d matrix(Mat4d::identity());
1059  matrix.setTranslation(mTranslation);
1060 
1061  AffineMap::Ptr affineMap(new AffineMap(matrix));
1062  return affineMap;
1063  }
1064 
1065 
1067 
1068 
1069  MapBase::Ptr preRotate(double radians, Axis axis) const
1070  {
1071  AffineMap::Ptr affineMap = getAffineMap();
1072  affineMap->accumPreRotation(axis, radians);
1073  return simplify(affineMap);
1074 
1075  }
1076  MapBase::Ptr preTranslate(const Vec3d& t) const
1077  {
1078  return MapBase::Ptr(new TranslationMap(t + mTranslation));
1079  }
1080 
1081  MapBase::Ptr preScale(const Vec3d& v) const;
1082 
1083  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
1084  {
1085  AffineMap::Ptr affineMap = getAffineMap();
1086  affineMap->accumPreShear(axis0, axis1, shear);
1087  return simplify(affineMap);
1088  }
1090 
1092 
1093 
1094  MapBase::Ptr postRotate(double radians, Axis axis) const
1095  {
1096  AffineMap::Ptr affineMap = getAffineMap();
1097  affineMap->accumPostRotation(axis, radians);
1098  return simplify(affineMap);
1099 
1100  }
1101  MapBase::Ptr postTranslate(const Vec3d& t) const
1102  { // post and pre are the same for this
1103  return MapBase::Ptr(new TranslationMap(t + mTranslation));
1104  }
1105 
1106  MapBase::Ptr postScale(const Vec3d& v) const;
1107 
1108  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
1109  {
1110  AffineMap::Ptr affineMap = getAffineMap();
1111  affineMap->accumPostShear(axis0, axis1, shear);
1112  return simplify(affineMap);
1113  }
1115 
1117 
1119 OPENVDB_DEPRECATED MapBase::Ptr rotate(double radians, Axis axis) const
1120  { return preRotate(radians, axis); }
1122 OPENVDB_DEPRECATED MapBase::Ptr translate(const Vec3d& t) const
1123  { return preTranslate(t); }
1125 OPENVDB_DEPRECATED MapBase::Ptr scale(const Vec3d& v) const;
1128  { return preShear(shear, axis0, axis1); }
1129 
1130 private:
1131  Vec3d mTranslation;
1132 }; // class TranslationMap
1133 
1134 
1136 
1137 
1142 {
1143 public:
1144  typedef boost::shared_ptr<ScaleTranslateMap> Ptr;
1145  typedef boost::shared_ptr<const ScaleTranslateMap> ConstPtr;
1146 
1148  MapBase(),
1149  mTranslation(Vec3d(0,0,0)),
1150  mScaleValues(Vec3d(1,1,1)),
1151  mVoxelSize(Vec3d(1,1,1)),
1152  mScaleValuesInverse(Vec3d(1,1,1)),
1153  mInvScaleSqr(1,1,1),
1154  mInvTwiceScale(0.5,0.5,0.5)
1155  {
1156  }
1157 
1158  ScaleTranslateMap(const Vec3d& scale, const Vec3d& translate):
1159  MapBase(),
1160  mTranslation(translate),
1161  mScaleValues(scale),
1162  mVoxelSize(std::abs(scale(0)), std::abs(scale(1)), std::abs(scale(2)))
1163  {
1164  if (scale.eq(Vec3d(0.0), 1.0e-10)) {
1165  OPENVDB_THROW(ArithmeticError, "Non-zero scale values required");
1166  }
1167  mScaleValuesInverse = 1.0 / mScaleValues;
1168  mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1169  mInvTwiceScale = mScaleValuesInverse / 2;
1170  }
1171 
1172  ScaleTranslateMap(const ScaleMap& scale, const TranslationMap& translate):
1173  MapBase(),
1174  mTranslation(translate.getTranslation()),
1175  mScaleValues(scale.getScale()),
1176  mVoxelSize(std::abs(mScaleValues(0)),
1177  std::abs(mScaleValues(1)),
1178  std::abs(mScaleValues(2))),
1179  mScaleValuesInverse(1.0 / scale.getScale())
1180  {
1181  mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1182  mInvTwiceScale = mScaleValuesInverse / 2;
1183  }
1184 
1186  MapBase(),
1187  mTranslation(other.mTranslation),
1188  mScaleValues(other.mScaleValues),
1189  mVoxelSize(other.mVoxelSize),
1190  mScaleValuesInverse(other.mScaleValuesInverse),
1191  mInvScaleSqr(other.mInvScaleSqr),
1192  mInvTwiceScale(other.mInvTwiceScale)
1193  {}
1194 
1196 
1200  MapBase::Ptr copy() const { return MapBase::Ptr(new ScaleTranslateMap(*this)); }
1201 
1202  static bool isRegistered() { return MapRegistry::isRegistered(ScaleTranslateMap::mapType()); }
1203 
1204  static void registerMap()
1205  {
1206  MapRegistry::registerMap(
1207  ScaleTranslateMap::mapType(),
1208  ScaleTranslateMap::create);
1209  }
1210 
1211  Name type() const { return mapType(); }
1212  static Name mapType() { return Name("ScaleTranslateMap"); }
1213 
1215  bool isLinear() const { return true; }
1216 
1219  bool hasUniformScale() const {
1220  bool value;
1221  value = isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.y()));
1222  value = value && isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.z()));
1223  return value;
1224  }
1225 
1227  Vec3d applyMap(const Vec3d& in) const
1228  {
1229  return Vec3d(
1230  in.x() * mScaleValues.x() + mTranslation.x(),
1231  in.y() * mScaleValues.y() + mTranslation.y(),
1232  in.z() * mScaleValues.z() + mTranslation.z());
1233  }
1235  Vec3d applyInverseMap(const Vec3d& in) const
1236  {
1237  return Vec3d(
1238  (in.x() - mTranslation.x() ) / mScaleValues.x(),
1239  (in.y() - mTranslation.y() ) / mScaleValues.y(),
1240  (in.z() - mTranslation.z() ) / mScaleValues.z());
1241  }
1244  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
1246  Vec3d applyIJT(const Vec3d& in) const
1247  {
1248  return Vec3d(
1249  in.x() / mScaleValues.x(),
1250  in.y() / mScaleValues.y(),
1251  in.z() / mScaleValues.z());
1252  }
1254  Mat3d applyIJC(const Mat3d& in) const {
1255  Mat3d tmp;
1256  for (int i=0; i<3; i++){
1257  tmp.setRow(i, in.row(i)*mScaleValuesInverse(i));
1258  }
1259  for (int i=0; i<3; i++){
1260  tmp.setCol(i, tmp.col(i)*mScaleValuesInverse(i));
1261  }
1262  return tmp;
1263  }
1264  Mat3d applyIJC(const Mat3d& in, const Vec3d&, const Vec3d& ) const {
1265  return applyIJC(in);
1266  }
1267 
1269  double determinant(const Vec3d& ) const { return determinant(); }
1271  double determinant() const { return mScaleValues.x()*mScaleValues.y()*mScaleValues.z(); }
1273  Vec3d voxelSize() const { return mVoxelSize;}
1275 OPENVDB_DEPRECATED Vec3d voxelDimensions() const { return voxelSize();}
1278  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
1280 OPENVDB_DEPRECATED Vec3d voxelDimensions(const Vec3d&) const { return voxelSize();}
1281 
1283  const Vec3d& getScale() const { return mScaleValues; }
1285  const Vec3d& getTranslation() const { return mTranslation; }
1286 
1288  const Vec3d& getInvScaleSqr() const {return mInvScaleSqr;}
1290  const Vec3d& getInvTwiceScale() const {return mInvTwiceScale;}
1292  const Vec3d& getInvScale() const {return mScaleValuesInverse; }
1293 
1295  void read(std::istream& is)
1296  {
1297  mTranslation.read(is);
1298  mScaleValues.read(is);
1299  mVoxelSize.read(is);
1300  mScaleValuesInverse.read(is);
1301  mInvScaleSqr.read(is);
1302  mInvTwiceScale.read(is);
1303  }
1305  void write(std::ostream& os) const
1306  {
1307  mTranslation.write(os);
1308  mScaleValues.write(os);
1309  mVoxelSize.write(os);
1310  mScaleValuesInverse.write(os);
1311  mInvScaleSqr.write(os);
1312  mInvTwiceScale.write(os);
1313  }
1315  std::string str() const
1316  {
1317  std::ostringstream buffer;
1318  buffer << " - translation: " << mTranslation << std::endl;
1319  buffer << " - scale: " << mScaleValues << std::endl;
1320  buffer << " - voxel dimensions: " << mVoxelSize << std::endl;
1321  return buffer.str();
1322  }
1323 
1324  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1325 
1326  bool operator==(const ScaleTranslateMap& other) const
1327  {
1328  // ::eq() uses a tolerance
1329  if (!mScaleValues.eq(other.mScaleValues)) { return false; }
1330  if (!mTranslation.eq(other.mTranslation)) { return false; }
1331  return true;
1332  }
1333 
1334  bool operator!=(const ScaleTranslateMap& other) const { return !(*this == other); }
1335 
1337  AffineMap::Ptr getAffineMap() const
1338  {
1339  AffineMap::Ptr affineMap(new AffineMap(math::scale<Mat4d>(mScaleValues)));
1340  affineMap->accumPostTranslation(mTranslation);
1341  return affineMap;
1342  }
1343 
1345 
1346 
1347  MapBase::Ptr preRotate(double radians, Axis axis) const
1348  {
1349  AffineMap::Ptr affineMap = getAffineMap();
1350  affineMap->accumPreRotation(axis, radians);
1351  return simplify(affineMap);
1352  }
1353  MapBase::Ptr preTranslate(const Vec3d& t) const
1354  {
1355  const Vec3d& s = mScaleValues;
1356  const Vec3d scaled_trans( t.x() * s.x(),
1357  t.y() * s.y(),
1358  t.z() * s.z() );
1359  return MapBase::Ptr( new ScaleTranslateMap(mScaleValues, mTranslation + scaled_trans));
1360  }
1361 
1362  MapBase::Ptr preScale(const Vec3d& v) const;
1363 
1364  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
1365  {
1366  AffineMap::Ptr affineMap = getAffineMap();
1367  affineMap->accumPreShear(axis0, axis1, shear);
1368  return simplify(affineMap);
1369  }
1371 
1373 
1374 
1375  MapBase::Ptr postRotate(double radians, Axis axis) const
1376  {
1377  AffineMap::Ptr affineMap = getAffineMap();
1378  affineMap->accumPostRotation(axis, radians);
1379  return simplify(affineMap);
1380  }
1381  MapBase::Ptr postTranslate(const Vec3d& t) const
1382  {
1383  return MapBase::Ptr( new ScaleTranslateMap(mScaleValues, mTranslation + t));
1384  }
1385 
1386  MapBase::Ptr postScale(const Vec3d& v) const;
1387 
1388  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
1389  {
1390  AffineMap::Ptr affineMap = getAffineMap();
1391  affineMap->accumPostShear(axis0, axis1, shear);
1392  return simplify(affineMap);
1393  }
1395 
1397 OPENVDB_DEPRECATED MapBase::Ptr rotate(double radians, Axis axis) const
1398  { return preRotate(radians, axis); }
1400 OPENVDB_DEPRECATED MapBase::Ptr translate(const Vec3d& t) const
1401  { return postTranslate(t); }
1403 OPENVDB_DEPRECATED MapBase::Ptr scale(const Vec3d& v) const;
1406  { return preShear(shear, axis0, axis1); }
1407 
1408 private:
1409  Vec3d mTranslation, mScaleValues, mVoxelSize, mScaleValuesInverse,
1410  mInvScaleSqr, mInvTwiceScale;
1411 }; // class ScaleTanslateMap
1412 
1413 
1414 inline MapBase::Ptr
1415 ScaleMap::postTranslate(const Vec3d& t) const
1416 {
1417  return MapBase::Ptr(new ScaleTranslateMap(mScaleValues, t));
1418 }
1419 
1420 
1421 inline MapBase::Ptr
1422 ScaleMap::preTranslate(const Vec3d& t) const
1423 {
1424 
1425  const Vec3d& s = mScaleValues;
1426  const Vec3d scaled_trans( t.x() * s.x(),
1427  t.y() * s.y(),
1428  t.z() * s.z() );
1429  return MapBase::Ptr(new ScaleTranslateMap(mScaleValues, scaled_trans));
1430 }
1431 
1433 ScaleMap::translate(const Vec3d& t) const
1434 { return postTranslate(t); }
1435 
1436 
1440 {
1441 public:
1442  typedef boost::shared_ptr<UniformScaleTranslateMap> Ptr;
1443  typedef boost::shared_ptr<const UniformScaleTranslateMap> ConstPtr;
1444 
1446  UniformScaleTranslateMap(double scale, const Vec3d& translate):
1447  ScaleTranslateMap(Vec3d(scale,scale,scale), translate) {}
1449  ScaleTranslateMap(scale.getScale(), translate.getTranslation()) {}
1450 
1453 
1457  MapBase::Ptr copy() const { return MapBase::Ptr(new UniformScaleTranslateMap(*this)); }
1458 
1459  static bool isRegistered()
1460  {
1461  return MapRegistry::isRegistered(UniformScaleTranslateMap::mapType());
1462  }
1463 
1464  static void registerMap()
1465  {
1466  MapRegistry::registerMap(
1467  UniformScaleTranslateMap::mapType(),
1468  UniformScaleTranslateMap::create);
1469  }
1470 
1471  Name type() const { return mapType(); }
1472  static Name mapType() { return Name("UniformScaleTranslateMap"); }
1473 
1474  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1475 
1476  bool operator==(const UniformScaleTranslateMap& other) const
1477  {
1478  return ScaleTranslateMap::operator==(other);
1479  }
1480  bool operator!=(const UniformScaleTranslateMap& other) const { return !(*this == other); }
1481 
1482 
1485  MapBase::Ptr preTranslate(const Vec3d& t) const
1486  {
1487  const double scale = this->getScale().x();
1488  const Vec3d new_trans = this->getTranslation() + scale * t;
1489  return MapBase::Ptr( new UniformScaleTranslateMap(scale, new_trans));
1490  }
1491 
1494  MapBase::Ptr postTranslate(const Vec3d& t) const
1495  {
1496  const double scale = this->getScale().x();
1497  return MapBase::Ptr( new UniformScaleTranslateMap(scale, this->getTranslation() + t));
1498  }
1500 OPENVDB_DEPRECATED MapBase::Ptr translate(const Vec3d& t) const
1501  { return postTranslate(t); }
1502 
1503 }; // class UniformScaleTanslateMap
1504 
1505 
1506 
1507 
1508 inline MapBase::Ptr
1509 UniformScaleMap::postTranslate(const Vec3d& t) const
1510 {
1511  const double scale = this->getScale().x();
1512  return MapBase::Ptr(new UniformScaleTranslateMap(scale, t));
1513 }
1514 
1515 
1516 inline MapBase::Ptr
1517 UniformScaleMap::preTranslate(const Vec3d& t) const
1518 {
1519  const double scale = this->getScale().x();
1520  return MapBase::Ptr(new UniformScaleTranslateMap(scale, scale*t));
1521 }
1522 
1523 inline MapBase::Ptr
1524 TranslationMap::preScale(const Vec3d& v) const
1525 {
1526  if (isApproxEqual(v[0],v[1]) && isApproxEqual(v[0],v[2])) {
1527  return MapBase::Ptr(new UniformScaleTranslateMap(v[0], mTranslation));
1528  } else {
1529  return MapBase::Ptr(new ScaleTranslateMap(v, mTranslation));
1530  }
1531 }
1532 
1533 inline MapBase::Ptr
1534 TranslationMap::postScale(const Vec3d& v) const
1535 {
1536  if (isApproxEqual(v[0],v[1]) && isApproxEqual(v[0],v[2])) {
1537  return MapBase::Ptr(new UniformScaleTranslateMap(v[0], v[0]*mTranslation));
1538  } else {
1539  const Vec3d trans(mTranslation.x()*v.x(),
1540  mTranslation.y()*v.y(),
1541  mTranslation.z()*v.z());
1542  return MapBase::Ptr(new ScaleTranslateMap(v, trans));
1543  }
1544 }
1545 
1546 inline MapBase::Ptr
1547 ScaleTranslateMap::preScale(const Vec3d& v) const
1548 {
1549  Vec3d new_scale( v * mScaleValues );
1550  if (isApproxEqual(new_scale[0],new_scale[1]) && isApproxEqual(new_scale[0],new_scale[2])) {
1551  return MapBase::Ptr( new UniformScaleTranslateMap(new_scale[0], mTranslation));
1552  } else {
1553  return MapBase::Ptr( new ScaleTranslateMap(new_scale, mTranslation));
1554  }
1555 }
1556 
1557 inline MapBase::Ptr
1558 ScaleTranslateMap::postScale(const Vec3d& v) const
1559 {
1560  const Vec3d new_scale( v * mScaleValues );
1561  const Vec3d new_trans( mTranslation.x()*v.x(),
1562  mTranslation.y()*v.y(),
1563  mTranslation.z()*v.z() );
1564 
1565  if (isApproxEqual(new_scale[0],new_scale[1]) && isApproxEqual(new_scale[0],new_scale[2])) {
1566  return MapBase::Ptr( new UniformScaleTranslateMap(new_scale[0], new_trans));
1567  } else {
1568  return MapBase::Ptr( new ScaleTranslateMap(new_scale, new_trans));
1569  }
1570 }
1571 
1572 
1574 UniformScaleMap::translate(const Vec3d& t) const
1575 { return postTranslate(t); }
1576 
1579 { return preScale(v); }
1580 
1583 { return preScale(v); }
1584 
1585 
1587 
1591 {
1592 public:
1593  typedef boost::shared_ptr<UnitaryMap> Ptr;
1594  typedef boost::shared_ptr<const UnitaryMap> ConstPtr;
1595 
1597  UnitaryMap(): mAffineMap(Mat4d::identity())
1598  {
1599  }
1600 
1601  UnitaryMap(const Vec3d& axis, double radians)
1602  {
1603  Mat3d matrix;
1604  matrix.setToRotation(axis, radians);
1605  mAffineMap = AffineMap(matrix);
1606  }
1607 
1608  UnitaryMap(Axis axis, double radians)
1609  {
1610  Mat4d matrix;
1611  matrix.setToRotation(axis, radians);
1612  mAffineMap = AffineMap(matrix);
1613  }
1614 
1615  UnitaryMap(const Mat3d& m)
1616  {
1617  // test that the mat3 is a rotation || reflection
1618  if (!isUnitary(m)) {
1619  OPENVDB_THROW(ArithmeticError, "Matrix initializing unitary map was not unitary");
1620  }
1621 
1622  Mat4d matrix(Mat4d::identity());
1623  matrix.setMat3(m);
1624  mAffineMap = AffineMap(matrix);
1625  }
1626 
1627  UnitaryMap(const Mat4d& m)
1628  {
1629  if (!isInvertible(m)) {
1631  "4x4 Matrix initializing unitary map was not unitary: not invertible");
1632  }
1633 
1634  if (!isAffine(m)) {
1636  "4x4 Matrix initializing unitary map was not unitary: not affine");
1637  }
1638 
1639  if (hasTranslation(m)) {
1641  "4x4 Matrix initializing unitary map was not unitary: had translation");
1642  }
1643 
1644  if (!isUnitary(m.getMat3())) {
1646  "4x4 Matrix initializing unitary map was not unitary");
1647  }
1648 
1649  mAffineMap = AffineMap(m);
1650  }
1651 
1652  UnitaryMap(const UnitaryMap& other):
1653  MapBase(other),
1654  mAffineMap(other.mAffineMap)
1655  {
1656  }
1657 
1658  UnitaryMap(const UnitaryMap& first, const UnitaryMap& second):
1659  mAffineMap(*(first.getAffineMap()), *(second.getAffineMap()))
1660  {
1661  }
1662 
1665  static MapBase::Ptr create() { return MapBase::Ptr(new UnitaryMap()); }
1667  MapBase::Ptr copy() const { return MapBase::Ptr(new UnitaryMap(*this)); }
1668 
1669  static bool isRegistered() { return MapRegistry::isRegistered(UnitaryMap::mapType()); }
1670 
1671  static void registerMap()
1672  {
1673  MapRegistry::registerMap(
1674  UnitaryMap::mapType(),
1675  UnitaryMap::create);
1676  }
1677 
1679  Name type() const { return mapType(); }
1681  static Name mapType() { return Name("UnitaryMap"); }
1682 
1684  bool isLinear() const { return true; }
1685 
1687  bool hasUniformScale() const { return true; }
1688 
1689  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1690 
1691  bool operator==(const UnitaryMap& other) const
1692  {
1693  // compare underlying linear map.
1694  if (mAffineMap!=other.mAffineMap) return false;
1695  return true;
1696  }
1697 
1698  bool operator!=(const UnitaryMap& other) const { return !(*this == other); }
1700  Vec3d applyMap(const Vec3d& in) const { return mAffineMap.applyMap(in); }
1702  Vec3d applyInverseMap(const Vec3d& in) const { return mAffineMap.applyInverseMap(in); }
1705  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
1707  Vec3d applyIJT(const Vec3d& in) const { return mAffineMap.applyIJT(in); }
1709  Mat3d applyIJC(const Mat3d& in) const {
1710  return mAffineMap.applyIJC(in);
1711  }
1712  Mat3d applyIJC(const Mat3d& in, const Vec3d&, const Vec3d& ) const {
1713  return applyIJC(in);
1714  }
1716  double determinant(const Vec3d& ) const { return determinant(); }
1718  double determinant() const { return mAffineMap.determinant(); }
1719 
1720 
1725  Vec3d voxelSize() const { return mAffineMap.voxelSize();}
1726  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
1728 OPENVDB_DEPRECATED Vec3d voxelDimensions() const { return voxelSize();}
1730 OPENVDB_DEPRECATED Vec3d voxelDimensions(const Vec3d&) const { return voxelSize();}
1731 
1733  void read(std::istream& is)
1734  {
1735  mAffineMap.read(is);
1736  }
1737 
1739  void write(std::ostream& os) const
1740  {
1741  mAffineMap.write(os);
1742  }
1744  std::string str() const
1745  {
1746  std::ostringstream buffer;
1747  buffer << mAffineMap.str();
1748  return buffer.str();
1749  }
1751  AffineMap::Ptr getAffineMap() const { return AffineMap::Ptr(new AffineMap(mAffineMap)); }
1752 
1754 OPENVDB_DEPRECATED MapBase::Ptr rotate(double radians, Axis axis) const
1755  { return postRotate(radians, axis); }
1757 OPENVDB_DEPRECATED MapBase::Ptr translate(const Vec3d& t) const
1758  { return postTranslate(t); }
1761  { return preScale(v); }
1764  { return preShear(shear, axis0, axis1); }
1765 
1767 
1768 
1769  MapBase::Ptr preRotate(double radians, Axis axis) const
1770  {
1771  UnitaryMap first(axis, radians);
1772  UnitaryMap::Ptr unitaryMap(new UnitaryMap(first, *this));
1773  return boost::static_pointer_cast<MapBase, UnitaryMap>(unitaryMap);
1774  }
1775  MapBase::Ptr preTranslate(const Vec3d& t) const
1776  {
1777  AffineMap::Ptr affineMap = getAffineMap();
1778  affineMap->accumPreTranslation(t);
1779  return simplify(affineMap);
1780  }
1781  MapBase::Ptr preScale(const Vec3d& v) const
1782  {
1783  AffineMap::Ptr affineMap = getAffineMap();
1784  affineMap->accumPreScale(v);
1785  return simplify(affineMap);
1786  }
1787  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
1788  {
1789  AffineMap::Ptr affineMap = getAffineMap();
1790  affineMap->accumPreShear(axis0, axis1, shear);
1791  return simplify(affineMap);
1792  }
1794 
1795 
1797 
1798 
1799  MapBase::Ptr postRotate(double radians, Axis axis) const
1800  {
1801  UnitaryMap second(axis, radians);
1802  UnitaryMap::Ptr unitaryMap(new UnitaryMap(*this, second));
1803  return boost::static_pointer_cast<MapBase, UnitaryMap>(unitaryMap);
1804  }
1805  MapBase::Ptr postTranslate(const Vec3d& t) const
1806  {
1807  AffineMap::Ptr affineMap = getAffineMap();
1808  affineMap->accumPostTranslation(t);
1809  return simplify(affineMap);
1810  }
1811  MapBase::Ptr postScale(const Vec3d& v) const
1812  {
1813  AffineMap::Ptr affineMap = getAffineMap();
1814  affineMap->accumPostScale(v);
1815  return simplify(affineMap);
1816  }
1817  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
1818  {
1819  AffineMap::Ptr affineMap = getAffineMap();
1820  affineMap->accumPostShear(axis0, axis1, shear);
1821  return simplify(affineMap);
1822  }
1824 
1825 private:
1826  AffineMap mAffineMap;
1827 }; // class UnitaryMap
1828 
1829 
1831 
1832 
1840 {
1841 public:
1842  typedef boost::shared_ptr<NonlinearFrustumMap> Ptr;
1843  typedef boost::shared_ptr<const NonlinearFrustumMap> ConstPtr;
1844 
1846  MapBase(),
1847  mBBox(Vec3d(0), Vec3d(1)),
1848  mTaper(1),
1849  mDepth(1)
1850  {
1851  init();
1852  }
1853 
1857 OPENVDB_DEPRECATED NonlinearFrustumMap(const CoordBBox& bb, double taper, double depth):
1858  MapBase(),mBBox(bb.min().asVec3d(), bb.max().asVec3d()), mTaper(taper), mDepth(depth)
1859  {
1860  init();
1861  }
1862 
1863  NonlinearFrustumMap(const BBoxd& bb, double taper, double depth):
1864  MapBase(),mBBox(bb), mTaper(taper), mDepth(depth)
1865  {
1866  init();
1867  }
1868 
1869 
1875 OPENVDB_DEPRECATED NonlinearFrustumMap(const CoordBBox& bb, double taper, double depth,
1876  const MapBase::Ptr& secondMap):
1877  mBBox(bb.min().asVec3d(),bb.max().asVec3d()), mTaper(taper), mDepth(depth)
1878  {
1879  if (!secondMap->isLinear() ) {
1881  "The second map in the Frustum transfrom must be linear");
1882  }
1883  mSecondMap = *( secondMap->getAffineMap() );
1884  init();
1885  }
1886 
1887  NonlinearFrustumMap(const BBoxd& bb, double taper, double depth,
1888  const MapBase::Ptr& secondMap):
1889  mBBox(bb), mTaper(taper), mDepth(depth)
1890  {
1891  if (!secondMap->isLinear() ) {
1893  "The second map in the Frustum transfrom must be linear");
1894  }
1895  mSecondMap = *( secondMap->getAffineMap() );
1896  init();
1897  }
1898 
1900  MapBase(),
1901  mBBox(other.mBBox),
1902  mTaper(other.mTaper),
1903  mDepth(other.mDepth),
1904  mSecondMap(other.mSecondMap),
1905  mHasSimpleAffine(other.mHasSimpleAffine)
1906  {
1907  init();
1908  }
1909 
1910 
1926  NonlinearFrustumMap(const Vec3d& position,
1927  const Vec3d& direction,
1928  const Vec3d& up,
1929  double aspect /* width / height */,
1930  double z_near, double depth,
1931  Coord::ValueType x_count, Coord::ValueType z_count) {
1932 
1936  if (!(depth > 0)) {
1938  "The frustum depth must be non-zero and positive");
1939  }
1940  if (!(up.length() > 0)) {
1942  "The frustum height must be non-zero and positive");
1943  }
1944  if (!(aspect > 0)) {
1946  "The frustum aspect ratio must be non-zero and positive");
1947  }
1948  if (!(isApproxEqual(up.dot(direction), 0.))) {
1950  "The frustum up orientation must be perpendicular to into-frustum direction");
1951  }
1952 
1953  double near_plane_height = 2 * up.length();
1954  double near_plane_width = aspect * near_plane_height;
1955 
1956  Coord::ValueType y_count = static_cast<int>(Round(x_count / aspect));
1957 
1958  mBBox = BBoxd(Vec3d(0,0,0), Vec3d(x_count, y_count, z_count));
1959  mDepth = depth / near_plane_width; // depth non-dimensionalized on width
1960  double gamma = near_plane_width / z_near;
1961  mTaper = 1./(mDepth*gamma + 1.);
1962 
1963  Vec3d direction_unit = direction;
1964  direction_unit.normalize();
1965 
1966  Mat4d r1(Mat4d::identity());
1967  r1.setToRotation(/*from*/Vec3d(0,0,1), /*to */direction_unit);
1968  Mat4d r2(Mat4d::identity());
1969  Vec3d temp = r1.inverse().transform(up);
1970  r2.setToRotation(/*from*/Vec3d(0,1,0), /*to*/temp );
1971  Mat4d scale = math::scale<Mat4d>(
1972  Vec3d(near_plane_width, near_plane_width, near_plane_width));
1973 
1974  // move the near plane to origin, rotate to align with axis, and scale down
1975  // T_inv * R1_inv * R2_inv * scale_inv
1976  Mat4d mat = scale * r2 * r1;
1977  mat.setTranslation(position + z_near*direction_unit);
1978 
1979  mSecondMap = AffineMap(mat);
1980 
1981  init();
1982  }
1983 
1988  MapBase::Ptr copy() const { return MapBase::Ptr(new NonlinearFrustumMap(*this)); }
1989 
1990  static bool isRegistered() { return MapRegistry::isRegistered(NonlinearFrustumMap::mapType()); }
1991 
1992  static void registerMap()
1993  {
1994  MapRegistry::registerMap(
1995  NonlinearFrustumMap::mapType(),
1996  NonlinearFrustumMap::create);
1997  }
1999  Name type() const { return mapType(); }
2001  static Name mapType() { return Name("NonlinearFrustumMap"); }
2002 
2004  bool isLinear() const { return false; }
2005 
2007  bool hasUniformScale() const { return false; }
2008 
2009  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
2010 
2011  bool operator==(const NonlinearFrustumMap& other) const
2012  {
2013  if (mBBox!=other.mBBox) return false;
2014  if (!isApproxEqual(mTaper, other.mTaper)) return false;
2015  if (!isApproxEqual(mDepth, other.mDepth)) return false;
2016 
2017  // Two linear transforms are equivalent iff they have the same translation
2018  // and have the same affects on orthongal spanning basis check translation
2019  Vec3d e(0,0,0);
2020  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
2022  e(0) = 1;
2023  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
2024  e(0) = 0;
2025  e(1) = 1;
2026  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
2027  e(1) = 0;
2028  e(2) = 1;
2029  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
2030  return true;
2031  }
2032 
2033  bool operator!=(const NonlinearFrustumMap& other) const { return !(*this == other); }
2034 
2036  Vec3d applyMap(const Vec3d& in) const
2037  {
2038  return mSecondMap.applyMap(applyFrustumMap(in));
2039  }
2040 
2042  Vec3d applyInverseMap(const Vec3d& in) const
2043  {
2044  return applyFrustumInverseMap(mSecondMap.applyInverseMap(in));
2045  }
2046 
2048  Vec3d applyIJT(const Vec3d& in) const { return mSecondMap.applyIJT(in); }
2049 
2050  // the Jacobian of the nonlinear part of the transform is a sparse matrix
2051  // Jacobian^(-T) =
2052  //
2053  // (Lx)( 1/s 0 0 )
2054  // ( 0 1/s 0 )
2055  // ( -(x-xo)g/(sLx) -(y-yo)g/(sLx) Lz/(Depth Lx) )
2058  Vec3d applyIJT(const Vec3d& d1_is, const Vec3d& ijk) const
2059  {
2060 
2061  const Vec3d loc = applyFrustumMap(ijk);
2062  const double s = mGamma * loc.z() + 1.;
2063 
2064  // verify that we aren't at the singularity
2065  if (isApproxEqual(s, 0.)) {
2067  "Tried to evaluate the frustum transform at the singular focal point (e.g. camera)");
2068  }
2069 
2070  const double sinv = 1.0/s; // 1/(z*gamma + 1)
2071  const double pt0 = mLx * sinv; // Lx / (z*gamma +1)
2072  const double pt1 = mGamma * pt0; // gamma * Lx / ( z*gamma +1)
2073  const double pt2 = pt1 * sinv; // gamma * Lx / ( z*gamma +1)**2
2074 
2075  const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
2076 
2077  // compute \frac{\partial E_i}{\partial x_j}
2078  Mat3d gradE(Mat3d::zero());
2079  for (int j = 0; j < 3; ++j ) {
2080  gradE(0,j) = pt0 * jacinv(0,j) - pt2 * loc.x()*jacinv(2,j);
2081  gradE(1,j) = pt0 * jacinv(1,j) - pt2 * loc.y()*jacinv(2,j);
2082  gradE(2,j) = (1./mDepthOnLz) * jacinv(2,j);
2083  }
2084 
2085  Vec3d result;
2086  for (int i = 0; i < 3; ++i) {
2087  result(0) = d1_is(0) * gradE(0,i) + d1_is(1) * gradE(1,i) + d1_is(2) * gradE(2,i);
2088  }
2089 
2090  return result;
2091 
2092  }
2093 
2095  Mat3d applyIJC(const Mat3d& in) const { return mSecondMap.applyIJC(in); }
2100  Mat3d applyIJC(const Mat3d& d2_is, const Vec3d& d1_is, const Vec3d& ijk) const
2101  {
2102  const Vec3d loc = applyFrustumMap(ijk);
2103 
2104  const double s = mGamma * loc.z() + 1.;
2105 
2106  // verify that we aren't at the singularity
2107  if (isApproxEqual(s, 0.)) {
2109  "Tried to evaluate the frustum transform at the singular focal point (e.g. camera)");
2110  }
2111 
2112  // precompute
2113  const double sinv = 1.0/s; // 1/(z*gamma + 1)
2114  const double pt0 = mLx * sinv; // Lx / (z*gamma +1)
2115  const double pt1 = mGamma * pt0; // gamma * Lx / ( z*gamma +1)
2116  const double pt2 = pt1 * sinv; // gamma * Lx / ( z*gamma +1)**2
2117  const double pt3 = pt2 * sinv; // gamma * Lx / ( z*gamma +1)**3
2118 
2119 
2120  const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
2121 
2122  // compute \frac{\partial^2 E_i}{\partial x_j \partial x_k}
2123 
2124  Mat3d matE0(Mat3d::zero());
2125  Mat3d matE1(Mat3d::zero()); // matE2 = 0
2126  for(int j = 0; j < 3; j++) {
2127  for (int k = 0; k < 3; k++) {
2128 
2129  const double pt4 = 2. * jacinv(2,j) * jacinv(2,k) * pt3;
2130 
2131  matE0(j,k) = -(jacinv(0,j) * jacinv(2,k) + jacinv(2,j) * jacinv(0,k)) * pt2 +
2132  pt4 * loc.x();
2133 
2134  matE1(j,k) = -(jacinv(1,j) * jacinv(2,k) + jacinv(2,j) * jacinv(1,k)) * pt2 +
2135  pt4 * loc.y();
2136  }
2137  }
2138 
2139  // compute \frac{\partial E_i}{\partial x_j}
2140  Mat3d gradE(Mat3d::zero());
2141  for (int j = 0; j < 3; ++j ) {
2142  gradE(0,j) = pt0 * jacinv(0,j) - pt2 * loc.x()*jacinv(2,j);
2143  gradE(1,j) = pt0 * jacinv(1,j) - pt2 * loc.y()*jacinv(2,j);
2144  gradE(2,j) = (1./mDepthOnLz) * jacinv(2,j);
2145  }
2146 
2147  Mat3d result(Mat3d::zero());
2148  // compute \fac{\partial E_j}{\partial x_m} \fac{\partial E_i}{\partial x_n}
2149  // \frac{\partial^2 input}{\partial E_i \partial E_j}
2150  for (int m = 0; m < 3; ++m ) {
2151  for ( int n = 0; n < 3; ++n) {
2152  for (int i = 0; i < 3; ++i ) {
2153  for (int j = 0; j < 3; ++j) {
2154  result(m, n) += gradE(j, m) * gradE(i, n) * d2_is(i, j);
2155  }
2156  }
2157  }
2158  }
2159 
2160  for (int m = 0; m < 3; ++m ) {
2161  for ( int n = 0; n < 3; ++n) {
2162  result(m, n) += matE0(m, n) * d1_is(0) + matE1(m, n) * d1_is(1);// + matE2(m, n) * d1_is(2);
2163  }
2164  }
2165 
2166  return result;
2167  }
2168 
2170  double determinant() const {return mSecondMap.determinant();} // no implementation
2171 
2174  double determinant(const Vec3d& loc) const
2175  {
2176  double s = mGamma * loc.z() + 1.0;
2177  double frustum_determinant = s * s * mDepthOnLzLxLx;
2178  return mSecondMap.determinant() * frustum_determinant;
2179  }
2180 
2182  Vec3d voxelSize() const {
2183  const Vec3d loc( 0.5*(mBBox.min().x() + mBBox.max().x()),
2184  0.5*(mBBox.min().y() + mBBox.max().y()),
2185  mBBox.min().z());
2186 
2187  return voxelSize(loc);
2188 
2189  }
2190 
2195  Vec3d voxelSize(const Vec3d& loc) const {
2196  Vec3d out, pos = applyMap(loc);
2197  out(0) = (applyMap(loc + Vec3d(1,0,0)) - pos).length();
2198  out(1) = (applyMap(loc + Vec3d(0,1,0)) - pos).length();
2199  out(2) = (applyMap(loc + Vec3d(0,0,1)) - pos).length();
2200  return out;
2201  }
2202 
2204  OPENVDB_DEPRECATED Vec3d voxelDimensions() const { return voxelSize(); }
2206  OPENVDB_DEPRECATED Vec3d voxelDimensions(const Vec3d& loc) const { return voxelSize(loc); }
2207 
2208  AffineMap::Ptr getAffineMap() const
2209  {
2210  return mSecondMap.getAffineMap();
2211  }
2212 
2214  void setTaper(double t) { mTaper = t; init();}
2216  double getTaper() const { return mTaper; }
2218  void setDepth(double d) { mDepth = d; init();}
2220  double getDepth() const { return mDepth; }
2221  // gamma a non-dimensional number: nearplane x-width / camera to near plane distance
2222  double getGamma() const { return mGamma; }
2223 
2226 OPENVDB_DEPRECATED void setBBox(const CoordBBox& bb) { mBBox = BBoxd(bb.min().asVec3d(), bb.max().asVec3d()); }
2228  const BBoxd& getBBox() const { return mBBox; }
2229 
2231 OPENVDB_DEPRECATED MapBase::Ptr linearMap() const { return getAffineMap(); }
2232  const AffineMap& secondMap() const { return mSecondMap; }
2235  bool isValid() const { return !mBBox.empty();}
2236 
2238  bool hasSimpleAffine() const { return mHasSimpleAffine; }
2239 
2241  void read(std::istream& is)
2242  {
2243  // for backward compatibility with earlier version
2245  CoordBBox bb;
2246  bb.read(is);
2247  mBBox = BBoxd(bb.min().asVec3d(), bb.max().asVec3d());
2248  } else {
2249  mBBox.read(is);
2250  }
2251 
2252  is.read(reinterpret_cast<char*>(&mTaper), sizeof(double));
2253  is.read(reinterpret_cast<char*>(&mDepth), sizeof(double));
2254 
2255  // Read the second maps type.
2256  Name type = readString(is);
2257 
2258  // Check if the map has been registered.
2259  if(!MapRegistry::isRegistered(type)) {
2260  OPENVDB_THROW(KeyError, "Map " << type << " is not registered");
2261  }
2262 
2263  // Create the second map of the type and then read it in.
2264  MapBase::Ptr proxy = math::MapRegistry::createMap(type);
2265  proxy->read(is);
2266  mSecondMap = *(proxy->getAffineMap());
2267  init();
2268  }
2269 
2271  void write(std::ostream& os) const
2272  {
2273  mBBox.write(os);
2274  os.write(reinterpret_cast<const char*>(&mTaper), sizeof(double));
2275  os.write(reinterpret_cast<const char*>(&mDepth), sizeof(double));
2276 
2277  writeString(os, mSecondMap.type());
2278  mSecondMap.write(os);
2279  }
2280 
2282  std::string str() const
2283  {
2284  std::ostringstream buffer;
2285  buffer << " - taper: " << mTaper << std::endl;
2286  buffer << " - depth: " << mDepth << std::endl;
2287  buffer << " SecondMap: "<< mSecondMap.type() << std::endl;
2288  buffer << mSecondMap.str() << std::endl;
2289  return buffer.str();
2290  }
2291 
2293 
2294 
2295  MapBase::Ptr preRotate(double radians, Axis axis = X_AXIS) const
2296  {
2297  return MapBase::Ptr(
2298  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.preRotate(radians, axis)));
2299  }
2300  MapBase::Ptr preTranslate(const Vec3d& t) const
2301  {
2302  return MapBase::Ptr(
2303  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.preTranslate(t)));
2304  }
2305  MapBase::Ptr preScale(const Vec3d& s) const
2306  {
2307  return MapBase::Ptr(
2308  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.preScale(s)));
2309  }
2310  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
2311  {
2312  return MapBase::Ptr(new NonlinearFrustumMap(
2313  mBBox, mTaper, mDepth, mSecondMap.preShear(shear, axis0, axis1)));
2314  }
2316 
2318 
2319 
2320  MapBase::Ptr postRotate(double radians, Axis axis = X_AXIS) const
2321  {
2322  return MapBase::Ptr(
2323  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.postRotate(radians, axis)));
2324  }
2325  MapBase::Ptr postTranslate(const Vec3d& t) const
2326  {
2327  return MapBase::Ptr(
2328  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.postTranslate(t)));
2329  }
2330  MapBase::Ptr postScale(const Vec3d& s) const
2331  {
2332  return MapBase::Ptr(
2333  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.postScale(s)));
2334  }
2335  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
2336  {
2337  return MapBase::Ptr(new NonlinearFrustumMap(
2338  mBBox, mTaper, mDepth, mSecondMap.postShear(shear, axis0, axis1)));
2339  }
2341 
2343 OPENVDB_DEPRECATED MapBase::Ptr rotate(double radians, Axis axis = X_AXIS) const
2344  { return preRotate(radians, axis); }
2346 OPENVDB_DEPRECATED MapBase::Ptr translate(const Vec3d& t) const
2347  { return postTranslate(t); }
2350  { return preScale(s); }
2353  { return preShear(shear, axis0, axis1); }
2354 
2355 
2356 private:
2357  void init() {
2358  // set up as a frustum
2359  mLx = mBBox.extents().x();
2360  mLy = mBBox.extents().y();
2361  mLz = mBBox.extents().z();
2362 
2363  if (isApproxEqual(mLx,0.) || isApproxEqual(mLy,0.) || isApproxEqual(mLz,0.) ) {
2365  "The index space bounding box must have at least two index points in each direction.");
2366  }
2367 
2368  mXo = 0.5* mLx;
2369  mYo = 0.5* mLy;
2370 
2371  // mDepth is non-dimensionalized on near
2372  mGamma = (1./mTaper - 1) / mDepth;
2373 
2374  mDepthOnLz = mDepth/mLz;
2375  mDepthOnLzLxLx = mDepthOnLz/(mLx * mLx);
2376 
2378  mHasSimpleAffine = true;
2379  Vec3d tmp = mSecondMap.voxelSize();
2380 
2382  if (!isApproxEqual(tmp(0), tmp(1))) { mHasSimpleAffine = false; return; }
2383  if (!isApproxEqual(tmp(0), tmp(2))) { mHasSimpleAffine = false; return; }
2384 
2385  Vec3d trans = mSecondMap.applyMap(Vec3d(0,0,0));
2387  Vec3d tmp1 = mSecondMap.applyMap(Vec3d(1,0,0)) - trans;
2388  Vec3d tmp2 = mSecondMap.applyMap(Vec3d(0,1,0)) - trans;
2389  Vec3d tmp3 = mSecondMap.applyMap(Vec3d(0,0,1)) - trans;
2390 
2392  if (!isApproxEqual(tmp1.dot(tmp2), 0., 1.e-7)) { mHasSimpleAffine = false; return; }
2393  if (!isApproxEqual(tmp2.dot(tmp3), 0., 1.e-7)) { mHasSimpleAffine = false; return; }
2394  if (!isApproxEqual(tmp3.dot(tmp1), 0., 1.e-7)) { mHasSimpleAffine = false; return; }
2395  }
2396 
2397  Vec3d applyFrustumMap(const Vec3d& in) const
2398  {
2399 
2400  // Move the center of the x-face of the bbox
2401  // to the origin in index space.
2402  Vec3d out(in);
2403  out = out - mBBox.min();
2404  out.x() -= mXo;
2405  out.y() -= mYo;
2406 
2407  // scale the z-direction on depth / K count
2408  out.z() *= mDepthOnLz;
2409 
2410  double scale = (mGamma * out.z() + 1.)/ mLx;
2411 
2412  // scale the x-y on the length I count and apply tapper
2413  out.x() *= scale ;
2414  out.y() *= scale ;
2415 
2416  return out;
2417  }
2418 
2419  Vec3d applyFrustumInverseMap(const Vec3d& in) const
2420  {
2421  // invert taper and resize: scale = 1/( (z+1)/2 (mt-1) + 1)
2422  Vec3d out(in);
2423  double invScale = mLx / (mGamma * out.z() + 1.);
2424  out.x() *= invScale;
2425  out.y() *= invScale;
2426 
2427  out.x() += mXo;
2428  out.y() += mYo;
2429 
2430  out.z() /= mDepthOnLz;
2431 
2432  // move back
2433  out = out + mBBox.min();
2434  return out;
2435  }
2436 
2437  // bounding box in index space used in Frustum transforms.
2438  BBoxd mBBox;
2439 
2440  // taper value used in constructing Frustums.
2441  double mTaper;
2442  double mDepth;
2443 
2444  // defines the second map
2445  AffineMap mSecondMap;
2446 
2447  // these are derived from the above.
2448  double mLx, mLy, mLz;
2449  double mXo, mYo, mGamma, mDepthOnLz, mDepthOnLzLxLx;
2450 
2451  // true: if the mSecondMap is linear and has no shear, and has no non-uniform scale
2452  bool mHasSimpleAffine;
2453 }; // class NonlinearFrustumMap
2454 
2455 
2457 
2458 
2462 template<typename FirstMapType, typename SecondMapType>
2463 class CompoundMap
2464 {
2465 public:
2467 
2468  typedef boost::shared_ptr<MyType> Ptr;
2469  typedef boost::shared_ptr<const MyType> ConstPtr;
2470 
2471 
2472  CompoundMap() { updateAffineMatrix(); }
2473 
2474  CompoundMap(const FirstMapType& f, const SecondMapType& s): mFirstMap(f), mSecondMap(s)
2475  {
2476  updateAffineMatrix();
2477  }
2478 
2479  CompoundMap(const MyType& other):
2480  mFirstMap(other.mFirstMap),
2481  mSecondMap(other.mSecondMap),
2482  mAffineMap(other.mAffineMap)
2483  {}
2484 
2485  Name type() const { return mapType(); }
2486  static Name mapType()
2487  {
2488  return (FirstMapType::mapType() + Name(":") + SecondMapType::mapType());
2489  }
2490 
2491  bool operator==(const MyType& other) const
2492  {
2493  if (mFirstMap != other.mFirstMap) return false;
2494  if (mSecondMap != other.mSecondMap) return false;
2495  if (mAffineMap != other.mAffineMap) return false;
2496  return true;
2497  }
2498 
2499  bool operator!=(const MyType& other) const { return !(*this == other); }
2500 
2501  MyType& operator=(const MyType& other)
2502  {
2503  mFirstMap = other.mFirstMap;
2504  mSecondMap = other.mSecondMap;
2505  mAffineMap = other.mAffineMap;
2506  return *this;
2507  }
2508 
2509  bool isIdentity() const {
2511  return mAffineMap.isIdentity();
2512  } else {
2513  return mFirstMap.isIdentity()&&mSecondMap.isIdentity();
2514  }
2515  }
2516 
2517  bool isDiagonal() const {
2519  return mAffineMap.isDiagonal();
2520  } else {
2521  return mFirstMap.isDiagonal()&&mSecondMap.isDiagonal();
2522  }
2523  }
2524  AffineMap::Ptr getAffineMap() const
2525  {
2527  AffineMap::Ptr affine(new AffineMap(mAffineMap));
2528  return affine;
2529  } else {
2531  "Constant affine matrix representation not possible for this nonlinear map");
2532  }
2533  }
2534 
2535  // direct decompotion
2536  const FirstMapType& firstMap() const { return mFirstMap; }
2537  const SecondMapType& secondMap() const {return mSecondMap; }
2538 
2539  void setFirstMap(const FirstMapType& first) { mFirstMap = first; updateAffineMatrix(); }
2540  void setSecondMap(const SecondMapType& second) { mSecondMap = second; updateAffineMatrix(); }
2541 
2542  void read(std::istream& is)
2543  {
2544  mAffineMap.read(is);
2545  mFirstMap.read(is);
2546  mSecondMap.read(is);
2547  }
2548  void write(std::ostream& os) const
2549  {
2550  mAffineMap.write(os);
2551  mFirstMap.write(os);
2552  mSecondMap.write(os);
2553  }
2554 
2555 private:
2556  void updateAffineMatrix()
2557  {
2559  // both maps need to be linear, these methods are only defined for linear maps
2560  AffineMap::Ptr first = mFirstMap.getAffineMap();
2561  AffineMap::Ptr second= mSecondMap.getAffineMap();
2562  mAffineMap = AffineMap(*first, *second);
2563  }
2564  }
2565 
2566  FirstMapType mFirstMap;
2567  SecondMapType mSecondMap;
2568  // used for acceleration
2569  AffineMap mAffineMap;
2570 }; // class CompoundMap
2571 
2572 
2573 } // namespace math
2574 } // namespace OPENVDB_VERSION_NAME
2575 } // namespace openvdb
2576 
2577 #endif // OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED
2578 
2579 // Copyright (c) 2012 DreamWorks Animation LLC
2580 // All rights reserved. This software is distributed under the
2581 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )