OpenVDB  1.2.0
Maps.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 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 "Math.h"
37 #include "Mat4.h"
38 #include "Vec3.h"
39 #include "BBox.h"
40 #include "Coord.h"
41 #include <openvdb/util/Name.h>
42 #include <openvdb/Types.h>
43 #include <boost/shared_ptr.hpp>
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 
152 
153 
155 
156 
159 {
160 public:
161  typedef boost::shared_ptr<MapBase> Ptr;
162  typedef boost::shared_ptr<const MapBase> ConstPtr;
163  typedef Ptr (*MapFactory)();
164 
165  virtual ~MapBase(){}
166 
167  virtual boost::shared_ptr<AffineMap> getAffineMap() const = 0;
168 
170  virtual Name type() const = 0;
171 
173  template<typename MapT> bool isType() const { return this->type() == MapT::mapType(); }
174 
176  virtual bool isEqual(const MapBase& other) const = 0;
177 
179  virtual bool isLinear() const = 0;
181  virtual bool hasUniformScale() const = 0;
182 
183  virtual Vec3d applyMap(const Vec3d& in) const = 0;
184  virtual Vec3d applyInverseMap(const Vec3d& in) const = 0;
185 
186  virtual Vec3d applyIJT(const Vec3d& in) const = 0;
187  virtual Vec3d applyIJT(const Vec3d& in, const Vec3d& pos) const = 0;
188  virtual Mat3d applyIJC(const Mat3d& m) const = 0;
189  virtual Mat3d applyIJC(const Mat3d& m, const Vec3d& v, const Vec3d& pos) const = 0;
190 
191 
192  virtual double determinant() const = 0;
193  virtual double determinant(const Vec3d&) const = 0;
194 
195 
197  virtual Vec3d voxelSize() const = 0;
199  virtual Vec3d voxelSize(const Vec3d&) const = 0;
201 
202  virtual void read(std::istream&) = 0;
203  virtual void write(std::ostream&) const = 0;
204 
205  virtual std::string str() const = 0;
206 
207  virtual MapBase::Ptr copy() const = 0;
208 
210  virtual MapBase::Ptr preRotate(double radians, Axis axis = X_AXIS) const = 0;
212  virtual MapBase::Ptr preTranslate(const Vec3d&) const = 0;
213  virtual MapBase::Ptr preScale(const Vec3d&) const = 0;
214  virtual MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const = 0;
215 
216  virtual MapBase::Ptr postRotate(double radians, Axis axis = X_AXIS) const = 0;
217  virtual MapBase::Ptr postTranslate(const Vec3d&) const = 0;
218  virtual MapBase::Ptr postScale(const Vec3d&) const = 0;
219  virtual MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const = 0;
221 
222 protected:
223  MapBase() {}
224 
225  template<typename MapT>
226  static bool isEqualBase(const MapT& self, const MapBase& other)
227  {
228  return other.isType<MapT>() && (self == *static_cast<const MapT*>(&other));
229  }
230 };
231 
232 
234 
235 
239 {
240 public:
241  typedef std::map<Name, MapBase::MapFactory> MapDictionary;
242 
243  static MapRegistry* instance();
244 
246  static MapBase::Ptr createMap(const Name&);
247 
249  static bool isRegistered(const Name&);
250 
252  static void registerMap(const Name&, MapBase::MapFactory);
253 
255  static void unregisterMap(const Name&);
256 
258  static void clear();
259 
260 private:
261  MapRegistry() {}
262 
263  static MapRegistry* staticInstance();
264 
265  static MapRegistry* mInstance;
266 
267  MapDictionary mMap;
268 };
269 
270 
272 
273 
277 {
278 public:
279  typedef boost::shared_ptr<AffineMap> Ptr;
280  typedef boost::shared_ptr<const AffineMap> ConstPtr;
281 
283  mMatrix(Mat4d::identity()),
284  mMatrixInv(Mat4d::identity()),
285  mJacobianInv(Mat3d::identity()),
286  mDeterminant(1),
287  mVoxelSize(Vec3d(1,1,1)),
288  mIsDiagonal(true),
289  mIsIdentity(true)
290  // the default constructor for translation is zero
291  {
292  }
293 
294  AffineMap(const Mat3d& m)
295  {
296  Mat4d mat4(Mat4d::identity());
297  mat4.setMat3(m);
298  mMatrix = mat4;
299  updateAcceleration();
300  }
301 
302  AffineMap(const Mat4d& m): mMatrix(m)
303  {
304  if (!isAffine(m)) {
306  "Tried to initialize an affine transform from a non-affine 4x4 matrix");
307  }
308  updateAcceleration();
309  }
310 
311  AffineMap(const AffineMap& other):
312  MapBase(other),
313  mMatrix(other.mMatrix),
314  mMatrixInv(other.mMatrixInv),
315  mJacobianInv(other.mJacobianInv),
316  mDeterminant(other.mDeterminant),
317  mVoxelSize(other.mVoxelSize),
318  mIsDiagonal(other.mIsDiagonal),
319  mIsIdentity(other.mIsIdentity)
320  {
321  }
322 
324  AffineMap(const AffineMap& first, const AffineMap& second):
325  mMatrix(first.mMatrix * second.mMatrix)
326  {
327  updateAcceleration();
328  }
329 
331 
333  static MapBase::Ptr create() { return MapBase::Ptr(new AffineMap()); }
335  MapBase::Ptr copy() const { return MapBase::Ptr(new AffineMap(*this)); }
336 
337  static bool isRegistered() { return MapRegistry::isRegistered(AffineMap::mapType()); }
338 
339  static void registerMap()
340  {
341  MapRegistry::registerMap(
342  AffineMap::mapType(),
343  AffineMap::create);
344  }
345 
346  Name type() const { return mapType(); }
347  static Name mapType() { return Name("AffineMap"); }
348 
350  bool isLinear() const { return true; }
351 
353  bool hasUniformScale() const
354  {
355  Mat3d mat = mMatrix.getMat3();
356  const double det = mat.det();
357  if (isApproxEqual(det, double(0))) {
358  return false;
359  } else {
360  mat *= (1.f / pow(std::abs(det),1./3.));
361  return isUnitary(mat);
362  }
363  }
364 
365  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
366 
367  bool operator==(const AffineMap& other) const
368  {
369  // the Mat.eq() is approximate
370  if (!mMatrix.eq(other.mMatrix)) { return false; }
371  if (!mMatrixInv.eq(other.mMatrixInv)) { return false; }
372  return true;
373  }
374 
375  bool operator!=(const AffineMap& other) const { return !(*this == other); }
376 
378  {
379  mMatrix = other.mMatrix;
380  mMatrixInv = other.mMatrixInv;
381 
382  mJacobianInv = other.mJacobianInv;
383  mDeterminant = other.mDeterminant;
384  mVoxelSize = other.mVoxelSize;
385  mIsDiagonal = other.mIsDiagonal;
386  mIsIdentity = other.mIsIdentity;
387  return *this;
388  }
390  Vec3d applyMap(const Vec3d& in) const { return in * mMatrix; }
392  Vec3d applyInverseMap(const Vec3d& in) const {return in * mMatrixInv; }
393 
395  Vec3d applyJacobian(const Vec3d& in, const Vec3d&) const { return applyJacobian(in); }
397  Vec3d applyJacobian(const Vec3d& in) const {
398  const double* m = mMatrix.asPointer();
399  return Vec3d( m[ 0] * in[0] + m[ 4] * in[1] + m[ 8] * in[2],
400  m[ 1] * in[0] + m[ 5] * in[1] + m[ 9] * in[2],
401  m[ 2] * in[0] + m[ 6] * in[1] + m[10] * in[2] );
402  }
403 
405  Vec3d applyIJT(const Vec3d& in, const Vec3d&) const { return applyIJT(in); }
407  Vec3d applyIJT(const Vec3d& in) const { return in * mJacobianInv; }
409  Mat3d applyIJC(const Mat3d& m) const {
410  return mJacobianInv.transpose()* m * mJacobianInv;
411  }
412  Mat3d applyIJC(const Mat3d& in, const Vec3d& , const Vec3d& ) const {
413  return applyIJC(in);
414  }
416  double determinant(const Vec3d& ) const { return determinant(); }
418  double determinant() const { return mDeterminant; }
419 
421  Vec3d voxelSize() const { return mVoxelSize; }
424  Vec3d voxelSize(const Vec3d&) const { return voxelSize(); }
426 
428  bool isIdentity() const { return mIsIdentity; }
430  bool isDiagonal() const { return mIsDiagonal; }
432  bool isScale() const { return isDiagonal(); }
434  bool isScaleTranslate() const { return math::isDiagonal(mMatrix.getMat3()); }
435 
436 
437  // Methods that modify the existing affine map
438 
440  void accumPreRotation(Axis axis, double radians)
442  {
443  mMatrix.preRotate(axis, radians);
444  updateAcceleration();
445  }
446  void accumPreScale(const Vec3d& v)
447  {
448  mMatrix.preScale(v);
449  updateAcceleration();
450  }
451  void accumPreTranslation(const Vec3d& v)
452  {
453  mMatrix.preTranslate(v);
454  updateAcceleration();
455  }
456  void accumPreShear(Axis axis0, Axis axis1, double shear)
457  {
458  mMatrix.preShear(axis0, axis1, shear);
459  updateAcceleration();
460  }
462 
463 
465  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  }
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  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  }
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  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  }
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 
584  Mat4d getMat4() const { return mMatrix;}
585  const Mat4d& getConstMat4() const {return mMatrix;}
586  const Mat3d& getConstJacobianInv() const {return mJacobianInv;}
587 
588 private:
589  void updateAcceleration() {
590  mDeterminant = mMatrix.getMat3().det();
591 
592  if (std::abs(mDeterminant) < (3.0 * tolerance<double>::value())) {
594  "Tried to initialize an affine transform from a nearly singular matrix");
595  }
596  mMatrixInv = mMatrix.inverse();
597  mJacobianInv = mMatrixInv.getMat3().transpose();
598  mIsDiagonal = math::isDiagonal(mMatrix);
599  mIsIdentity = math::isIdentity(mMatrix);
600  Vec3d pos = applyMap(Vec3d(0,0,0));
601  mVoxelSize(0) = (applyMap(Vec3d(1,0,0)) - pos).length();
602  mVoxelSize(1) = (applyMap(Vec3d(0,1,0)) - pos).length();
603  mVoxelSize(2) = (applyMap(Vec3d(0,0,1)) - pos).length();
604  }
605 
606  // the underlying matrix
607  Mat4d mMatrix;
608 
609  // stored for acceleration
610  Mat4d mMatrixInv;
611  Mat3d mJacobianInv;
612  double mDeterminant;
613  Vec3d mVoxelSize;
614  bool mIsDiagonal, mIsIdentity;
615 }; // class AffineMap
616 
617 
619 
620 
624 {
625 public:
626  typedef boost::shared_ptr<ScaleMap> Ptr;
627  typedef boost::shared_ptr<const ScaleMap> ConstPtr;
628 
629  ScaleMap(): MapBase(), mScaleValues(Vec3d(1,1,1)), mVoxelSize(Vec3d(1,1,1)),
630  mInvScaleSqr(1,1,1), mInvTwiceScale(0.5,0.5,0.5){}
631 
633  MapBase(),
634  mScaleValues(scale),
635  mVoxelSize(Vec3d(std::abs(scale(0)),std::abs(scale(1)), std::abs(scale(2))))
636  {
637  double determinant = scale[0]* scale[1] * scale[2];
638  if (std::abs(determinant) < 3.0 * tolerance<double>::value()) {
639  OPENVDB_THROW(ArithmeticError, "Non-zero scale values required");
640  }
641  mScaleValuesInverse = 1.0 / mScaleValues;
642  mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
643  mInvTwiceScale = mScaleValuesInverse / 2;
644  }
645 
646  ScaleMap(const ScaleMap& other):
647  MapBase(),
648  mScaleValues(other.mScaleValues),
649  mVoxelSize(other.mVoxelSize),
650  mScaleValuesInverse(other.mScaleValuesInverse),
651  mInvScaleSqr(other.mInvScaleSqr),
652  mInvTwiceScale(other.mInvTwiceScale)
653  {
654  }
655 
657 
659  static MapBase::Ptr create() { return MapBase::Ptr(new ScaleMap()); }
661  MapBase::Ptr copy() const { return MapBase::Ptr(new ScaleMap(*this)); }
662 
663  static bool isRegistered() { return MapRegistry::isRegistered(ScaleMap::mapType()); }
664 
665  static void registerMap()
666  {
667  MapRegistry::registerMap(
668  ScaleMap::mapType(),
669  ScaleMap::create);
670  }
671 
672  Name type() const { return mapType(); }
673  static Name mapType() { return Name("ScaleMap"); }
674 
676  bool isLinear() const { return true; }
677 
679  bool hasUniformScale() const
680  {
681  bool value = isApproxEqual(
682  std::abs(mScaleValues.x()), std::abs(mScaleValues.y()), double(5e-7));
683  value = value && isApproxEqual(
684  std::abs(mScaleValues.x()), std::abs(mScaleValues.z()), double(5e-7));
685  return value;
686  }
687 
689  Vec3d applyMap(const Vec3d& in) const
690  {
691  return Vec3d(
692  in.x() * mScaleValues.x(),
693  in.y() * mScaleValues.y(),
694  in.z() * mScaleValues.z());
695  }
697  Vec3d applyInverseMap(const Vec3d& in) const
698  {
699  return Vec3d(
700  in.x() * mScaleValuesInverse.x(),
701  in.y() * mScaleValuesInverse.y(),
702  in.z() * mScaleValuesInverse.z());
703  }
705  Vec3d applyJacobian(const Vec3d& in, const Vec3d&) const { return applyJacobian(in); }
707  Vec3d applyJacobian(const Vec3d& in) const { return applyMap(in); }
708 
709 
712  Vec3d applyIJT(const Vec3d& in, const Vec3d&) const { return applyIJT(in);}
714  Vec3d applyIJT(const Vec3d& in) const { return applyInverseMap(in); }
716  Mat3d applyIJC(const Mat3d& in) const
717  {
718  Mat3d tmp;
719  for (int i = 0; i < 3; i++) {
720  tmp.setRow(i, in.row(i) * mScaleValuesInverse(i));
721  }
722  for (int i = 0; i < 3; i++) {
723  tmp.setCol(i, tmp.col(i) * mScaleValuesInverse(i));
724  }
725  return tmp;
726  }
727  Mat3d applyIJC(const Mat3d& in, const Vec3d&, const Vec3d&) const { return applyIJC(in); }
729  double determinant(const Vec3d&) const { return determinant(); }
731  double determinant() const { return mScaleValues.x() * mScaleValues.y() * mScaleValues.z(); }
732 
734  const Vec3d& getScale() const {return mScaleValues;}
735 
737  const Vec3d& getInvScaleSqr() const { return mInvScaleSqr; }
739  const Vec3d& getInvTwiceScale() const { return mInvTwiceScale; }
741  const Vec3d& getInvScale() const { return mScaleValuesInverse; }
742 
744  Vec3d voxelSize() const { return mVoxelSize; }
749  Vec3d voxelSize(const Vec3d&) const { return voxelSize(); }
751 
753  void read(std::istream& is)
754  {
755  mScaleValues.read(is);
756  mVoxelSize.read(is);
757  mScaleValuesInverse.read(is);
758  mInvScaleSqr.read(is);
759  mInvTwiceScale.read(is);
760  }
762  void write(std::ostream& os) const
763  {
764  mScaleValues.write(os);
765  mVoxelSize.write(os);
766  mScaleValuesInverse.write(os);
767  mInvScaleSqr.write(os);
768  mInvTwiceScale.write(os);
769  }
771  std::string str() const
772  {
773  std::ostringstream buffer;
774  buffer << " - scale: " << mScaleValues << std::endl;
775  buffer << " - voxel dimensions: " << mVoxelSize << std::endl;
776  return buffer.str();
777  }
778 
779  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
780 
781  bool operator==(const ScaleMap& other) const
782  {
783  // ::eq() uses a tolerance
784  if (!mScaleValues.eq(other.mScaleValues)) { return false; }
785  return true;
786  }
787 
788  bool operator!=(const ScaleMap& other) const { return !(*this == other); }
789 
792  {
793  return AffineMap::Ptr(new AffineMap(math::scale<Mat4d>(mScaleValues)));
794  }
795 
796 
797 
799  MapBase::Ptr preRotate(double radians, Axis axis) const
802  {
803  AffineMap::Ptr affineMap = getAffineMap();
804  affineMap->accumPreRotation(axis, radians);
805  return simplify(affineMap);
806  }
807 
808  MapBase::Ptr preTranslate(const Vec3d& tr) const;
809 
810  MapBase::Ptr preScale(const Vec3d& v) const;
811 
812  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
813  {
814  AffineMap::Ptr affineMap = getAffineMap();
815  affineMap->accumPreShear(axis0, axis1, shear);
816  return simplify(affineMap);
817  }
819 
820 
822  MapBase::Ptr postRotate(double radians, Axis axis) const
825  {
826  AffineMap::Ptr affineMap = getAffineMap();
827  affineMap->accumPostRotation(axis, radians);
828  return simplify(affineMap);
829  }
830 
831  MapBase::Ptr postTranslate(const Vec3d& tr) const;
832 
833  MapBase::Ptr postScale(const Vec3d& v) const;
834 
835  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
836  {
837  AffineMap::Ptr affineMap = getAffineMap();
838  affineMap->accumPostShear(axis0, axis1, shear);
839  return simplify(affineMap);
840  }
842 
843 private:
844  Vec3d mScaleValues, mVoxelSize, mScaleValuesInverse, mInvScaleSqr, mInvTwiceScale;
845 }; // class ScaleMap
846 
847 
851 {
852 public:
853  typedef boost::shared_ptr<UniformScaleMap> Ptr;
854  typedef boost::shared_ptr<const UniformScaleMap> ConstPtr;
855 
857  UniformScaleMap(double scale): ScaleMap(Vec3d(scale, scale, scale)) {}
858  UniformScaleMap(const UniformScaleMap& other): ScaleMap(other) {}
860 
862  static MapBase::Ptr create() { return MapBase::Ptr(new UniformScaleMap()); }
864  MapBase::Ptr copy() const { return MapBase::Ptr(new UniformScaleMap(*this)); }
865 
866  static bool isRegistered() { return MapRegistry::isRegistered(UniformScaleMap::mapType()); }
867  static void registerMap()
868  {
869  MapRegistry::registerMap(
870  UniformScaleMap::mapType(),
871  UniformScaleMap::create);
872  }
873 
874  Name type() const { return mapType(); }
875  static Name mapType() { return Name("UniformScaleMap"); }
876 
877  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
878 
879  bool operator==(const UniformScaleMap& other) const { return ScaleMap::operator==(other); }
880  bool operator!=(const UniformScaleMap& other) const { return !(*this == other); }
881 
884  MapBase::Ptr preTranslate(const Vec3d& tr) const;
885 
888  MapBase::Ptr postTranslate(const Vec3d& tr) const;
889 
890 }; // class UniformScaleMap
891 
892 
894 
895 
896 inline MapBase::Ptr
897 ScaleMap::preScale(const Vec3d& v) const
898 {
899  const Vec3d new_scale(v * mScaleValues);
900  if (isApproxEqual(new_scale[0],new_scale[1]) && isApproxEqual(new_scale[0],new_scale[2])) {
901  return MapBase::Ptr(new UniformScaleMap(new_scale[0]));
902  } else {
903  return MapBase::Ptr(new ScaleMap(new_scale));
904  }
905 }
906 
907 
908 inline MapBase::Ptr
909 ScaleMap::postScale(const Vec3d& v) const
910 { // pre-post Scale are the same for a scale map
911  return preScale(v);
912 }
913 
914 
917 {
918 public:
919  typedef boost::shared_ptr<TranslationMap> Ptr;
920  typedef boost::shared_ptr<const TranslationMap> ConstPtr;
921 
922  // default constructor is a translation by zero.
923  TranslationMap(): MapBase(), mTranslation(Vec3d(0,0,0)) {}
924  TranslationMap(const Vec3d& t): MapBase(), mTranslation(t) {}
925  TranslationMap(const TranslationMap& other): MapBase(), mTranslation(other.mTranslation) {}
926 
928 
930  static MapBase::Ptr create() { return MapBase::Ptr(new TranslationMap()); }
932  MapBase::Ptr copy() const { return MapBase::Ptr(new TranslationMap(*this)); }
933 
934  static bool isRegistered() { return MapRegistry::isRegistered(TranslationMap::mapType()); }
935 
936  static void registerMap()
937  {
938  MapRegistry::registerMap(
939  TranslationMap::mapType(),
940  TranslationMap::create);
941  }
942 
943  Name type() const { return mapType(); }
944  static Name mapType() { return Name("TranslationMap"); }
945 
947  bool isLinear() const { return true; }
948 
950  bool hasUniformScale() const { return true; }
951 
953  Vec3d applyMap(const Vec3d& in) const { return in + mTranslation; }
955  Vec3d applyInverseMap(const Vec3d& in) const { return in - mTranslation; }
957  Vec3d applyJacobian(const Vec3d& in, const Vec3d&) const { return applyJacobian(in); }
959  Vec3d applyJacobian(const Vec3d& in) const { return in; }
960 
963  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
966  Vec3d applyIJT(const Vec3d& in) const {return in;}
968  Mat3d applyIJC(const Mat3d& mat) const {return mat;}
969  Mat3d applyIJC(const Mat3d& mat, const Vec3d&, const Vec3d&) const { return applyIJC(mat); }
970 
972  double determinant(const Vec3d& ) const { return determinant(); }
974  double determinant() const { return 1.0; }
975 
977  Vec3d voxelSize() const { return Vec3d(1,1,1);}
979  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
980 
982  const Vec3d& getTranslation() const { return mTranslation; }
984  void read(std::istream& is) { mTranslation.read(is); }
986  void write(std::ostream& os) const { mTranslation.write(os); }
987 
989  std::string str() const
990  {
991  std::ostringstream buffer;
992  buffer << " - translation: " << mTranslation << std::endl;
993  return buffer.str();
994  }
995 
996  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
997 
998  bool operator==(const TranslationMap& other) const
999  {
1000  // ::eq() uses a tolerance
1001  return mTranslation.eq(other.mTranslation);
1002  }
1003 
1004  bool operator!=(const TranslationMap& other) const { return !(*this == other); }
1005 
1008  {
1009  Mat4d matrix(Mat4d::identity());
1010  matrix.setTranslation(mTranslation);
1011 
1012  AffineMap::Ptr affineMap(new AffineMap(matrix));
1013  return affineMap;
1014  }
1015 
1017  MapBase::Ptr preRotate(double radians, Axis axis) const
1020  {
1021  AffineMap::Ptr affineMap = getAffineMap();
1022  affineMap->accumPreRotation(axis, radians);
1023  return simplify(affineMap);
1024 
1025  }
1027  {
1028  return MapBase::Ptr(new TranslationMap(t + mTranslation));
1029  }
1030 
1031  MapBase::Ptr preScale(const Vec3d& v) const;
1032 
1033  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
1034  {
1035  AffineMap::Ptr affineMap = getAffineMap();
1036  affineMap->accumPreShear(axis0, axis1, shear);
1037  return simplify(affineMap);
1038  }
1040 
1042  MapBase::Ptr postRotate(double radians, Axis axis) const
1045  {
1046  AffineMap::Ptr affineMap = getAffineMap();
1047  affineMap->accumPostRotation(axis, radians);
1048  return simplify(affineMap);
1049 
1050  }
1052  { // post and pre are the same for this
1053  return MapBase::Ptr(new TranslationMap(t + mTranslation));
1054  }
1055 
1056  MapBase::Ptr postScale(const Vec3d& v) const;
1057 
1058  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
1059  {
1060  AffineMap::Ptr affineMap = getAffineMap();
1061  affineMap->accumPostShear(axis0, axis1, shear);
1062  return simplify(affineMap);
1063  }
1065 
1066 private:
1067  Vec3d mTranslation;
1068 }; // class TranslationMap
1069 
1070 
1072 
1073 
1078 {
1079 public:
1080  typedef boost::shared_ptr<ScaleTranslateMap> Ptr;
1081  typedef boost::shared_ptr<const ScaleTranslateMap> ConstPtr;
1082 
1084  MapBase(),
1085  mTranslation(Vec3d(0,0,0)),
1086  mScaleValues(Vec3d(1,1,1)),
1087  mVoxelSize(Vec3d(1,1,1)),
1088  mScaleValuesInverse(Vec3d(1,1,1)),
1089  mInvScaleSqr(1,1,1),
1090  mInvTwiceScale(0.5,0.5,0.5)
1091  {
1092  }
1093 
1094  ScaleTranslateMap(const Vec3d& scale, const Vec3d& translate):
1095  MapBase(),
1096  mTranslation(translate),
1097  mScaleValues(scale),
1098  mVoxelSize(std::abs(scale(0)), std::abs(scale(1)), std::abs(scale(2)))
1099  {
1100  const double determinant = scale[0]* scale[1] * scale[2];
1101  if (std::abs(determinant) < 3.0 * tolerance<double>::value()) {
1102  OPENVDB_THROW(ArithmeticError, "Non-zero scale values required");
1103  }
1104  mScaleValuesInverse = 1.0 / mScaleValues;
1105  mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1106  mInvTwiceScale = mScaleValuesInverse / 2;
1107  }
1108 
1109  ScaleTranslateMap(const ScaleMap& scale, const TranslationMap& translate):
1110  MapBase(),
1111  mTranslation(translate.getTranslation()),
1112  mScaleValues(scale.getScale()),
1113  mVoxelSize(std::abs(mScaleValues(0)),
1114  std::abs(mScaleValues(1)),
1115  std::abs(mScaleValues(2))),
1116  mScaleValuesInverse(1.0 / scale.getScale())
1117  {
1118  mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1119  mInvTwiceScale = mScaleValuesInverse / 2;
1120  }
1121 
1123  MapBase(),
1124  mTranslation(other.mTranslation),
1125  mScaleValues(other.mScaleValues),
1126  mVoxelSize(other.mVoxelSize),
1127  mScaleValuesInverse(other.mScaleValuesInverse),
1128  mInvScaleSqr(other.mInvScaleSqr),
1129  mInvTwiceScale(other.mInvTwiceScale)
1130  {}
1131 
1133 
1137  MapBase::Ptr copy() const { return MapBase::Ptr(new ScaleTranslateMap(*this)); }
1138 
1139  static bool isRegistered() { return MapRegistry::isRegistered(ScaleTranslateMap::mapType()); }
1140 
1141  static void registerMap()
1142  {
1143  MapRegistry::registerMap(
1144  ScaleTranslateMap::mapType(),
1145  ScaleTranslateMap::create);
1146  }
1147 
1148  Name type() const { return mapType(); }
1149  static Name mapType() { return Name("ScaleTranslateMap"); }
1150 
1152  bool isLinear() const { return true; }
1153 
1156  bool hasUniformScale() const
1157  {
1158  bool value = isApproxEqual(
1159  std::abs(mScaleValues.x()), std::abs(mScaleValues.y()), double(5e-7));
1160  value = value && isApproxEqual(
1161  std::abs(mScaleValues.x()), std::abs(mScaleValues.z()), double(5e-7));
1162  return value;
1163  }
1164 
1166  Vec3d applyMap(const Vec3d& in) const
1167  {
1168  return Vec3d(
1169  in.x() * mScaleValues.x() + mTranslation.x(),
1170  in.y() * mScaleValues.y() + mTranslation.y(),
1171  in.z() * mScaleValues.z() + mTranslation.z());
1172  }
1174  Vec3d applyInverseMap(const Vec3d& in) const
1175  {
1176  return Vec3d(
1177  (in.x() - mTranslation.x() ) / mScaleValues.x(),
1178  (in.y() - mTranslation.y() ) / mScaleValues.y(),
1179  (in.z() - mTranslation.z() ) / mScaleValues.z());
1180  }
1181 
1183  Vec3d applyJacobian(const Vec3d& in, const Vec3d&) const { return applyJacobian(in); }
1185  Vec3d applyJacobian(const Vec3d& in) const { return in * mScaleValues; }
1186 
1189  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
1191  Vec3d applyIJT(const Vec3d& in) const
1192  {
1193  return Vec3d(
1194  in.x() / mScaleValues.x(),
1195  in.y() / mScaleValues.y(),
1196  in.z() / mScaleValues.z());
1197  }
1199  Mat3d applyIJC(const Mat3d& in) const
1200  {
1201  Mat3d tmp;
1202  for (int i=0; i<3; i++){
1203  tmp.setRow(i, in.row(i)*mScaleValuesInverse(i));
1204  }
1205  for (int i=0; i<3; i++){
1206  tmp.setCol(i, tmp.col(i)*mScaleValuesInverse(i));
1207  }
1208  return tmp;
1209  }
1210  Mat3d applyIJC(const Mat3d& in, const Vec3d&, const Vec3d& ) const { return applyIJC(in); }
1211 
1213  double determinant(const Vec3d& ) const { return determinant(); }
1215  double determinant() const { return mScaleValues.x()*mScaleValues.y()*mScaleValues.z(); }
1217  Vec3d voxelSize() const { return mVoxelSize;}
1220  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
1221 
1223  const Vec3d& getScale() const { return mScaleValues; }
1225  const Vec3d& getTranslation() const { return mTranslation; }
1226 
1228  const Vec3d& getInvScaleSqr() const {return mInvScaleSqr;}
1230  const Vec3d& getInvTwiceScale() const {return mInvTwiceScale;}
1232  const Vec3d& getInvScale() const {return mScaleValuesInverse; }
1233 
1235  void read(std::istream& is)
1236  {
1237  mTranslation.read(is);
1238  mScaleValues.read(is);
1239  mVoxelSize.read(is);
1240  mScaleValuesInverse.read(is);
1241  mInvScaleSqr.read(is);
1242  mInvTwiceScale.read(is);
1243  }
1245  void write(std::ostream& os) const
1246  {
1247  mTranslation.write(os);
1248  mScaleValues.write(os);
1249  mVoxelSize.write(os);
1250  mScaleValuesInverse.write(os);
1251  mInvScaleSqr.write(os);
1252  mInvTwiceScale.write(os);
1253  }
1255  std::string str() const
1256  {
1257  std::ostringstream buffer;
1258  buffer << " - translation: " << mTranslation << std::endl;
1259  buffer << " - scale: " << mScaleValues << std::endl;
1260  buffer << " - voxel dimensions: " << mVoxelSize << std::endl;
1261  return buffer.str();
1262  }
1263 
1264  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1265 
1266  bool operator==(const ScaleTranslateMap& other) const
1267  {
1268  // ::eq() uses a tolerance
1269  if (!mScaleValues.eq(other.mScaleValues)) { return false; }
1270  if (!mTranslation.eq(other.mTranslation)) { return false; }
1271  return true;
1272  }
1273 
1274  bool operator!=(const ScaleTranslateMap& other) const { return !(*this == other); }
1275 
1278  {
1279  AffineMap::Ptr affineMap(new AffineMap(math::scale<Mat4d>(mScaleValues)));
1280  affineMap->accumPostTranslation(mTranslation);
1281  return affineMap;
1282  }
1283 
1285  MapBase::Ptr preRotate(double radians, Axis axis) const
1288  {
1289  AffineMap::Ptr affineMap = getAffineMap();
1290  affineMap->accumPreRotation(axis, radians);
1291  return simplify(affineMap);
1292  }
1294  {
1295  const Vec3d& s = mScaleValues;
1296  const Vec3d scaled_trans( t.x() * s.x(),
1297  t.y() * s.y(),
1298  t.z() * s.z() );
1299  return MapBase::Ptr( new ScaleTranslateMap(mScaleValues, mTranslation + scaled_trans));
1300  }
1301 
1302  MapBase::Ptr preScale(const Vec3d& v) const;
1303 
1304  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
1305  {
1306  AffineMap::Ptr affineMap = getAffineMap();
1307  affineMap->accumPreShear(axis0, axis1, shear);
1308  return simplify(affineMap);
1309  }
1311 
1313  MapBase::Ptr postRotate(double radians, Axis axis) const
1316  {
1317  AffineMap::Ptr affineMap = getAffineMap();
1318  affineMap->accumPostRotation(axis, radians);
1319  return simplify(affineMap);
1320  }
1322  {
1323  return MapBase::Ptr( new ScaleTranslateMap(mScaleValues, mTranslation + t));
1324  }
1325 
1326  MapBase::Ptr postScale(const Vec3d& v) const;
1327 
1328  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
1329  {
1330  AffineMap::Ptr affineMap = getAffineMap();
1331  affineMap->accumPostShear(axis0, axis1, shear);
1332  return simplify(affineMap);
1333  }
1335 
1336 private:
1337  Vec3d mTranslation, mScaleValues, mVoxelSize, mScaleValuesInverse,
1338  mInvScaleSqr, mInvTwiceScale;
1339 }; // class ScaleTanslateMap
1340 
1341 
1342 inline MapBase::Ptr
1343 ScaleMap::postTranslate(const Vec3d& t) const
1344 {
1345  return MapBase::Ptr(new ScaleTranslateMap(mScaleValues, t));
1346 }
1347 
1348 
1349 inline MapBase::Ptr
1350 ScaleMap::preTranslate(const Vec3d& t) const
1351 {
1352 
1353  const Vec3d& s = mScaleValues;
1354  const Vec3d scaled_trans( t.x() * s.x(),
1355  t.y() * s.y(),
1356  t.z() * s.z() );
1357  return MapBase::Ptr(new ScaleTranslateMap(mScaleValues, scaled_trans));
1358 }
1359 
1360 
1364 {
1365 public:
1366  typedef boost::shared_ptr<UniformScaleTranslateMap> Ptr;
1367  typedef boost::shared_ptr<const UniformScaleTranslateMap> ConstPtr;
1368 
1370  UniformScaleTranslateMap(double scale, const Vec3d& translate):
1371  ScaleTranslateMap(Vec3d(scale,scale,scale), translate) {}
1373  ScaleTranslateMap(scale.getScale(), translate.getTranslation()) {}
1374 
1377 
1381  MapBase::Ptr copy() const { return MapBase::Ptr(new UniformScaleTranslateMap(*this)); }
1382 
1383  static bool isRegistered()
1384  {
1385  return MapRegistry::isRegistered(UniformScaleTranslateMap::mapType());
1386  }
1387 
1388  static void registerMap()
1389  {
1390  MapRegistry::registerMap(
1391  UniformScaleTranslateMap::mapType(),
1392  UniformScaleTranslateMap::create);
1393  }
1394 
1395  Name type() const { return mapType(); }
1396  static Name mapType() { return Name("UniformScaleTranslateMap"); }
1397 
1398  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1399 
1400  bool operator==(const UniformScaleTranslateMap& other) const
1401  {
1402  return ScaleTranslateMap::operator==(other);
1403  }
1404  bool operator!=(const UniformScaleTranslateMap& other) const { return !(*this == other); }
1405 
1409  {
1410  const double scale = this->getScale().x();
1411  const Vec3d new_trans = this->getTranslation() + scale * t;
1412  return MapBase::Ptr( new UniformScaleTranslateMap(scale, new_trans));
1413  }
1414 
1418  {
1419  const double scale = this->getScale().x();
1420  return MapBase::Ptr( new UniformScaleTranslateMap(scale, this->getTranslation() + t));
1421  }
1422 }; // class UniformScaleTanslateMap
1423 
1424 
1425 inline MapBase::Ptr
1426 UniformScaleMap::postTranslate(const Vec3d& t) const
1427 {
1428  const double scale = this->getScale().x();
1429  return MapBase::Ptr(new UniformScaleTranslateMap(scale, t));
1430 }
1431 
1432 
1433 inline MapBase::Ptr
1434 UniformScaleMap::preTranslate(const Vec3d& t) const
1435 {
1436  const double scale = this->getScale().x();
1437  return MapBase::Ptr(new UniformScaleTranslateMap(scale, scale*t));
1438 }
1439 
1440 
1441 inline MapBase::Ptr
1442 TranslationMap::preScale(const Vec3d& v) const
1443 {
1444  if (isApproxEqual(v[0],v[1]) && isApproxEqual(v[0],v[2])) {
1445  return MapBase::Ptr(new UniformScaleTranslateMap(v[0], mTranslation));
1446  } else {
1447  return MapBase::Ptr(new ScaleTranslateMap(v, mTranslation));
1448  }
1449 }
1450 
1451 
1452 inline MapBase::Ptr
1453 TranslationMap::postScale(const Vec3d& v) const
1454 {
1455  if (isApproxEqual(v[0],v[1]) && isApproxEqual(v[0],v[2])) {
1456  return MapBase::Ptr(new UniformScaleTranslateMap(v[0], v[0]*mTranslation));
1457  } else {
1458  const Vec3d trans(mTranslation.x()*v.x(),
1459  mTranslation.y()*v.y(),
1460  mTranslation.z()*v.z());
1461  return MapBase::Ptr(new ScaleTranslateMap(v, trans));
1462  }
1463 }
1464 
1465 
1466 inline MapBase::Ptr
1467 ScaleTranslateMap::preScale(const Vec3d& v) const
1468 {
1469  const Vec3d new_scale( v * mScaleValues );
1470  if (isApproxEqual(new_scale[0],new_scale[1]) && isApproxEqual(new_scale[0],new_scale[2])) {
1471  return MapBase::Ptr( new UniformScaleTranslateMap(new_scale[0], mTranslation));
1472  } else {
1473  return MapBase::Ptr( new ScaleTranslateMap(new_scale, mTranslation));
1474  }
1475 }
1476 
1477 
1478 inline MapBase::Ptr
1479 ScaleTranslateMap::postScale(const Vec3d& v) const
1480 {
1481  const Vec3d new_scale( v * mScaleValues );
1482  const Vec3d new_trans( mTranslation.x()*v.x(),
1483  mTranslation.y()*v.y(),
1484  mTranslation.z()*v.z() );
1485 
1486  if (isApproxEqual(new_scale[0],new_scale[1]) && isApproxEqual(new_scale[0],new_scale[2])) {
1487  return MapBase::Ptr( new UniformScaleTranslateMap(new_scale[0], new_trans));
1488  } else {
1489  return MapBase::Ptr( new ScaleTranslateMap(new_scale, new_trans));
1490  }
1491 }
1492 
1493 
1495 
1496 
1500 {
1501 public:
1502  typedef boost::shared_ptr<UnitaryMap> Ptr;
1503  typedef boost::shared_ptr<const UnitaryMap> ConstPtr;
1504 
1506  UnitaryMap(): mAffineMap(Mat4d::identity())
1507  {
1508  }
1509 
1510  UnitaryMap(const Vec3d& axis, double radians)
1511  {
1512  Mat3d matrix;
1513  matrix.setToRotation(axis, radians);
1514  mAffineMap = AffineMap(matrix);
1515  }
1516 
1517  UnitaryMap(Axis axis, double radians)
1518  {
1519  Mat4d matrix;
1520  matrix.setToRotation(axis, radians);
1521  mAffineMap = AffineMap(matrix);
1522  }
1523 
1524  UnitaryMap(const Mat3d& m)
1525  {
1526  // test that the mat3 is a rotation || reflection
1527  if (!isUnitary(m)) {
1528  OPENVDB_THROW(ArithmeticError, "Matrix initializing unitary map was not unitary");
1529  }
1530 
1531  Mat4d matrix(Mat4d::identity());
1532  matrix.setMat3(m);
1533  mAffineMap = AffineMap(matrix);
1534  }
1535 
1536  UnitaryMap(const Mat4d& m)
1537  {
1538  if (!isInvertible(m)) {
1540  "4x4 Matrix initializing unitary map was not unitary: not invertible");
1541  }
1542 
1543  if (!isAffine(m)) {
1545  "4x4 Matrix initializing unitary map was not unitary: not affine");
1546  }
1547 
1548  if (hasTranslation(m)) {
1550  "4x4 Matrix initializing unitary map was not unitary: had translation");
1551  }
1552 
1553  if (!isUnitary(m.getMat3())) {
1555  "4x4 Matrix initializing unitary map was not unitary");
1556  }
1557 
1558  mAffineMap = AffineMap(m);
1559  }
1560 
1561  UnitaryMap(const UnitaryMap& other):
1562  MapBase(other),
1563  mAffineMap(other.mAffineMap)
1564  {
1565  }
1566 
1567  UnitaryMap(const UnitaryMap& first, const UnitaryMap& second):
1568  mAffineMap(*(first.getAffineMap()), *(second.getAffineMap()))
1569  {
1570  }
1571 
1574  static MapBase::Ptr create() { return MapBase::Ptr(new UnitaryMap()); }
1576  MapBase::Ptr copy() const { return MapBase::Ptr(new UnitaryMap(*this)); }
1577 
1578  static bool isRegistered() { return MapRegistry::isRegistered(UnitaryMap::mapType()); }
1579 
1580  static void registerMap()
1581  {
1582  MapRegistry::registerMap(
1583  UnitaryMap::mapType(),
1584  UnitaryMap::create);
1585  }
1586 
1588  Name type() const { return mapType(); }
1590  static Name mapType() { return Name("UnitaryMap"); }
1591 
1593  bool isLinear() const { return true; }
1594 
1596  bool hasUniformScale() const { return true; }
1597 
1598  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1599 
1600  bool operator==(const UnitaryMap& other) const
1601  {
1602  // compare underlying linear map.
1603  if (mAffineMap!=other.mAffineMap) return false;
1604  return true;
1605  }
1606 
1607  bool operator!=(const UnitaryMap& other) const { return !(*this == other); }
1609  Vec3d applyMap(const Vec3d& in) const { return mAffineMap.applyMap(in); }
1611  Vec3d applyInverseMap(const Vec3d& in) const { return mAffineMap.applyInverseMap(in); }
1612 
1613  Vec3d applyJacobian(const Vec3d& in, const Vec3d&) const { return applyJacobian(in); }
1615  Vec3d applyJacobian(const Vec3d& in) const { return applyMap(in); }
1616 
1619  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
1621  Vec3d applyIJT(const Vec3d& in) const { return mAffineMap.applyIJT(in); }
1623  Mat3d applyIJC(const Mat3d& in) const { return mAffineMap.applyIJC(in); }
1624  Mat3d applyIJC(const Mat3d& in, const Vec3d&, const Vec3d& ) const { return applyIJC(in); }
1626  double determinant(const Vec3d& ) const { return determinant(); }
1628  double determinant() const { return mAffineMap.determinant(); }
1629 
1630 
1635  Vec3d voxelSize() const { return mAffineMap.voxelSize();}
1636  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
1637 
1639  void read(std::istream& is)
1640  {
1641  mAffineMap.read(is);
1642  }
1643 
1645  void write(std::ostream& os) const
1646  {
1647  mAffineMap.write(os);
1648  }
1650  std::string str() const
1651  {
1652  std::ostringstream buffer;
1653  buffer << mAffineMap.str();
1654  return buffer.str();
1655  }
1657  AffineMap::Ptr getAffineMap() const { return AffineMap::Ptr(new AffineMap(mAffineMap)); }
1658 
1660  MapBase::Ptr preRotate(double radians, Axis axis) const
1663  {
1664  UnitaryMap first(axis, radians);
1665  UnitaryMap::Ptr unitaryMap(new UnitaryMap(first, *this));
1666  return boost::static_pointer_cast<MapBase, UnitaryMap>(unitaryMap);
1667  }
1669  {
1670  AffineMap::Ptr affineMap = getAffineMap();
1671  affineMap->accumPreTranslation(t);
1672  return simplify(affineMap);
1673  }
1674  MapBase::Ptr preScale(const Vec3d& v) const
1675  {
1676  AffineMap::Ptr affineMap = getAffineMap();
1677  affineMap->accumPreScale(v);
1678  return simplify(affineMap);
1679  }
1680  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
1681  {
1682  AffineMap::Ptr affineMap = getAffineMap();
1683  affineMap->accumPreShear(axis0, axis1, shear);
1684  return simplify(affineMap);
1685  }
1687 
1688 
1690  MapBase::Ptr postRotate(double radians, Axis axis) const
1693  {
1694  UnitaryMap second(axis, radians);
1695  UnitaryMap::Ptr unitaryMap(new UnitaryMap(*this, second));
1696  return boost::static_pointer_cast<MapBase, UnitaryMap>(unitaryMap);
1697  }
1699  {
1700  AffineMap::Ptr affineMap = getAffineMap();
1701  affineMap->accumPostTranslation(t);
1702  return simplify(affineMap);
1703  }
1704  MapBase::Ptr postScale(const Vec3d& v) const
1705  {
1706  AffineMap::Ptr affineMap = getAffineMap();
1707  affineMap->accumPostScale(v);
1708  return simplify(affineMap);
1709  }
1710  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
1711  {
1712  AffineMap::Ptr affineMap = getAffineMap();
1713  affineMap->accumPostShear(axis0, axis1, shear);
1714  return simplify(affineMap);
1715  }
1717 
1718 private:
1719  AffineMap mAffineMap;
1720 }; // class UnitaryMap
1721 
1722 
1724 
1725 
1733 {
1734 public:
1735  typedef boost::shared_ptr<NonlinearFrustumMap> Ptr;
1736  typedef boost::shared_ptr<const NonlinearFrustumMap> ConstPtr;
1737 
1739  MapBase(),
1740  mBBox(Vec3d(0), Vec3d(1)),
1741  mTaper(1),
1742  mDepth(1)
1743  {
1744  init();
1745  }
1746 
1750  NonlinearFrustumMap(const BBoxd& bb, double taper, double depth):
1751  MapBase(),mBBox(bb), mTaper(taper), mDepth(depth)
1752  {
1753  init();
1754  }
1755 
1761  NonlinearFrustumMap(const BBoxd& bb, double taper, double depth,
1762  const MapBase::Ptr& secondMap):
1763  mBBox(bb), mTaper(taper), mDepth(depth)
1764  {
1765  if (!secondMap->isLinear() ) {
1767  "The second map in the Frustum transfrom must be linear");
1768  }
1769  mSecondMap = *( secondMap->getAffineMap() );
1770  init();
1771  }
1772 
1774  MapBase(),
1775  mBBox(other.mBBox),
1776  mTaper(other.mTaper),
1777  mDepth(other.mDepth),
1778  mSecondMap(other.mSecondMap),
1779  mHasSimpleAffine(other.mHasSimpleAffine)
1780  {
1781  init();
1782  }
1783 
1799  NonlinearFrustumMap(const Vec3d& position,
1800  const Vec3d& direction,
1801  const Vec3d& up,
1802  double aspect /* width / height */,
1803  double z_near, double depth,
1804  Coord::ValueType x_count, Coord::ValueType z_count) {
1805 
1809  if (!(depth > 0)) {
1811  "The frustum depth must be non-zero and positive");
1812  }
1813  if (!(up.length() > 0)) {
1815  "The frustum height must be non-zero and positive");
1816  }
1817  if (!(aspect > 0)) {
1819  "The frustum aspect ratio must be non-zero and positive");
1820  }
1821  if (!(isApproxEqual(up.dot(direction), 0.))) {
1823  "The frustum up orientation must be perpendicular to into-frustum direction");
1824  }
1825 
1826  double near_plane_height = 2 * up.length();
1827  double near_plane_width = aspect * near_plane_height;
1828 
1829  Coord::ValueType y_count = static_cast<int>(Round(x_count / aspect));
1830 
1831  mBBox = BBoxd(Vec3d(0,0,0), Vec3d(x_count, y_count, z_count));
1832  mDepth = depth / near_plane_width; // depth non-dimensionalized on width
1833  double gamma = near_plane_width / z_near;
1834  mTaper = 1./(mDepth*gamma + 1.);
1835 
1836  Vec3d direction_unit = direction;
1837  direction_unit.normalize();
1838 
1839  Mat4d r1(Mat4d::identity());
1840  r1.setToRotation(/*from*/Vec3d(0,0,1), /*to */direction_unit);
1841  Mat4d r2(Mat4d::identity());
1842  Vec3d temp = r1.inverse().transform(up);
1843  r2.setToRotation(/*from*/Vec3d(0,1,0), /*to*/temp );
1844  Mat4d scale = math::scale<Mat4d>(
1845  Vec3d(near_plane_width, near_plane_width, near_plane_width));
1846 
1847  // move the near plane to origin, rotate to align with axis, and scale down
1848  // T_inv * R1_inv * R2_inv * scale_inv
1849  Mat4d mat = scale * r2 * r1;
1850  mat.setTranslation(position + z_near*direction_unit);
1851 
1852  mSecondMap = AffineMap(mat);
1853 
1854  init();
1855  }
1856 
1861  MapBase::Ptr copy() const { return MapBase::Ptr(new NonlinearFrustumMap(*this)); }
1862 
1863  static bool isRegistered() { return MapRegistry::isRegistered(NonlinearFrustumMap::mapType()); }
1864 
1865  static void registerMap()
1866  {
1867  MapRegistry::registerMap(
1868  NonlinearFrustumMap::mapType(),
1869  NonlinearFrustumMap::create);
1870  }
1872  Name type() const { return mapType(); }
1874  static Name mapType() { return Name("NonlinearFrustumMap"); }
1875 
1877  bool isLinear() const { return false; }
1878 
1880  bool hasUniformScale() const { return false; }
1881 
1883  bool isIdentity() const
1884  {
1885  // The frustum can only be consistent with a linear map if the taper value is 1
1886  if (!isApproxEqual(mTaper, double(1)) ) return false;
1887 
1888  // There are various ways an identity can decomposed between the two parts of the
1889  // map. Best to just check that the principle vectors are stationary.
1890  const Vec3d e1(1,0,0);
1891  if (!applyMap(e1).eq(e1)) return false;
1892 
1893  const Vec3d e2(0,1,0);
1894  if (!applyMap(e2).eq(e2)) return false;
1895 
1896  const Vec3d e3(0,0,1);
1897  if (!applyMap(e3).eq(e3)) return false;
1898 
1899  return true;
1900  }
1901 
1902  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1903 
1904  bool operator==(const NonlinearFrustumMap& other) const
1905  {
1906  if (mBBox!=other.mBBox) return false;
1907  if (!isApproxEqual(mTaper, other.mTaper)) return false;
1908  if (!isApproxEqual(mDepth, other.mDepth)) return false;
1909 
1910  // Two linear transforms are equivalent iff they have the same translation
1911  // and have the same affects on orthongal spanning basis check translation
1912  Vec3d e(0,0,0);
1913  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
1915  e(0) = 1;
1916  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
1917  e(0) = 0;
1918  e(1) = 1;
1919  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
1920  e(1) = 0;
1921  e(2) = 1;
1922  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
1923  return true;
1924  }
1925 
1926  bool operator!=(const NonlinearFrustumMap& other) const { return !(*this == other); }
1927 
1929  Vec3d applyMap(const Vec3d& in) const
1930  {
1931  return mSecondMap.applyMap(applyFrustumMap(in));
1932  }
1933 
1935  Vec3d applyInverseMap(const Vec3d& in) const
1936  {
1937  return applyFrustumInverseMap(mSecondMap.applyInverseMap(in));
1938  }
1940  Vec3d applyJacobian(const Vec3d& in) const { return mSecondMap.applyJacobian(in); }
1942  Vec3d applyJacobian(const Vec3d& in, const Vec3d& loc) const
1943  {
1944  // Move the center of the x-face of the bbox
1945  // to the origin in index space.
1946  Vec3d centered(loc);
1947  centered = centered - mBBox.min();
1948  centered.x() -= mXo;
1949  centered.y() -= mYo;
1950 
1951  // scale the z-direction on depth / K count
1952  const double zprime = centered.z()*mDepthOnLz;
1953 
1954  const double scale = (mGamma * zprime + 1.) / mLx;
1955  const double scale2 = mGamma * mDepthOnLz / mLx;
1956 
1957  const Vec3d tmp(scale * in.x() + scale2 * centered.x()* in.z(),
1958  scale * in.y() + scale2 * centered.y()* in.z(),
1959  mDepthOnLz * in.z());
1960 
1961  return mSecondMap.applyJacobian(tmp);
1962  }
1963 
1965  Vec3d applyIJT(const Vec3d& in) const { return mSecondMap.applyIJT(in); }
1966 
1967  // the Jacobian of the nonlinear part of the transform is a sparse matrix
1968  // Jacobian^(-T) =
1969  //
1970  // (Lx)( 1/s 0 0 )
1971  // ( 0 1/s 0 )
1972  // ( -(x-xo)g/(sLx) -(y-yo)g/(sLx) Lz/(Depth Lx) )
1975  Vec3d applyIJT(const Vec3d& d1_is, const Vec3d& ijk) const
1976  {
1977  const Vec3d loc = applyFrustumMap(ijk);
1978  const double s = mGamma * loc.z() + 1.;
1979 
1980  // verify that we aren't at the singularity
1981  if (isApproxEqual(s, 0.)) {
1982  OPENVDB_THROW(ArithmeticError, "Tried to evaluate the frustum transform"
1983  " at the singular focal point (e.g. camera)");
1984  }
1985 
1986  const double sinv = 1.0/s; // 1/(z*gamma + 1)
1987  const double pt0 = mLx * sinv; // Lx / (z*gamma +1)
1988  const double pt1 = mGamma * pt0; // gamma * Lx / ( z*gamma +1)
1989  const double pt2 = pt1 * sinv; // gamma * Lx / ( z*gamma +1)**2
1990 
1991  const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
1992 
1993  // compute \frac{\partial E_i}{\partial x_j}
1994  Mat3d gradE(Mat3d::zero());
1995  for (int j = 0; j < 3; ++j ) {
1996  gradE(0,j) = pt0 * jacinv(0,j) - pt2 * loc.x()*jacinv(2,j);
1997  gradE(1,j) = pt0 * jacinv(1,j) - pt2 * loc.y()*jacinv(2,j);
1998  gradE(2,j) = (1./mDepthOnLz) * jacinv(2,j);
1999  }
2000 
2001  Vec3d result;
2002  for (int i = 0; i < 3; ++i) {
2003  result(0) = d1_is(0) * gradE(0,i) + d1_is(1) * gradE(1,i) + d1_is(2) * gradE(2,i);
2004  }
2005 
2006  return result;
2007 
2008  }
2009 
2011  Mat3d applyIJC(const Mat3d& in) const { return mSecondMap.applyIJC(in); }
2016  Mat3d applyIJC(const Mat3d& d2_is, const Vec3d& d1_is, const Vec3d& ijk) const
2017  {
2018  const Vec3d loc = applyFrustumMap(ijk);
2019 
2020  const double s = mGamma * loc.z() + 1.;
2021 
2022  // verify that we aren't at the singularity
2023  if (isApproxEqual(s, 0.)) {
2024  OPENVDB_THROW(ArithmeticError, "Tried to evaluate the frustum transform"
2025  " at the singular focal point (e.g. camera)");
2026  }
2027 
2028  // precompute
2029  const double sinv = 1.0/s; // 1/(z*gamma + 1)
2030  const double pt0 = mLx * sinv; // Lx / (z*gamma +1)
2031  const double pt1 = mGamma * pt0; // gamma * Lx / ( z*gamma +1)
2032  const double pt2 = pt1 * sinv; // gamma * Lx / ( z*gamma +1)**2
2033  const double pt3 = pt2 * sinv; // gamma * Lx / ( z*gamma +1)**3
2034 
2035  const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
2036 
2037  // compute \frac{\partial^2 E_i}{\partial x_j \partial x_k}
2038 
2039  Mat3d matE0(Mat3d::zero());
2040  Mat3d matE1(Mat3d::zero()); // matE2 = 0
2041  for(int j = 0; j < 3; j++) {
2042  for (int k = 0; k < 3; k++) {
2043 
2044  const double pt4 = 2. * jacinv(2,j) * jacinv(2,k) * pt3;
2045 
2046  matE0(j,k) = -(jacinv(0,j) * jacinv(2,k) + jacinv(2,j) * jacinv(0,k)) * pt2 +
2047  pt4 * loc.x();
2048 
2049  matE1(j,k) = -(jacinv(1,j) * jacinv(2,k) + jacinv(2,j) * jacinv(1,k)) * pt2 +
2050  pt4 * loc.y();
2051  }
2052  }
2053 
2054  // compute \frac{\partial E_i}{\partial x_j}
2055  Mat3d gradE(Mat3d::zero());
2056  for (int j = 0; j < 3; ++j ) {
2057  gradE(0,j) = pt0 * jacinv(0,j) - pt2 * loc.x()*jacinv(2,j);
2058  gradE(1,j) = pt0 * jacinv(1,j) - pt2 * loc.y()*jacinv(2,j);
2059  gradE(2,j) = (1./mDepthOnLz) * jacinv(2,j);
2060  }
2061 
2062  Mat3d result(Mat3d::zero());
2063  // compute \fac{\partial E_j}{\partial x_m} \fac{\partial E_i}{\partial x_n}
2064  // \frac{\partial^2 input}{\partial E_i \partial E_j}
2065  for (int m = 0; m < 3; ++m ) {
2066  for ( int n = 0; n < 3; ++n) {
2067  for (int i = 0; i < 3; ++i ) {
2068  for (int j = 0; j < 3; ++j) {
2069  result(m, n) += gradE(j, m) * gradE(i, n) * d2_is(i, j);
2070  }
2071  }
2072  }
2073  }
2074 
2075  for (int m = 0; m < 3; ++m ) {
2076  for ( int n = 0; n < 3; ++n) {
2077  result(m, n) +=
2078  matE0(m, n) * d1_is(0) + matE1(m, n) * d1_is(1);// + matE2(m, n) * d1_is(2);
2079  }
2080  }
2081 
2082  return result;
2083  }
2084 
2086  double determinant() const {return mSecondMap.determinant();} // no implementation
2087 
2090  double determinant(const Vec3d& loc) const
2091  {
2092  double s = mGamma * loc.z() + 1.0;
2093  double frustum_determinant = s * s * mDepthOnLzLxLx;
2094  return mSecondMap.determinant() * frustum_determinant;
2095  }
2096 
2099  {
2100  const Vec3d loc( 0.5*(mBBox.min().x() + mBBox.max().x()),
2101  0.5*(mBBox.min().y() + mBBox.max().y()),
2102  mBBox.min().z());
2103 
2104  return voxelSize(loc);
2105 
2106  }
2107 
2112  Vec3d voxelSize(const Vec3d& loc) const
2113  {
2114  Vec3d out, pos = applyMap(loc);
2115  out(0) = (applyMap(loc + Vec3d(1,0,0)) - pos).length();
2116  out(1) = (applyMap(loc + Vec3d(0,1,0)) - pos).length();
2117  out(2) = (applyMap(loc + Vec3d(0,0,1)) - pos).length();
2118  return out;
2119  }
2120 
2121  AffineMap::Ptr getAffineMap() const { return mSecondMap.getAffineMap(); }
2122 
2124  void setTaper(double t) { mTaper = t; init();}
2126  double getTaper() const { return mTaper; }
2128  void setDepth(double d) { mDepth = d; init();}
2130  double getDepth() const { return mDepth; }
2131  // gamma a non-dimensional number: nearplane x-width / camera to near plane distance
2132  double getGamma() const { return mGamma; }
2133 
2135  const BBoxd& getBBox() const { return mBBox; }
2136 
2138  const AffineMap& secondMap() const { return mSecondMap; }
2141  bool isValid() const { return !mBBox.empty();}
2142 
2144  bool hasSimpleAffine() const { return mHasSimpleAffine; }
2145 
2147  void read(std::istream& is)
2148  {
2149  // for backward compatibility with earlier version
2151  CoordBBox bb;
2152  bb.read(is);
2153  mBBox = BBoxd(bb.min().asVec3d(), bb.max().asVec3d());
2154  } else {
2155  mBBox.read(is);
2156  }
2157 
2158  is.read(reinterpret_cast<char*>(&mTaper), sizeof(double));
2159  is.read(reinterpret_cast<char*>(&mDepth), sizeof(double));
2160 
2161  // Read the second maps type.
2162  Name type = readString(is);
2163 
2164  // Check if the map has been registered.
2165  if(!MapRegistry::isRegistered(type)) {
2166  OPENVDB_THROW(KeyError, "Map " << type << " is not registered");
2167  }
2168 
2169  // Create the second map of the type and then read it in.
2170  MapBase::Ptr proxy = math::MapRegistry::createMap(type);
2171  proxy->read(is);
2172  mSecondMap = *(proxy->getAffineMap());
2173  init();
2174  }
2175 
2177  void write(std::ostream& os) const
2178  {
2179  mBBox.write(os);
2180  os.write(reinterpret_cast<const char*>(&mTaper), sizeof(double));
2181  os.write(reinterpret_cast<const char*>(&mDepth), sizeof(double));
2182 
2183  writeString(os, mSecondMap.type());
2184  mSecondMap.write(os);
2185  }
2186 
2188  std::string str() const
2189  {
2190  std::ostringstream buffer;
2191  buffer << " - taper: " << mTaper << std::endl;
2192  buffer << " - depth: " << mDepth << std::endl;
2193  buffer << " SecondMap: "<< mSecondMap.type() << std::endl;
2194  buffer << mSecondMap.str() << std::endl;
2195  return buffer.str();
2196  }
2197 
2199  MapBase::Ptr preRotate(double radians, Axis axis = X_AXIS) const
2202  {
2203  return MapBase::Ptr(
2204  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.preRotate(radians, axis)));
2205  }
2207  {
2208  return MapBase::Ptr(
2209  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.preTranslate(t)));
2210  }
2211  MapBase::Ptr preScale(const Vec3d& s) const
2212  {
2213  return MapBase::Ptr(
2214  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.preScale(s)));
2215  }
2216  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
2217  {
2218  return MapBase::Ptr(new NonlinearFrustumMap(
2219  mBBox, mTaper, mDepth, mSecondMap.preShear(shear, axis0, axis1)));
2220  }
2222 
2224  MapBase::Ptr postRotate(double radians, Axis axis = X_AXIS) const
2227  {
2228  return MapBase::Ptr(
2229  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.postRotate(radians, axis)));
2230  }
2232  {
2233  return MapBase::Ptr(
2234  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.postTranslate(t)));
2235  }
2236  MapBase::Ptr postScale(const Vec3d& s) const
2237  {
2238  return MapBase::Ptr(
2239  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.postScale(s)));
2240  }
2241  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
2242  {
2243  return MapBase::Ptr(new NonlinearFrustumMap(
2244  mBBox, mTaper, mDepth, mSecondMap.postShear(shear, axis0, axis1)));
2245  }
2247 
2248 private:
2249  void init()
2250  {
2251  // set up as a frustum
2252  mLx = mBBox.extents().x();
2253  mLy = mBBox.extents().y();
2254  mLz = mBBox.extents().z();
2255 
2256  if (isApproxEqual(mLx,0.) || isApproxEqual(mLy,0.) || isApproxEqual(mLz,0.) ) {
2257  OPENVDB_THROW(ArithmeticError, "The index space bounding box"
2258  " must have at least two index points in each direction.");
2259  }
2260 
2261  mXo = 0.5* mLx;
2262  mYo = 0.5* mLy;
2263 
2264  // mDepth is non-dimensionalized on near
2265  mGamma = (1./mTaper - 1) / mDepth;
2266 
2267  mDepthOnLz = mDepth/mLz;
2268  mDepthOnLzLxLx = mDepthOnLz/(mLx * mLx);
2269 
2271  mHasSimpleAffine = true;
2272  Vec3d tmp = mSecondMap.voxelSize();
2273 
2275  if (!isApproxEqual(tmp(0), tmp(1))) { mHasSimpleAffine = false; return; }
2276  if (!isApproxEqual(tmp(0), tmp(2))) { mHasSimpleAffine = false; return; }
2277 
2278  Vec3d trans = mSecondMap.applyMap(Vec3d(0,0,0));
2280  Vec3d tmp1 = mSecondMap.applyMap(Vec3d(1,0,0)) - trans;
2281  Vec3d tmp2 = mSecondMap.applyMap(Vec3d(0,1,0)) - trans;
2282  Vec3d tmp3 = mSecondMap.applyMap(Vec3d(0,0,1)) - trans;
2283 
2285  if (!isApproxEqual(tmp1.dot(tmp2), 0., 1.e-7)) { mHasSimpleAffine = false; return; }
2286  if (!isApproxEqual(tmp2.dot(tmp3), 0., 1.e-7)) { mHasSimpleAffine = false; return; }
2287  if (!isApproxEqual(tmp3.dot(tmp1), 0., 1.e-7)) { mHasSimpleAffine = false; return; }
2288  }
2289 
2290  Vec3d applyFrustumMap(const Vec3d& in) const
2291  {
2292 
2293  // Move the center of the x-face of the bbox
2294  // to the origin in index space.
2295  Vec3d out(in);
2296  out = out - mBBox.min();
2297  out.x() -= mXo;
2298  out.y() -= mYo;
2299 
2300  // scale the z-direction on depth / K count
2301  out.z() *= mDepthOnLz;
2302 
2303  double scale = (mGamma * out.z() + 1.)/ mLx;
2304 
2305  // scale the x-y on the length I count and apply tapper
2306  out.x() *= scale ;
2307  out.y() *= scale ;
2308 
2309  return out;
2310  }
2311 
2312  Vec3d applyFrustumInverseMap(const Vec3d& in) const
2313  {
2314  // invert taper and resize: scale = 1/( (z+1)/2 (mt-1) + 1)
2315  Vec3d out(in);
2316  double invScale = mLx / (mGamma * out.z() + 1.);
2317  out.x() *= invScale;
2318  out.y() *= invScale;
2319 
2320  out.x() += mXo;
2321  out.y() += mYo;
2322 
2323  out.z() /= mDepthOnLz;
2324 
2325  // move back
2326  out = out + mBBox.min();
2327  return out;
2328  }
2329 
2330  // bounding box in index space used in Frustum transforms.
2331  BBoxd mBBox;
2332 
2333  // taper value used in constructing Frustums.
2334  double mTaper;
2335  double mDepth;
2336 
2337  // defines the second map
2338  AffineMap mSecondMap;
2339 
2340  // these are derived from the above.
2341  double mLx, mLy, mLz;
2342  double mXo, mYo, mGamma, mDepthOnLz, mDepthOnLzLxLx;
2343 
2344  // true: if the mSecondMap is linear and has no shear, and has no non-uniform scale
2345  bool mHasSimpleAffine;
2346 }; // class NonlinearFrustumMap
2347 
2348 
2350 
2351 
2355 template<typename FirstMapType, typename SecondMapType>
2356 class CompoundMap
2357 {
2358 public:
2360 
2361  typedef boost::shared_ptr<MyType> Ptr;
2362  typedef boost::shared_ptr<const MyType> ConstPtr;
2363 
2364 
2365  CompoundMap() { updateAffineMatrix(); }
2366 
2367  CompoundMap(const FirstMapType& f, const SecondMapType& s): mFirstMap(f), mSecondMap(s)
2368  {
2369  updateAffineMatrix();
2370  }
2371 
2372  CompoundMap(const MyType& other):
2373  mFirstMap(other.mFirstMap),
2374  mSecondMap(other.mSecondMap),
2375  mAffineMap(other.mAffineMap)
2376  {}
2377 
2378  Name type() const { return mapType(); }
2379  static Name mapType()
2380  {
2381  return (FirstMapType::mapType() + Name(":") + SecondMapType::mapType());
2382  }
2383 
2384  bool operator==(const MyType& other) const
2385  {
2386  if (mFirstMap != other.mFirstMap) return false;
2387  if (mSecondMap != other.mSecondMap) return false;
2388  if (mAffineMap != other.mAffineMap) return false;
2389  return true;
2390  }
2391 
2392  bool operator!=(const MyType& other) const { return !(*this == other); }
2393 
2394  MyType& operator=(const MyType& other)
2395  {
2396  mFirstMap = other.mFirstMap;
2397  mSecondMap = other.mSecondMap;
2398  mAffineMap = other.mAffineMap;
2399  return *this;
2400  }
2401 
2402  bool isIdentity() const
2403  {
2405  return mAffineMap.isIdentity();
2406  } else {
2407  return mFirstMap.isIdentity()&&mSecondMap.isIdentity();
2408  }
2409  }
2410 
2411  bool isDiagonal() const {
2413  return mAffineMap.isDiagonal();
2414  } else {
2415  return mFirstMap.isDiagonal()&&mSecondMap.isDiagonal();
2416  }
2417  }
2418 
2420  {
2422  AffineMap::Ptr affine(new AffineMap(mAffineMap));
2423  return affine;
2424  } else {
2426  "Constant affine matrix representation not possible for this nonlinear map");
2427  }
2428  }
2429 
2430  // direct decompotion
2431  const FirstMapType& firstMap() const { return mFirstMap; }
2432  const SecondMapType& secondMap() const {return mSecondMap; }
2433 
2434  void setFirstMap(const FirstMapType& first) { mFirstMap = first; updateAffineMatrix(); }
2435  void setSecondMap(const SecondMapType& second) { mSecondMap = second; updateAffineMatrix(); }
2436 
2437  void read(std::istream& is)
2438  {
2439  mAffineMap.read(is);
2440  mFirstMap.read(is);
2441  mSecondMap.read(is);
2442  }
2443  void write(std::ostream& os) const
2444  {
2445  mAffineMap.write(os);
2446  mFirstMap.write(os);
2447  mSecondMap.write(os);
2448  }
2449 
2450 private:
2451  void updateAffineMatrix()
2452  {
2454  // both maps need to be linear, these methods are only defined for linear maps
2455  AffineMap::Ptr first = mFirstMap.getAffineMap();
2456  AffineMap::Ptr second= mSecondMap.getAffineMap();
2457  mAffineMap = AffineMap(*first, *second);
2458  }
2459  }
2460 
2461  FirstMapType mFirstMap;
2462  SecondMapType mSecondMap;
2463  // used for acceleration
2464  AffineMap mAffineMap;
2465 }; // class CompoundMap
2466 
2467 } // namespace math
2468 } // namespace OPENVDB_VERSION_NAME
2469 } // namespace openvdb
2470 
2471 #endif // OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED
2472 
2473 // Copyright (c) 2012-2013 DreamWorks Animation LLC
2474 // All rights reserved. This software is distributed under the
2475 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )