33 #ifndef OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED
34 #define OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED
41 #include <openvdb/util/Name.h>
42 #include <openvdb/Types.h>
43 #include <boost/shared_ptr.hpp>
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; };
161 typedef boost::shared_ptr<MapBase>
Ptr;
163 typedef Ptr (*MapFactory)();
167 virtual boost::shared_ptr<AffineMap> getAffineMap()
const = 0;
170 virtual Name type()
const = 0;
173 template<
typename MapT>
bool isType()
const {
return this->type() == MapT::mapType(); }
176 virtual bool isEqual(
const MapBase& other)
const = 0;
179 virtual bool isLinear()
const = 0;
181 virtual bool hasUniformScale()
const = 0;
183 virtual Vec3d applyMap(
const Vec3d& in)
const = 0;
184 virtual Vec3d applyInverseMap(
const Vec3d& in)
const = 0;
186 virtual Vec3d applyIJT(
const Vec3d& in)
const = 0;
188 virtual Mat3d applyIJC(
const Mat3d& m)
const = 0;
192 virtual double determinant()
const = 0;
193 virtual double determinant(
const Vec3d&)
const = 0;
197 virtual Vec3d voxelSize()
const = 0;
199 virtual Vec3d voxelSize(
const Vec3d&)
const = 0;
202 virtual void read(std::istream&) = 0;
203 virtual void write(std::ostream&)
const = 0;
205 virtual std::string str()
const = 0;
225 template<
typename MapT>
228 return other.
isType<MapT>() && (
self == *static_cast<const MapT*>(&other));
249 static bool isRegistered(
const Name&);
252 static void registerMap(
const Name&, MapBase::MapFactory);
255 static void unregisterMap(
const Name&);
279 typedef boost::shared_ptr<AffineMap>
Ptr;
280 typedef boost::shared_ptr<const AffineMap>
ConstPtr;
283 mMatrix(
Mat4d::identity()),
284 mMatrixInv(
Mat4d::identity()),
285 mJacobianInv(
Mat3d::identity()),
287 mVoxelSize(
Vec3d(1,1,1)),
299 updateAcceleration();
306 "Tried to initialize an affine transform from a non-affine 4x4 matrix");
308 updateAcceleration();
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)
325 mMatrix(first.mMatrix * second.mMatrix)
327 updateAcceleration();
337 static bool isRegistered() {
return MapRegistry::isRegistered(AffineMap::mapType()); }
341 MapRegistry::registerMap(
342 AffineMap::mapType(),
355 Mat3d mat = mMatrix.getMat3();
356 const double det = mat.
det();
360 mat *= (1.f / pow(std::abs(det),1./3.));
365 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
370 if (!mMatrix.eq(other.mMatrix)) {
return false; }
371 if (!mMatrixInv.eq(other.mMatrixInv)) {
return false; }
379 mMatrix = other.mMatrix;
380 mMatrixInv = other.mMatrixInv;
382 mJacobianInv = other.mJacobianInv;
383 mDeterminant = other.mDeterminant;
384 mVoxelSize = other.mVoxelSize;
385 mIsDiagonal = other.mIsDiagonal;
386 mIsIdentity = other.mIsIdentity;
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] );
410 return mJacobianInv.
transpose()* m * mJacobianInv;
421 Vec3d voxelSize()
const {
return mVoxelSize; }
440 void accumPreRotation(
Axis axis,
double radians)
443 mMatrix.preRotate(axis, radians);
444 updateAcceleration();
449 updateAcceleration();
453 mMatrix.preTranslate(v);
454 updateAcceleration();
458 mMatrix.preShear(axis0, axis1, shear);
459 updateAcceleration();
465 void accumPostRotation(
Axis axis,
double radians)
468 mMatrix.postRotate(axis, radians);
469 updateAcceleration();
473 mMatrix.postScale(v);
474 updateAcceleration();
478 mMatrix.postTranslate(v);
479 updateAcceleration();
483 mMatrix.postShear(axis0, axis1, shear);
484 updateAcceleration();
493 updateAcceleration();
505 std::ostringstream buffer;
506 buffer <<
" - mat4:\n" << mMatrix.str() << std::endl;
507 buffer <<
" - voxel dimensions: " << mVoxelSize << std::endl;
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);
589 void updateAcceleration() {
590 mDeterminant = mMatrix.getMat3().
det();
594 "Tried to initialize an affine transform from a nearly singular matrix");
596 mMatrixInv = mMatrix.inverse();
597 mJacobianInv = mMatrixInv.getMat3().transpose();
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();
614 bool mIsDiagonal, mIsIdentity;
626 typedef boost::shared_ptr<ScaleMap>
Ptr;
627 typedef boost::shared_ptr<const ScaleMap>
ConstPtr;
630 mInvScaleSqr(1,1,1), mInvTwiceScale(0.5,0.5,0.5){}
635 mVoxelSize(
Vec3d(std::abs(scale(0)),std::abs(scale(1)), std::abs(scale(2))))
637 double determinant = scale[0]* scale[1] * scale[2];
641 mScaleValuesInverse = 1.0 / mScaleValues;
642 mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
643 mInvTwiceScale = mScaleValuesInverse / 2;
648 mScaleValues(other.mScaleValues),
649 mVoxelSize(other.mVoxelSize),
650 mScaleValuesInverse(other.mScaleValuesInverse),
651 mInvScaleSqr(other.mInvScaleSqr),
652 mInvTwiceScale(other.mInvTwiceScale)
663 static bool isRegistered() {
return MapRegistry::isRegistered(ScaleMap::mapType()); }
667 MapRegistry::registerMap(
682 std::abs(mScaleValues.x()), std::abs(mScaleValues.y()),
double(5e-7));
684 std::abs(mScaleValues.x()), std::abs(mScaleValues.z()),
double(5e-7));
692 in.
x() * mScaleValues.x(),
693 in.
y() * mScaleValues.y(),
694 in.
z() * mScaleValues.z());
700 in.
x() * mScaleValuesInverse.x(),
701 in.
y() * mScaleValuesInverse.y(),
702 in.
z() * mScaleValuesInverse.z());
719 for (
int i = 0; i < 3; i++) {
720 tmp.
setRow(i, in.
row(i) * mScaleValuesInverse(i));
722 for (
int i = 0; i < 3; i++) {
723 tmp.
setCol(i, tmp.
col(i) * mScaleValuesInverse(i));
731 double determinant()
const {
return mScaleValues.x() * mScaleValues.y() * mScaleValues.z(); }
744 Vec3d voxelSize()
const {
return mVoxelSize; }
755 mScaleValues.read(is);
757 mScaleValuesInverse.read(is);
758 mInvScaleSqr.read(is);
759 mInvTwiceScale.read(is);
764 mScaleValues.write(os);
765 mVoxelSize.write(os);
766 mScaleValuesInverse.write(os);
767 mInvScaleSqr.write(os);
768 mInvTwiceScale.write(os);
773 std::ostringstream buffer;
774 buffer <<
" - scale: " << mScaleValues << std::endl;
775 buffer <<
" - voxel dimensions: " << mVoxelSize << std::endl;
779 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
784 if (!mScaleValues.eq(other.mScaleValues)) {
return false; }
804 affineMap->accumPreRotation(axis, radians);
808 MapBase::Ptr preTranslate(
const Vec3d& tr)
const;
810 MapBase::Ptr preScale(
const Vec3d& v)
const;
815 affineMap->accumPreShear(axis0, axis1, shear);
827 affineMap->accumPostRotation(axis, radians);
831 MapBase::Ptr postTranslate(
const Vec3d& tr)
const;
833 MapBase::Ptr postScale(
const Vec3d& v)
const;
838 affineMap->accumPostShear(axis0, axis1, shear);
844 Vec3d mScaleValues, mVoxelSize, mScaleValuesInverse, mInvScaleSqr, mInvTwiceScale;
853 typedef boost::shared_ptr<UniformScaleMap>
Ptr;
854 typedef boost::shared_ptr<const UniformScaleMap>
ConstPtr;
866 static bool isRegistered() {
return MapRegistry::isRegistered(UniformScaleMap::mapType()); }
869 MapRegistry::registerMap(
870 UniformScaleMap::mapType(),
871 UniformScaleMap::create);
877 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
897 ScaleMap::preScale(
const Vec3d& v)
const
899 const Vec3d new_scale(v * mScaleValues);
909 ScaleMap::postScale(
const Vec3d& v)
const
919 typedef boost::shared_ptr<TranslationMap>
Ptr;
920 typedef boost::shared_ptr<const TranslationMap>
ConstPtr;
934 static bool isRegistered() {
return MapRegistry::isRegistered(TranslationMap::mapType()); }
938 MapRegistry::registerMap(
939 TranslationMap::mapType(),
940 TranslationMap::create);
984 void read(std::istream& is) { mTranslation.read(is); }
986 void write(std::ostream& os)
const { mTranslation.write(os); }
991 std::ostringstream buffer;
992 buffer <<
" - translation: " << mTranslation << std::endl;
996 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1001 return mTranslation.eq(other.mTranslation);
1022 affineMap->accumPreRotation(axis, radians);
1036 affineMap->accumPreShear(axis0, axis1, shear);
1047 affineMap->accumPostRotation(axis, radians);
1061 affineMap->accumPostShear(axis0, axis1, shear);
1080 typedef boost::shared_ptr<ScaleTranslateMap>
Ptr;
1081 typedef boost::shared_ptr<const ScaleTranslateMap>
ConstPtr;
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)
1096 mTranslation(translate),
1097 mScaleValues(scale),
1098 mVoxelSize(std::abs(scale(0)), std::abs(scale(1)), std::abs(scale(2)))
1100 const double determinant = scale[0]* scale[1] * scale[2];
1104 mScaleValuesInverse = 1.0 / mScaleValues;
1105 mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1106 mInvTwiceScale = mScaleValuesInverse / 2;
1111 mTranslation(translate.getTranslation()),
1113 mVoxelSize(std::abs(mScaleValues(0)),
1114 std::abs(mScaleValues(1)),
1115 std::abs(mScaleValues(2))),
1116 mScaleValuesInverse(1.0 / scale.
getScale())
1118 mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1119 mInvTwiceScale = mScaleValuesInverse / 2;
1124 mTranslation(other.mTranslation),
1125 mScaleValues(other.mScaleValues),
1126 mVoxelSize(other.mVoxelSize),
1127 mScaleValuesInverse(other.mScaleValuesInverse),
1128 mInvScaleSqr(other.mInvScaleSqr),
1129 mInvTwiceScale(other.mInvTwiceScale)
1139 static bool isRegistered() {
return MapRegistry::isRegistered(ScaleTranslateMap::mapType()); }
1143 MapRegistry::registerMap(
1144 ScaleTranslateMap::mapType(),
1145 ScaleTranslateMap::create);
1159 std::abs(mScaleValues.x()), std::abs(mScaleValues.y()),
double(5e-7));
1161 std::abs(mScaleValues.x()), std::abs(mScaleValues.z()),
double(5e-7));
1169 in.
x() * mScaleValues.x() + mTranslation.x(),
1170 in.
y() * mScaleValues.y() + mTranslation.y(),
1171 in.
z() * mScaleValues.z() + mTranslation.z());
1177 (in.
x() - mTranslation.x() ) / mScaleValues.x(),
1178 (in.
y() - mTranslation.y() ) / mScaleValues.y(),
1179 (in.
z() - mTranslation.z() ) / mScaleValues.z());
1194 in.
x() / mScaleValues.x(),
1195 in.
y() / mScaleValues.y(),
1196 in.
z() / mScaleValues.z());
1202 for (
int i=0; i<3; i++){
1203 tmp.
setRow(i, in.
row(i)*mScaleValuesInverse(i));
1205 for (
int i=0; i<3; i++){
1206 tmp.
setCol(i, tmp.
col(i)*mScaleValuesInverse(i));
1215 double determinant()
const {
return mScaleValues.x()*mScaleValues.y()*mScaleValues.z(); }
1237 mTranslation.read(is);
1238 mScaleValues.read(is);
1239 mVoxelSize.read(is);
1240 mScaleValuesInverse.read(is);
1241 mInvScaleSqr.read(is);
1242 mInvTwiceScale.read(is);
1247 mTranslation.write(os);
1248 mScaleValues.write(os);
1249 mVoxelSize.write(os);
1250 mScaleValuesInverse.write(os);
1251 mInvScaleSqr.write(os);
1252 mInvTwiceScale.write(os);
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();
1264 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1269 if (!mScaleValues.eq(other.mScaleValues)) {
return false; }
1270 if (!mTranslation.eq(other.mTranslation)) {
return false; }
1280 affineMap->accumPostTranslation(mTranslation);
1290 affineMap->accumPreRotation(axis, radians);
1295 const Vec3d& s = mScaleValues;
1296 const Vec3d scaled_trans( t.
x() * s.
x(),
1307 affineMap->accumPreShear(axis0, axis1, shear);
1318 affineMap->accumPostRotation(axis, radians);
1331 affineMap->accumPostShear(axis0, axis1, shear);
1337 Vec3d mTranslation, mScaleValues, mVoxelSize, mScaleValuesInverse,
1338 mInvScaleSqr, mInvTwiceScale;
1343 ScaleMap::postTranslate(
const Vec3d& t)
const
1350 ScaleMap::preTranslate(
const Vec3d& t)
const
1353 const Vec3d& s = mScaleValues;
1354 const Vec3d scaled_trans( t.
x() * s.
x(),
1366 typedef boost::shared_ptr<UniformScaleTranslateMap>
Ptr;
1367 typedef boost::shared_ptr<const UniformScaleTranslateMap>
ConstPtr;
1385 return MapRegistry::isRegistered(UniformScaleTranslateMap::mapType());
1390 MapRegistry::registerMap(
1391 UniformScaleTranslateMap::mapType(),
1392 UniformScaleTranslateMap::create);
1398 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1411 const Vec3d new_trans = this->getTranslation() + scale * t;
1426 UniformScaleMap::postTranslate(
const Vec3d& t)
const
1434 UniformScaleMap::preTranslate(
const Vec3d& t)
const
1442 TranslationMap::preScale(
const Vec3d& v)
const
1453 TranslationMap::postScale(
const Vec3d& v)
const
1458 const Vec3d trans(mTranslation.x()*v.
x(),
1459 mTranslation.y()*v.
y(),
1460 mTranslation.z()*v.
z());
1467 ScaleTranslateMap::preScale(
const Vec3d& v)
const
1469 const Vec3d new_scale( v * mScaleValues );
1479 ScaleTranslateMap::postScale(
const Vec3d& v)
const
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() );
1502 typedef boost::shared_ptr<UnitaryMap>
Ptr;
1540 "4x4 Matrix initializing unitary map was not unitary: not invertible");
1545 "4x4 Matrix initializing unitary map was not unitary: not affine");
1550 "4x4 Matrix initializing unitary map was not unitary: had translation");
1555 "4x4 Matrix initializing unitary map was not unitary");
1563 mAffineMap(other.mAffineMap)
1568 mAffineMap(*(first.getAffineMap()), *(second.getAffineMap()))
1578 static bool isRegistered() {
return MapRegistry::isRegistered(UnitaryMap::mapType()); }
1582 MapRegistry::registerMap(
1583 UnitaryMap::mapType(),
1584 UnitaryMap::create);
1598 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1603 if (mAffineMap!=other.mAffineMap)
return false;
1641 mAffineMap.read(is);
1647 mAffineMap.write(os);
1652 std::ostringstream buffer;
1653 buffer << mAffineMap.str();
1654 return buffer.str();
1671 affineMap->accumPreTranslation(t);
1677 affineMap->accumPreScale(v);
1683 affineMap->accumPreShear(axis0, axis1, shear);
1701 affineMap->accumPostTranslation(t);
1707 affineMap->accumPostScale(v);
1713 affineMap->accumPostShear(axis0, axis1, shear);
1735 typedef boost::shared_ptr<NonlinearFrustumMap>
Ptr;
1736 typedef boost::shared_ptr<const NonlinearFrustumMap>
ConstPtr;
1751 MapBase(),mBBox(bb), mTaper(taper), mDepth(depth)
1763 mBBox(bb), mTaper(taper), mDepth(depth)
1765 if (!secondMap->isLinear() ) {
1767 "The second map in the Frustum transfrom must be linear");
1769 mSecondMap = *( secondMap->getAffineMap() );
1776 mTaper(other.mTaper),
1777 mDepth(other.mDepth),
1778 mSecondMap(other.mSecondMap),
1779 mHasSimpleAffine(other.mHasSimpleAffine)
1800 const Vec3d& direction,
1803 double z_near,
double depth,
1811 "The frustum depth must be non-zero and positive");
1813 if (!(up.
length() > 0)) {
1815 "The frustum height must be non-zero and positive");
1817 if (!(aspect > 0)) {
1819 "The frustum aspect ratio must be non-zero and positive");
1823 "The frustum up orientation must be perpendicular to into-frustum direction");
1826 double near_plane_height = 2 * up.
length();
1827 double near_plane_width = aspect * near_plane_height;
1832 mDepth = depth / near_plane_width;
1833 double gamma = near_plane_width / z_near;
1834 mTaper = 1./(mDepth*gamma + 1.);
1836 Vec3d direction_unit = direction;
1845 Vec3d(near_plane_width, near_plane_width, near_plane_width));
1849 Mat4d mat = scale * r2 * r1;
1863 static bool isRegistered() {
return MapRegistry::isRegistered(NonlinearFrustumMap::mapType()); }
1867 MapRegistry::registerMap(
1868 NonlinearFrustumMap::mapType(),
1869 NonlinearFrustumMap::create);
1890 const Vec3d e1(1,0,0);
1891 if (!applyMap(e1).eq(e1))
return false;
1893 const Vec3d e2(0,1,0);
1894 if (!applyMap(e2).eq(e2))
return false;
1896 const Vec3d e3(0,0,1);
1897 if (!applyMap(e3).eq(e3))
return false;
1902 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1906 if (mBBox!=other.mBBox)
return false;
1913 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
1916 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
1919 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
1922 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
1931 return mSecondMap.applyMap(applyFrustumMap(in));
1937 return applyFrustumInverseMap(mSecondMap.applyInverseMap(in));
1946 Vec3d centered(loc);
1947 centered = centered - mBBox.min();
1948 centered.
x() -= mXo;
1949 centered.
y() -= mYo;
1952 const double zprime = centered.
z()*mDepthOnLz;
1954 const double scale = (mGamma * zprime + 1.) / mLx;
1955 const double scale2 = mGamma * mDepthOnLz / mLx;
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());
1961 return mSecondMap.applyJacobian(tmp);
1977 const Vec3d loc = applyFrustumMap(ijk);
1978 const double s = mGamma * loc.
z() + 1.;
1983 " at the singular focal point (e.g. camera)");
1986 const double sinv = 1.0/s;
1987 const double pt0 = mLx * sinv;
1988 const double pt1 = mGamma * pt0;
1989 const double pt2 = pt1 * sinv;
1991 const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
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);
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);
2018 const Vec3d loc = applyFrustumMap(ijk);
2020 const double s = mGamma * loc.
z() + 1.;
2025 " at the singular focal point (e.g. camera)");
2029 const double sinv = 1.0/s;
2030 const double pt0 = mLx * sinv;
2031 const double pt1 = mGamma * pt0;
2032 const double pt2 = pt1 * sinv;
2033 const double pt3 = pt2 * sinv;
2035 const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
2039 Mat3d matE0(Mat3d::zero());
2040 Mat3d matE1(Mat3d::zero());
2041 for(
int j = 0; j < 3; j++) {
2042 for (
int k = 0; k < 3; k++) {
2044 const double pt4 = 2. * jacinv(2,j) * jacinv(2,k) * pt3;
2046 matE0(j,k) = -(jacinv(0,j) * jacinv(2,k) + jacinv(2,j) * jacinv(0,k)) * pt2 +
2049 matE1(j,k) = -(jacinv(1,j) * jacinv(2,k) + jacinv(2,j) * jacinv(1,k)) * pt2 +
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);
2062 Mat3d result(Mat3d::zero());
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);
2075 for (
int m = 0; m < 3; ++m ) {
2076 for (
int n = 0; n < 3; ++n) {
2078 matE0(m, n) * d1_is(0) + matE1(m, n) * d1_is(1);
2092 double s = mGamma * loc.
z() + 1.0;
2093 double frustum_determinant = s * s * mDepthOnLzLxLx;
2094 return mSecondMap.determinant() * frustum_determinant;
2100 const Vec3d loc( 0.5*(mBBox.min().x() + mBBox.max().x()),
2101 0.5*(mBBox.min().y() + mBBox.max().y()),
2104 return voxelSize(loc);
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();
2158 is.read(reinterpret_cast<char*>(&mTaper),
sizeof(
double));
2159 is.read(reinterpret_cast<char*>(&mDepth),
sizeof(
double));
2165 if(!MapRegistry::isRegistered(type)) {
2170 MapBase::Ptr proxy = math::MapRegistry::createMap(type);
2172 mSecondMap = *(proxy->getAffineMap());
2180 os.write(reinterpret_cast<const char*>(&mTaper),
sizeof(
double));
2181 os.write(reinterpret_cast<const char*>(&mDepth),
sizeof(
double));
2184 mSecondMap.write(os);
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();
2219 mBBox, mTaper, mDepth, mSecondMap.preShear(shear, axis0, axis1)));
2244 mBBox, mTaper, mDepth, mSecondMap.postShear(shear, axis0, axis1)));
2252 mLx = mBBox.extents().x();
2253 mLy = mBBox.extents().y();
2254 mLz = mBBox.extents().z();
2258 " must have at least two index points in each direction.");
2265 mGamma = (1./mTaper - 1) / mDepth;
2267 mDepthOnLz = mDepth/mLz;
2268 mDepthOnLzLxLx = mDepthOnLz/(mLx * mLx);
2271 mHasSimpleAffine =
true;
2272 Vec3d tmp = mSecondMap.voxelSize();
2275 if (!
isApproxEqual(tmp(0), tmp(1))) { mHasSimpleAffine =
false;
return; }
2276 if (!
isApproxEqual(tmp(0), tmp(2))) { mHasSimpleAffine =
false;
return; }
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;
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; }
2296 out = out - mBBox.min();
2301 out.z() *= mDepthOnLz;
2303 double scale = (mGamma * out.z() + 1.)/ mLx;
2312 Vec3d applyFrustumInverseMap(
const Vec3d& in)
const
2316 double invScale = mLx / (mGamma * out.z() + 1.);
2317 out.x() *= invScale;
2318 out.y() *= invScale;
2323 out.z() /= mDepthOnLz;
2326 out = out + mBBox.min();
2338 AffineMap mSecondMap;
2341 double mLx, mLy, mLz;
2342 double mXo, mYo, mGamma, mDepthOnLz, mDepthOnLzLxLx;
2345 bool mHasSimpleAffine;
2355 template<
typename FirstMapType,
typename SecondMapType>
2361 typedef boost::shared_ptr<MyType>
Ptr;
2367 CompoundMap(
const FirstMapType& f,
const SecondMapType& s): mFirstMap(f), mSecondMap(s)
2369 updateAffineMatrix();
2373 mFirstMap(other.mFirstMap),
2374 mSecondMap(other.mSecondMap),
2375 mAffineMap(other.mAffineMap)
2381 return (FirstMapType::mapType() +
Name(
":") + SecondMapType::mapType());
2386 if (mFirstMap != other.mFirstMap)
return false;
2387 if (mSecondMap != other.mSecondMap)
return false;
2388 if (mAffineMap != other.mAffineMap)
return false;
2396 mFirstMap = other.mFirstMap;
2397 mSecondMap = other.mSecondMap;
2398 mAffineMap = other.mAffineMap;
2405 return mAffineMap.isIdentity();
2407 return mFirstMap.isIdentity()&&mSecondMap.isIdentity();
2413 return mAffineMap.isDiagonal();
2415 return mFirstMap.isDiagonal()&&mSecondMap.isDiagonal();
2426 "Constant affine matrix representation not possible for this nonlinear map");
2431 const FirstMapType&
firstMap()
const {
return mFirstMap; }
2432 const SecondMapType&
secondMap()
const {
return mSecondMap; }
2434 void setFirstMap(
const FirstMapType& first) { mFirstMap = first; updateAffineMatrix(); }
2435 void setSecondMap(
const SecondMapType& second) { mSecondMap = second; updateAffineMatrix(); }
2439 mAffineMap.read(is);
2441 mSecondMap.read(is);
2445 mAffineMap.write(os);
2446 mFirstMap.write(os);
2447 mSecondMap.write(os);
2451 void updateAffineMatrix()
2457 mAffineMap =
AffineMap(*first, *second);
2461 FirstMapType mFirstMap;
2462 SecondMapType mSecondMap;
2464 AffineMap mAffineMap;
2471 #endif // OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED