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;
198 virtual void read(std::istream&) = 0;
199 virtual void write(std::ostream&)
const = 0;
201 virtual std::string str()
const = 0;
221 template<
typename MapT>
222 static bool isEqualBase(
const MapT&
self,
const MapBase& other)
224 return other.
isType<MapT>() && (
self == *static_cast<const MapT*>(&other));
239 typedef Mutex::scoped_lock
Lock;
247 static bool isRegistered(
const Name&);
250 static void registerMap(
const Name&, MapBase::MapFactory);
253 static void unregisterMap(
const Name&);
262 mutable Mutex mMutex;
275 typedef boost::shared_ptr<AffineMap>
Ptr;
276 typedef boost::shared_ptr<const AffineMap>
ConstPtr;
279 mMatrix(
Mat4d::identity()),
280 mMatrixInv(
Mat4d::identity()),
281 mJacobianInv(
Mat3d::identity()),
283 mVoxelSize(
Vec3d(1,1,1)),
295 updateAcceleration();
302 "Tried to initialize an affine transform from a non-affine 4x4 matrix");
304 updateAcceleration();
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)
321 mMatrix(first.mMatrix * second.mMatrix)
323 updateAcceleration();
333 static bool isRegistered() {
return MapRegistry::isRegistered(AffineMap::mapType()); }
335 static void registerMap()
337 MapRegistry::registerMap(
338 AffineMap::mapType(),
351 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
356 if (!mMatrix.eq(other.mMatrix)) {
return false; }
357 if (!mMatrixInv.eq(other.mMatrixInv)) {
return false; }
365 mMatrix = other.mMatrix;
366 mMatrixInv = other.mMatrixInv;
368 mJacobianInv = other.mJacobianInv;
369 mDeterminant = other.mDeterminant;
370 mVoxelSize = other.mVoxelSize;
371 mIsDiagonal = other.mIsDiagonal;
372 mIsIdentity = other.mIsIdentity;
386 return mJacobianInv.
transpose()* m * mJacobianInv;
418 void accumPreRotation(
Axis axis,
double radians)
420 mMatrix.preRotate(axis, radians);
421 updateAcceleration();
426 updateAcceleration();
428 void accumPreTranslation(
const Vec3d& v)
430 mMatrix.preTranslate(v);
431 updateAcceleration();
435 mMatrix.preShear(axis0, axis1, shear);
436 updateAcceleration();
444 void accumPostRotation(
Axis axis,
double radians)
446 mMatrix.postRotate(axis, radians);
447 updateAcceleration();
451 mMatrix.postScale(v);
452 updateAcceleration();
454 void accumPostTranslation(
const Vec3d& v)
456 mMatrix.postTranslate(v);
457 updateAcceleration();
461 mMatrix.postShear(axis0, axis1, shear);
462 updateAcceleration();
468 void read(std::istream& is)
471 updateAcceleration();
475 void write(std::ostream& os)
const
481 std::string str()
const
483 std::ostringstream buffer;
484 buffer <<
" - mat4:\n" << mMatrix.str() << std::endl;
485 buffer <<
" - voxel dimensions: " << mVoxelSize << std::endl;
490 boost::shared_ptr<FullyDecomposedMap> createDecomposedMap()
508 affineMap->accumPreRotation(axis, radians);
514 affineMap->accumPreTranslation(t);
520 affineMap->accumPreScale(s);
526 affineMap->accumPreShear(axis0, axis1, shear);
538 affineMap->accumPostRotation(axis, radians);
544 affineMap->accumPostTranslation(t);
550 affineMap->accumPostScale(s);
556 affineMap->accumPostShear(axis0, axis1, shear);
567 void updateAcceleration() {
568 mDeterminant = mMatrix.getMat3().
det();
572 "Tried to initialize an affine transform from a nearly signular matrix");
574 mMatrixInv = mMatrix.inverse();
575 mJacobianInv = mMatrixInv.getMat3().transpose();
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();
592 bool mIsDiagonal, mIsIdentity;
604 typedef boost::shared_ptr<ScaleMap>
Ptr;
605 typedef boost::shared_ptr<const ScaleMap>
ConstPtr;
608 mInvScaleSqr(1,1,1), mInvTwiceScale(0.5,0.5,0.5){}
613 mVoxelSize(
Vec3d(std::abs(scale(0)),std::abs(scale(1)), std::abs(scale(2))))
615 if (scale.
eq(
Vec3d(0.0), 1.0e-10)) {
618 mScaleValuesInverse = 1.0 / mScaleValues;
619 mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
620 mInvTwiceScale = mScaleValuesInverse / 2;
625 mScaleValues(other.mScaleValues),
626 mVoxelSize(other.mVoxelSize),
627 mScaleValuesInverse(other.mScaleValuesInverse),
628 mInvScaleSqr(other.mInvScaleSqr),
629 mInvTwiceScale(other.mInvTwiceScale)
640 static bool isRegistered() {
return MapRegistry::isRegistered(ScaleMap::mapType()); }
642 static void registerMap()
644 MapRegistry::registerMap(
658 value =
isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.y()));
659 value = value &&
isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.z()));
666 in.
x() * mScaleValues.x(),
667 in.
y() * mScaleValues.y(),
668 in.
z() * mScaleValues.z());
674 in.
x() * mScaleValuesInverse.x(),
675 in.
y() * mScaleValuesInverse.y(),
676 in.
z() * mScaleValuesInverse.z());
686 for (
int i=0; i<3; i++){
687 tmp.
setRow(i, in.
row(i)*mScaleValuesInverse(i));
689 for (
int i=0; i<3; i++){
690 tmp.
setCol(i, tmp.
col(i)*mScaleValuesInverse(i));
700 double determinant()
const {
return mScaleValues.x()*mScaleValues.y()*mScaleValues.z(); }
722 void read(std::istream& is)
724 mScaleValues.read(is);
726 mScaleValuesInverse.read(is);
727 mInvScaleSqr.read(is);
728 mInvTwiceScale.read(is);
731 void write(std::ostream& os)
const
733 mScaleValues.write(os);
734 mVoxelSize.write(os);
735 mScaleValuesInverse.write(os);
736 mInvScaleSqr.write(os);
737 mInvTwiceScale.write(os);
740 std::string str()
const
742 std::ostringstream buffer;
743 buffer <<
" - scale: " << mScaleValues << std::endl;
744 buffer <<
" - voxel dimensions: " << mVoxelSize << std::endl;
748 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
753 if (!mScaleValues.eq(other.mScaleValues)) {
return false; }
773 affineMap->accumPreRotation(axis, radians);
784 affineMap->accumPreShear(axis0, axis1, shear);
796 affineMap->accumPostRotation(axis, radians);
807 affineMap->accumPostShear(axis0, axis1, shear);
813 Vec3d mScaleValues, mVoxelSize, mScaleValuesInverse, mInvScaleSqr, mInvTwiceScale;
822 typedef boost::shared_ptr<UniformScaleMap>
Ptr;
823 typedef boost::shared_ptr<const UniformScaleMap>
ConstPtr;
835 static bool isRegistered() {
return MapRegistry::isRegistered(UniformScaleMap::mapType()); }
836 static void registerMap()
838 MapRegistry::registerMap(
839 UniformScaleMap::mapType(),
840 UniformScaleMap::create);
846 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
866 ScaleMap::preScale(
const Vec3d& v)
const
868 Vec3d new_scale(v * mScaleValues);
878 ScaleMap::postScale(
const Vec3d& v)
const
888 typedef boost::shared_ptr<TranslationMap>
Ptr;
889 typedef boost::shared_ptr<const TranslationMap>
ConstPtr;
903 static bool isRegistered() {
return MapRegistry::isRegistered(TranslationMap::mapType()); }
905 static void registerMap()
907 MapRegistry::registerMap(
908 TranslationMap::mapType(),
909 TranslationMap::create);
935 return applyIJC(mat);
951 void read(std::istream& is) { mTranslation.read(is); }
953 void write(std::ostream& os)
const { mTranslation.write(os); }
956 std::string str()
const
958 std::ostringstream buffer;
959 buffer <<
" - translation: " << mTranslation << std::endl;
963 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
968 if (!mTranslation.eq(other.mTranslation)) {
return false; }
990 affineMap->accumPreRotation(axis, radians);
1004 affineMap->accumPreShear(axis0, axis1, shear);
1015 affineMap->accumPostRotation(axis, radians);
1029 affineMap->accumPostShear(axis0, axis1, shear);
1048 typedef boost::shared_ptr<ScaleTranslateMap>
Ptr;
1049 typedef boost::shared_ptr<const ScaleTranslateMap>
ConstPtr;
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)
1064 mTranslation(translate),
1065 mScaleValues(scale),
1066 mVoxelSize(std::abs(scale(0)), std::abs(scale(1)), std::abs(scale(2)))
1068 if (scale.
eq(
Vec3d(0.0), 1.0e-10)) {
1071 mScaleValuesInverse = 1.0 / mScaleValues;
1072 mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1073 mInvTwiceScale = mScaleValuesInverse / 2;
1078 mTranslation(translate.getTranslation()),
1080 mVoxelSize(std::abs(mScaleValues(0)),
1081 std::abs(mScaleValues(1)),
1082 std::abs(mScaleValues(2))),
1083 mScaleValuesInverse(1.0 / scale.
getScale())
1085 mInvScaleSqr = mScaleValuesInverse * mScaleValuesInverse;
1086 mInvTwiceScale = mScaleValuesInverse / 2;
1091 mTranslation(other.mTranslation),
1092 mScaleValues(other.mScaleValues),
1093 mVoxelSize(other.mVoxelSize),
1094 mScaleValuesInverse(other.mScaleValuesInverse),
1095 mInvScaleSqr(other.mInvScaleSqr),
1096 mInvTwiceScale(other.mInvTwiceScale)
1106 static bool isRegistered() {
return MapRegistry::isRegistered(ScaleTranslateMap::mapType()); }
1108 static void registerMap()
1110 MapRegistry::registerMap(
1111 ScaleTranslateMap::mapType(),
1112 ScaleTranslateMap::create);
1125 value =
isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.y()));
1126 value = value &&
isApproxEqual(std::abs(mScaleValues.x()), std::abs(mScaleValues.z()));
1134 in.
x() * mScaleValues.x() + mTranslation.x(),
1135 in.
y() * mScaleValues.y() + mTranslation.y(),
1136 in.
z() * mScaleValues.z() + mTranslation.z());
1142 (in.
x() - mTranslation.x() ) / mScaleValues.x(),
1143 (in.
y() - mTranslation.y() ) / mScaleValues.y(),
1144 (in.
z() - mTranslation.z() ) / mScaleValues.z());
1153 in.
x() / mScaleValues.x(),
1154 in.
y() / mScaleValues.y(),
1155 in.
z() / mScaleValues.z());
1160 for (
int i=0; i<3; i++){
1161 tmp.
setRow(i, in.
row(i)*mScaleValuesInverse(i));
1163 for (
int i=0; i<3; i++){
1164 tmp.
setCol(i, tmp.
col(i)*mScaleValuesInverse(i));
1169 return applyIJC(in);
1175 double determinant()
const {
return mScaleValues.x()*mScaleValues.y()*mScaleValues.z(); }
1195 void read(std::istream& is)
1197 mTranslation.
read(is);
1198 mScaleValues.read(is);
1199 mVoxelSize.read(is);
1200 mScaleValuesInverse.read(is);
1201 mInvScaleSqr.read(is);
1202 mInvTwiceScale.read(is);
1205 void write(std::ostream& os)
const
1207 mTranslation.write(os);
1208 mScaleValues.write(os);
1209 mVoxelSize.write(os);
1210 mScaleValuesInverse.write(os);
1211 mInvScaleSqr.write(os);
1212 mInvTwiceScale.write(os);
1215 std::string str()
const
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();
1224 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1229 if (!mScaleValues.eq(other.mScaleValues)) {
return false; }
1230 if (!mTranslation.eq(other.mTranslation)) {
return false; }
1240 affineMap->accumPostTranslation(mTranslation);
1250 affineMap->accumPreRotation(axis, radians);
1255 const Vec3d& s = mScaleValues;
1256 const Vec3d scaled_trans( t.
x() * s.
x(),
1267 affineMap->accumPreShear(axis0, axis1, shear);
1278 affineMap->accumPostRotation(axis, radians);
1291 affineMap->accumPostShear(axis0, axis1, shear);
1297 Vec3d mTranslation, mScaleValues, mVoxelSize, mScaleValuesInverse,
1298 mInvScaleSqr, mInvTwiceScale;
1303 ScaleMap::postTranslate(
const Vec3d& t)
const
1310 ScaleMap::preTranslate(
const Vec3d& t)
const
1313 const Vec3d& s = mScaleValues;
1314 const Vec3d scaled_trans( t.
x() * s.
x(),
1326 typedef boost::shared_ptr<UniformScaleTranslateMap>
Ptr;
1327 typedef boost::shared_ptr<const UniformScaleTranslateMap>
ConstPtr;
1343 static bool isRegistered()
1345 return MapRegistry::isRegistered(UniformScaleTranslateMap::mapType());
1348 static void registerMap()
1350 MapRegistry::registerMap(
1351 UniformScaleTranslateMap::mapType(),
1352 UniformScaleTranslateMap::create);
1358 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1371 const Vec3d new_trans = this->getTranslation() + scale * t;
1386 UniformScaleMap::postTranslate(
const Vec3d& t)
const
1394 UniformScaleMap::preTranslate(
const Vec3d& t)
const
1401 TranslationMap::preScale(
const Vec3d& v)
const
1411 TranslationMap::postScale(
const Vec3d& v)
const
1416 const Vec3d trans(mTranslation.x()*v.
x(),
1417 mTranslation.y()*v.
y(),
1418 mTranslation.z()*v.
z());
1424 ScaleTranslateMap::preScale(
const Vec3d& v)
const
1426 Vec3d new_scale( v * mScaleValues );
1435 ScaleTranslateMap::postScale(
const Vec3d& v)
const
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() );
1458 typedef boost::shared_ptr<UnitaryMap>
Ptr;
1496 "4x4 Matrix initializing unitary map was not unitary: not invertible");
1501 "4x4 Matrix initializing unitary map was not unitary: not affine");
1506 "4x4 Matrix initializing unitary map was not unitary: had translation");
1511 "4x4 Matrix initializing unitary map was not unitary");
1519 mAffineMap(other.mAffineMap)
1524 mAffineMap(*(first.getAffineMap()), *(second.getAffineMap()))
1534 static bool isRegistered() {
return MapRegistry::isRegistered(UnitaryMap::mapType()); }
1536 static void registerMap()
1538 MapRegistry::registerMap(
1539 UnitaryMap::mapType(),
1540 UnitaryMap::create);
1554 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1559 if (mAffineMap!=other.mAffineMap)
return false;
1575 return mAffineMap.applyIJC(in);
1578 return applyIJC(in);
1594 void read(std::istream& is)
1596 mAffineMap.read(is);
1600 void write(std::ostream& os)
const
1602 mAffineMap.write(os);
1605 std::string str()
const
1607 std::ostringstream buffer;
1608 buffer << mAffineMap.str();
1609 return buffer.str();
1626 affineMap->accumPreTranslation(t);
1632 affineMap->accumPreScale(v);
1638 affineMap->accumPreShear(axis0, axis1, shear);
1656 affineMap->accumPostTranslation(t);
1662 affineMap->accumPostScale(v);
1668 affineMap->accumPostShear(axis0, axis1, shear);
1690 typedef boost::shared_ptr<NonlinearFrustumMap>
Ptr;
1691 typedef boost::shared_ptr<const NonlinearFrustumMap>
ConstPtr;
1706 MapBase(),mBBox(bb), mTaper(taper), mDepth(depth)
1718 mBBox(bb), mTaper(taper), mDepth(depth)
1720 if (!secondMap->isLinear() ) {
1722 "The second map in the Frustum transfrom must be linear");
1724 mSecondMap = *( secondMap->getAffineMap() );
1731 mTaper(other.mTaper),
1732 mDepth(other.mDepth),
1733 mSecondMap(other.mSecondMap),
1734 mHasSimpleAffine(other.mHasSimpleAffine)
1755 const Vec3d& direction,
1758 double z_near,
double depth,
1766 "The frustum depth must be non-zero and positive");
1768 if (!(up.
length() > 0)) {
1770 "The frustum height must be non-zero and positive");
1772 if (!(aspect > 0)) {
1774 "The frustum aspect ratio must be non-zero and positive");
1778 "The frustum up orientation must be perpendicular to into-frustum direction");
1781 double near_plane_height = 2 * up.
length();
1782 double near_plane_width = aspect * near_plane_height;
1787 mDepth = depth / near_plane_width;
1788 double gamma = near_plane_width / z_near;
1789 mTaper = 1./(mDepth*gamma + 1.);
1791 Vec3d direction_unit = direction;
1800 Vec3d(near_plane_width, near_plane_width, near_plane_width));
1804 Mat4d mat = scale * r2 * r1;
1818 static bool isRegistered() {
return MapRegistry::isRegistered(NonlinearFrustumMap::mapType()); }
1820 static void registerMap()
1822 MapRegistry::registerMap(
1823 NonlinearFrustumMap::mapType(),
1824 NonlinearFrustumMap::create);
1845 const Vec3d e1(1,0,0);
1846 if (!applyMap(e1).eq(e1))
return false;
1848 const Vec3d e2(0,1,0);
1849 if (!applyMap(e2).eq(e2))
return false;
1851 const Vec3d e3(0,0,1);
1852 if (!applyMap(e3).eq(e3))
return false;
1857 virtual bool isEqual(
const MapBase& other)
const {
return isEqualBase(*
this, other); }
1861 if (mBBox!=other.mBBox)
return false;
1868 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
1871 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
1874 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
1877 if (!mSecondMap.applyMap(e).eq(other.mSecondMap.
applyMap(e)))
return false;
1886 return mSecondMap.applyMap(applyFrustumMap(in));
1892 return applyFrustumInverseMap(mSecondMap.applyInverseMap(in));
1908 const Vec3d loc = applyFrustumMap(ijk);
1909 const double s = mGamma * loc.
z() + 1.;
1914 " at the singular focal point (e.g. camera)");
1917 const double sinv = 1.0/s;
1918 const double pt0 = mLx * sinv;
1919 const double pt1 = mGamma * pt0;
1920 const double pt2 = pt1 * sinv;
1922 const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
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);
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);
1949 const Vec3d loc = applyFrustumMap(ijk);
1951 const double s = mGamma * loc.
z() + 1.;
1956 " at the singular focal point (e.g. camera)");
1960 const double sinv = 1.0/s;
1961 const double pt0 = mLx * sinv;
1962 const double pt1 = mGamma * pt0;
1963 const double pt2 = pt1 * sinv;
1964 const double pt3 = pt2 * sinv;
1967 const Mat3d& jacinv = mSecondMap.getConstJacobianInv();
1971 Mat3d matE0(Mat3d::zero());
1972 Mat3d matE1(Mat3d::zero());
1973 for(
int j = 0; j < 3; j++) {
1974 for (
int k = 0; k < 3; k++) {
1976 const double pt4 = 2. * jacinv(2,j) * jacinv(2,k) * pt3;
1978 matE0(j,k) = -(jacinv(0,j) * jacinv(2,k) + jacinv(2,j) * jacinv(0,k)) * pt2 +
1981 matE1(j,k) = -(jacinv(1,j) * jacinv(2,k) + jacinv(2,j) * jacinv(1,k)) * pt2 +
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);
1994 Mat3d result(Mat3d::zero());
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);
2007 for (
int m = 0; m < 3; ++m ) {
2008 for (
int n = 0; n < 3; ++n) {
2010 matE0(m, n) * d1_is(0) + matE1(m, n) * d1_is(1);
2022 double determinant(
const Vec3d& loc)
const
2024 double s = mGamma * loc.
z() + 1.0;
2025 double frustum_determinant = s * s * mDepthOnLzLxLx;
2026 return mSecondMap.determinant() * frustum_determinant;
2031 const Vec3d loc( 0.5*(mBBox.min().x() + mBBox.max().x()),
2032 0.5*(mBBox.min().y() + mBBox.max().y()),
2035 return voxelSize(loc);
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();
2053 return mSecondMap.getAffineMap();
2080 void read(std::istream& is)
2091 is.read(reinterpret_cast<char*>(&mTaper),
sizeof(
double));
2092 is.read(reinterpret_cast<char*>(&mDepth),
sizeof(
double));
2098 if(!MapRegistry::isRegistered(type)) {
2103 MapBase::Ptr proxy = math::MapRegistry::createMap(type);
2105 mSecondMap = *(proxy->getAffineMap());
2110 void write(std::ostream& os)
const
2113 os.write(reinterpret_cast<const char*>(&mTaper),
sizeof(
double));
2114 os.write(reinterpret_cast<const char*>(&mDepth),
sizeof(
double));
2117 mSecondMap.write(os);
2121 std::string str()
const
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();
2152 mBBox, mTaper, mDepth, mSecondMap.preShear(shear, axis0, axis1)));
2177 mBBox, mTaper, mDepth, mSecondMap.postShear(shear, axis0, axis1)));
2185 mLx = mBBox.extents().x();
2186 mLy = mBBox.extents().y();
2187 mLz = mBBox.extents().z();
2191 " must have at least two index points in each direction.");
2198 mGamma = (1./mTaper - 1) / mDepth;
2200 mDepthOnLz = mDepth/mLz;
2201 mDepthOnLzLxLx = mDepthOnLz/(mLx * mLx);
2204 mHasSimpleAffine =
true;
2205 Vec3d tmp = mSecondMap.voxelSize();
2208 if (!
isApproxEqual(tmp(0), tmp(1))) { mHasSimpleAffine =
false;
return; }
2209 if (!
isApproxEqual(tmp(0), tmp(2))) { mHasSimpleAffine =
false;
return; }
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;
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; }
2229 out = out - mBBox.min();
2234 out.z() *= mDepthOnLz;
2236 double scale = (mGamma * out.z() + 1.)/ mLx;
2245 Vec3d applyFrustumInverseMap(
const Vec3d& in)
const
2249 double invScale = mLx / (mGamma * out.z() + 1.);
2250 out.x() *= invScale;
2251 out.y() *= invScale;
2256 out.z() /= mDepthOnLz;
2259 out = out + mBBox.min();
2271 AffineMap mSecondMap;
2274 double mLx, mLy, mLz;
2275 double mXo, mYo, mGamma, mDepthOnLz, mDepthOnLzLxLx;
2278 bool mHasSimpleAffine;
2288 template<
typename FirstMapType,
typename SecondMapType>
2294 typedef boost::shared_ptr<MyType>
Ptr;
2300 CompoundMap(
const FirstMapType& f,
const SecondMapType& s): mFirstMap(f), mSecondMap(s)
2302 updateAffineMatrix();
2306 mFirstMap(other.mFirstMap),
2307 mSecondMap(other.mSecondMap),
2308 mAffineMap(other.mAffineMap)
2314 return (FirstMapType::mapType() +
Name(
":") + SecondMapType::mapType());
2319 if (mFirstMap != other.mFirstMap)
return false;
2320 if (mSecondMap != other.mSecondMap)
return false;
2321 if (mAffineMap != other.mAffineMap)
return false;
2329 mFirstMap = other.mFirstMap;
2330 mSecondMap = other.mSecondMap;
2331 mAffineMap = other.mAffineMap;
2337 return mAffineMap.isIdentity();
2339 return mFirstMap.isIdentity()&&mSecondMap.isIdentity();
2345 return mAffineMap.isDiagonal();
2347 return mFirstMap.isDiagonal()&&mSecondMap.isDiagonal();
2357 "Constant affine matrix representation not possible for this nonlinear map");
2362 const FirstMapType&
firstMap()
const {
return mFirstMap; }
2363 const SecondMapType&
secondMap()
const {
return mSecondMap; }
2365 void setFirstMap(
const FirstMapType& first) { mFirstMap = first; updateAffineMatrix(); }
2366 void setSecondMap(
const SecondMapType& second) { mSecondMap = second; updateAffineMatrix(); }
2368 void read(std::istream& is)
2370 mAffineMap.read(is);
2372 mSecondMap.read(is);
2374 void write(std::ostream& os)
const
2376 mAffineMap.write(os);
2377 mFirstMap.write(os);
2378 mSecondMap.write(os);
2382 void updateAffineMatrix()
2388 mAffineMap =
AffineMap(*first, *second);
2392 FirstMapType mFirstMap;
2393 SecondMapType mSecondMap;
2395 AffineMap mAffineMap;
2402 #endif // OPENVDB_MATH_MAPS_HAS_BEEN_INCLUDED