OpenVDB  1.1.0
Maps.h
Go to the documentation of this file.
1 
2 //
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 "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 
198  virtual void read(std::istream&) = 0;
199  virtual void write(std::ostream&) const = 0;
200 
201  virtual std::string str() const = 0;
202 
203  virtual MapBase::Ptr copy() const = 0;
204 
206 
207  virtual MapBase::Ptr preRotate(double radians, Axis axis = X_AXIS) const = 0;
208  virtual MapBase::Ptr preTranslate(const Vec3d&) const = 0;
209  virtual MapBase::Ptr preScale(const Vec3d&) const = 0;
210  virtual MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const = 0;
211 
212  virtual MapBase::Ptr postRotate(double radians, Axis axis = X_AXIS) const = 0;
213  virtual MapBase::Ptr postTranslate(const Vec3d&) const = 0;
214  virtual MapBase::Ptr postScale(const Vec3d&) const = 0;
215  virtual MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const = 0;
217 
218 protected:
219  MapBase() {}
220 
221  template<typename MapT>
222  static bool isEqualBase(const MapT& self, const MapBase& other)
223  {
224  return other.isType<MapT>() && (self == *static_cast<const MapT*>(&other));
225  }
226 };
227 
228 
230 
231 
235 {
236 public:
237  typedef std::map<Name, MapBase::MapFactory> MapDictionary;
238  typedef tbb::mutex Mutex;
239  typedef Mutex::scoped_lock Lock;
240 
241  static MapRegistry* instance();
242 
244  static MapBase::Ptr createMap(const Name&);
245 
247  static bool isRegistered(const Name&);
248 
250  static void registerMap(const Name&, MapBase::MapFactory);
251 
253  static void unregisterMap(const Name&);
254 
256  static void clear();
257 
258 private:
259  MapRegistry() {}
260  static MapRegistry* mInstance;
261 
262  mutable Mutex mMutex;
263  MapDictionary mMap;
264 };
265 
266 
268 
269 
273 {
274 public:
275  typedef boost::shared_ptr<AffineMap> Ptr;
276  typedef boost::shared_ptr<const AffineMap> ConstPtr;
277 
279  mMatrix(Mat4d::identity()),
280  mMatrixInv(Mat4d::identity()),
281  mJacobianInv(Mat3d::identity()),
282  mDeterminant(1),
283  mVoxelSize(Vec3d(1,1,1)),
284  mIsDiagonal(true),
285  mIsIdentity(true)
286  // the default constructor for translation is zero
287  {
288  }
289 
290  AffineMap(const Mat3d& m)
291  {
292  Mat4d mat4(Mat4d::identity());
293  mat4.setMat3(m);
294  mMatrix = mat4;
295  updateAcceleration();
296  }
297 
298  AffineMap(const Mat4d& m): mMatrix(m)
299  {
300  if (!isAffine(m)) {
302  "Tried to initialize an affine transform from a non-affine 4x4 matrix");
303  }
304  updateAcceleration();
305  }
306 
307  AffineMap(const AffineMap& other):
308  MapBase(other),
309  mMatrix(other.mMatrix),
310  mMatrixInv(other.mMatrixInv),
311  mJacobianInv(other.mJacobianInv),
312  mDeterminant(other.mDeterminant),
313  mVoxelSize(other.mVoxelSize),
314  mIsDiagonal(other.mIsDiagonal),
315  mIsIdentity(other.mIsIdentity)
316  {
317  }
318 
320  AffineMap(const AffineMap& first, const AffineMap& second):
321  mMatrix(first.mMatrix * second.mMatrix)
322  {
323  updateAcceleration();
324  }
325 
327 
329  static MapBase::Ptr create() { return MapBase::Ptr(new AffineMap()); }
331  MapBase::Ptr copy() const { return MapBase::Ptr(new AffineMap(*this)); }
332 
333  static bool isRegistered() { return MapRegistry::isRegistered(AffineMap::mapType()); }
334 
335  static void registerMap()
336  {
337  MapRegistry::registerMap(
338  AffineMap::mapType(),
339  AffineMap::create);
340  }
341 
342  Name type() const { return mapType(); }
343  static Name mapType() { return Name("AffineMap"); }
344 
346  bool isLinear() const { return true; }
347 
349  bool hasUniformScale() const { return isUnitary(mMatrix.getMat3());}
350 
351  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
352 
353  bool operator==(const AffineMap& other) const
354  {
355  // the Mat.eq() is approximate
356  if (!mMatrix.eq(other.mMatrix)) { return false; }
357  if (!mMatrixInv.eq(other.mMatrixInv)) { return false; }
358  return true;
359  }
360 
361  bool operator!=(const AffineMap& other) const { return !(*this == other); }
362 
363  AffineMap& operator=(const AffineMap& other)
364  {
365  mMatrix = other.mMatrix;
366  mMatrixInv = other.mMatrixInv;
367 
368  mJacobianInv = other.mJacobianInv;
369  mDeterminant = other.mDeterminant;
370  mVoxelSize = other.mVoxelSize;
371  mIsDiagonal = other.mIsDiagonal;
372  mIsIdentity = other.mIsIdentity;
373  return *this;
374  }
376  Vec3d applyMap(const Vec3d& in) const { return in * mMatrix; }
378  Vec3d applyInverseMap(const Vec3d& in) const {return in * mMatrixInv; }
379 
381  Vec3d applyIJT(const Vec3d& in, const Vec3d&) const { return applyIJT(in); }
383  Vec3d applyIJT(const Vec3d& in) const { return in * mJacobianInv; }
385  Mat3d applyIJC(const Mat3d& m) const {
386  return mJacobianInv.transpose()* m * mJacobianInv;
387  }
388  Mat3d applyIJC(const Mat3d& in, const Vec3d& , const Vec3d& ) const {
389  return applyIJC(in);
390  }
392  double determinant(const Vec3d& ) const { return determinant(); }
394  double determinant() const { return mDeterminant; }
395 
397 
398 
399  Vec3d voxelSize() const { return mVoxelSize; }
400  Vec3d voxelSize(const Vec3d&) const { return voxelSize(); }
402 
404  bool isIdentity() const { return mIsIdentity; }
406  bool isDiagonal() const { return mIsDiagonal; }
408  bool isScale() const { return isDiagonal(); }
410  bool isScaleTranslate() const { return math::isDiagonal(mMatrix.getMat3()); }
411 
412 
413  // Methods that modify the existing affine map
414 
416 
417 
418  void accumPreRotation(Axis axis, double radians)
419  {
420  mMatrix.preRotate(axis, radians);
421  updateAcceleration();
422  }
423  void accumPreScale(const Vec3d& v)
424  {
425  mMatrix.preScale(v);
426  updateAcceleration();
427  }
428  void accumPreTranslation(const Vec3d& v)
429  {
430  mMatrix.preTranslate(v);
431  updateAcceleration();
432  }
433  void accumPreShear(Axis axis0, Axis axis1, double shear)
434  {
435  mMatrix.preShear(axis0, axis1, shear);
436  updateAcceleration();
437  }
439 
440 
442 
443 
444  void accumPostRotation(Axis axis, double radians)
445  {
446  mMatrix.postRotate(axis, radians);
447  updateAcceleration();
448  }
449  void accumPostScale(const Vec3d& v)
450  {
451  mMatrix.postScale(v);
452  updateAcceleration();
453  }
454  void accumPostTranslation(const Vec3d& v)
455  {
456  mMatrix.postTranslate(v);
457  updateAcceleration();
458  }
459  void accumPostShear(Axis axis0, Axis axis1, double shear)
460  {
461  mMatrix.postShear(axis0, axis1, shear);
462  updateAcceleration();
463  }
465 
466 
468  void read(std::istream& is)
469  {
470  mMatrix.read(is);
471  updateAcceleration();
472  }
473 
475  void write(std::ostream& os) const
476  {
477  mMatrix.write(os);
478  }
479 
481  std::string str() const
482  {
483  std::ostringstream buffer;
484  buffer << " - mat4:\n" << mMatrix.str() << std::endl;
485  buffer << " - voxel dimensions: " << mVoxelSize << std::endl;
486  return buffer.str();
487  }
488 
490  boost::shared_ptr<FullyDecomposedMap> createDecomposedMap()
491  {
492  return createFullyDecomposedMap(mMatrix);
493  }
494 
496  AffineMap::Ptr getAffineMap() const { return AffineMap::Ptr(new AffineMap(*this)); }
497 
499  AffineMap::Ptr inverse() const { return AffineMap::Ptr(new AffineMap(mMatrixInv)); }
500 
501 
503 
504 
505  MapBase::Ptr preRotate(double radians, Axis axis = X_AXIS) const
506  {
507  AffineMap::Ptr affineMap = getAffineMap();
508  affineMap->accumPreRotation(axis, radians);
509  return simplify(affineMap);
510  }
511  MapBase::Ptr preTranslate(const Vec3d& t) const
512  {
513  AffineMap::Ptr affineMap = getAffineMap();
514  affineMap->accumPreTranslation(t);
515  return boost::static_pointer_cast<MapBase, AffineMap>(affineMap);
516  }
517  MapBase::Ptr preScale(const Vec3d& s) const
518  {
519  AffineMap::Ptr affineMap = getAffineMap();
520  affineMap->accumPreScale(s);
521  return boost::static_pointer_cast<MapBase, AffineMap>(affineMap);
522  }
523  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
524  {
525  AffineMap::Ptr affineMap = getAffineMap();
526  affineMap->accumPreShear(axis0, axis1, shear);
527  return simplify(affineMap);
528  }
530 
531 
533 
534 
535  MapBase::Ptr postRotate(double radians, Axis axis = X_AXIS) const
536  {
537  AffineMap::Ptr affineMap = getAffineMap();
538  affineMap->accumPostRotation(axis, radians);
539  return simplify(affineMap);
540  }
541  MapBase::Ptr postTranslate(const Vec3d& t) const
542  {
543  AffineMap::Ptr affineMap = getAffineMap();
544  affineMap->accumPostTranslation(t);
545  return boost::static_pointer_cast<MapBase, AffineMap>(affineMap);
546  }
547  MapBase::Ptr postScale(const Vec3d& s) const
548  {
549  AffineMap::Ptr affineMap = getAffineMap();
550  affineMap->accumPostScale(s);
551  return boost::static_pointer_cast<MapBase, AffineMap>(affineMap);
552  }
553  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
554  {
555  AffineMap::Ptr affineMap = getAffineMap();
556  affineMap->accumPostShear(axis0, axis1, shear);
557  return simplify(affineMap);
558  }
560 
562  Mat4d getMat4() const { return mMatrix;}
563  const Mat4d& getConstMat4() const {return mMatrix;}
564  const Mat3d& getConstJacobianInv() const {return mJacobianInv;}
565 
566 private:
567  void updateAcceleration() {
568  mDeterminant = mMatrix.getMat3().det();
569 
570  if (std::abs(mDeterminant) < (3.0 * tolerance<double>::value())) {
572  "Tried to initialize an affine transform from a nearly signular matrix");
573  }
574  mMatrixInv = mMatrix.inverse();
575  mJacobianInv = mMatrixInv.getMat3().transpose();
576  mIsDiagonal = math::isDiagonal(mMatrix);
577  mIsIdentity = math::isIdentity(mMatrix);
578  Vec3d pos = applyMap(Vec3d(0,0,0));
579  mVoxelSize(0) = (applyMap(Vec3d(1,0,0)) - pos).length();
580  mVoxelSize(1) = (applyMap(Vec3d(0,1,0)) - pos).length();
581  mVoxelSize(2) = (applyMap(Vec3d(0,0,1)) - pos).length();
582  }
583 
584  // the underlying matrix
585  Mat4d mMatrix;
586 
587  // stored for acceleration
588  Mat4d mMatrixInv;
589  Mat3d mJacobianInv;
590  double mDeterminant;
591  Vec3d mVoxelSize;
592  bool mIsDiagonal, mIsIdentity;
593 }; // class AffineMap
594 
595 
597 
598 
602 {
603 public:
604  typedef boost::shared_ptr<ScaleMap> Ptr;
605  typedef boost::shared_ptr<const ScaleMap> ConstPtr;
606 
607  ScaleMap(): MapBase(), mScaleValues(Vec3d(1,1,1)), mVoxelSize(Vec3d(1,1,1)),
608  mInvScaleSqr(1,1,1), mInvTwiceScale(0.5,0.5,0.5){}
609 
611  MapBase(),
612  mScaleValues(scale),
613  mVoxelSize(Vec3d(std::abs(scale(0)),std::abs(scale(1)), std::abs(scale(2))))
614  {
615  if (scale.eq(Vec3d(0.0), 1.0e-10)) {
616  OPENVDB_THROW(ArithmeticError, "Non-zero scale values required");
617  }
618  mScaleValuesInverse = 1.0 / mScaleValues;
619  mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
620  mInvTwiceScale = mScaleValuesInverse / 2;
621  }
622 
623  ScaleMap(const ScaleMap& other):
624  MapBase(),
625  mScaleValues(other.mScaleValues),
626  mVoxelSize(other.mVoxelSize),
627  mScaleValuesInverse(other.mScaleValuesInverse),
628  mInvScaleSqr(other.mInvScaleSqr),
629  mInvTwiceScale(other.mInvTwiceScale)
630  {
631  }
632 
634 
636  static MapBase::Ptr create() { return MapBase::Ptr(new ScaleMap()); }
638  MapBase::Ptr copy() const { return MapBase::Ptr(new ScaleMap(*this)); }
639 
640  static bool isRegistered() { return MapRegistry::isRegistered(ScaleMap::mapType()); }
641 
642  static void registerMap()
643  {
644  MapRegistry::registerMap(
645  ScaleMap::mapType(),
646  ScaleMap::create);
647  }
648 
649  Name type() const { return mapType(); }
650  static Name mapType() { return Name("ScaleMap"); }
651 
653  bool isLinear() const { return true; }
654 
656  bool hasUniformScale() const {
657  bool value;
658  value = isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.y()));
659  value = value && isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.z()));
660  return value; }
661 
663  Vec3d applyMap(const Vec3d& in) const
664  {
665  return Vec3d(
666  in.x() * mScaleValues.x(),
667  in.y() * mScaleValues.y(),
668  in.z() * mScaleValues.z());
669  }
671  Vec3d applyInverseMap(const Vec3d& in) const
672  {
673  return Vec3d(
674  in.x() * mScaleValuesInverse.x(),
675  in.y() * mScaleValuesInverse.y(),
676  in.z() * mScaleValuesInverse.z());
677  }
680  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
682  Vec3d applyIJT(const Vec3d& in) const { return applyInverseMap(in); }
684  Mat3d applyIJC(const Mat3d& in) const {
685  Mat3d tmp;
686  for (int i=0; i<3; i++){
687  tmp.setRow(i, in.row(i)*mScaleValuesInverse(i));
688  }
689  for (int i=0; i<3; i++){
690  tmp.setCol(i, tmp.col(i)*mScaleValuesInverse(i));
691  }
692  return tmp;
693  }
694  Mat3d applyIJC(const Mat3d& in, const Vec3d&, const Vec3d& ) const {
695  return applyIJC(in);
696  }
698  double determinant(const Vec3d& ) const { return determinant(); }
700  double determinant() const { return mScaleValues.x()*mScaleValues.y()*mScaleValues.z(); }
701 
703  const Vec3d& getScale() const {return mScaleValues;}
704 
706  const Vec3d& getInvScaleSqr() const {return mInvScaleSqr;}
708  const Vec3d& getInvTwiceScale() const {return mInvTwiceScale;}
710  const Vec3d& getInvScale() const {return mScaleValuesInverse; }
711 
713 
714 
715 
716 
717  Vec3d voxelSize() const { return mVoxelSize; }
718  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
720 
722  void read(std::istream& is)
723  {
724  mScaleValues.read(is);
725  mVoxelSize.read(is);
726  mScaleValuesInverse.read(is);
727  mInvScaleSqr.read(is);
728  mInvTwiceScale.read(is);
729  }
731  void write(std::ostream& os) const
732  {
733  mScaleValues.write(os);
734  mVoxelSize.write(os);
735  mScaleValuesInverse.write(os);
736  mInvScaleSqr.write(os);
737  mInvTwiceScale.write(os);
738  }
740  std::string str() const
741  {
742  std::ostringstream buffer;
743  buffer << " - scale: " << mScaleValues << std::endl;
744  buffer << " - voxel dimensions: " << mVoxelSize << std::endl;
745  return buffer.str();
746  }
747 
748  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
749 
750  bool operator==(const ScaleMap& other) const
751  {
752  // ::eq() uses a tolerance
753  if (!mScaleValues.eq(other.mScaleValues)) { return false; }
754  return true;
755  }
756 
757  bool operator!=(const ScaleMap& other) const { return !(*this == other); }
758 
760  AffineMap::Ptr getAffineMap() const
761  {
762  return AffineMap::Ptr(new AffineMap(math::scale<Mat4d>(mScaleValues)));
763  }
764 
765 
766 
768 
769 
770  MapBase::Ptr preRotate(double radians, Axis axis) const
771  {
772  AffineMap::Ptr affineMap = getAffineMap();
773  affineMap->accumPreRotation(axis, radians);
774  return simplify(affineMap);
775  }
776 
777  MapBase::Ptr preTranslate(const Vec3d& tr) const;
778 
779  MapBase::Ptr preScale(const Vec3d& v) const;
780 
781  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
782  {
783  AffineMap::Ptr affineMap = getAffineMap();
784  affineMap->accumPreShear(axis0, axis1, shear);
785  return simplify(affineMap);
786  }
788 
789 
791 
792 
793  MapBase::Ptr postRotate(double radians, Axis axis) const
794  {
795  AffineMap::Ptr affineMap = getAffineMap();
796  affineMap->accumPostRotation(axis, radians);
797  return simplify(affineMap);
798  }
799 
800  MapBase::Ptr postTranslate(const Vec3d& tr) const;
801 
802  MapBase::Ptr postScale(const Vec3d& v) const;
803 
804  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
805  {
806  AffineMap::Ptr affineMap = getAffineMap();
807  affineMap->accumPostShear(axis0, axis1, shear);
808  return simplify(affineMap);
809  }
811 
812 private:
813  Vec3d mScaleValues, mVoxelSize, mScaleValuesInverse, mInvScaleSqr, mInvTwiceScale;
814 }; // class ScaleMap
815 
816 
820 {
821 public:
822  typedef boost::shared_ptr<UniformScaleMap> Ptr;
823  typedef boost::shared_ptr<const UniformScaleMap> ConstPtr;
824 
826  UniformScaleMap(double scale): ScaleMap(Vec3d(scale, scale, scale)) {}
827  UniformScaleMap(const UniformScaleMap& other): ScaleMap(other) {}
829 
831  static MapBase::Ptr create() { return MapBase::Ptr(new UniformScaleMap()); }
833  MapBase::Ptr copy() const { return MapBase::Ptr(new UniformScaleMap(*this)); }
834 
835  static bool isRegistered() { return MapRegistry::isRegistered(UniformScaleMap::mapType()); }
836  static void registerMap()
837  {
838  MapRegistry::registerMap(
839  UniformScaleMap::mapType(),
840  UniformScaleMap::create);
841  }
842 
843  Name type() const { return mapType(); }
844  static Name mapType() { return Name("UniformScaleMap"); }
845 
846  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
847 
848  bool operator==(const UniformScaleMap& other) const { return ScaleMap::operator==(other); }
849  bool operator!=(const UniformScaleMap& other) const { return !(*this == other); }
850 
853  MapBase::Ptr preTranslate(const Vec3d& tr) const;
854 
857  MapBase::Ptr postTranslate(const Vec3d& tr) const;
858 
859 }; // class UniformScaleMap
860 
861 
863 
864 
865 inline MapBase::Ptr
866 ScaleMap::preScale(const Vec3d& v) const
867 {
868  Vec3d new_scale(v * mScaleValues);
869  if (isApproxEqual(new_scale[0],new_scale[1]) && isApproxEqual(new_scale[0],new_scale[2])) {
870  return MapBase::Ptr(new UniformScaleMap(new_scale[0]));
871  } else {
872  return MapBase::Ptr(new ScaleMap(new_scale));
873  }
874 }
875 
876 
877 inline MapBase::Ptr
878 ScaleMap::postScale(const Vec3d& v) const
879 { // pre-post Scale are the same for a scale map
880  return preScale(v);
881 }
882 
883 
886 {
887 public:
888  typedef boost::shared_ptr<TranslationMap> Ptr;
889  typedef boost::shared_ptr<const TranslationMap> ConstPtr;
890 
891  // default constructor is a translation by zero.
892  TranslationMap(): MapBase(), mTranslation(Vec3d(0,0,0)) {}
893  TranslationMap(const Vec3d& t): MapBase(), mTranslation(t) {}
894  TranslationMap(const TranslationMap& other): MapBase(), mTranslation(other.mTranslation) {}
895 
897 
899  static MapBase::Ptr create() { return MapBase::Ptr(new TranslationMap()); }
901  MapBase::Ptr copy() const { return MapBase::Ptr(new TranslationMap(*this)); }
902 
903  static bool isRegistered() { return MapRegistry::isRegistered(TranslationMap::mapType()); }
904 
905  static void registerMap()
906  {
907  MapRegistry::registerMap(
908  TranslationMap::mapType(),
909  TranslationMap::create);
910  }
911 
912  Name type() const { return mapType(); }
913  static Name mapType() { return Name("TranslationMap"); }
914 
916  bool isLinear() const { return true; }
917 
919  bool hasUniformScale() const { return true; }
920 
922  Vec3d applyMap(const Vec3d& in) const { return in + mTranslation; }
924  Vec3d applyInverseMap(const Vec3d& in) const { return in - mTranslation; }
925 
928  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
931  Vec3d applyIJT(const Vec3d& in) const {return in;}
933  Mat3d applyIJC(const Mat3d& mat) const {return mat;}
934  Mat3d applyIJC(const Mat3d& mat, const Vec3d&, const Vec3d&) const {
935  return applyIJC(mat);
936  }
937 
939  double determinant(const Vec3d& ) const { return determinant(); }
941  double determinant() const { return 1.0; }
942 
944  Vec3d voxelSize() const { return Vec3d(1,1,1);}
946  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
947 
949  const Vec3d& getTranslation() const { return mTranslation; }
951  void read(std::istream& is) { mTranslation.read(is); }
953  void write(std::ostream& os) const { mTranslation.write(os); }
954 
956  std::string str() const
957  {
958  std::ostringstream buffer;
959  buffer << " - translation: " << mTranslation << std::endl;
960  return buffer.str();
961  }
962 
963  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
964 
965  bool operator==(const TranslationMap& other) const
966  {
967  // ::eq() uses a tolerance
968  if (!mTranslation.eq(other.mTranslation)) { return false; }
969  return true;
970  }
971 
972  bool operator!=(const TranslationMap& other) const { return !(*this == other); }
973 
975  AffineMap::Ptr getAffineMap() const
976  {
977  Mat4d matrix(Mat4d::identity());
978  matrix.setTranslation(mTranslation);
979 
980  AffineMap::Ptr affineMap(new AffineMap(matrix));
981  return affineMap;
982  }
983 
985 
986 
987  MapBase::Ptr preRotate(double radians, Axis axis) const
988  {
989  AffineMap::Ptr affineMap = getAffineMap();
990  affineMap->accumPreRotation(axis, radians);
991  return simplify(affineMap);
992 
993  }
994  MapBase::Ptr preTranslate(const Vec3d& t) const
995  {
996  return MapBase::Ptr(new TranslationMap(t + mTranslation));
997  }
998 
999  MapBase::Ptr preScale(const Vec3d& v) const;
1000 
1001  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
1002  {
1003  AffineMap::Ptr affineMap = getAffineMap();
1004  affineMap->accumPreShear(axis0, axis1, shear);
1005  return simplify(affineMap);
1006  }
1008 
1010 
1011 
1012  MapBase::Ptr postRotate(double radians, Axis axis) const
1013  {
1014  AffineMap::Ptr affineMap = getAffineMap();
1015  affineMap->accumPostRotation(axis, radians);
1016  return simplify(affineMap);
1017 
1018  }
1019  MapBase::Ptr postTranslate(const Vec3d& t) const
1020  { // post and pre are the same for this
1021  return MapBase::Ptr(new TranslationMap(t + mTranslation));
1022  }
1023 
1024  MapBase::Ptr postScale(const Vec3d& v) const;
1025 
1026  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
1027  {
1028  AffineMap::Ptr affineMap = getAffineMap();
1029  affineMap->accumPostShear(axis0, axis1, shear);
1030  return simplify(affineMap);
1031  }
1033 
1034 private:
1035  Vec3d mTranslation;
1036 }; // class TranslationMap
1037 
1038 
1040 
1041 
1046 {
1047 public:
1048  typedef boost::shared_ptr<ScaleTranslateMap> Ptr;
1049  typedef boost::shared_ptr<const ScaleTranslateMap> ConstPtr;
1050 
1052  MapBase(),
1053  mTranslation(Vec3d(0,0,0)),
1054  mScaleValues(Vec3d(1,1,1)),
1055  mVoxelSize(Vec3d(1,1,1)),
1056  mScaleValuesInverse(Vec3d(1,1,1)),
1057  mInvScaleSqr(1,1,1),
1058  mInvTwiceScale(0.5,0.5,0.5)
1059  {
1060  }
1061 
1062  ScaleTranslateMap(const Vec3d& scale, const Vec3d& translate):
1063  MapBase(),
1064  mTranslation(translate),
1065  mScaleValues(scale),
1066  mVoxelSize(std::abs(scale(0)), std::abs(scale(1)), std::abs(scale(2)))
1067  {
1068  if (scale.eq(Vec3d(0.0), 1.0e-10)) {
1069  OPENVDB_THROW(ArithmeticError, "Non-zero scale values required");
1070  }
1071  mScaleValuesInverse = 1.0 / mScaleValues;
1072  mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1073  mInvTwiceScale = mScaleValuesInverse / 2;
1074  }
1075 
1076  ScaleTranslateMap(const ScaleMap& scale, const TranslationMap& translate):
1077  MapBase(),
1078  mTranslation(translate.getTranslation()),
1079  mScaleValues(scale.getScale()),
1080  mVoxelSize(std::abs(mScaleValues(0)),
1081  std::abs(mScaleValues(1)),
1082  std::abs(mScaleValues(2))),
1083  mScaleValuesInverse(1.0 / scale.getScale())
1084  {
1085  mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1086  mInvTwiceScale = mScaleValuesInverse / 2;
1087  }
1088 
1090  MapBase(),
1091  mTranslation(other.mTranslation),
1092  mScaleValues(other.mScaleValues),
1093  mVoxelSize(other.mVoxelSize),
1094  mScaleValuesInverse(other.mScaleValuesInverse),
1095  mInvScaleSqr(other.mInvScaleSqr),
1096  mInvTwiceScale(other.mInvTwiceScale)
1097  {}
1098 
1100 
1104  MapBase::Ptr copy() const { return MapBase::Ptr(new ScaleTranslateMap(*this)); }
1105 
1106  static bool isRegistered() { return MapRegistry::isRegistered(ScaleTranslateMap::mapType()); }
1107 
1108  static void registerMap()
1109  {
1110  MapRegistry::registerMap(
1111  ScaleTranslateMap::mapType(),
1112  ScaleTranslateMap::create);
1113  }
1114 
1115  Name type() const { return mapType(); }
1116  static Name mapType() { return Name("ScaleTranslateMap"); }
1117 
1119  bool isLinear() const { return true; }
1120 
1123  bool hasUniformScale() const {
1124  bool value;
1125  value = isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.y()));
1126  value = value && isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.z()));
1127  return value;
1128  }
1129 
1131  Vec3d applyMap(const Vec3d& in) const
1132  {
1133  return Vec3d(
1134  in.x() * mScaleValues.x() + mTranslation.x(),
1135  in.y() * mScaleValues.y() + mTranslation.y(),
1136  in.z() * mScaleValues.z() + mTranslation.z());
1137  }
1139  Vec3d applyInverseMap(const Vec3d& in) const
1140  {
1141  return Vec3d(
1142  (in.x() - mTranslation.x() ) / mScaleValues.x(),
1143  (in.y() - mTranslation.y() ) / mScaleValues.y(),
1144  (in.z() - mTranslation.z() ) / mScaleValues.z());
1145  }
1148  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
1150  Vec3d applyIJT(const Vec3d& in) const
1151  {
1152  return Vec3d(
1153  in.x() / mScaleValues.x(),
1154  in.y() / mScaleValues.y(),
1155  in.z() / mScaleValues.z());
1156  }
1158  Mat3d applyIJC(const Mat3d& in) const {
1159  Mat3d tmp;
1160  for (int i=0; i<3; i++){
1161  tmp.setRow(i, in.row(i)*mScaleValuesInverse(i));
1162  }
1163  for (int i=0; i<3; i++){
1164  tmp.setCol(i, tmp.col(i)*mScaleValuesInverse(i));
1165  }
1166  return tmp;
1167  }
1168  Mat3d applyIJC(const Mat3d& in, const Vec3d&, const Vec3d& ) const {
1169  return applyIJC(in);
1170  }
1171 
1173  double determinant(const Vec3d& ) const { return determinant(); }
1175  double determinant() const { return mScaleValues.x()*mScaleValues.y()*mScaleValues.z(); }
1177  Vec3d voxelSize() const { return mVoxelSize;}
1180  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
1181 
1183  const Vec3d& getScale() const { return mScaleValues; }
1185  const Vec3d& getTranslation() const { return mTranslation; }
1186 
1188  const Vec3d& getInvScaleSqr() const {return mInvScaleSqr;}
1190  const Vec3d& getInvTwiceScale() const {return mInvTwiceScale;}
1192  const Vec3d& getInvScale() const {return mScaleValuesInverse; }
1193 
1195  void read(std::istream& is)
1196  {
1197  mTranslation.read(is);
1198  mScaleValues.read(is);
1199  mVoxelSize.read(is);
1200  mScaleValuesInverse.read(is);
1201  mInvScaleSqr.read(is);
1202  mInvTwiceScale.read(is);
1203  }
1205  void write(std::ostream& os) const
1206  {
1207  mTranslation.write(os);
1208  mScaleValues.write(os);
1209  mVoxelSize.write(os);
1210  mScaleValuesInverse.write(os);
1211  mInvScaleSqr.write(os);
1212  mInvTwiceScale.write(os);
1213  }
1215  std::string str() const
1216  {
1217  std::ostringstream buffer;
1218  buffer << " - translation: " << mTranslation << std::endl;
1219  buffer << " - scale: " << mScaleValues << std::endl;
1220  buffer << " - voxel dimensions: " << mVoxelSize << std::endl;
1221  return buffer.str();
1222  }
1223 
1224  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1225 
1226  bool operator==(const ScaleTranslateMap& other) const
1227  {
1228  // ::eq() uses a tolerance
1229  if (!mScaleValues.eq(other.mScaleValues)) { return false; }
1230  if (!mTranslation.eq(other.mTranslation)) { return false; }
1231  return true;
1232  }
1233 
1234  bool operator!=(const ScaleTranslateMap& other) const { return !(*this == other); }
1235 
1237  AffineMap::Ptr getAffineMap() const
1238  {
1239  AffineMap::Ptr affineMap(new AffineMap(math::scale<Mat4d>(mScaleValues)));
1240  affineMap->accumPostTranslation(mTranslation);
1241  return affineMap;
1242  }
1243 
1245 
1246 
1247  MapBase::Ptr preRotate(double radians, Axis axis) const
1248  {
1249  AffineMap::Ptr affineMap = getAffineMap();
1250  affineMap->accumPreRotation(axis, radians);
1251  return simplify(affineMap);
1252  }
1253  MapBase::Ptr preTranslate(const Vec3d& t) const
1254  {
1255  const Vec3d& s = mScaleValues;
1256  const Vec3d scaled_trans( t.x() * s.x(),
1257  t.y() * s.y(),
1258  t.z() * s.z() );
1259  return MapBase::Ptr( new ScaleTranslateMap(mScaleValues, mTranslation + scaled_trans));
1260  }
1261 
1262  MapBase::Ptr preScale(const Vec3d& v) const;
1263 
1264  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
1265  {
1266  AffineMap::Ptr affineMap = getAffineMap();
1267  affineMap->accumPreShear(axis0, axis1, shear);
1268  return simplify(affineMap);
1269  }
1271 
1273 
1274 
1275  MapBase::Ptr postRotate(double radians, Axis axis) const
1276  {
1277  AffineMap::Ptr affineMap = getAffineMap();
1278  affineMap->accumPostRotation(axis, radians);
1279  return simplify(affineMap);
1280  }
1281  MapBase::Ptr postTranslate(const Vec3d& t) const
1282  {
1283  return MapBase::Ptr( new ScaleTranslateMap(mScaleValues, mTranslation + t));
1284  }
1285 
1286  MapBase::Ptr postScale(const Vec3d& v) const;
1287 
1288  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
1289  {
1290  AffineMap::Ptr affineMap = getAffineMap();
1291  affineMap->accumPostShear(axis0, axis1, shear);
1292  return simplify(affineMap);
1293  }
1295 
1296 private:
1297  Vec3d mTranslation, mScaleValues, mVoxelSize, mScaleValuesInverse,
1298  mInvScaleSqr, mInvTwiceScale;
1299 }; // class ScaleTanslateMap
1300 
1301 
1302 inline MapBase::Ptr
1303 ScaleMap::postTranslate(const Vec3d& t) const
1304 {
1305  return MapBase::Ptr(new ScaleTranslateMap(mScaleValues, t));
1306 }
1307 
1308 
1309 inline MapBase::Ptr
1310 ScaleMap::preTranslate(const Vec3d& t) const
1311 {
1312 
1313  const Vec3d& s = mScaleValues;
1314  const Vec3d scaled_trans( t.x() * s.x(),
1315  t.y() * s.y(),
1316  t.z() * s.z() );
1317  return MapBase::Ptr(new ScaleTranslateMap(mScaleValues, scaled_trans));
1318 }
1319 
1320 
1324 {
1325 public:
1326  typedef boost::shared_ptr<UniformScaleTranslateMap> Ptr;
1327  typedef boost::shared_ptr<const UniformScaleTranslateMap> ConstPtr;
1328 
1330  UniformScaleTranslateMap(double scale, const Vec3d& translate):
1331  ScaleTranslateMap(Vec3d(scale,scale,scale), translate) {}
1333  ScaleTranslateMap(scale.getScale(), translate.getTranslation()) {}
1334 
1337 
1341  MapBase::Ptr copy() const { return MapBase::Ptr(new UniformScaleTranslateMap(*this)); }
1342 
1343  static bool isRegistered()
1344  {
1345  return MapRegistry::isRegistered(UniformScaleTranslateMap::mapType());
1346  }
1347 
1348  static void registerMap()
1349  {
1350  MapRegistry::registerMap(
1351  UniformScaleTranslateMap::mapType(),
1352  UniformScaleTranslateMap::create);
1353  }
1354 
1355  Name type() const { return mapType(); }
1356  static Name mapType() { return Name("UniformScaleTranslateMap"); }
1357 
1358  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1359 
1360  bool operator==(const UniformScaleTranslateMap& other) const
1361  {
1362  return ScaleTranslateMap::operator==(other);
1363  }
1364  bool operator!=(const UniformScaleTranslateMap& other) const { return !(*this == other); }
1365 
1368  MapBase::Ptr preTranslate(const Vec3d& t) const
1369  {
1370  const double scale = this->getScale().x();
1371  const Vec3d new_trans = this->getTranslation() + scale * t;
1372  return MapBase::Ptr( new UniformScaleTranslateMap(scale, new_trans));
1373  }
1374 
1377  MapBase::Ptr postTranslate(const Vec3d& t) const
1378  {
1379  const double scale = this->getScale().x();
1380  return MapBase::Ptr( new UniformScaleTranslateMap(scale, this->getTranslation() + t));
1381  }
1382 }; // class UniformScaleTanslateMap
1383 
1384 
1385 inline MapBase::Ptr
1386 UniformScaleMap::postTranslate(const Vec3d& t) const
1387 {
1388  const double scale = this->getScale().x();
1389  return MapBase::Ptr(new UniformScaleTranslateMap(scale, t));
1390 }
1391 
1392 
1393 inline MapBase::Ptr
1394 UniformScaleMap::preTranslate(const Vec3d& t) const
1395 {
1396  const double scale = this->getScale().x();
1397  return MapBase::Ptr(new UniformScaleTranslateMap(scale, scale*t));
1398 }
1399 
1400 inline MapBase::Ptr
1401 TranslationMap::preScale(const Vec3d& v) const
1402 {
1403  if (isApproxEqual(v[0],v[1]) && isApproxEqual(v[0],v[2])) {
1404  return MapBase::Ptr(new UniformScaleTranslateMap(v[0], mTranslation));
1405  } else {
1406  return MapBase::Ptr(new ScaleTranslateMap(v, mTranslation));
1407  }
1408 }
1409 
1410 inline MapBase::Ptr
1411 TranslationMap::postScale(const Vec3d& v) const
1412 {
1413  if (isApproxEqual(v[0],v[1]) && isApproxEqual(v[0],v[2])) {
1414  return MapBase::Ptr(new UniformScaleTranslateMap(v[0], v[0]*mTranslation));
1415  } else {
1416  const Vec3d trans(mTranslation.x()*v.x(),
1417  mTranslation.y()*v.y(),
1418  mTranslation.z()*v.z());
1419  return MapBase::Ptr(new ScaleTranslateMap(v, trans));
1420  }
1421 }
1422 
1423 inline MapBase::Ptr
1424 ScaleTranslateMap::preScale(const Vec3d& v) const
1425 {
1426  Vec3d new_scale( v * mScaleValues );
1427  if (isApproxEqual(new_scale[0],new_scale[1]) && isApproxEqual(new_scale[0],new_scale[2])) {
1428  return MapBase::Ptr( new UniformScaleTranslateMap(new_scale[0], mTranslation));
1429  } else {
1430  return MapBase::Ptr( new ScaleTranslateMap(new_scale, mTranslation));
1431  }
1432 }
1433 
1434 inline MapBase::Ptr
1435 ScaleTranslateMap::postScale(const Vec3d& v) const
1436 {
1437  const Vec3d new_scale( v * mScaleValues );
1438  const Vec3d new_trans( mTranslation.x()*v.x(),
1439  mTranslation.y()*v.y(),
1440  mTranslation.z()*v.z() );
1441 
1442  if (isApproxEqual(new_scale[0],new_scale[1]) && isApproxEqual(new_scale[0],new_scale[2])) {
1443  return MapBase::Ptr( new UniformScaleTranslateMap(new_scale[0], new_trans));
1444  } else {
1445  return MapBase::Ptr( new ScaleTranslateMap(new_scale, new_trans));
1446  }
1447 }
1448 
1449 
1451 
1452 
1456 {
1457 public:
1458  typedef boost::shared_ptr<UnitaryMap> Ptr;
1459  typedef boost::shared_ptr<const UnitaryMap> ConstPtr;
1460 
1462  UnitaryMap(): mAffineMap(Mat4d::identity())
1463  {
1464  }
1465 
1466  UnitaryMap(const Vec3d& axis, double radians)
1467  {
1468  Mat3d matrix;
1469  matrix.setToRotation(axis, radians);
1470  mAffineMap = AffineMap(matrix);
1471  }
1472 
1473  UnitaryMap(Axis axis, double radians)
1474  {
1475  Mat4d matrix;
1476  matrix.setToRotation(axis, radians);
1477  mAffineMap = AffineMap(matrix);
1478  }
1479 
1480  UnitaryMap(const Mat3d& m)
1481  {
1482  // test that the mat3 is a rotation || reflection
1483  if (!isUnitary(m)) {
1484  OPENVDB_THROW(ArithmeticError, "Matrix initializing unitary map was not unitary");
1485  }
1486 
1487  Mat4d matrix(Mat4d::identity());
1488  matrix.setMat3(m);
1489  mAffineMap = AffineMap(matrix);
1490  }
1491 
1492  UnitaryMap(const Mat4d& m)
1493  {
1494  if (!isInvertible(m)) {
1496  "4x4 Matrix initializing unitary map was not unitary: not invertible");
1497  }
1498 
1499  if (!isAffine(m)) {
1501  "4x4 Matrix initializing unitary map was not unitary: not affine");
1502  }
1503 
1504  if (hasTranslation(m)) {
1506  "4x4 Matrix initializing unitary map was not unitary: had translation");
1507  }
1508 
1509  if (!isUnitary(m.getMat3())) {
1511  "4x4 Matrix initializing unitary map was not unitary");
1512  }
1513 
1514  mAffineMap = AffineMap(m);
1515  }
1516 
1517  UnitaryMap(const UnitaryMap& other):
1518  MapBase(other),
1519  mAffineMap(other.mAffineMap)
1520  {
1521  }
1522 
1523  UnitaryMap(const UnitaryMap& first, const UnitaryMap& second):
1524  mAffineMap(*(first.getAffineMap()), *(second.getAffineMap()))
1525  {
1526  }
1527 
1530  static MapBase::Ptr create() { return MapBase::Ptr(new UnitaryMap()); }
1532  MapBase::Ptr copy() const { return MapBase::Ptr(new UnitaryMap(*this)); }
1533 
1534  static bool isRegistered() { return MapRegistry::isRegistered(UnitaryMap::mapType()); }
1535 
1536  static void registerMap()
1537  {
1538  MapRegistry::registerMap(
1539  UnitaryMap::mapType(),
1540  UnitaryMap::create);
1541  }
1542 
1544  Name type() const { return mapType(); }
1546  static Name mapType() { return Name("UnitaryMap"); }
1547 
1549  bool isLinear() const { return true; }
1550 
1552  bool hasUniformScale() const { return true; }
1553 
1554  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1555 
1556  bool operator==(const UnitaryMap& other) const
1557  {
1558  // compare underlying linear map.
1559  if (mAffineMap!=other.mAffineMap) return false;
1560  return true;
1561  }
1562 
1563  bool operator!=(const UnitaryMap& other) const { return !(*this == other); }
1565  Vec3d applyMap(const Vec3d& in) const { return mAffineMap.applyMap(in); }
1567  Vec3d applyInverseMap(const Vec3d& in) const { return mAffineMap.applyInverseMap(in); }
1570  Vec3d applyIJT(const Vec3d& in, const Vec3d& ) const { return applyIJT(in);}
1572  Vec3d applyIJT(const Vec3d& in) const { return mAffineMap.applyIJT(in); }
1574  Mat3d applyIJC(const Mat3d& in) const {
1575  return mAffineMap.applyIJC(in);
1576  }
1577  Mat3d applyIJC(const Mat3d& in, const Vec3d&, const Vec3d& ) const {
1578  return applyIJC(in);
1579  }
1581  double determinant(const Vec3d& ) const { return determinant(); }
1583  double determinant() const { return mAffineMap.determinant(); }
1584 
1585 
1590  Vec3d voxelSize() const { return mAffineMap.voxelSize();}
1591  Vec3d voxelSize(const Vec3d&) const { return voxelSize();}
1592 
1594  void read(std::istream& is)
1595  {
1596  mAffineMap.read(is);
1597  }
1598 
1600  void write(std::ostream& os) const
1601  {
1602  mAffineMap.write(os);
1603  }
1605  std::string str() const
1606  {
1607  std::ostringstream buffer;
1608  buffer << mAffineMap.str();
1609  return buffer.str();
1610  }
1612  AffineMap::Ptr getAffineMap() const { return AffineMap::Ptr(new AffineMap(mAffineMap)); }
1613 
1615 
1616 
1617  MapBase::Ptr preRotate(double radians, Axis axis) const
1618  {
1619  UnitaryMap first(axis, radians);
1620  UnitaryMap::Ptr unitaryMap(new UnitaryMap(first, *this));
1621  return boost::static_pointer_cast<MapBase, UnitaryMap>(unitaryMap);
1622  }
1623  MapBase::Ptr preTranslate(const Vec3d& t) const
1624  {
1625  AffineMap::Ptr affineMap = getAffineMap();
1626  affineMap->accumPreTranslation(t);
1627  return simplify(affineMap);
1628  }
1629  MapBase::Ptr preScale(const Vec3d& v) const
1630  {
1631  AffineMap::Ptr affineMap = getAffineMap();
1632  affineMap->accumPreScale(v);
1633  return simplify(affineMap);
1634  }
1635  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
1636  {
1637  AffineMap::Ptr affineMap = getAffineMap();
1638  affineMap->accumPreShear(axis0, axis1, shear);
1639  return simplify(affineMap);
1640  }
1642 
1643 
1645 
1646 
1647  MapBase::Ptr postRotate(double radians, Axis axis) const
1648  {
1649  UnitaryMap second(axis, radians);
1650  UnitaryMap::Ptr unitaryMap(new UnitaryMap(*this, second));
1651  return boost::static_pointer_cast<MapBase, UnitaryMap>(unitaryMap);
1652  }
1653  MapBase::Ptr postTranslate(const Vec3d& t) const
1654  {
1655  AffineMap::Ptr affineMap = getAffineMap();
1656  affineMap->accumPostTranslation(t);
1657  return simplify(affineMap);
1658  }
1659  MapBase::Ptr postScale(const Vec3d& v) const
1660  {
1661  AffineMap::Ptr affineMap = getAffineMap();
1662  affineMap->accumPostScale(v);
1663  return simplify(affineMap);
1664  }
1665  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
1666  {
1667  AffineMap::Ptr affineMap = getAffineMap();
1668  affineMap->accumPostShear(axis0, axis1, shear);
1669  return simplify(affineMap);
1670  }
1672 
1673 private:
1674  AffineMap mAffineMap;
1675 }; // class UnitaryMap
1676 
1677 
1679 
1680 
1688 {
1689 public:
1690  typedef boost::shared_ptr<NonlinearFrustumMap> Ptr;
1691  typedef boost::shared_ptr<const NonlinearFrustumMap> ConstPtr;
1692 
1694  MapBase(),
1695  mBBox(Vec3d(0), Vec3d(1)),
1696  mTaper(1),
1697  mDepth(1)
1698  {
1699  init();
1700  }
1701 
1705  NonlinearFrustumMap(const BBoxd& bb, double taper, double depth):
1706  MapBase(),mBBox(bb), mTaper(taper), mDepth(depth)
1707  {
1708  init();
1709  }
1710 
1716  NonlinearFrustumMap(const BBoxd& bb, double taper, double depth,
1717  const MapBase::Ptr& secondMap):
1718  mBBox(bb), mTaper(taper), mDepth(depth)
1719  {
1720  if (!secondMap->isLinear() ) {
1722  "The second map in the Frustum transfrom must be linear");
1723  }
1724  mSecondMap = *( secondMap->getAffineMap() );
1725  init();
1726  }
1727 
1729  MapBase(),
1730  mBBox(other.mBBox),
1731  mTaper(other.mTaper),
1732  mDepth(other.mDepth),
1733  mSecondMap(other.mSecondMap),
1734  mHasSimpleAffine(other.mHasSimpleAffine)
1735  {
1736  init();
1737  }
1738 
1754  NonlinearFrustumMap(const Vec3d& position,
1755  const Vec3d& direction,
1756  const Vec3d& up,
1757  double aspect /* width / height */,
1758  double z_near, double depth,
1759  Coord::ValueType x_count, Coord::ValueType z_count) {
1760 
1764  if (!(depth > 0)) {
1766  "The frustum depth must be non-zero and positive");
1767  }
1768  if (!(up.length() > 0)) {
1770  "The frustum height must be non-zero and positive");
1771  }
1772  if (!(aspect > 0)) {
1774  "The frustum aspect ratio must be non-zero and positive");
1775  }
1776  if (!(isApproxEqual(up.dot(direction), 0.))) {
1778  "The frustum up orientation must be perpendicular to into-frustum direction");
1779  }
1780 
1781  double near_plane_height = 2 * up.length();
1782  double near_plane_width = aspect * near_plane_height;
1783 
1784  Coord::ValueType y_count = static_cast<int>(Round(x_count / aspect));
1785 
1786  mBBox = BBoxd(Vec3d(0,0,0), Vec3d(x_count, y_count, z_count));
1787  mDepth = depth / near_plane_width; // depth non-dimensionalized on width
1788  double gamma = near_plane_width / z_near;
1789  mTaper = 1./(mDepth*gamma + 1.);
1790 
1791  Vec3d direction_unit = direction;
1792  direction_unit.normalize();
1793 
1794  Mat4d r1(Mat4d::identity());
1795  r1.setToRotation(/*from*/Vec3d(0,0,1), /*to */direction_unit);
1796  Mat4d r2(Mat4d::identity());
1797  Vec3d temp = r1.inverse().transform(up);
1798  r2.setToRotation(/*from*/Vec3d(0,1,0), /*to*/temp );
1799  Mat4d scale = math::scale<Mat4d>(
1800  Vec3d(near_plane_width, near_plane_width, near_plane_width));
1801 
1802  // move the near plane to origin, rotate to align with axis, and scale down
1803  // T_inv * R1_inv * R2_inv * scale_inv
1804  Mat4d mat = scale * r2 * r1;
1805  mat.setTranslation(position + z_near*direction_unit);
1806 
1807  mSecondMap = AffineMap(mat);
1808 
1809  init();
1810  }
1811 
1816  MapBase::Ptr copy() const { return MapBase::Ptr(new NonlinearFrustumMap(*this)); }
1817 
1818  static bool isRegistered() { return MapRegistry::isRegistered(NonlinearFrustumMap::mapType()); }
1819 
1820  static void registerMap()
1821  {
1822  MapRegistry::registerMap(
1823  NonlinearFrustumMap::mapType(),
1824  NonlinearFrustumMap::create);
1825  }
1827  Name type() const { return mapType(); }
1829  static Name mapType() { return Name("NonlinearFrustumMap"); }
1830 
1832  bool isLinear() const { return false; }
1833 
1835  bool hasUniformScale() const { return false; }
1836 
1838  bool isIdentity() const {
1839 
1840  // The frustum can only be consistent with a linear map if the taper value is 1
1841  if (!isApproxEqual(mTaper, double(1)) ) return false;
1842 
1843  // There are various ways an identity can decomposed between the two parts of the
1844  // map. Best to just check that the principle vectors are stationary.
1845  const Vec3d e1(1,0,0);
1846  if (!applyMap(e1).eq(e1)) return false;
1847 
1848  const Vec3d e2(0,1,0);
1849  if (!applyMap(e2).eq(e2)) return false;
1850 
1851  const Vec3d e3(0,0,1);
1852  if (!applyMap(e3).eq(e3)) return false;
1853 
1854  return true;
1855  }
1856 
1857  virtual bool isEqual(const MapBase& other) const { return isEqualBase(*this, other); }
1858 
1859  bool operator==(const NonlinearFrustumMap& other) const
1860  {
1861  if (mBBox!=other.mBBox) return false;
1862  if (!isApproxEqual(mTaper, other.mTaper)) return false;
1863  if (!isApproxEqual(mDepth, other.mDepth)) return false;
1864 
1865  // Two linear transforms are equivalent iff they have the same translation
1866  // and have the same affects on orthongal spanning basis check translation
1867  Vec3d e(0,0,0);
1868  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
1870  e(0) = 1;
1871  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
1872  e(0) = 0;
1873  e(1) = 1;
1874  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
1875  e(1) = 0;
1876  e(2) = 1;
1877  if (!mSecondMap.applyMap(e).eq(other.mSecondMap.applyMap(e))) return false;
1878  return true;
1879  }
1880 
1881  bool operator!=(const NonlinearFrustumMap& other) const { return !(*this == other); }
1882 
1884  Vec3d applyMap(const Vec3d& in) const
1885  {
1886  return mSecondMap.applyMap(applyFrustumMap(in));
1887  }
1888 
1890  Vec3d applyInverseMap(const Vec3d& in) const
1891  {
1892  return applyFrustumInverseMap(mSecondMap.applyInverseMap(in));
1893  }
1894 
1896  Vec3d applyIJT(const Vec3d& in) const { return mSecondMap.applyIJT(in); }
1897 
1898  // the Jacobian of the nonlinear part of the transform is a sparse matrix
1899  // Jacobian^(-T) =
1900  //
1901  // (Lx)( 1/s 0 0 )
1902  // ( 0 1/s 0 )
1903  // ( -(x-xo)g/(sLx) -(y-yo)g/(sLx) Lz/(Depth Lx) )
1906  Vec3d applyIJT(const Vec3d& d1_is, const Vec3d& ijk) const
1907  {
1908  const Vec3d loc = applyFrustumMap(ijk);
1909  const double s = mGamma * loc.z() + 1.;
1910 
1911  // verify that we aren't at the singularity
1912  if (isApproxEqual(s, 0.)) {
1913  OPENVDB_THROW(ArithmeticError, "Tried to evaluate the frustum transform"
1914  " at the singular focal point (e.g. camera)");
1915  }
1916 
1917  const double sinv = 1.0/s; // 1/(z*gamma + 1)
1918  const double pt0 = mLx * sinv; // Lx / (z*gamma +1)
1919  const double pt1 = mGamma * pt0; // gamma * Lx / ( z*gamma +1)
1920  const double pt2 = pt1 * sinv; // gamma * Lx / ( z*gamma +1)**2
1921 
1922  const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
1923 
1924  // compute \frac{\partial E_i}{\partial x_j}
1925  Mat3d gradE(Mat3d::zero());
1926  for (int j = 0; j < 3; ++j ) {
1927  gradE(0,j) = pt0 * jacinv(0,j) - pt2 * loc.x()*jacinv(2,j);
1928  gradE(1,j) = pt0 * jacinv(1,j) - pt2 * loc.y()*jacinv(2,j);
1929  gradE(2,j) = (1./mDepthOnLz) * jacinv(2,j);
1930  }
1931 
1932  Vec3d result;
1933  for (int i = 0; i < 3; ++i) {
1934  result(0) = d1_is(0) * gradE(0,i) + d1_is(1) * gradE(1,i) + d1_is(2) * gradE(2,i);
1935  }
1936 
1937  return result;
1938 
1939  }
1940 
1942  Mat3d applyIJC(const Mat3d& in) const { return mSecondMap.applyIJC(in); }
1947  Mat3d applyIJC(const Mat3d& d2_is, const Vec3d& d1_is, const Vec3d& ijk) const
1948  {
1949  const Vec3d loc = applyFrustumMap(ijk);
1950 
1951  const double s = mGamma * loc.z() + 1.;
1952 
1953  // verify that we aren't at the singularity
1954  if (isApproxEqual(s, 0.)) {
1955  OPENVDB_THROW(ArithmeticError, "Tried to evaluate the frustum transform"
1956  " at the singular focal point (e.g. camera)");
1957  }
1958 
1959  // precompute
1960  const double sinv = 1.0/s; // 1/(z*gamma + 1)
1961  const double pt0 = mLx * sinv; // Lx / (z*gamma +1)
1962  const double pt1 = mGamma * pt0; // gamma * Lx / ( z*gamma +1)
1963  const double pt2 = pt1 * sinv; // gamma * Lx / ( z*gamma +1)**2
1964  const double pt3 = pt2 * sinv; // gamma * Lx / ( z*gamma +1)**3
1965 
1966 
1967  const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
1968 
1969  // compute \frac{\partial^2 E_i}{\partial x_j \partial x_k}
1970 
1971  Mat3d matE0(Mat3d::zero());
1972  Mat3d matE1(Mat3d::zero()); // matE2 = 0
1973  for(int j = 0; j < 3; j++) {
1974  for (int k = 0; k < 3; k++) {
1975 
1976  const double pt4 = 2. * jacinv(2,j) * jacinv(2,k) * pt3;
1977 
1978  matE0(j,k) = -(jacinv(0,j) * jacinv(2,k) + jacinv(2,j) * jacinv(0,k)) * pt2 +
1979  pt4 * loc.x();
1980 
1981  matE1(j,k) = -(jacinv(1,j) * jacinv(2,k) + jacinv(2,j) * jacinv(1,k)) * pt2 +
1982  pt4 * loc.y();
1983  }
1984  }
1985 
1986  // compute \frac{\partial E_i}{\partial x_j}
1987  Mat3d gradE(Mat3d::zero());
1988  for (int j = 0; j < 3; ++j ) {
1989  gradE(0,j) = pt0 * jacinv(0,j) - pt2 * loc.x()*jacinv(2,j);
1990  gradE(1,j) = pt0 * jacinv(1,j) - pt2 * loc.y()*jacinv(2,j);
1991  gradE(2,j) = (1./mDepthOnLz) * jacinv(2,j);
1992  }
1993 
1994  Mat3d result(Mat3d::zero());
1995  // compute \fac{\partial E_j}{\partial x_m} \fac{\partial E_i}{\partial x_n}
1996  // \frac{\partial^2 input}{\partial E_i \partial E_j}
1997  for (int m = 0; m < 3; ++m ) {
1998  for ( int n = 0; n < 3; ++n) {
1999  for (int i = 0; i < 3; ++i ) {
2000  for (int j = 0; j < 3; ++j) {
2001  result(m, n) += gradE(j, m) * gradE(i, n) * d2_is(i, j);
2002  }
2003  }
2004  }
2005  }
2006 
2007  for (int m = 0; m < 3; ++m ) {
2008  for ( int n = 0; n < 3; ++n) {
2009  result(m, n) +=
2010  matE0(m, n) * d1_is(0) + matE1(m, n) * d1_is(1);// + matE2(m, n) * d1_is(2);
2011  }
2012  }
2013 
2014  return result;
2015  }
2016 
2018  double determinant() const {return mSecondMap.determinant();} // no implementation
2019 
2022  double determinant(const Vec3d& loc) const
2023  {
2024  double s = mGamma * loc.z() + 1.0;
2025  double frustum_determinant = s * s * mDepthOnLzLxLx;
2026  return mSecondMap.determinant() * frustum_determinant;
2027  }
2028 
2030  Vec3d voxelSize() const {
2031  const Vec3d loc( 0.5*(mBBox.min().x() + mBBox.max().x()),
2032  0.5*(mBBox.min().y() + mBBox.max().y()),
2033  mBBox.min().z());
2034 
2035  return voxelSize(loc);
2036 
2037  }
2038 
2043  Vec3d voxelSize(const Vec3d& loc) const {
2044  Vec3d out, pos = applyMap(loc);
2045  out(0) = (applyMap(loc + Vec3d(1,0,0)) - pos).length();
2046  out(1) = (applyMap(loc + Vec3d(0,1,0)) - pos).length();
2047  out(2) = (applyMap(loc + Vec3d(0,0,1)) - pos).length();
2048  return out;
2049  }
2050 
2051  AffineMap::Ptr getAffineMap() const
2052  {
2053  return mSecondMap.getAffineMap();
2054  }
2055 
2057  void setTaper(double t) { mTaper = t; init();}
2059  double getTaper() const { return mTaper; }
2061  void setDepth(double d) { mDepth = d; init();}
2063  double getDepth() const { return mDepth; }
2064  // gamma a non-dimensional number: nearplane x-width / camera to near plane distance
2065  double getGamma() const { return mGamma; }
2066 
2068  const BBoxd& getBBox() const { return mBBox; }
2069 
2071  const AffineMap& secondMap() const { return mSecondMap; }
2074  bool isValid() const { return !mBBox.empty();}
2075 
2077  bool hasSimpleAffine() const { return mHasSimpleAffine; }
2078 
2080  void read(std::istream& is)
2081  {
2082  // for backward compatibility with earlier version
2084  CoordBBox bb;
2085  bb.read(is);
2086  mBBox = BBoxd(bb.min().asVec3d(), bb.max().asVec3d());
2087  } else {
2088  mBBox.read(is);
2089  }
2090 
2091  is.read(reinterpret_cast<char*>(&mTaper), sizeof(double));
2092  is.read(reinterpret_cast<char*>(&mDepth), sizeof(double));
2093 
2094  // Read the second maps type.
2095  Name type = readString(is);
2096 
2097  // Check if the map has been registered.
2098  if(!MapRegistry::isRegistered(type)) {
2099  OPENVDB_THROW(KeyError, "Map " << type << " is not registered");
2100  }
2101 
2102  // Create the second map of the type and then read it in.
2103  MapBase::Ptr proxy = math::MapRegistry::createMap(type);
2104  proxy->read(is);
2105  mSecondMap = *(proxy->getAffineMap());
2106  init();
2107  }
2108 
2110  void write(std::ostream& os) const
2111  {
2112  mBBox.write(os);
2113  os.write(reinterpret_cast<const char*>(&mTaper), sizeof(double));
2114  os.write(reinterpret_cast<const char*>(&mDepth), sizeof(double));
2115 
2116  writeString(os, mSecondMap.type());
2117  mSecondMap.write(os);
2118  }
2119 
2121  std::string str() const
2122  {
2123  std::ostringstream buffer;
2124  buffer << " - taper: " << mTaper << std::endl;
2125  buffer << " - depth: " << mDepth << std::endl;
2126  buffer << " SecondMap: "<< mSecondMap.type() << std::endl;
2127  buffer << mSecondMap.str() << std::endl;
2128  return buffer.str();
2129  }
2130 
2132 
2133 
2134  MapBase::Ptr preRotate(double radians, Axis axis = X_AXIS) const
2135  {
2136  return MapBase::Ptr(
2137  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.preRotate(radians, axis)));
2138  }
2139  MapBase::Ptr preTranslate(const Vec3d& t) const
2140  {
2141  return MapBase::Ptr(
2142  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.preTranslate(t)));
2143  }
2144  MapBase::Ptr preScale(const Vec3d& s) const
2145  {
2146  return MapBase::Ptr(
2147  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.preScale(s)));
2148  }
2149  MapBase::Ptr preShear(double shear, Axis axis0, Axis axis1) const
2150  {
2151  return MapBase::Ptr(new NonlinearFrustumMap(
2152  mBBox, mTaper, mDepth, mSecondMap.preShear(shear, axis0, axis1)));
2153  }
2155 
2157 
2158 
2159  MapBase::Ptr postRotate(double radians, Axis axis = X_AXIS) const
2160  {
2161  return MapBase::Ptr(
2162  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.postRotate(radians, axis)));
2163  }
2164  MapBase::Ptr postTranslate(const Vec3d& t) const
2165  {
2166  return MapBase::Ptr(
2167  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.postTranslate(t)));
2168  }
2169  MapBase::Ptr postScale(const Vec3d& s) const
2170  {
2171  return MapBase::Ptr(
2172  new NonlinearFrustumMap(mBBox, mTaper, mDepth, mSecondMap.postScale(s)));
2173  }
2174  MapBase::Ptr postShear(double shear, Axis axis0, Axis axis1) const
2175  {
2176  return MapBase::Ptr(new NonlinearFrustumMap(
2177  mBBox, mTaper, mDepth, mSecondMap.postShear(shear, axis0, axis1)));
2178  }
2180 
2181 
2182 private:
2183  void init() {
2184  // set up as a frustum
2185  mLx = mBBox.extents().x();
2186  mLy = mBBox.extents().y();
2187  mLz = mBBox.extents().z();
2188 
2189  if (isApproxEqual(mLx,0.) || isApproxEqual(mLy,0.) || isApproxEqual(mLz,0.) ) {
2190  OPENVDB_THROW(ArithmeticError, "The index space bounding box"
2191  " must have at least two index points in each direction.");
2192  }
2193 
2194  mXo = 0.5* mLx;
2195  mYo = 0.5* mLy;
2196 
2197  // mDepth is non-dimensionalized on near
2198  mGamma = (1./mTaper - 1) / mDepth;
2199 
2200  mDepthOnLz = mDepth/mLz;
2201  mDepthOnLzLxLx = mDepthOnLz/(mLx * mLx);
2202 
2204  mHasSimpleAffine = true;
2205  Vec3d tmp = mSecondMap.voxelSize();
2206 
2208  if (!isApproxEqual(tmp(0), tmp(1))) { mHasSimpleAffine = false; return; }
2209  if (!isApproxEqual(tmp(0), tmp(2))) { mHasSimpleAffine = false; return; }
2210 
2211  Vec3d trans = mSecondMap.applyMap(Vec3d(0,0,0));
2213  Vec3d tmp1 = mSecondMap.applyMap(Vec3d(1,0,0)) - trans;
2214  Vec3d tmp2 = mSecondMap.applyMap(Vec3d(0,1,0)) - trans;
2215  Vec3d tmp3 = mSecondMap.applyMap(Vec3d(0,0,1)) - trans;
2216 
2218  if (!isApproxEqual(tmp1.dot(tmp2), 0., 1.e-7)) { mHasSimpleAffine = false; return; }
2219  if (!isApproxEqual(tmp2.dot(tmp3), 0., 1.e-7)) { mHasSimpleAffine = false; return; }
2220  if (!isApproxEqual(tmp3.dot(tmp1), 0., 1.e-7)) { mHasSimpleAffine = false; return; }
2221  }
2222 
2223  Vec3d applyFrustumMap(const Vec3d& in) const
2224  {
2225 
2226  // Move the center of the x-face of the bbox
2227  // to the origin in index space.
2228  Vec3d out(in);
2229  out = out - mBBox.min();
2230  out.x() -= mXo;
2231  out.y() -= mYo;
2232 
2233  // scale the z-direction on depth / K count
2234  out.z() *= mDepthOnLz;
2235 
2236  double scale = (mGamma * out.z() + 1.)/ mLx;
2237 
2238  // scale the x-y on the length I count and apply tapper
2239  out.x() *= scale ;
2240  out.y() *= scale ;
2241 
2242  return out;
2243  }
2244 
2245  Vec3d applyFrustumInverseMap(const Vec3d& in) const
2246  {
2247  // invert taper and resize: scale = 1/( (z+1)/2 (mt-1) + 1)
2248  Vec3d out(in);
2249  double invScale = mLx / (mGamma * out.z() + 1.);
2250  out.x() *= invScale;
2251  out.y() *= invScale;
2252 
2253  out.x() += mXo;
2254  out.y() += mYo;
2255 
2256  out.z() /= mDepthOnLz;
2257 
2258  // move back
2259  out = out + mBBox.min();
2260  return out;
2261  }
2262 
2263  // bounding box in index space used in Frustum transforms.
2264  BBoxd mBBox;
2265 
2266  // taper value used in constructing Frustums.
2267  double mTaper;
2268  double mDepth;
2269 
2270  // defines the second map
2271  AffineMap mSecondMap;
2272 
2273  // these are derived from the above.
2274  double mLx, mLy, mLz;
2275  double mXo, mYo, mGamma, mDepthOnLz, mDepthOnLzLxLx;
2276 
2277  // true: if the mSecondMap is linear and has no shear, and has no non-uniform scale
2278  bool mHasSimpleAffine;
2279 }; // class NonlinearFrustumMap
2280 
2281 
2283 
2284 
2288 template<typename FirstMapType, typename SecondMapType>
2289 class CompoundMap
2290 {
2291 public:
2293 
2294  typedef boost::shared_ptr<MyType> Ptr;
2295  typedef boost::shared_ptr<const MyType> ConstPtr;
2296 
2297 
2298  CompoundMap() { updateAffineMatrix(); }
2299 
2300  CompoundMap(const FirstMapType& f, const SecondMapType& s): mFirstMap(f), mSecondMap(s)
2301  {
2302  updateAffineMatrix();
2303  }
2304 
2305  CompoundMap(const MyType& other):
2306  mFirstMap(other.mFirstMap),
2307  mSecondMap(other.mSecondMap),
2308  mAffineMap(other.mAffineMap)
2309  {}
2310 
2311  Name type() const { return mapType(); }
2312  static Name mapType()
2313  {
2314  return (FirstMapType::mapType() + Name(":") + SecondMapType::mapType());
2315  }
2316 
2317  bool operator==(const MyType& other) const
2318  {
2319  if (mFirstMap != other.mFirstMap) return false;
2320  if (mSecondMap != other.mSecondMap) return false;
2321  if (mAffineMap != other.mAffineMap) return false;
2322  return true;
2323  }
2324 
2325  bool operator!=(const MyType& other) const { return !(*this == other); }
2326 
2327  MyType& operator=(const MyType& other)
2328  {
2329  mFirstMap = other.mFirstMap;
2330  mSecondMap = other.mSecondMap;
2331  mAffineMap = other.mAffineMap;
2332  return *this;
2333  }
2334 
2335  bool isIdentity() const {
2337  return mAffineMap.isIdentity();
2338  } else {
2339  return mFirstMap.isIdentity()&&mSecondMap.isIdentity();
2340  }
2341  }
2342 
2343  bool isDiagonal() const {
2345  return mAffineMap.isDiagonal();
2346  } else {
2347  return mFirstMap.isDiagonal()&&mSecondMap.isDiagonal();
2348  }
2349  }
2350  AffineMap::Ptr getAffineMap() const
2351  {
2353  AffineMap::Ptr affine(new AffineMap(mAffineMap));
2354  return affine;
2355  } else {
2357  "Constant affine matrix representation not possible for this nonlinear map");
2358  }
2359  }
2360 
2361  // direct decompotion
2362  const FirstMapType& firstMap() const { return mFirstMap; }
2363  const SecondMapType& secondMap() const {return mSecondMap; }
2364 
2365  void setFirstMap(const FirstMapType& first) { mFirstMap = first; updateAffineMatrix(); }
2366  void setSecondMap(const SecondMapType& second) { mSecondMap = second; updateAffineMatrix(); }
2367 
2368  void read(std::istream& is)
2369  {
2370  mAffineMap.read(is);
2371  mFirstMap.read(is);
2372  mSecondMap.read(is);
2373  }
2374  void write(std::ostream& os) const
2375  {
2376  mAffineMap.write(os);
2377  mFirstMap.write(os);
2378  mSecondMap.write(os);
2379  }
2380 
2381 private:
2382  void updateAffineMatrix()
2383  {
2385  // both maps need to be linear, these methods are only defined for linear maps
2386  AffineMap::Ptr first = mFirstMap.getAffineMap();
2387  AffineMap::Ptr second= mSecondMap.getAffineMap();
2388  mAffineMap = AffineMap(*first, *second);
2389  }
2390  }
2391 
2392  FirstMapType mFirstMap;
2393  SecondMapType mSecondMap;
2394  // used for acceleration
2395  AffineMap mAffineMap;
2396 }; // class CompoundMap
2397 
2398 } // namespace math
2399 } // namespace OPENVDB_VERSION_NAME
2400 } // namespace openvdb
2401 
2402 #endif // OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED
2403 
2404 // Copyright (c) 2012-2013 DreamWorks Animation LLC
2405 // All rights reserved. This software is distributed under the
2406 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )