From 1d922d71f7b11e8ae611cafb047e00aef9d19d0b Mon Sep 17 00:00:00 2001 From: Nick Rasmussen Date: Fri, 8 Apr 2011 22:24:40 +0000 Subject: [PATCH] Updated Imath with the latest internal changes: + near() and far() in ImathFrustum are now updated to be nearPlane() and farPlane() to prevent conflicts with #defines on windows. The hither() and yon() functions remain but can be considered deprecated. + ImathBox now supports infinite functions analogous to the empty functions (isInfinite(), makeInfinite()), and the ImathBoxAlgo transform code has been updated accordingly. + ImathBoxAlgo has an additional transform() overload that takes the result as the last argument. + ImathMatrix now has matrix minor and detriminant functions. + ImathMatrixAlgo now has outerProduct(), addOffset() and computeLocalFrame() functions. --- IlmBase/Imath/ImathBox.h | 63 +++++ IlmBase/Imath/ImathBoxAlgo.h | 157 ++++++++++- IlmBase/Imath/ImathFrustum.h | 178 ++++++------ IlmBase/Imath/ImathMatrix.h | 178 +++++++++++- IlmBase/Imath/ImathMatrixAlgo.h | 245 ++++++++++++++-- IlmBase/ImathTest/Makefile.am | 3 +- IlmBase/ImathTest/main.cpp | 6 +- IlmBase/ImathTest/testBox.cpp | 145 ++++++++++ IlmBase/ImathTest/testBoxAlgo.cpp | 13 + IlmBase/ImathTest/testExtractEuler.cpp | 6 +- IlmBase/ImathTest/testFrustum.cpp | 4 +- IlmBase/ImathTest/testMatrix.cpp | 237 ++++++++++++++-- IlmBase/ImathTest/testMiscMatrixAlgo.cpp | 341 +++++++++++++++++++++++ IlmBase/ImathTest/testMiscMatrixAlgo.h | 39 +++ 14 files changed, 1447 insertions(+), 168 deletions(-) create mode 100644 IlmBase/ImathTest/testMiscMatrixAlgo.cpp create mode 100644 IlmBase/ImathTest/testMiscMatrixAlgo.h diff --git a/IlmBase/Imath/ImathBox.h b/IlmBase/Imath/ImathBox.h index 0e8af4f..7dc4eb7 100644 --- a/IlmBase/Imath/ImathBox.h +++ b/IlmBase/Imath/ImathBox.h @@ -101,6 +101,7 @@ class Box void makeEmpty (); void extendBy (const T &point); void extendBy (const Box &box); + void makeInfinite (); //--------------------------------------------------- // Query functions - these compute results each time @@ -119,6 +120,7 @@ class Box bool isEmpty () const; bool hasVolume () const; + bool isInfinite () const; }; @@ -186,6 +188,13 @@ inline void Box::makeEmpty() max = T(T::baseTypeMin()); } +template +inline void Box::makeInfinite() +{ + min = T(T::baseTypeMin()); + max = T(T::baseTypeMax()); +} + template inline void @@ -277,6 +286,19 @@ Box::isEmpty() const return false; } +template +inline bool +Box::isInfinite() const +{ + for (unsigned int i = 0; i < min.dimensions(); i++) + { + if (min[i] != T::baseTypeMin() || max[i] != T::baseTypeMax()) + return false; + } + + return true; +} + template inline bool @@ -350,6 +372,7 @@ class Box > void makeEmpty(); void extendBy (const Vec2 &point); void extendBy (const Box > &box); + void makeInfinite(); //--------------------------------------------------- // Query functions - these compute results each time @@ -368,6 +391,7 @@ class Box > bool isEmpty() const; bool hasVolume() const; + bool isInfinite() const; }; @@ -420,6 +444,13 @@ inline void Box >::makeEmpty() max = Vec2(Vec2::baseTypeMin()); } +template +inline void Box >::makeInfinite() +{ + min = Vec2(Vec2::baseTypeMin()); + max = Vec2(Vec2::baseTypeMax()); +} + template inline void @@ -511,6 +542,17 @@ Box >::isEmpty() const return false; } +template +inline bool +Box > ::isInfinite() const +{ + if (min[0] != limits::min() || max[0] != limits::max() || + min[1] != limits::min() || max[1] != limits::max()) + return false; + + return true; +} + template inline bool @@ -572,6 +614,7 @@ class Box > void makeEmpty(); void extendBy (const Vec3 &point); void extendBy (const Box > &box); + void makeInfinite (); //--------------------------------------------------- // Query functions - these compute results each time @@ -590,6 +633,7 @@ class Box > bool isEmpty() const; bool hasVolume() const; + bool isInfinite() const; }; @@ -643,6 +687,13 @@ inline void Box >::makeEmpty() max = Vec3(Vec3::baseTypeMin()); } +template +inline void Box >::makeInfinite() +{ + min = Vec3(Vec3::baseTypeMin()); + max = Vec3(Vec3::baseTypeMax()); +} + template inline void @@ -749,6 +800,18 @@ Box >::isEmpty() const return false; } +template +inline bool +Box >::isInfinite() const +{ + if (min[0] != limits::min() || max[0] != limits::max() || + min[1] != limits::min() || max[1] != limits::max() || + min[2] != limits::min() || max[2] != limits::max()) + return false; + + return true; +} + template inline bool diff --git a/IlmBase/Imath/ImathBoxAlgo.h b/IlmBase/Imath/ImathBoxAlgo.h index 4d3003f..54bd71a 100644 --- a/IlmBase/Imath/ImathBoxAlgo.h +++ b/IlmBase/Imath/ImathBoxAlgo.h @@ -1,6 +1,6 @@ /////////////////////////////////////////////////////////////////////////// // -// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Copyright (c) 2002-2010, Industrial Light & Magic, a division of Lucas // Digital Ltd. LLC // // All rights reserved. @@ -54,7 +54,13 @@ // // Vec3 closestPointInBox(const Vec3&, const Box>& ) // -// void transform(Box>&, const Matrix44&) +// Box< Vec3 > transform(const Box>&, const Matrix44&) +// Box< Vec3 > affineTransform(const Box>&, const Matrix44&) +// +// void transform(const Box>&, const Matrix44&, Box>&) +// void affineTransform(const Box>&, +// const Matrix44&, +// Box>&) // // bool findEntryAndExitPoints(const Line &line, // const Box< Vec3 > &box, @@ -167,10 +173,11 @@ transform (const Box< Vec3 > &box, const Matrix44 &m) // // - // A transformed empty box is still empty + // A transformed empty box is still empty, and a transformed infinite box + // is still infinite // - if (box.isEmpty()) + if (box.isEmpty() || box.isInfinite()) return box; // @@ -234,6 +241,86 @@ transform (const Box< Vec3 > &box, const Matrix44 &m) return newBox; } +template +void +transform (const Box< Vec3 > &box, + const Matrix44 &m, + Box< Vec3 > &result) +{ + // + // Transform a 3D box by a matrix, and compute a new box that + // tightly encloses the transformed box. + // + // If m is an affine transform, then we use James Arvo's fast + // method as described in "Graphics Gems", Academic Press, 1990, + // pp. 548-550. + // + + // + // A transformed empty box is still empty, and a transformed infinite + // box is still infinite + // + + if (box.isEmpty() || box.isInfinite()) + { + return; + } + + // + // If the last column of m is (0 0 0 1) then m is an affine + // transform, and we use the fast Graphics Gems trick. + // + + if (m[0][3] == 0 && m[1][3] == 0 && m[2][3] == 0 && m[3][3] == 1) + { + for (int i = 0; i < 3; i++) + { + result.min[i] = result.max[i] = (S) m[3][i]; + + for (int j = 0; j < 3; j++) + { + float a, b; + + a = (S) m[j][i] * box.min[j]; + b = (S) m[j][i] * box.max[j]; + + if (a < b) + { + result.min[i] += a; + result.max[i] += b; + } + else + { + result.min[i] += b; + result.max[i] += a; + } + } + } + + return; + } + + // + // M is a projection matrix. Do things the naive way: + // Transform the eight corners of the box, and find an + // axis-parallel box that encloses the transformed corners. + // + + Vec3 points[8]; + + points[0][0] = points[1][0] = points[2][0] = points[3][0] = box.min[0]; + points[4][0] = points[5][0] = points[6][0] = points[7][0] = box.max[0]; + + points[0][1] = points[1][1] = points[4][1] = points[5][1] = box.min[1]; + points[2][1] = points[3][1] = points[6][1] = points[7][1] = box.max[1]; + + points[0][2] = points[2][2] = points[4][2] = points[6][2] = box.min[2]; + points[1][2] = points[3][2] = points[5][2] = points[7][2] = box.max[2]; + + for (int i = 0; i < 8; i++) + result.extendBy (points[i] * m); +} + template Box< Vec3 > @@ -248,10 +335,10 @@ affineTransform (const Box< Vec3 > &box, const Matrix44 &m) // fast method. // - if (box.isEmpty()) + if (box.isEmpty() || box.isInfinite()) { // - // A transformed empty box is still empty + // A transformed empty or infinite box is still empty or infinite // return box; @@ -286,6 +373,64 @@ affineTransform (const Box< Vec3 > &box, const Matrix44 &m) return newBox; } +template +void +affineTransform (const Box< Vec3 > &box, + const Matrix44 &m, + Box > &result) +{ + // + // Transform a 3D box by a matrix whose rightmost column + // is (0 0 0 1), and compute a new box that tightly encloses + // the transformed box. + // + // As in the transform() function, above, we use James Arvo's + // fast method. + // + + if (box.isEmpty()) + { + // + // A transformed empty box is still empty + // + result.makeEmpty(); + return; + } + + if (box.isInfinite()) + { + // + // A transformed infinite box is still infinite + // + result.makeInfinite(); + return; + } + + for (int i = 0; i < 3; i++) + { + result.min[i] = result.max[i] = (S) m[3][i]; + + for (int j = 0; j < 3; j++) + { + float a, b; + + a = (S) m[j][i] * box.min[j]; + b = (S) m[j][i] * box.max[j]; + + if (a < b) + { + result.min[i] += a; + result.max[i] += b; + } + else + { + result.min[i] += b; + result.max[i] += a; + } + } + } +} + template bool diff --git a/IlmBase/Imath/ImathFrustum.h b/IlmBase/Imath/ImathFrustum.h index 160ffc0..cf5d2be 100644 --- a/IlmBase/Imath/ImathFrustum.h +++ b/IlmBase/Imath/ImathFrustum.h @@ -46,17 +46,6 @@ #include "ImathFun.h" #include "IexMathExc.h" -#if defined _WIN32 || defined _WIN64 - #ifdef near - #define _redef_near - #undef near - #endif - #ifdef far - #define _redef_far - #undef far - #endif -#endif - namespace Imath { // @@ -70,6 +59,9 @@ namespace Imath { // +Z. Additional functions are provided for conversion from // and from various camera coordinate spaces. // +// nearPlane/farPlane: near/far are keywords used by Microsoft's +// compiler, so we use nearPlane/farPlane instead to avoid +// issues. template @@ -78,8 +70,8 @@ class Frustum public: Frustum(); Frustum(const Frustum &); - Frustum(T near, T far, T left, T right, T top, T bottom, bool ortho=false); - Frustum(T near, T far, T fovx, T fovy, T aspect); + Frustum(T nearPlane, T farPlane, T left, T right, T top, T bottom, bool ortho=false); + Frustum(T nearPlane, T farPlane, T fovx, T fovy, T aspect); virtual ~Frustum(); //-------------------- @@ -99,29 +91,29 @@ class Frustum // Set functions change the entire state of the Frustum //-------------------------------------------------------- - void set(T near, T far, + void set(T nearPlane, T farPlane, T left, T right, T top, T bottom, bool ortho=false); - void set(T near, T far, T fovx, T fovy, T aspect); + void set(T nearPlane, T farPlane, T fovx, T fovy, T aspect); //------------------------------------------------------ // These functions modify an already valid frustum state //------------------------------------------------------ - void modifyNearAndFar(T near, T far); + void modifyNearAndFar(T nearPlane, T farPlane); void setOrthographic(bool); //-------------- // Access //-------------- - + bool orthographic() const { return _orthographic; } - T near() const { return _near; } - T hither() const { return _near; } - T far() const { return _far; } - T yon() const { return _far; } + T nearPlane() const { return _nearPlane; } + T hither() const { return _nearPlane; } + T farPlane() const { return _farPlane; } + T yon() const { return _farPlane; } T left() const { return _left; } T right() const { return _right; } T bottom() const { return _bottom; } @@ -177,8 +169,8 @@ class Frustum Vec2 localToScreen( const Vec2 & ) const; protected: - T _near; - T _far; + T _nearPlane; + T _farPlane; T _left; T _right; T _top; @@ -212,9 +204,9 @@ inline Frustum::Frustum(T n, T f, T l, T r, T t, T b, bool o) } template -inline Frustum::Frustum(T near, T far, T fovx, T fovy, T aspect) +inline Frustum::Frustum(T nearPlane, T farPlane, T fovx, T fovy, T aspect) { - set(near,far,fovx,fovy,aspect); + set(nearPlane,farPlane,fovx,fovy,aspect); } template @@ -226,8 +218,8 @@ template const Frustum & Frustum::operator = (const Frustum &f) { - _near = f._near; - _far = f._far; + _nearPlane = f._nearPlane; + _farPlane = f._farPlane; _left = f._left; _right = f._right; _top = f._top; @@ -242,8 +234,8 @@ bool Frustum::operator == (const Frustum &src) const { return - _near == src._near && - _far == src._far && + _nearPlane == src._nearPlane && + _farPlane == src._farPlane && _left == src._left && _right == src._right && _top == src._top && @@ -261,8 +253,8 @@ Frustum::operator != (const Frustum &src) const template void Frustum::set(T n, T f, T l, T r, T t, T b, bool o) { - _near = n; - _far = f; + _nearPlane = n; + _farPlane = f; _left = l; _right = r; _bottom = b; @@ -275,27 +267,27 @@ void Frustum::modifyNearAndFar(T n, T f) { if ( _orthographic ) { - _near = n; + _nearPlane = n; } else { - Line3 lowerLeft( Vec3(0,0,0), Vec3(_left,_bottom,-_near) ); - Line3 upperRight( Vec3(0,0,0), Vec3(_right,_top,-_near) ); + Line3 lowerLeft( Vec3(0,0,0), Vec3(_left,_bottom,-_nearPlane) ); + Line3 upperRight( Vec3(0,0,0), Vec3(_right,_top,-_nearPlane) ); Plane3 nearPlane( Vec3(0,0,-1), n ); Vec3 ll,ur; nearPlane.intersect(lowerLeft,ll); nearPlane.intersect(upperRight,ur); - _left = ll.x; - _right = ur.x; - _top = ur.y; - _bottom = ll.y; - _near = n; - _far = f; + _left = ll.x; + _right = ur.x; + _top = ur.y; + _bottom = ll.y; + _nearPlane = n; + _farPlane = f; } - _far = f; + _farPlane = f; } template @@ -305,40 +297,42 @@ void Frustum::setOrthographic(bool ortho) } template -void Frustum::set(T near, T far, T fovx, T fovy, T aspect) +void Frustum::set(T nearPlane, T farPlane, T fovx, T fovy, T aspect) { if (fovx != 0 && fovy != 0) throw Iex::ArgExc ("fovx and fovy cannot both be non-zero."); + const T two = static_cast(2); + if (fovx != 0) { - _right = near * Math::tan (fovx / 2); + _right = nearPlane * Math::tan(fovx / two); _left = -_right; - _top = ((_right - _left) / aspect) / 2; + _top = ((_right - _left) / aspect) / two; _bottom = -_top; } else { - _top = near * Math::tan (fovy / 2); + _top = nearPlane * Math::tan(fovy / two); _bottom = -_top; - _right = (_top - _bottom) * aspect / 2; + _right = (_top - _bottom) * aspect / two; _left = -_right; } - _near = near; - _far = far; + _nearPlane = nearPlane; + _farPlane = farPlane; _orthographic = false; } template T Frustum::fovx() const { - return Math::atan2(_right,_near) - Math::atan2(_left,_near); + return Math::atan2(_right,_nearPlane) - Math::atan2(_left,_nearPlane); } template T Frustum::fovy() const { - return Math::atan2(_top,_near) - Math::atan2(_bottom,_near); + return Math::atan2(_top,_nearPlane) - Math::atan2(_bottom,_nearPlane); } template @@ -366,8 +360,8 @@ Matrix44 Frustum::projectionMatrix() const T topPlusBottom = _top+_bottom; T topMinusBottom = _top-_bottom; - T farPlusNear = _far+_near; - T farMinusNear = _far-_near; + T farPlusNear = _farPlane+_nearPlane; + T farMinusNear = _farPlane-_nearPlane; if ((abs(rightMinusLeft) < 1 && abs(rightPlusLeft) > limits::max() * abs(rightMinusLeft)) || @@ -412,7 +406,7 @@ Matrix44 Frustum::projectionMatrix() const T B = topPlusBottom / topMinusBottom; T C = -farPlusNear / farMinusNear; - T farTimesNear = -2 * _far * _near; + T farTimesNear = -2 * _farPlane * _nearPlane; if (abs(farMinusNear) < 1 && abs(farTimesNear) > limits::max() * abs(farMinusNear)) { @@ -422,7 +416,7 @@ Matrix44 Frustum::projectionMatrix() const T D = farTimesNear / farMinusNear; - T twoTimesNear = 2 * _near; + T twoTimesNear = 2 * _nearPlane; if ((abs(rightMinusLeft) < 1 && abs(twoTimesNear) > limits::max() * abs(rightMinusLeft)) || @@ -451,7 +445,7 @@ Frustum Frustum::window(T l, T r, T t, T b) const Vec2 bl = screenToLocal( Vec2(l,b) ); Vec2 tr = screenToLocal( Vec2(r,t) ); - return Frustum(_near, _far, bl.x, tr.x, tr.y, bl.y, _orthographic); + return Frustum(_nearPlane, _farPlane, bl.x, tr.x, tr.y, bl.y, _orthographic); } @@ -490,9 +484,9 @@ Line3 Frustum::projectScreenToRay(const Vec2 &p) const Vec2 point = screenToLocal(p); if (orthographic()) return Line3( Vec3(point.x,point.y, 0.0), - Vec3(point.x,point.y,-_near)); + Vec3(point.x,point.y,-_nearPlane)); else - return Line3( Vec3(0, 0, 0), Vec3(point.x,point.y,-_near)); + return Line3( Vec3(0, 0, 0), Vec3(point.x,point.y,-_nearPlane)); } template @@ -501,8 +495,8 @@ Vec2 Frustum::projectPointToScreen(const Vec3 &point) const if (orthographic() || point.z == 0) return localToScreen( Vec2( point.x, point.y ) ); else - return localToScreen( Vec2( point.x * _near / -point.z, - point.y * _near / -point.z ) ); + return localToScreen( Vec2( point.x * _nearPlane / -point.z, + point.y * _nearPlane / -point.z ) ); } template @@ -529,12 +523,12 @@ T Frustum::normalizedZToDepth(T zval) const if ( _orthographic ) { - return -(Zp*(_far-_near) + (_far+_near))/2; + return -(Zp*(_farPlane-_nearPlane) + (_farPlane+_nearPlane))/2; } else { - T farTimesNear = 2 * _far * _near; - T farMinusNear = Zp * (_far - _near) - _far - _near; + T farTimesNear = 2 * _farPlane * _nearPlane; + T farMinusNear = Zp * (_farPlane - _nearPlane) - _farPlane - _nearPlane; if (abs(farMinusNear) < 1 && abs(farTimesNear) > limits::max() * abs(farMinusNear)) @@ -553,11 +547,11 @@ template long Frustum::DepthToZ(T depth,long zmin,long zmax) const { long zdiff = zmax - zmin; - T farMinusNear = _far-_near; + T farMinusNear = _farPlane-_nearPlane; if ( _orthographic ) { - T farPlusNear = 2*depth + _far + _near; + T farPlusNear = 2*depth + _farPlane + _nearPlane; if (abs(farMinusNear) < 1 && abs(farPlusNear) > limits::max() * abs(farMinusNear)) @@ -574,7 +568,7 @@ long Frustum::DepthToZ(T depth,long zmin,long zmax) const { // Perspective - T farTimesNear = 2*_far*_near; + T farTimesNear = 2*_farPlane*_nearPlane; if (abs(depth) < 1 && abs(farTimesNear) > limits::max() * abs(depth)) { @@ -583,7 +577,7 @@ long Frustum::DepthToZ(T depth,long zmin,long zmax) const "is too small"); } - T farPlusNear = farTimesNear/depth + _far + _near; + T farPlusNear = farTimesNear/depth + _farPlane + _nearPlane; if (abs(farMinusNear) < 1 && abs(farPlusNear) > limits::max() * abs(farMinusNear)) { @@ -602,17 +596,17 @@ T Frustum::screenRadius(const Vec3 &p, T radius) const { // Derivation: // Consider X-Z plane. - // X coord of projection of p = xp = p.x * (-_near / p.z) + // X coord of projection of p = xp = p.x * (-_nearPlane / p.z) // Let q be p + (radius, 0, 0). - // X coord of projection of q = xq = (p.x - radius) * (-_near / p.z) + // X coord of projection of q = xq = (p.x - radius) * (-_nearPlane / p.z) // X coord of projection of segment from p to q = r = xp - xq - // = radius * (-_near / p.z) + // = radius * (-_nearPlane / p.z) // A similar analysis holds in the Y-Z plane. // So r is the quantity we want to return. - if (abs(p.z) > 1 || abs(-_near) < limits::max() * abs(p.z)) + if (abs(p.z) > 1 || abs(-_nearPlane) < limits::max() * abs(p.z)) { - return radius * (-_near / p.z); + return radius * (-_nearPlane / p.z); } else { @@ -621,15 +615,15 @@ T Frustum::screenRadius(const Vec3 &p, T radius) const "is too small"); } - return radius * (-_near / p.z); + return radius * (-_nearPlane / p.z); } template T Frustum::worldRadius(const Vec3 &p, T radius) const { - if (abs(-_near) > 1 || abs(p.z) < limits::max() * abs(-_near)) + if (abs(-_nearPlane) > 1 || abs(p.z) < limits::max() * abs(-_nearPlane)) { - return radius * (p.z / -_near); + return radius * (p.z / -_nearPlane); } else { @@ -649,10 +643,10 @@ void Frustum::planes(Plane3 p[6]) if (! _orthographic) { - Vec3 a( _left, _bottom, -_near); - Vec3 b( _left, _top, -_near); - Vec3 c( _right, _top, -_near); - Vec3 d( _right, _bottom, -_near); + Vec3 a( _left, _bottom, -_nearPlane); + Vec3 b( _left, _top, -_nearPlane); + Vec3 c( _right, _top, -_nearPlane); + Vec3 d( _right, _bottom, -_nearPlane); Vec3 o(0,0,0); p[0].set( o, c, b ); @@ -667,8 +661,8 @@ void Frustum::planes(Plane3 p[6]) p[2].set( Vec3( 0,-1, 0),-_bottom ); p[3].set( Vec3(-1, 0, 0),-_left ); } - p[4].set( Vec3(0, 0, 1), -_near ); - p[5].set( Vec3(0, 0,-1), _far ); + p[4].set( Vec3(0, 0, 1), -_nearPlane ); + p[5].set( Vec3(0, 0,-1), _farPlane ); } @@ -680,20 +674,20 @@ void Frustum::planes(Plane3 p[6], const Matrix44 &M) // Normals point outwards. // - Vec3 a = Vec3( _left, _bottom, -_near) * M; - Vec3 b = Vec3( _left, _top, -_near) * M; - Vec3 c = Vec3( _right, _top, -_near) * M; - Vec3 d = Vec3( _right, _bottom, -_near) * M; + Vec3 a = Vec3( _left, _bottom, -_nearPlane) * M; + Vec3 b = Vec3( _left, _top, -_nearPlane) * M; + Vec3 c = Vec3( _right, _top, -_nearPlane) * M; + Vec3 d = Vec3( _right, _bottom, -_nearPlane) * M; if (! _orthographic) { - double s = _far / double(_near); + double s = _farPlane / double(_nearPlane); T farLeft = (T) (s * _left); T farRight = (T) (s * _right); T farTop = (T) (s * _top); T farBottom = (T) (s * _bottom); - Vec3 e = Vec3( farLeft, farBottom, -_far) * M; - Vec3 f = Vec3( farLeft, farTop, -_far) * M; - Vec3 g = Vec3( farRight, farTop, -_far) * M; + Vec3 e = Vec3( farLeft, farBottom, -_farPlane) * M; + Vec3 f = Vec3( farLeft, farTop, -_farPlane) * M; + Vec3 g = Vec3( farRight, farTop, -_farPlane) * M; Vec3 o = Vec3(0,0,0) * M; p[0].set( o, c, b ); p[1].set( o, d, c ); @@ -704,10 +698,10 @@ void Frustum::planes(Plane3 p[6], const Matrix44 &M) } else { - Vec3 e = Vec3( _left, _bottom, -_far) * M; - Vec3 f = Vec3( _left, _top, -_far) * M; - Vec3 g = Vec3( _right, _top, -_far) * M; - Vec3 h = Vec3( _right, _bottom, -_far) * M; + Vec3 e = Vec3( _left, _bottom, -_farPlane) * M; + Vec3 f = Vec3( _left, _top, -_farPlane) * M; + Vec3 g = Vec3( _right, _top, -_farPlane) * M; + Vec3 h = Vec3( _right, _bottom, -_farPlane) * M; p[0].set( c, g, f ); p[1].set( d, h, g ); p[2].set( a, e, h ); diff --git a/IlmBase/Imath/ImathMatrix.h b/IlmBase/Imath/ImathMatrix.h index 54939a8..7de0408 100644 --- a/IlmBase/Imath/ImathMatrix.h +++ b/IlmBase/Imath/ImathMatrix.h @@ -273,6 +273,25 @@ template class Matrix33 throw (Iex::MathExc); + //------------------------------------------------ + // Calculate the matrix minor of the (r,c) element + //------------------------------------------------ + + T minorOf (const int r, const int c) const; + + //--------------------------------------------------- + // Build a minor using the specified rows and columns + //--------------------------------------------------- + + T fastMinor (const int r0, const int r1, + const int c0, const int c1) const; + + //------------ + // Determinant + //------------ + + T determinant() const; + //----------------------------------------- // Set matrix to rotation by r (in radians) //----------------------------------------- @@ -628,6 +647,25 @@ template class Matrix44 throw (Iex::MathExc); + //------------------------------------------------ + // Calculate the matrix minor of the (r,c) element + //------------------------------------------------ + + T minorOf (const int r, const int c) const; + + //--------------------------------------------------- + // Build a minor using the specified rows and columns + //--------------------------------------------------- + + T fastMinor (const int r0, const int r1, const int r2, + const int c0, const int c1, const int c2) const; + + //------------ + // Determinant + //------------ + + T determinant() const; + //-------------------------------------------------------- // Set matrix to rotation by XYZ euler angles (in radians) //-------------------------------------------------------- @@ -734,6 +772,13 @@ template class Matrix44 template const Matrix44 & shear (const Vec3 &h); + //-------------------------------------------------------- + // Number of the row and column dimensions, since + // Matrix44 is a square matrix. + //-------------------------------------------------------- + + static unsigned int dimensions() {return 4;} + //------------------------------------------------------------ // Shear the matrix by the given factors. The composed matrix @@ -750,14 +795,6 @@ template class Matrix44 const Matrix44 & shear (const Shear6 &h); - //-------------------------------------------------------- - // Number of the row and column dimensions, since - // Matrix44 is a square matrix. - //-------------------------------------------------------- - - static unsigned int dimensions() {return 4;} - - //------------------------------------------------- // Limitations of type T (see also class limits) //------------------------------------------------- @@ -1620,6 +1657,56 @@ Matrix33::inverse (bool singExc) const throw (Iex::MathExc) } } +template +inline T +Matrix33::minorOf (const int r, const int c) const +{ + int copy[9] = {1,1,1, + 1,1,1, + 1,1,1}; + + copy[ c] = 0; + copy[ 3+c] = 0; + copy[ 6+c] = 0; + + int rowOffset = r*3; + copy[rowOffset++] = 0; + copy[rowOffset++] = 0; + copy[rowOffset ] = 0; + + int offset = 0; + T working[4] = { (T)0,(T)0,(T)0,(T)0 }; + T *elem = working; + *elem += copy[0]*x[0][0]; elem += copy[0]; + *elem += copy[1]*x[0][1]; elem += copy[1]; + *elem += copy[2]*x[0][2]; elem += copy[2]; + *elem += copy[3]*x[1][0]; elem += copy[3]; + *elem += copy[4]*x[1][1]; elem += copy[4]; + *elem += copy[5]*x[1][2]; elem += copy[5]; + *elem += copy[6]*x[2][0]; elem += copy[6]; + *elem += copy[7]*x[2][1]; elem += copy[7]; + *elem += copy[8]*x[2][2]; + + return working[0]*working[3] - working[1]*working[2]; +} + +template +inline T +Matrix33::fastMinor( const int r0, const int r1, + const int c0, const int c1) const +{ + return x[r0][c0]*x[r1][c1] - x[r0][c1]*x[r1][c0]; +} + +template +inline T +Matrix33::determinant () const +{ + return x[0][0]*(x[1][1]*x[2][2] - x[1][2]*x[2][1]) + + x[0][1]*(x[1][2]*x[2][0] - x[1][0]*x[2][2]) + + x[0][2]*(x[1][0]*x[2][1] - x[1][1]*x[2][0]); +} + template template const Matrix33 & @@ -2459,9 +2546,9 @@ Matrix44::multiply (const Matrix44 &a, const Matrix44 &b, Matrix44 &c) { - register const T * IMATH_RESTRICT ap = &a.x[0][0]; - register const T * IMATH_RESTRICT bp = &b.x[0][0]; - register T * IMATH_RESTRICT cp = &c.x[0][0]; + register const T * IMATH_RESTRICT ap = &a.x[0][0]; + register const T * IMATH_RESTRICT bp = &b.x[0][0]; + register T * IMATH_RESTRICT cp = &c.x[0][0]; register T a0, a1, a2, a3; @@ -2818,6 +2905,75 @@ Matrix44::inverse (bool singExc) const throw (Iex::MathExc) return s; } +template +inline T +Matrix44::fastMinor( const int r0, const int r1, const int r2, + const int c0, const int c1, const int c2) const +{ + return x[r0][c0] * (x[r1][c1]*x[r2][c2] - x[r1][c2]*x[r2][c1]) + + x[r0][c1] * (x[r1][c2]*x[r2][c0] - x[r1][c0]*x[r2][c2]) + + x[r0][c2] * (x[r1][c0]*x[r2][c1] - x[r1][c1]*x[r2][c0]); +} + +template +inline T +Matrix44::minorOf (const int r, const int c) const +{ + int copy[16] = {1,1,1,1, + 1,1,1,1, + 1,1,1,1, + 1,1,1,1}; + + copy[ c] = 0; + copy[ 4+c] = 0; + copy[ 8+c] = 0; + copy[12+c] = 0; + + int rowOffset = r*4; + copy[rowOffset++] = 0; + copy[rowOffset++] = 0; + copy[rowOffset++] = 0; + copy[rowOffset ] = 0; + + int offset = 0; + Matrix33 working ((T)0,(T)0,(T)0, + (T)0,(T)0,(T)0, + (T)0,(T)0,(T)0); + T *elem = working.getValue(); + *elem += copy[ 0]*x[0][0]; elem += copy[ 0]; + *elem += copy[ 1]*x[0][1]; elem += copy[ 1]; + *elem += copy[ 2]*x[0][2]; elem += copy[ 2]; + *elem += copy[ 3]*x[0][3]; elem += copy[ 3]; + *elem += copy[ 4]*x[1][0]; elem += copy[ 4]; + *elem += copy[ 5]*x[1][1]; elem += copy[ 5]; + *elem += copy[ 6]*x[1][2]; elem += copy[ 6]; + *elem += copy[ 7]*x[1][3]; elem += copy[ 7]; + *elem += copy[ 8]*x[2][0]; elem += copy[ 8]; + *elem += copy[ 9]*x[2][1]; elem += copy[ 9]; + *elem += copy[10]*x[2][2]; elem += copy[10]; + *elem += copy[11]*x[2][3]; elem += copy[11]; + *elem += copy[12]*x[3][0]; elem += copy[12]; + *elem += copy[13]*x[3][1]; elem += copy[13]; + *elem += copy[14]*x[3][2]; elem += copy[14]; + *elem = copy[15]*x[3][3]; + + return working.determinant(); +} + +template +inline T +Matrix44::determinant () const +{ + T sum = (T)0; + + if (x[0][3] != 0.) sum -= x[0][3] * fastMinor(1,2,3,0,1,2); + if (x[1][3] != 0.) sum += x[1][3] * fastMinor(0,2,3,0,1,2); + if (x[2][3] != 0.) sum -= x[2][3] * fastMinor(0,1,3,0,1,2); + if (x[3][3] != 0.) sum += x[3][3] * fastMinor(0,1,2,0,1,2); + + return sum; +} + template template const Matrix44 & diff --git a/IlmBase/Imath/ImathMatrixAlgo.h b/IlmBase/Imath/ImathMatrixAlgo.h index 63c5202..b1e060f 100644 --- a/IlmBase/Imath/ImathMatrixAlgo.h +++ b/IlmBase/Imath/ImathMatrixAlgo.h @@ -148,6 +148,11 @@ template Matrix44 sansScalingAndShear (const Matrix44 &mat, bool exc = true); +template void sansScalingAndShear + (Matrix44 &result, + const Matrix44 &mat, + bool exc = true); + template bool removeScalingAndShear (Matrix44 &mat, bool exc = true); @@ -226,7 +231,7 @@ template Matrix44 rotationMatrixWithUpDir // -// Returns a matrix that rotates the z-axis so that it +// Constructs a matrix that rotates the z-axis so that it // points towards "targetDir". You must also specify // that you want the up vector to be pointing in a // certain direction "upDir". @@ -238,9 +243,49 @@ template Matrix44 rotationMatrixWithUpDir // (b) when any of the given direction vectors have zero length // -template Matrix44 alignZAxisWithTargetDir - (Vec3 targetDir, - Vec3 upDir); +template void alignZAxisWithTargetDir + (Matrix44 &result, + Vec3 targetDir, + Vec3 upDir); + + +// Compute an orthonormal direct frame from : a position, an x axis direction and a normal to the y axis +// If the x axis and normal are perpendicular, then the normal will have the same direction as the z axis. +// Inputs are : +// -the position of the frame +// -the x axis direction of the frame +// -a normal to the y axis of the frame +// Return is the orthonormal frame +template Matrix44 computeLocalFrame( const Vec3& p, + const Vec3& xDir, + const Vec3& normal); + +// Add a translate/rotate/scale offset to an input frame +// and put it in another frame of reference +// Inputs are : +// - input frame +// - translate offset +// - rotate offset in degrees +// - scale offset +// - frame of reference +// Output is the offsetted frame +template Matrix44 addOffset( const Matrix44& inMat, + const Vec3& tOffset, + const Vec3& rOffset, + const Vec3& sOffset, + const Vec3& ref); + +// Compute Translate/Rotate/Scale matrix from matrix A with the Rotate/Scale of Matrix B +// Inputs are : +// -keepRotateA : if true keep rotate from matrix A, use B otherwise +// -keepScaleA : if true keep scale from matrix A, use B otherwise +// -Matrix A +// -Matrix B +// Return Matrix A with tweaked rotation/scale +template Matrix44 computeRSMatrix( bool keepRotateA, + bool keepScaleA, + const Matrix44& A, + const Matrix44& B); //---------------------------------------------------------------------- @@ -299,7 +344,9 @@ template bool checkForZeroScaleInRow const Vec2 &row, bool exc = true); - +template Matrix33 outerProduct + ( const Vec3 &a, + const Vec3 &b); //----------------------------------------------------------------------------- @@ -393,6 +440,18 @@ sansScalingAndShear (const Matrix44 &mat, bool exc) } +template +void +sansScalingAndShear (Matrix44 &result, const Matrix44 &mat, bool exc) +{ + Vec3 scl; + Vec3 shr; + + if (! extractAndRemoveScalingAndShear (result, scl, shr, exc)) + result = mat; +} + + template bool removeScalingAndShear (Matrix44 &mat, bool exc) @@ -543,7 +602,6 @@ extractEulerXYZ (const Matrix44 &mat, Vec3 &rot) k[0], k[1], k[2], 0, 0, 0, 0, 1); - // // Extract the first angle, rot.x. // @@ -555,18 +613,18 @@ extractEulerXYZ (const Matrix44 &mat, Vec3 &rot) // rotation, N, is only around two axes, and gimbal lock // cannot occur. // + Matrix44 N; N.rotate (Vec3 (-rot.x, 0, 0)); + N = N * M; - N = N * M; - + // // Extract the other two angles, rot.y and rot.z, from N. // T cy = Math::sqrt (N[0][0]*N[0][0] + N[0][1]*N[0][1]); rot.y = Math::atan2 (-N[0][2], cy); rot.z = Math::atan2 (-N[1][0], N[1][1]); - } @@ -633,9 +691,9 @@ extractQuat (const Matrix44 &mat) // check the diagonal if (tr > 0.0) { - s = Math::sqrt (tr + 1.0); - quat.r = s / 2.0; - s = 0.5 / s; + s = Math::sqrt (tr + T(1.0)); + quat.r = s / T(2.0); + s = T(0.5) / s; quat.v.x = (mat[1][2] - mat[2][1]) * s; quat.v.y = (mat[2][0] - mat[0][2]) * s; @@ -651,11 +709,11 @@ extractQuat (const Matrix44 &mat) j = nxt[i]; k = nxt[j]; - s = Math::sqrt ((mat[i][i] - (mat[j][j] + mat[k][k])) + 1.0); + s = Math::sqrt ((mat[i][i] - (mat[j][j] + mat[k][k])) + T(1.0)); - q[i] = s * 0.5; - if (s != 0.0) - s = 0.5 / s; + q[i] = s * T(0.5); + if (s != T(0.0)) + s = T(0.5) / s; q[3] = (mat[j][k] - mat[k][j]) * s; q[j] = (mat[i][j] + mat[j][i]) * s; @@ -779,12 +837,13 @@ rotationMatrixWithUpDir (const Vec3 &fromDir, else { - Matrix44 zAxis2FromDir = alignZAxisWithTargetDir - (fromDir, Vec3 (0, 1, 0)); + Matrix44 zAxis2FromDir( Imath::UNINITIALIZED ); + alignZAxisWithTargetDir (zAxis2FromDir, fromDir, Vec3 (0, 1, 0)); Matrix44 fromDir2zAxis = zAxis2FromDir.transposed (); - Matrix44 zAxis2ToDir = alignZAxisWithTargetDir (toDir, upDir); + Matrix44 zAxis2ToDir( Imath::UNINITIALIZED ); + alignZAxisWithTargetDir (zAxis2ToDir, toDir, upDir); return fromDir2zAxis * zAxis2ToDir; } @@ -792,8 +851,8 @@ rotationMatrixWithUpDir (const Vec3 &fromDir, template -Matrix44 -alignZAxisWithTargetDir (Vec3 targetDir, Vec3 upDir) +void +alignZAxisWithTargetDir (Matrix44 &result, Vec3 targetDir, Vec3 upDir) { // // Ensure that the target direction is non-zero. @@ -840,10 +899,136 @@ alignZAxisWithTargetDir (Vec3 targetDir, Vec3 upDir) row[1] = targetUpDir .normalized (); row[2] = targetDir .normalized (); - Matrix44 mat ( row[0][0], row[0][1], row[0][2], 0, - row[1][0], row[1][1], row[1][2], 0, - row[2][0], row[2][1], row[2][2], 0, - 0, 0, 0, 1 ); + result.x[0][0] = row[0][0]; + result.x[0][1] = row[0][1]; + result.x[0][2] = row[0][2]; + result.x[0][3] = (T)0; + + result.x[1][0] = row[1][0]; + result.x[1][1] = row[1][1]; + result.x[1][2] = row[1][2]; + result.x[1][3] = (T)0; + + result.x[2][0] = row[2][0]; + result.x[2][1] = row[2][1]; + result.x[2][2] = row[2][2]; + result.x[2][3] = (T)0; + + result.x[3][0] = (T)0; + result.x[3][1] = (T)0; + result.x[3][2] = (T)0; + result.x[3][3] = (T)1; +} + + +// Compute an orthonormal direct frame from : a position, an x axis direction and a normal to the y axis +// If the x axis and normal are perpendicular, then the normal will have the same direction as the z axis. +// Inputs are : +// -the position of the frame +// -the x axis direction of the frame +// -a normal to the y axis of the frame +// Return is the orthonormal frame +template +Matrix44 +computeLocalFrame( const Vec3& p, + const Vec3& xDir, + const Vec3& normal) +{ + Vec3 _xDir(xDir); + Vec3 x = _xDir.normalize(); + Vec3 y = (normal % x).normalize(); + Vec3 z = (x % y).normalize(); + + Matrix44 L; + L[0][0] = x[0]; + L[0][1] = x[1]; + L[0][2] = x[2]; + L[0][3] = 0.0; + + L[1][0] = y[0]; + L[1][1] = y[1]; + L[1][2] = y[2]; + L[1][3] = 0.0; + + L[2][0] = z[0]; + L[2][1] = z[1]; + L[2][2] = z[2]; + L[2][3] = 0.0; + + L[3][0] = p[0]; + L[3][1] = p[1]; + L[3][2] = p[2]; + L[3][3] = 1.0; + + return L; +} + +// Add a translate/rotate/scale offset to an input frame +// and put it in another frame of reference +// Inputs are : +// - input frame +// - translate offset +// - rotate offset in degrees +// - scale offset +// - frame of reference +// Output is the offsetted frame +template +Matrix44 +addOffset( const Matrix44& inMat, + const Vec3& tOffset, + const Vec3& rOffset, + const Vec3& sOffset, + const Matrix44& ref) +{ + Matrix44 O; + + Vec3 _rOffset(rOffset); + _rOffset *= M_PI / 180.0; + O.rotate (_rOffset); + + O[3][0] = tOffset[0]; + O[3][1] = tOffset[1]; + O[3][2] = tOffset[2]; + + Matrix44 S; + S.scale (sOffset); + + Matrix44 X = S * O * inMat * ref; + + return X; +} + +// Compute Translate/Rotate/Scale matrix from matrix A with the Rotate/Scale of Matrix B +// Inputs are : +// -keepRotateA : if true keep rotate from matrix A, use B otherwise +// -keepScaleA : if true keep scale from matrix A, use B otherwise +// -Matrix A +// -Matrix B +// Return Matrix A with tweaked rotation/scale +template +Matrix44 +computeRSMatrix( bool keepRotateA, + bool keepScaleA, + const Matrix44& A, + const Matrix44& B) +{ + Vec3 as, ah, ar, at; + extractSHRT (A, as, ah, ar, at); + + Vec3 bs, bh, br, bt; + extractSHRT (B, bs, bh, br, bt); + + if (!keepRotateA) + ar = br; + + if (!keepScaleA) + as = bs; + + Matrix44 mat; + mat.makeIdentity(); + mat.translate (at); + mat.rotate (ar); + mat.scale (as); return mat; } @@ -1110,6 +1295,16 @@ checkForZeroScaleInRow (const T& scl, } +template +Matrix33 +outerProduct (const Vec3 &a, const Vec3 &b ) +{ + return Matrix33 (a.x*b.x, a.x*b.y, a.x*b.z, + a.y*b.x, a.y*b.y, a.y*b.z, + a.z*b.x, a.z*b.y, a.z*b.z ); +} + + } // namespace Imath #endif diff --git a/IlmBase/ImathTest/Makefile.am b/IlmBase/ImathTest/Makefile.am index 899e174..e3630e1 100644 --- a/IlmBase/ImathTest/Makefile.am +++ b/IlmBase/ImathTest/Makefile.am @@ -12,7 +12,8 @@ ImathTest_SOURCES = main.cpp testExtractEuler.cpp testExtractSHRT.cpp \ testQuatSetRotation.h testLineAlgo.cpp testLineAlgo.h \ testQuatSlerp.cpp testQuatSlerp.h testQuat.cpp \ testQuat.h testBoxAlgo.cpp testBoxAlgo.h \ - testVec.cpp testVec.h testBox.cpp testBox.h + testVec.cpp testVec.h testBox.cpp testBox.h \ + testMiscMatrixAlgo.cpp testMiscMatrixAlgo.h INCLUDES = -I$(top_builddir) -I$(top_srcdir)/Imath -I$(top_srcdir)/Iex -I$(top_srcdir)/Half \ -I$(top_srcdir)/config diff --git a/IlmBase/ImathTest/main.cpp b/IlmBase/ImathTest/main.cpp index 6188db0..d43cf9d 100644 --- a/IlmBase/ImathTest/main.cpp +++ b/IlmBase/ImathTest/main.cpp @@ -43,6 +43,7 @@ #include #include #include +#include #include #include #include @@ -59,10 +60,11 @@ int main (int argc, char *argv[]) { - TEST (testVec); + TEST (testVec); TEST (testColor); TEST (testShear); TEST (testMatrix); + TEST (testMiscMatrixAlgo); TEST (testRoots); TEST (testFun); TEST (testInvert); @@ -76,6 +78,6 @@ main (int argc, char *argv[]) TEST (testLineAlgo); TEST (testBoxAlgo); TEST (testBox); - + return 0; } diff --git a/IlmBase/ImathTest/testBox.cpp b/IlmBase/ImathTest/testBox.cpp index 1bf3681..3d138f9 100644 --- a/IlmBase/ImathTest/testBox.cpp +++ b/IlmBase/ImathTest/testBox.cpp @@ -184,6 +184,51 @@ testMakeEmpty(const char *type) } } +template +void +testMakeInfinite(const char *type) +{ + cout << " makeInfinite() for type " << type << endl; + + // + // Infinite box + // + { + Imath::Box b; + b.makeInfinite(); + assert(b.min == T(T::baseTypeMin()) && + b.max == T(T::baseTypeMax())); + } + + // + // Non-empty, has volume + // + { + Imath::Box b(T(-1), T(1)); + b.makeInfinite(); + assert(b.min == T(T::baseTypeMin()) && + b.max == T(T::baseTypeMax())); + } + + // + // Non-empty, no volume + // Boxes are: + // 2D: [(0, 0), (0, 1) ] + // 3D: [(0, 0, 0), (0, 0, 1) ] + // 4D: [(0, 0, 0, 0), (0, 0, 0, 1)] + // + { + T min(0); + T max(0); + max[T::dimensions() - 1] = 1; + + Imath::Box b(min, max); + b.makeInfinite(); + assert(b.min == T(T::baseTypeMin()) && + b.max == T(T::baseTypeMax())); + } +} + template void testExtendByPoint(const char *type) @@ -736,6 +781,61 @@ testIsEmpty(const char *type) } +template +void +testIsInfinite(const char *type) +{ + cout << " isInfinite() for type " << type << endl; + + // + // Infinite box. + // + { + Imath::Box b; + b.makeInfinite(); + assert(b.isInfinite()); + } + + // + // Non-empty, has-volume box. + // 2D: [(-2, -4), ( 8, 2) ] + // 3D: [(-2, -4, -6), (12, 8, 2) ] + // 4D: [(-2, -4, -6, -8), (16, 12, 8, 4) ] + // + { + Imath::Box b0(T(-1), T(1)); + assert(!b0.isInfinite()); + + T p0; + T p1; + for (unsigned int i = 0; i < T::dimensions(); i++) + { + p0[i] = -pow(2, i + 1); + p1[i] = pow(2, T::dimensions() - i); + } + Imath::Box b1(p0, p1); + assert(!b1.isInfinite()); + } + + // + // Non-empty, no-volume box. + // Boxes are: + // 2D: [(0, 0), (0, 2) ] + // 3D: [(0, 0, 0), (0, 0, 2) ] + // 4D: [(0, 0, 0, 0), (0, 0, 0, 2)] + // + { + T min(0); + T max = min; + max[T::dimensions() - 1] = 2; + + Imath::Box b(min, max); + + assert(!b.isInfinite()); + } +} + + template void testHasVolume(const char *type) @@ -750,6 +850,15 @@ testHasVolume(const char *type) assert(!b.hasVolume()); } + // + // Infinite box. + // + { + Imath::Box b; + b.makeInfinite(); + assert(b.hasVolume()); + } + // // Non-empty, has-volume box. // 2D: [(-2, -4), ( 8, 2) ] @@ -897,6 +1006,24 @@ testBox() testMakeEmpty("V4f"); testMakeEmpty("V4d"); + // + // makeInfinite() + // + testMakeInfinite("V2s"); + testMakeInfinite("V2i"); + testMakeInfinite("V2f"); + testMakeInfinite("V2d"); + + testMakeInfinite("V3s"); + testMakeInfinite("V3i"); + testMakeInfinite("V3f"); + testMakeInfinite("V3d"); + + testMakeInfinite("V4s"); + testMakeInfinite("V4i"); + testMakeInfinite("V4f"); + testMakeInfinite("V4d"); + // // extendBy() (point) // @@ -1005,6 +1132,24 @@ testBox() testIsEmpty("V4f"); testIsEmpty("V4d"); + // + // isInfinite() + // + testIsInfinite("V2s"); + testIsInfinite("V2i"); + testIsInfinite("V2f"); + testIsInfinite("V2d"); + + testIsInfinite("V3s"); + testIsInfinite("V3i"); + testIsInfinite("V3f"); + testIsInfinite("V3d"); + + testIsInfinite("V4s"); + testIsInfinite("V4i"); + testIsInfinite("V4f"); + testIsInfinite("V4d"); + // // hasVolume() // diff --git a/IlmBase/ImathTest/testBoxAlgo.cpp b/IlmBase/ImathTest/testBoxAlgo.cpp index c82bcd9..1c4b577 100644 --- a/IlmBase/ImathTest/testBoxAlgo.cpp +++ b/IlmBase/ImathTest/testBoxAlgo.cpp @@ -878,12 +878,20 @@ boxMatrixTransform () Box3f b2 = transform (b1, M); Box3f b3 = affineTransform (b1, M); Box3f b4 = transformSimple (b1, M); + Box3f b21; + Box3f b31; + + transform (b1, M, b21); + affineTransform (b1, M, b31); assert (approximatelyEqual (b2.min, b4.min, e)); assert (approximatelyEqual (b2.max, b4.max, e)); assert (approximatelyEqual (b3.max, b4.max, e)); assert (approximatelyEqual (b3.max, b4.max, e)); + assert (b21 == b2); + assert (b31 == b3); + M[0][3] = 1; M[1][3] = 2; M[2][3] = 3; @@ -891,9 +899,14 @@ boxMatrixTransform () Box3f b5 = transform (b1, M); Box3f b6 = transformSimple (b1, M); + Box3f b51; + + transform (b1, M, b51); assert (approximatelyEqual (b5.min, b6.min, e)); assert (approximatelyEqual (b5.max, b6.max, e)); + assert (b51 == b5); + } diff --git a/IlmBase/ImathTest/testExtractEuler.cpp b/IlmBase/ImathTest/testExtractEuler.cpp index 6ca5936..378a676 100644 --- a/IlmBase/ImathTest/testExtractEuler.cpp +++ b/IlmBase/ImathTest/testExtractEuler.cpp @@ -122,7 +122,7 @@ testRandomAngles (M44f (*matrixEulerMatrix)(const M44f &, Eulerf::Order), { Rand48 r(0); - for (int i = 0; i < 100000; ++i) + for (int i = 0; i < 100000; ++i) { // // Create a rotation matrix, M @@ -196,11 +196,11 @@ test (M44f (*matrixEulerMatrix)(const M44f &, Eulerf::Order), { cout << "order = " << setbase (16) << int (order) << setbase (10) << endl; - cout << "random angles" << endl; + // cout << "random angles" << endl; testRandomAngles (matrixEulerMatrix, order); - cout << "special angles" << endl; + // cout << "special angles" << endl; for (int i = 0; i < 360; i += 90) for (int j = 0; j < 360; j += 90) diff --git a/IlmBase/ImathTest/testFrustum.cpp b/IlmBase/ImathTest/testFrustum.cpp index c60d300..9e731ac 100644 --- a/IlmBase/ImathTest/testFrustum.cpp +++ b/IlmBase/ImathTest/testFrustum.cpp @@ -137,7 +137,7 @@ testFrustumPlanes (Imath::Frustumf &frustum) mView.multDirMatrix (front, front); assert ((front ^ planes[4].normal) > 0.0); - pt = Imath::V3f (0.0f, 0.0f, -frustum.near()); + pt = Imath::V3f (0.0f, 0.0f, -frustum.nearPlane()); d = planes0[4].distanceTo (pt); assert (Imath::iszero (d, eps)); pt = pt * mView; @@ -149,7 +149,7 @@ testFrustumPlanes (Imath::Frustumf &frustum) mView.multDirMatrix (back, back); assert ((back ^ planes[5].normal) > 0.0); - pt = Imath::V3f (0.0f, 0.0f, -frustum.far()); + pt = Imath::V3f (0.0f, 0.0f, -frustum.farPlane()); d = planes0[5].distanceTo (pt); assert (Imath::iszero (d, eps)); pt = pt * mView; diff --git a/IlmBase/ImathTest/testMatrix.cpp b/IlmBase/ImathTest/testMatrix.cpp index 4a4f7e2..9d8627a 100644 --- a/IlmBase/ImathTest/testMatrix.cpp +++ b/IlmBase/ImathTest/testMatrix.cpp @@ -36,10 +36,12 @@ #include #include "ImathMatrix.h" +#include "ImathMatrixAlgo.h" #include "ImathVec.h" #include "ImathLimits.h" #include "ImathMath.h" #include "ImathInt64.h" +#include "ImathRandom.h" #include #include @@ -241,32 +243,215 @@ testMatrix () assert(m2[0][0] == (float)m1[0][0]); } - // VC 2005 64 bits compiler has a bug with __restrict keword. - // Pointers with __restrict should not alias the same symbol. - // But, with optimization on, VC removes intermediate temp variable - // and ignores __restrict. - { - cout << "M44 multiplicaftion test" << endl; - Imath::M44f M ( 1.0f, 2.0f, 3.0f, 4.0f, - 5.0f, 6.0f, 7.0f, 8.0f, - 9.0f, 10.0f, 11.0f, 12.0f, - 13.0f, 14.0f, 15.0f, 16.0f); - - Imath::M44f N; N.makeIdentity(); - - // N should be equal to M - // This typical test fails - // when __restrict is used for pointers in "multiply" function. - N = N * M; - - assert(N == M); - - if (N != M) { - cout << "M44 multiplication test has failed, error." << endl - << "M" << endl << M << endl - << "N" << endl << N << endl; - } - } + // Matrix minors + { + cout << "3x3 Matrix minors" << endl; + + Imath::M33f a(1,2,3,4,5,6,7,8,9); + + assert (a.minorOf(0,0) == a.fastMinor(1,2,1,2)); + assert (a.minorOf(0,1) == a.fastMinor(1,2,0,2)); + assert (a.minorOf(0,2) == a.fastMinor(1,2,0,1)); + assert (a.minorOf(1,0) == a.fastMinor(0,2,1,2)); + assert (a.minorOf(1,1) == a.fastMinor(0,2,0,2)); + assert (a.minorOf(1,2) == a.fastMinor(0,2,0,1)); + assert (a.minorOf(2,0) == a.fastMinor(0,1,1,2)); + assert (a.minorOf(2,1) == a.fastMinor(0,1,0,2)); + assert (a.minorOf(2,2) == a.fastMinor(0,1,0,1)); + } + { + Imath::M33d a(1,2,3,4,5,6,7,8,9); + + assert (a.minorOf(0,0) == a.fastMinor(1,2,1,2)); + assert (a.minorOf(0,1) == a.fastMinor(1,2,0,2)); + assert (a.minorOf(0,2) == a.fastMinor(1,2,0,1)); + assert (a.minorOf(1,0) == a.fastMinor(0,2,1,2)); + assert (a.minorOf(1,1) == a.fastMinor(0,2,0,2)); + assert (a.minorOf(1,2) == a.fastMinor(0,2,0,1)); + assert (a.minorOf(2,0) == a.fastMinor(0,1,1,2)); + assert (a.minorOf(2,1) == a.fastMinor(0,1,0,2)); + assert (a.minorOf(2,2) == a.fastMinor(0,1,0,1)); + } + + // Determinants (by building a random singular value decomposition) + { + cout << "3x3 determinant" << endl; + + Imath::Rand32 random; + + Imath::M33f u; + Imath::M33f v; + Imath::M33f s; + + u.setRotation( random.nextf() ); + v.setRotation( random.nextf() ); + s[0][0] = random.nextf(); + s[1][1] = random.nextf(); + s[2][2] = random.nextf(); + + Imath::M33f c = u * s * v.transpose(); + assert (fabsf(c.determinant() - s[0][0]*s[1][1]*s[2][2]) <= u.baseTypeEpsilon()); + } + { + Imath::Rand32 random; + + Imath::M33d u; + Imath::M33d v; + Imath::M33d s; + + u.setRotation( (double)random.nextf() ); + v.setRotation( (double)random.nextf() ); + s[0][0] = (double)random.nextf(); + s[1][1] = (double)random.nextf(); + s[2][2] = (double)random.nextf(); + + Imath::M33d c = u * s * v.transpose(); + assert (fabs(c.determinant() - s[0][0]*s[1][1]*s[2][2]) <= u.baseTypeEpsilon()); + } + + // Outer product of two 3D vectors + { + cout << "Outer product of two 3D vectors" << endl; + + Imath::V3f a(1,2,3); + Imath::V3f b(4,5,6); + Imath::M33f p = Imath::outerProduct(a,b); + + for (int i=0; i<3; i++ ) + { + for (int j=0; j<3; j++) + { + assert (p[i][j] == a[i]*b[j]); + } + } + } + { + Imath::V3d a(1,2,3); + Imath::V3d b(4,5,6); + Imath::M33d p = Imath::outerProduct(a,b); + + for (int i=0; i<3; i++ ) + { + for (int j=0; j<3; j++) + { + assert (p[i][j] == a[i]*b[j]); + } + } + } + + + // Determinants (by building a random singular value decomposition) + { + cout << "4x4 determinants" << endl; + + Imath::Rand32 random; + + Imath::M44f u = Imath::rotationMatrix + ( Imath::V3f(random.nextf(),random.nextf(),random.nextf()).normalize(), + Imath::V3f(random.nextf(),random.nextf(),random.nextf()).normalize() ); + Imath::M44f v = Imath::rotationMatrix + ( Imath::V3f(random.nextf(),random.nextf(),random.nextf()).normalize(), + Imath::V3f(random.nextf(),random.nextf(),random.nextf()).normalize() ); + Imath::M44f s; + + s[0][0] = random.nextf(); + s[1][1] = random.nextf(); + s[2][2] = random.nextf(); + s[3][3] = random.nextf(); + + Imath::M44f c = u * s * v.transpose(); + assert (fabsf(c.determinant() - s[0][0]*s[1][1]*s[2][2]*s[3][3]) <= u.baseTypeEpsilon()); + } + { + Imath::Rand32 random; + + Imath::M44d u = Imath::rotationMatrix + ( Imath::V3d(random.nextf(),random.nextf(),random.nextf()).normalize(), + Imath::V3d(random.nextf(),random.nextf(),random.nextf()).normalize() ); + Imath::M44d v = Imath::rotationMatrix + ( Imath::V3d(random.nextf(),random.nextf(),random.nextf()).normalize(), + Imath::V3d(random.nextf(),random.nextf(),random.nextf()).normalize() ); + Imath::M44d s; + + s[0][0] = random.nextf(); + s[1][1] = random.nextf(); + s[2][2] = random.nextf(); + s[3][3] = random.nextf(); + + Imath::M44d c = u * s * v.transpose(); + assert (fabs(c.determinant() - s[0][0]*s[1][1]*s[2][2]*s[3][3]) <= u.baseTypeEpsilon()); + } + + // Matrix minors + { + cout << "4x4 matrix minors" << endl; + + Imath::M44d a(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16); + + assert (a.minorOf(0,0) == a.fastMinor(1,2,3,1,2,3)); + assert (a.minorOf(0,1) == a.fastMinor(1,2,3,0,2,3)); + assert (a.minorOf(0,2) == a.fastMinor(1,2,3,0,1,3)); + assert (a.minorOf(0,3) == a.fastMinor(1,2,3,0,1,2)); + assert (a.minorOf(1,0) == a.fastMinor(0,2,3,1,2,3)); + assert (a.minorOf(1,1) == a.fastMinor(0,2,3,0,2,3)); + assert (a.minorOf(1,2) == a.fastMinor(0,2,3,0,1,3)); + assert (a.minorOf(1,3) == a.fastMinor(0,2,3,0,1,2)); + assert (a.minorOf(2,0) == a.fastMinor(0,1,3,1,2,3)); + assert (a.minorOf(2,1) == a.fastMinor(0,1,3,0,2,3)); + assert (a.minorOf(2,2) == a.fastMinor(0,1,3,0,1,3)); + assert (a.minorOf(2,3) == a.fastMinor(0,1,3,0,1,2)); + assert (a.minorOf(3,0) == a.fastMinor(0,1,2,1,2,3)); + assert (a.minorOf(3,1) == a.fastMinor(0,1,2,0,2,3)); + assert (a.minorOf(3,2) == a.fastMinor(0,1,2,0,1,3)); + assert (a.minorOf(3,3) == a.fastMinor(0,1,2,0,1,2)); + } + { + Imath::M44f a(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16); + + assert (a.minorOf(0,0) == a.fastMinor(1,2,3,1,2,3)); + assert (a.minorOf(0,1) == a.fastMinor(1,2,3,0,2,3)); + assert (a.minorOf(0,2) == a.fastMinor(1,2,3,0,1,3)); + assert (a.minorOf(0,3) == a.fastMinor(1,2,3,0,1,2)); + assert (a.minorOf(1,0) == a.fastMinor(0,2,3,1,2,3)); + assert (a.minorOf(1,1) == a.fastMinor(0,2,3,0,2,3)); + assert (a.minorOf(1,2) == a.fastMinor(0,2,3,0,1,3)); + assert (a.minorOf(1,3) == a.fastMinor(0,2,3,0,1,2)); + assert (a.minorOf(2,0) == a.fastMinor(0,1,3,1,2,3)); + assert (a.minorOf(2,1) == a.fastMinor(0,1,3,0,2,3)); + assert (a.minorOf(2,2) == a.fastMinor(0,1,3,0,1,3)); + assert (a.minorOf(2,3) == a.fastMinor(0,1,3,0,1,2)); + assert (a.minorOf(3,0) == a.fastMinor(0,1,2,1,2,3)); + assert (a.minorOf(3,1) == a.fastMinor(0,1,2,0,2,3)); + assert (a.minorOf(3,2) == a.fastMinor(0,1,2,0,1,3)); + assert (a.minorOf(3,3) == a.fastMinor(0,1,2,0,1,2)); + } + + // VC 2005 64 bits compiler has a bug with __restrict keword. + // Pointers with __restrict should not alias the same symbol. + // But, with optimization on, VC removes intermediate temp variable + // and ignores __restrict. + { + cout << "M44 multiplicaftion test" << endl; + Imath::M44f M ( 1.0f, 2.0f, 3.0f, 4.0f, + 5.0f, 6.0f, 7.0f, 8.0f, + 9.0f, 10.0f, 11.0f, 12.0f, + 13.0f, 14.0f, 15.0f, 16.0f); + + Imath::M44f N; N.makeIdentity(); + + // N should be equal to M + // This typical test fails + // when __restrict is used for pointers in "multiply" function. + N = N * M; + + assert(N == M); + + if (N != M) { + cout << "M44 multiplication test has failed, error." << endl + << "M" << endl << M << endl + << "N" << endl << N << endl; + } + } cout << "ok\n" << endl; } diff --git a/IlmBase/ImathTest/testMiscMatrixAlgo.cpp b/IlmBase/ImathTest/testMiscMatrixAlgo.cpp new file mode 100644 index 0000000..d0f1863 --- /dev/null +++ b/IlmBase/ImathTest/testMiscMatrixAlgo.cpp @@ -0,0 +1,341 @@ +/////////////////////////////////////////////////////////////////////////// +// +// Copyright (c) 2009, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +/////////////////////////////////////////////////////////////////////////// + + +#include +#include "ImathMatrixAlgo.h" +#include "ImathRandom.h" +#include +#include +#include +#include + + +#if 0 + #define debug(x) (printf x, fflush (stdout)) +#else + #define debug(x) +#endif + + +using namespace std; +using namespace Imath; + +namespace { + +float rad (float deg) {return deg * (M_PI / 180);} + +void +testComputeLocalFrame () +{ + float eps = 0.00005; + Rand48 random(0); + for (int i = 0; i < 100000; ++i) + { + debug (("iteration: %d\n", i)); + + // Random pos + V3f p (random.nextf (-10, 10), + random.nextf (-10, 10), + random.nextf (-10, 10)); + + // Random xDir + V3f xDir (random.nextf (-10, 10), + random.nextf (-10, 10), + random.nextf (-10, 10)); + + // Random normalDir + V3f normalDir (random.nextf (-10, 10), + random.nextf (-10, 10), + random.nextf (-10, 10)); + + // Run computeLocalFrame we want to test + M44f L = computeLocalFrame(p, xDir, normalDir); + + // test position + for (int j=0; j<3; j++) + { + if ( abs(L[3][j] - p[j])>eps ) + assert(false); + } + if (abs (L[3][3] - 1.0)>eps) + assert(false ); + + // check that xAxis has the same dir as xDir and that is is normalized + V3f x( L[0][0], L[0][1], L[0][2]); + assert( (x%xDir).length() < eps ); + if (abs (L[0][3])>eps) + assert(false); + assert((abs(x.length()-1.f)eps) + assert(false ); + assert(abs(x^y)eps) + assert(false ); + assert((abs(z.length()-1.f)0); + } +} + +void +getRandTRS(Rand48& random, V3f& trans, V3f& rot, V3f& scale) +{ + // Translate + trans = V3f (random.nextf (-10, 10), + random.nextf (-10, 10), + random.nextf (-10, 10)); + // Rotate + rot = V3f (rad (random.nextf (-180, 180)), + rad (random.nextf (-180, 180)), + rad (random.nextf (-180, 180))); + + // Scale + V3f s(random.nextf (0.000001, 2.0), + random.nextf (0.000001, 2.0), + random.nextf (0.000001, 2.0)); + for (int j=0; j < 3; j++) + if (random.nextf (0.0, 1.0) >= 0.5) + s[j] *= -1; + scale = s; +} + +M44f +createRandomMat(Rand48& random, V3f& trans, V3f& rot, V3f& scale) +{ + + M44f M; + V3f t, r, s; + getRandTRS(random, t, r, s); + + M.translate (t); + M.rotate (r); + + // Shear M. + V3f h (random.nextf (0.000001, 2.0), + random.nextf (0.000001, 2.0), + random.nextf (0.000001, 2.0)); + + for (int j=0; j < 3; j++) + if (random.nextf (0.0, 1.0) >= 0.5) + h[j] *= -1; + M.shear (h); + + M.scale (s); + + // + // Add a small random error to the elements of M + // + for (int j = 0; j < 4; ++j) + for (int k = 0; k < 3; ++k) + M[j][k] += random.nextf (-1e-7, 1e-7); + + V3f sh; + extractSHRT (M, scale, sh, rot, trans); + + debug (("Scale : %f %f %f\n", s[0], s[1], s[2])); + debug (("Shear : %f %f %f\n", h[0], h[1], h[2])); + debug (("Rot : %f %f %f\n", r[0], r[1], r[2])); + debug (("Trans : %f %f %f\n", t[0], t[1], t[2])); + + return M; +} + +void +compareMat(M44f& M, M44f& N) +{ + float eps = 0.0001; + + /// Verify that the entries in M and N do not + // differ too much. + + M44f D (M - N); + + for (int j = 0; j < 4; ++j) + { + for (int k = 0; k < 4; ++k) + { + //cout << "diff="< eps) + { + cout << "unexpectedly diff "<< + D[j][k] << endl; + + cout << j << " " << k << endl; + + cout << "M\n" << M << endl; + cout << "N\n" << N << endl; + cout << "D\n" << D << endl; + + assert (false); + } + } + } +} + +void +testAddOffset() +{ + Rand48 random(0); + + for (int i = 0; i < 100000; ++i) + { + debug (("iteration: %d\n", i)); + + V3f transA, transB, rotA, rotB, scaleA, scaleB; + V3f tOffset, rOffset, sOffset; + M44f inMat = createRandomMat(random, transA, rotA, scaleA); + M44f refMat = createRandomMat(random, transB, rotB, scaleB); + getRandTRS(random, tOffset, rOffset, sOffset); + + // addOffset : function to test + M44f outMat = addOffset( inMat, tOffset, rOffset, sOffset, refMat); + + // add the inverse offset + M44f invO; + invO.rotate (V3f(rad(rOffset[0]), rad(rOffset[1]), rad(rOffset[2]))); + invO[3][0] = tOffset[0]; + invO[3][1] = tOffset[1]; + invO[3][2] = tOffset[2]; + invO.invert(); + + M44f invS; + invS.scale (sOffset); + invS.invert(); // zero scale is avoided in getRandTRS + + // in ref mat from the function result + M44f outInRefMat = invO*invS*outMat; + + // in ref mat from the inputs + M44f inRefMat = inMat*refMat; + + // compare the mat + compareMat(outInRefMat, inRefMat); + } +} + +void +testRSMatrix(M44f& M, V3f& t, V3f& r, V3f& s) +{ + M44f N; + N.makeIdentity(); + N.translate (t); // ... matrix compositions + N.rotate (r); + N.scale (s); + + compareMat(M, N); +} + +void +testComputeRSMatrix () +{ + Rand48 random(0); + + for (int i = 0; i < 100000; ++i) + { + debug (("iteration: %d\n", i)); + + V3f transA, transB, rotA, rotB, scaleA, scaleB; + + M44f A = createRandomMat(random, transA, rotA, scaleA); + M44f B = createRandomMat(random, transB, rotB, scaleB); + + M44f ArAsA = computeRSMatrix( true, true, A, B); + M44f ArBsB = computeRSMatrix( false, false, A, B); + M44f ArAsB = computeRSMatrix( true, false, A, B); + M44f ArBsA = computeRSMatrix( false, true, A, B); + + testRSMatrix(ArAsA, transA, rotA, scaleA); + testRSMatrix(ArBsB, transA, rotB, scaleB); + testRSMatrix(ArAsB, transA, rotA, scaleB); + testRSMatrix(ArBsA, transA, rotB, scaleA); + + debug (("\n")); + } +} + +} // namespace + + +void +testMiscMatrixAlgo () +{ + try + { + cout << "Testing misc functions in ImathMatrixAlgo.h" << endl; + + cout << "Testing the building of an orthonormal direct frame from : a position, " + << "an x axis direction and a normal to the y axis" << endl; + cout << "Imath::computeLocalFrame()" << endl; + + testComputeLocalFrame (); + + cout << "ok\n" << endl; + + cout << "Add a translate/rotate/scale offset to an input frame " + << "and put it in another frame of reference" << endl; + cout << "Imath::addOffset()" << endl; + + testAddOffset (); + + cout << "ok\n" << endl; + + cout << "Compute Translate/Rotate/Scale matrix from matrix A "<