33 #ifndef OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED
34 #define OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED
40 #include <openvdb/util/Name.h>
41 #include <openvdb/Types.h>
42 #include <boost/shared_ptr.hpp>
43 #include <tbb/mutex.h>
59 class ScaleTranslateMap;
60 class UniformScaleMap;
61 class UniformScaleTranslateMap;
64 class NonlinearFrustumMap;
79 template<
typename T>
struct is_linear {
static const bool value =
false; };
99 static const bool value =
true;
103 template<
typename T>
struct is_scale {
static const bool value =
false; };
157 typedef boost::shared_ptr<MapBase>
Ptr;
159 typedef Ptr (*MapFactory)();
163 virtual boost::shared_ptr<AffineMap> getAffineMap()
const = 0;
166 virtual Name type()
const = 0;
169 template<
typename MapT>
bool isType()
const {
return this->type() == MapT::mapType(); }
172 virtual bool isEqual(
const MapBase& other)
const = 0;
175 virtual bool isLinear()
const = 0;
177 virtual bool hasUniformScale()
const = 0;
179 virtual Vec3d applyMap(
const Vec3d& in)
const = 0;
180 virtual Vec3d applyInverseMap(
const Vec3d& in)
const = 0;
182 virtual Vec3d applyIJT(
const Vec3d& in)
const = 0;
184 virtual Mat3d applyIJC(
const Mat3d& m)
const = 0;
188 virtual double determinant()
const = 0;
189 virtual double determinant(
const Vec3d&)
const = 0;
194 virtual Vec3d voxelSize()
const = 0;
195 virtual Vec3d voxelSize(
const Vec3d&)
const = 0;
201 virtual Vec3d voxelDimensions()
const = 0;
202 virtual Vec3d voxelDimensions(
const Vec3d&)
const = 0;
205 virtual void read(std::istream&) = 0;
206 virtual void write(std::ostream&)
const = 0;
208 virtual std::string str()
const = 0;
238 template<
typename MapT>
239 static bool isEqualBase(
const MapT&
self,
const MapBase& other)
241 return other.
isType<MapT>() && (
self == *static_cast<const MapT*>(&other));
256 typedef Mutex::scoped_lock
Lock;
264 static bool isRegistered(
const Name&);
267 static void registerMap(
const Name&, MapBase::MapFactory);
270 static void unregisterMap(
const Name&);
292 typedef boost::shared_ptr<AffineMap>
Ptr;
293 typedef boost::shared_ptr<const AffineMap>
ConstPtr;
296 mMatrix(
Mat4d::identity()),
297 mMatrixInv(
Mat4d::identity()),
298 mJacobianInv(
Mat3d::identity()),
300 mVoxelSize(
Vec3d(1,1,1)),
312 updateAcceleration();
319 "Tried to initialize an affine transform from a non-affine 4x4 matrix");
321 updateAcceleration();
326 mMatrix(other.mMatrix),
327 mMatrixInv(other.mMatrixInv),
328 mJacobianInv(other.mJacobianInv),
329 mDeterminant(other.mDeterminant),
330 mVoxelSize(other.mVoxelSize),
331 mIsDiagonal(other.mIsDiagonal),
332 mIsIdentity(other.mIsIdentity)
338 mMatrix(first.mMatrix * second.mMatrix)
340 updateAcceleration();
350 static bool isRegistered() {
return MapRegistry::isRegistered(AffineMap::mapType()); }
352 static void registerMap()
354 MapRegistry::registerMap(
355 AffineMap::mapType(),
368 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
373 if (!mMatrix.eq(other.mMatrix)) {
return false; }
374 if (!mMatrixInv.eq(other.mMatrixInv)) {
return false; }
382 mMatrix = other.mMatrix;
383 mMatrixInv = other.mMatrixInv;
385 mJacobianInv = other.mJacobianInv;
386 mDeterminant = other.mDeterminant;
387 mVoxelSize = other.mVoxelSize;
388 mIsDiagonal = other.mIsDiagonal;
389 mIsIdentity = other.mIsIdentity;
403 return mJacobianInv.
transpose()* m * mJacobianInv;
440 void accumPreRotation(
Axis axis,
double radians)
442 mMatrix.preRotate(axis, radians);
443 updateAcceleration();
448 updateAcceleration();
450 void accumPreTranslation(
const Vec3d& v)
452 mMatrix.preTranslate(v);
453 updateAcceleration();
457 mMatrix.preShear(axis0, axis1, shear);
458 updateAcceleration();
466 void accumPostRotation(
Axis axis,
double radians)
468 mMatrix.postRotate(axis, radians);
469 updateAcceleration();
473 mMatrix.postScale(v);
474 updateAcceleration();
476 void accumPostTranslation(
const Vec3d& v)
478 mMatrix.postTranslate(v);
479 updateAcceleration();
483 mMatrix.postShear(axis0, axis1, shear);
484 updateAcceleration();
490 void read(std::istream& is)
493 updateAcceleration();
497 void write(std::ostream& os)
const
503 std::string str()
const
505 std::ostringstream buffer;
506 buffer <<
" - mat4:\n" << mMatrix.str() << std::endl;
507 buffer <<
" - voxel dimensions: " << mVoxelSize << std::endl;
512 boost::shared_ptr<FullyDecomposedMap> createDecomposedMap()
530 affineMap->accumPreRotation(axis, radians);
536 affineMap->accumPreTranslation(t);
542 affineMap->accumPreScale(s);
548 affineMap->accumPreShear(axis0, axis1, shear);
560 affineMap->accumPostRotation(axis, radians);
566 affineMap->accumPostTranslation(t);
572 affineMap->accumPostScale(s);
578 affineMap->accumPostShear(axis0, axis1, shear);
586 { accumPreRotation(axis, radians); }
589 { accumPreScale(v); }
592 { accumPostTranslation(v); }
595 { accumPreShear(axis0, axis1, shear); }
599 {
return preRotate(radians, axis); }
602 {
return postTranslate(t); }
605 {
return preScale(s); }
608 {
return preShear(shear, axis0, axis1); }
615 void updateAcceleration() {
616 mDeterminant = mMatrix.getMat3().
det();
620 "Tried to initialize an affine transform from a nearly signular matrix");
622 mMatrixInv = mMatrix.inverse();
623 mJacobianInv = mMatrixInv.getMat3().transpose();
627 mVoxelSize(0) = (applyMap(
Vec3d(1,0,0)) - pos).length();
628 mVoxelSize(1) = (applyMap(
Vec3d(0,1,0)) - pos).length();
629 mVoxelSize(2) = (applyMap(
Vec3d(0,0,1)) - pos).length();
640 bool mIsDiagonal, mIsIdentity;
652 typedef boost::shared_ptr<ScaleMap>
Ptr;
653 typedef boost::shared_ptr<const ScaleMap>
ConstPtr;
656 mInvScaleSqr(1,1,1), mInvTwiceScale(0.5,0.5,0.5){}
661 mVoxelSize(
Vec3d(std::abs(scale(0)),std::abs(scale(1)), std::abs(scale(2))))
663 if (scale.
eq(
Vec3d(0.0), 1.0e-10)) {
666 mScaleValuesInverse = 1.0 / mScaleValues;
667 mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
668 mInvTwiceScale = mScaleValuesInverse / 2;
673 mScaleValues(other.mScaleValues),
674 mVoxelSize(other.mVoxelSize),
675 mScaleValuesInverse(other.mScaleValuesInverse),
676 mInvScaleSqr(other.mInvScaleSqr),
677 mInvTwiceScale(other.mInvTwiceScale)
688 static bool isRegistered() {
return MapRegistry::isRegistered(ScaleMap::mapType()); }
690 static void registerMap()
692 MapRegistry::registerMap(
706 value =
isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.y()));
707 value = value &&
isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.z()));
714 in.
x() * mScaleValues.x(),
715 in.
y() * mScaleValues.y(),
716 in.
z() * mScaleValues.z());
722 in.
x() * mScaleValuesInverse.x(),
723 in.
y() * mScaleValuesInverse.y(),
724 in.
z() * mScaleValuesInverse.z());
734 for (
int i=0; i<3; i++){
735 tmp.
setRow(i, in.
row(i)*mScaleValuesInverse(i));
737 for (
int i=0; i<3; i++){
738 tmp.
setCol(i, tmp.
col(i)*mScaleValuesInverse(i));
748 double determinant()
const {
return mScaleValues.x()*mScaleValues.y()*mScaleValues.z(); }
776 void read(std::istream& is)
778 mScaleValues.read(is);
780 mScaleValuesInverse.read(is);
781 mInvScaleSqr.read(is);
782 mInvTwiceScale.read(is);
785 void write(std::ostream& os)
const
787 mScaleValues.write(os);
788 mVoxelSize.write(os);
789 mScaleValuesInverse.write(os);
790 mInvScaleSqr.write(os);
791 mInvTwiceScale.write(os);
794 std::string str()
const
796 std::ostringstream buffer;
797 buffer <<
" - scale: " << mScaleValues << std::endl;
798 buffer <<
" - voxel dimensions: " << mVoxelSize << std::endl;
802 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
807 if (!mScaleValues.eq(other.mScaleValues)) {
return false; }
827 affineMap->accumPreRotation(axis, radians);
838 affineMap->accumPreShear(axis0, axis1, shear);
850 affineMap->accumPostRotation(axis, radians);
861 affineMap->accumPostShear(axis0, axis1, shear);
868 {
return preRotate(radians, axis); }
875 {
return preShear(shear, axis0, axis1); }
878 Vec3d mScaleValues, mVoxelSize, mScaleValuesInverse, mInvScaleSqr, mInvTwiceScale;
887 typedef boost::shared_ptr<UniformScaleMap>
Ptr;
888 typedef boost::shared_ptr<const UniformScaleMap>
ConstPtr;
900 static bool isRegistered() {
return MapRegistry::isRegistered(UniformScaleMap::mapType()); }
901 static void registerMap()
903 MapRegistry::registerMap(
904 UniformScaleMap::mapType(),
905 UniformScaleMap::create);
911 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
933 ScaleMap::preScale(
const Vec3d& v)
const
935 Vec3d new_scale(v * mScaleValues);
945 ScaleMap::postScale(
const Vec3d& v)
const
953 {
return preScale(v); }
959 typedef boost::shared_ptr<TranslationMap>
Ptr;
960 typedef boost::shared_ptr<const TranslationMap>
ConstPtr;
980 static bool isRegistered() {
return MapRegistry::isRegistered(TranslationMap::mapType()); }
982 static void registerMap()
984 MapRegistry::registerMap(
985 TranslationMap::mapType(),
986 TranslationMap::create);
1012 return applyIJC(mat);
1032 void read(std::istream& is) { mTranslation.read(is); }
1034 void write(std::ostream& os)
const { mTranslation.write(os); }
1037 std::string str()
const
1039 std::ostringstream buffer;
1040 buffer <<
" - translation: " << mTranslation << std::endl;
1041 return buffer.str();
1044 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1049 if (!mTranslation.eq(other.mTranslation)) {
return false; }
1072 affineMap->accumPreRotation(axis, radians);
1086 affineMap->accumPreShear(axis0, axis1, shear);
1097 affineMap->accumPostRotation(axis, radians);
1111 affineMap->accumPostShear(axis0, axis1, shear);
1120 {
return preRotate(radians, axis); }
1123 {
return preTranslate(t); }
1128 {
return preShear(shear, axis0, axis1); }
1144 typedef boost::shared_ptr<ScaleTranslateMap>
Ptr;
1145 typedef boost::shared_ptr<const ScaleTranslateMap>
ConstPtr;
1149 mTranslation(
Vec3d(0,0,0)),
1150 mScaleValues(
Vec3d(1,1,1)),
1151 mVoxelSize(
Vec3d(1,1,1)),
1152 mScaleValuesInverse(
Vec3d(1,1,1)),
1153 mInvScaleSqr(1,1,1),
1154 mInvTwiceScale(0.5,0.5,0.5)
1160 mTranslation(translate),
1161 mScaleValues(scale),
1162 mVoxelSize(std::abs(scale(0)), std::abs(scale(1)), std::abs(scale(2)))
1164 if (scale.
eq(
Vec3d(0.0), 1.0e-10)) {
1167 mScaleValuesInverse = 1.0 / mScaleValues;
1168 mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1169 mInvTwiceScale = mScaleValuesInverse / 2;
1174 mTranslation(translate.getTranslation()),
1176 mVoxelSize(std::abs(mScaleValues(0)),
1177 std::abs(mScaleValues(1)),
1178 std::abs(mScaleValues(2))),
1179 mScaleValuesInverse(1.0 / scale.
getScale())
1181 mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1182 mInvTwiceScale = mScaleValuesInverse / 2;
1187 mTranslation(other.mTranslation),
1188 mScaleValues(other.mScaleValues),
1189 mVoxelSize(other.mVoxelSize),
1190 mScaleValuesInverse(other.mScaleValuesInverse),
1191 mInvScaleSqr(other.mInvScaleSqr),
1192 mInvTwiceScale(other.mInvTwiceScale)
1202 static bool isRegistered() {
return MapRegistry::isRegistered(ScaleTranslateMap::mapType()); }
1204 static void registerMap()
1206 MapRegistry::registerMap(
1207 ScaleTranslateMap::mapType(),
1208 ScaleTranslateMap::create);
1221 value =
isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.y()));
1222 value = value &&
isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.z()));
1230 in.
x() * mScaleValues.x() + mTranslation.x(),
1231 in.
y() * mScaleValues.y() + mTranslation.y(),
1232 in.
z() * mScaleValues.z() + mTranslation.z());
1238 (in.
x() - mTranslation.x() ) / mScaleValues.x(),
1239 (in.
y() - mTranslation.y() ) / mScaleValues.y(),
1240 (in.
z() - mTranslation.z() ) / mScaleValues.z());
1249 in.
x() / mScaleValues.x(),
1250 in.
y() / mScaleValues.y(),
1251 in.
z() / mScaleValues.z());
1256 for (
int i=0; i<3; i++){
1257 tmp.
setRow(i, in.
row(i)*mScaleValuesInverse(i));
1259 for (
int i=0; i<3; i++){
1260 tmp.
setCol(i, tmp.
col(i)*mScaleValuesInverse(i));
1265 return applyIJC(in);
1271 double determinant()
const {
return mScaleValues.x()*mScaleValues.y()*mScaleValues.z(); }
1295 void read(std::istream& is)
1297 mTranslation.
read(is);
1298 mScaleValues.read(is);
1299 mVoxelSize.read(is);
1300 mScaleValuesInverse.read(is);
1301 mInvScaleSqr.read(is);
1302 mInvTwiceScale.read(is);
1305 void write(std::ostream& os)
const
1307 mTranslation.write(os);
1308 mScaleValues.write(os);
1309 mVoxelSize.write(os);
1310 mScaleValuesInverse.write(os);
1311 mInvScaleSqr.write(os);
1312 mInvTwiceScale.write(os);
1315 std::string str()
const
1317 std::ostringstream buffer;
1318 buffer <<
" - translation: " << mTranslation << std::endl;
1319 buffer <<
" - scale: " << mScaleValues << std::endl;
1320 buffer <<
" - voxel dimensions: " << mVoxelSize << std::endl;
1321 return buffer.str();
1324 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1329 if (!mScaleValues.eq(other.mScaleValues)) {
return false; }
1330 if (!mTranslation.eq(other.mTranslation)) {
return false; }
1340 affineMap->accumPostTranslation(mTranslation);
1350 affineMap->accumPreRotation(axis, radians);
1355 const Vec3d& s = mScaleValues;
1356 const Vec3d scaled_trans( t.
x() * s.
x(),
1367 affineMap->accumPreShear(axis0, axis1, shear);
1378 affineMap->accumPostRotation(axis, radians);
1391 affineMap->accumPostShear(axis0, axis1, shear);
1398 {
return preRotate(radians, axis); }
1401 {
return postTranslate(t); }
1406 {
return preShear(shear, axis0, axis1); }
1409 Vec3d mTranslation, mScaleValues, mVoxelSize, mScaleValuesInverse,
1410 mInvScaleSqr, mInvTwiceScale;
1415 ScaleMap::postTranslate(
const Vec3d& t)
const
1422 ScaleMap::preTranslate(
const Vec3d& t)
const
1425 const Vec3d& s = mScaleValues;
1426 const Vec3d scaled_trans( t.
x() * s.
x(),
1433 ScaleMap::translate(
const Vec3d& t)
const
1434 {
return postTranslate(t); }
1442 typedef boost::shared_ptr<UniformScaleTranslateMap>
Ptr;
1443 typedef boost::shared_ptr<const UniformScaleTranslateMap>
ConstPtr;
1459 static bool isRegistered()
1461 return MapRegistry::isRegistered(UniformScaleTranslateMap::mapType());
1464 static void registerMap()
1466 MapRegistry::registerMap(
1467 UniformScaleTranslateMap::mapType(),
1468 UniformScaleTranslateMap::create);
1474 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1488 const Vec3d new_trans = this->getTranslation() + scale * t;
1501 {
return postTranslate(t); }
1509 UniformScaleMap::postTranslate(
const Vec3d& t)
const
1517 UniformScaleMap::preTranslate(
const Vec3d& t)
const
1524 TranslationMap::preScale(
const Vec3d& v)
const
1534 TranslationMap::postScale(
const Vec3d& v)
const
1539 const Vec3d trans(mTranslation.x()*v.
x(),
1540 mTranslation.y()*v.
y(),
1541 mTranslation.z()*v.
z());
1547 ScaleTranslateMap::preScale(
const Vec3d& v)
const
1549 Vec3d new_scale( v * mScaleValues );
1558 ScaleTranslateMap::postScale(
const Vec3d& v)
const
1560 const Vec3d new_scale( v * mScaleValues );
1561 const Vec3d new_trans( mTranslation.x()*v.
x(),
1562 mTranslation.y()*v.
y(),
1563 mTranslation.z()*v.
z() );
1574 UniformScaleMap::translate(
const Vec3d& t)
const
1575 {
return postTranslate(t); }
1579 {
return preScale(v); }
1583 {
return preScale(v); }
1593 typedef boost::shared_ptr<UnitaryMap>
Ptr;
1631 "4x4 Matrix initializing unitary map was not unitary: not invertible");
1636 "4x4 Matrix initializing unitary map was not unitary: not affine");
1641 "4x4 Matrix initializing unitary map was not unitary: had translation");
1646 "4x4 Matrix initializing unitary map was not unitary");
1654 mAffineMap(other.mAffineMap)
1659 mAffineMap(*(first.getAffineMap()), *(second.getAffineMap()))
1669 static bool isRegistered() {
return MapRegistry::isRegistered(UnitaryMap::mapType()); }
1671 static void registerMap()
1673 MapRegistry::registerMap(
1674 UnitaryMap::mapType(),
1675 UnitaryMap::create);
1689 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1694 if (mAffineMap!=other.mAffineMap)
return false;
1710 return mAffineMap.applyIJC(in);
1713 return applyIJC(in);
1733 void read(std::istream& is)
1735 mAffineMap.read(is);
1739 void write(std::ostream& os)
const
1741 mAffineMap.write(os);
1744 std::string str()
const
1746 std::ostringstream buffer;
1747 buffer << mAffineMap.str();
1748 return buffer.str();
1755 {
return postRotate(radians, axis); }
1758 {
return postTranslate(t); }
1761 {
return preScale(v); }
1764 {
return preShear(shear, axis0, axis1); }
1778 affineMap->accumPreTranslation(t);
1784 affineMap->accumPreScale(v);
1790 affineMap->accumPreShear(axis0, axis1, shear);
1808 affineMap->accumPostTranslation(t);
1814 affineMap->accumPostScale(v);
1820 affineMap->accumPostShear(axis0, axis1, shear);
1842 typedef boost::shared_ptr<NonlinearFrustumMap>
Ptr;
1843 typedef boost::shared_ptr<const NonlinearFrustumMap>
ConstPtr;
1858 MapBase(),mBBox(bb.
min().asVec3d(), bb.
max().asVec3d()), mTaper(taper), mDepth(depth)
1864 MapBase(),mBBox(bb), mTaper(taper), mDepth(depth)
1877 mBBox(bb.
min().asVec3d(),bb.
max().asVec3d()), mTaper(taper), mDepth(depth)
1879 if (!secondMap->isLinear() ) {
1881 "The second map in the Frustum transfrom must be linear");
1883 mSecondMap = *( secondMap->getAffineMap() );
1889 mBBox(bb), mTaper(taper), mDepth(depth)
1891 if (!secondMap->isLinear() ) {
1893 "The second map in the Frustum transfrom must be linear");
1895 mSecondMap = *( secondMap->getAffineMap() );
1902 mTaper(other.mTaper),
1903 mDepth(other.mDepth),
1904 mSecondMap(other.mSecondMap),
1905 mHasSimpleAffine(other.mHasSimpleAffine)
1927 const Vec3d& direction,
1930 double z_near,
double depth,
1938 "The frustum depth must be non-zero and positive");
1940 if (!(up.
length() > 0)) {
1942 "The frustum height must be non-zero and positive");
1944 if (!(aspect > 0)) {
1946 "The frustum aspect ratio must be non-zero and positive");
1950 "The frustum up orientation must be perpendicular to into-frustum direction");
1953 double near_plane_height = 2 * up.
length();
1954 double near_plane_width = aspect * near_plane_height;
1959 mDepth = depth / near_plane_width;
1960 double gamma = near_plane_width / z_near;
1961 mTaper = 1./(mDepth*gamma + 1.);
1963 Vec3d direction_unit = direction;
1972 Vec3d(near_plane_width, near_plane_width, near_plane_width));
1976 Mat4d mat = scale * r2 * r1;
1990 static bool isRegistered() {
return MapRegistry::isRegistered(NonlinearFrustumMap::mapType()); }
1992 static void registerMap()
1994 MapRegistry::registerMap(
1995 NonlinearFrustumMap::mapType(),
1996 NonlinearFrustumMap::create);
2009 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
2013 if (mBBox!=other.mBBox)
return false;
2020 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
2023 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
2026 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
2029 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
2038 return mSecondMap.applyMap(applyFrustumMap(in));
2044 return applyFrustumInverseMap(mSecondMap.applyInverseMap(in));
2061 const Vec3d loc = applyFrustumMap(ijk);
2062 const double s = mGamma * loc.
z() + 1.;
2067 "Tried to evaluate the frustum transform at the singular focal point (e.g. camera)");
2070 const double sinv = 1.0/s;
2071 const double pt0 = mLx * sinv;
2072 const double pt1 = mGamma * pt0;
2073 const double pt2 = pt1 * sinv;
2075 const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
2078 Mat3d gradE(Mat3d::zero());
2079 for (
int j = 0; j < 3; ++j ) {
2080 gradE(0,j) = pt0 * jacinv(0,j) - pt2 * loc.
x()*jacinv(2,j);
2081 gradE(1,j) = pt0 * jacinv(1,j) - pt2 * loc.
y()*jacinv(2,j);
2082 gradE(2,j) = (1./mDepthOnLz) * jacinv(2,j);
2086 for (
int i = 0; i < 3; ++i) {
2087 result(0) = d1_is(0) * gradE(0,i) + d1_is(1) * gradE(1,i) + d1_is(2) * gradE(2,i);
2102 const Vec3d loc = applyFrustumMap(ijk);
2104 const double s = mGamma * loc.
z() + 1.;
2109 "Tried to evaluate the frustum transform at the singular focal point (e.g. camera)");
2113 const double sinv = 1.0/s;
2114 const double pt0 = mLx * sinv;
2115 const double pt1 = mGamma * pt0;
2116 const double pt2 = pt1 * sinv;
2117 const double pt3 = pt2 * sinv;
2120 const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
2124 Mat3d matE0(Mat3d::zero());
2125 Mat3d matE1(Mat3d::zero());
2126 for(
int j = 0; j < 3; j++) {
2127 for (
int k = 0; k < 3; k++) {
2129 const double pt4 = 2. * jacinv(2,j) * jacinv(2,k) * pt3;
2131 matE0(j,k) = -(jacinv(0,j) * jacinv(2,k) + jacinv(2,j) * jacinv(0,k)) * pt2 +
2134 matE1(j,k) = -(jacinv(1,j) * jacinv(2,k) + jacinv(2,j) * jacinv(1,k)) * pt2 +
2140 Mat3d gradE(Mat3d::zero());
2141 for (
int j = 0; j < 3; ++j ) {
2142 gradE(0,j) = pt0 * jacinv(0,j) - pt2 * loc.
x()*jacinv(2,j);
2143 gradE(1,j) = pt0 * jacinv(1,j) - pt2 * loc.
y()*jacinv(2,j);
2144 gradE(2,j) = (1./mDepthOnLz) * jacinv(2,j);
2147 Mat3d result(Mat3d::zero());
2150 for (
int m = 0; m < 3; ++m ) {
2151 for (
int n = 0; n < 3; ++n) {
2152 for (
int i = 0; i < 3; ++i ) {
2153 for (
int j = 0; j < 3; ++j) {
2154 result(m, n) += gradE(j, m) * gradE(i, n) * d2_is(i, j);
2160 for (
int m = 0; m < 3; ++m ) {
2161 for (
int n = 0; n < 3; ++n) {
2162 result(m, n) += matE0(m, n) * d1_is(0) + matE1(m, n) * d1_is(1);
2174 double determinant(
const Vec3d& loc)
const
2176 double s = mGamma * loc.
z() + 1.0;
2177 double frustum_determinant = s * s * mDepthOnLzLxLx;
2178 return mSecondMap.determinant() * frustum_determinant;
2183 const Vec3d loc( 0.5*(mBBox.min().x() + mBBox.max().x()),
2184 0.5*(mBBox.min().y() + mBBox.max().y()),
2187 return voxelSize(loc);
2196 Vec3d out, pos = applyMap(loc);
2197 out(0) = (applyMap(loc +
Vec3d(1,0,0)) - pos).length();
2198 out(1) = (applyMap(loc +
Vec3d(0,1,0)) - pos).length();
2199 out(2) = (applyMap(loc +
Vec3d(0,0,1)) - pos).length();
2210 return mSecondMap.getAffineMap();
2241 void read(std::istream& is)
2252 is.read(reinterpret_cast<char*>(&mTaper),
sizeof(
double));
2253 is.read(reinterpret_cast<char*>(&mDepth),
sizeof(
double));
2259 if(!MapRegistry::isRegistered(type)) {
2264 MapBase::Ptr proxy = math::MapRegistry::createMap(type);
2266 mSecondMap = *(proxy->getAffineMap());
2271 void write(std::ostream& os)
const
2274 os.write(reinterpret_cast<const char*>(&mTaper),
sizeof(
double));
2275 os.write(reinterpret_cast<const char*>(&mDepth),
sizeof(
double));
2278 mSecondMap.write(os);
2282 std::string str()
const
2284 std::ostringstream buffer;
2285 buffer <<
" - taper: " << mTaper << std::endl;
2286 buffer <<
" - depth: " << mDepth << std::endl;
2287 buffer <<
" SecondMap: "<< mSecondMap.type() << std::endl;
2288 buffer << mSecondMap.str() << std::endl;
2289 return buffer.str();
2313 mBBox, mTaper, mDepth, mSecondMap.preShear(shear, axis0, axis1)));
2338 mBBox, mTaper, mDepth, mSecondMap.postShear(shear, axis0, axis1)));
2344 {
return preRotate(radians, axis); }
2347 {
return postTranslate(t); }
2350 {
return preScale(s); }
2353 {
return preShear(shear, axis0, axis1); }
2359 mLx = mBBox.extents().x();
2360 mLy = mBBox.extents().y();
2361 mLz = mBBox.extents().z();
2365 "The index space bounding box must have at least two index points in each direction.");
2372 mGamma = (1./mTaper - 1) / mDepth;
2374 mDepthOnLz = mDepth/mLz;
2375 mDepthOnLzLxLx = mDepthOnLz/(mLx * mLx);
2378 mHasSimpleAffine =
true;
2379 Vec3d tmp = mSecondMap.voxelSize();
2382 if (!
isApproxEqual(tmp(0), tmp(1))) { mHasSimpleAffine =
false;
return; }
2383 if (!
isApproxEqual(tmp(0), tmp(2))) { mHasSimpleAffine =
false;
return; }
2385 Vec3d trans = mSecondMap.applyMap(
Vec3d(0,0,0));
2387 Vec3d tmp1 = mSecondMap.applyMap(
Vec3d(1,0,0)) - trans;
2388 Vec3d tmp2 = mSecondMap.applyMap(
Vec3d(0,1,0)) - trans;
2389 Vec3d tmp3 = mSecondMap.applyMap(
Vec3d(0,0,1)) - trans;
2392 if (!
isApproxEqual(tmp1.dot(tmp2), 0., 1.e-7)) { mHasSimpleAffine =
false;
return; }
2393 if (!
isApproxEqual(tmp2.dot(tmp3), 0., 1.e-7)) { mHasSimpleAffine =
false;
return; }
2394 if (!
isApproxEqual(tmp3.dot(tmp1), 0., 1.e-7)) { mHasSimpleAffine =
false;
return; }
2403 out = out - mBBox.min();
2408 out.z() *= mDepthOnLz;
2410 double scale = (mGamma * out.z() + 1.)/ mLx;
2419 Vec3d applyFrustumInverseMap(
const Vec3d& in)
const
2423 double invScale = mLx / (mGamma * out.z() + 1.);
2424 out.x() *= invScale;
2425 out.y() *= invScale;
2430 out.z() /= mDepthOnLz;
2433 out = out + mBBox.min();
2445 AffineMap mSecondMap;
2448 double mLx, mLy, mLz;
2449 double mXo, mYo, mGamma, mDepthOnLz, mDepthOnLzLxLx;
2452 bool mHasSimpleAffine;
2462 template<
typename FirstMapType,
typename SecondMapType>
2468 typedef boost::shared_ptr<MyType>
Ptr;
2474 CompoundMap(
const FirstMapType& f,
const SecondMapType& s): mFirstMap(f), mSecondMap(s)
2476 updateAffineMatrix();
2480 mFirstMap(other.mFirstMap),
2481 mSecondMap(other.mSecondMap),
2482 mAffineMap(other.mAffineMap)
2488 return (FirstMapType::mapType() +
Name(
":") + SecondMapType::mapType());
2493 if (mFirstMap != other.mFirstMap)
return false;
2494 if (mSecondMap != other.mSecondMap)
return false;
2495 if (mAffineMap != other.mAffineMap)
return false;
2503 mFirstMap = other.mFirstMap;
2504 mSecondMap = other.mSecondMap;
2505 mAffineMap = other.mAffineMap;
2511 return mAffineMap.isIdentity();
2513 return mFirstMap.isIdentity()&&mSecondMap.isIdentity();
2519 return mAffineMap.isDiagonal();
2521 return mFirstMap.isDiagonal()&&mSecondMap.isDiagonal();
2531 "Constant affine matrix representation not possible for this nonlinear map");
2536 const FirstMapType&
firstMap()
const {
return mFirstMap; }
2537 const SecondMapType&
secondMap()
const {
return mSecondMap; }
2539 void setFirstMap(
const FirstMapType& first) { mFirstMap = first; updateAffineMatrix(); }
2540 void setSecondMap(
const SecondMapType& second) { mSecondMap = second; updateAffineMatrix(); }
2542 void read(std::istream& is)
2544 mAffineMap.read(is);
2546 mSecondMap.read(is);
2548 void write(std::ostream& os)
const
2550 mAffineMap.write(os);
2551 mFirstMap.write(os);
2552 mSecondMap.write(os);
2556 void updateAffineMatrix()
2562 mAffineMap =
AffineMap(*first, *second);
2566 FirstMapType mFirstMap;
2567 SecondMapType mSecondMap;
2569 AffineMap mAffineMap;
2577 #endif // OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED