FGFS-4.1/OpenSceneGraph
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Various investigations into culling errors w.r.t matrix inversion resulted in

Browse Source 
the conclusion that the osg::Matrix::inverse was broken, have lifted a new
implementation from sgl and it seems to work fine.  Will need further testing
but looks good.
This commit is contained in:
Robert Osfield
2002-02-07 01:15:15 +00:00
parent 7c7f195d6e
commit d12a726d5b
8 changed files with 177 additions and 50 deletions
Download Patch File Download Diff File Expand all files Collapse all files

160
src/osg/Matrix.cpp
View File

@@ -21,6 +21,105 @@ using namespace osg;
+((a)._mat[r][3] * (b)._mat[3][c])
template <class T>
inline T SGL_ABS(T a)
{
return (a >= 0 ? a : -a);
}
#ifndef SGL_SWAP
#define SGL_SWAP(a,b,temp) ((temp)=(a),(a)=(b),(b)=(temp))
#endif
bool inverse(const Matrix& mat,Matrix& m_matrix)
{
unsigned int indxc[4], indxr[4], ipiv[4];
unsigned int i,j,k,l,ll;
unsigned int icol = 0;
unsigned int irow = 0;
float temp, pivinv, dum, big;
// copy in place this may be unnecessary
m_matrix = mat;
for (j=0; j<4; j++) ipiv[j]=0;
for(i=0;i<4;i++)
{
big=(float)0.0;
for (j=0; j<4; j++)
if (ipiv[j] != 1)
for (k=0; k<4; k++)
{
if (ipiv[k] == 0)
{
if (SGL_ABS(m_matrix(j,k)) >= big)
{
big = SGL_ABS(m_matrix(j,k));
irow=j;
icol=k;
}
}
else if (ipiv[k] > 1)
return false;
}
++(ipiv[icol]);
if (irow != icol)
for (l=0; l<4; l++) SGL_SWAP(m_matrix(irow,l),
m_matrix(icol,l),
temp);
indxr[i]=irow;
indxc[i]=icol;
if (m_matrix(icol,icol) == 0)
return false;
pivinv = 1.0/m_matrix(icol,icol);
m_matrix(icol,icol) = 1;
for (l=0; l<4; l++) m_matrix(icol,l) *= pivinv;
for (ll=0; ll<4; ll++)
if (ll != icol)
{
dum=m_matrix(ll,icol);
m_matrix(ll,icol) = 0;
for (l=0; l<4; l++) m_matrix(ll,l) -= m_matrix(icol,l)*dum;
}
}
for (int lx=4; lx>0; --lx)
{
if (indxr[lx-1] != indxc[lx-1])
for (k=0; k<4; k++) SGL_SWAP(m_matrix(k,indxr[lx-1]),
m_matrix(k,indxc[lx-1]),temp);
}
return true;
}
bool inverseAffine(const Matrix& mat,Matrix& m_matrix)
{
// | R p |' | R' -R'p |'
// | | -> | |
// | 0 0 0 1 | | 0 0 0 1 |
for (unsigned int i=0; i<3; i++)
{
m_matrix(i,3) = 0;
m_matrix(3,i) = -(mat(i,0)*mat(3,0) +
mat(i,1)*mat(3,1) +
mat(i,2)*mat(3,2));
for (unsigned int j=0; j<3; j++)
{
m_matrix(i,j) = mat(j,i);
}
}
m_matrix(3,3) = 1;
return true;
}
Matrix::Matrix() : Object()
{
makeIdentity();
@@ -47,7 +146,6 @@ Matrix::Matrix( float a00, float a01, float a02, float a03,
SET_ROW(3, a30, a31, a32, a33 )
}
Matrix& Matrix::operator = (const Matrix& other )
{
if( &other == this ) return *this;
@@ -169,6 +267,15 @@ void Matrix::makeRotate( float yaw, float pitch, float roll)
void Matrix::mult( const Matrix& lhs, const Matrix& rhs )
{
if (&lhs==this)
{
postMult(rhs);
return;
}
if (&rhs==this)
{
preMult(lhs);
}
// PRECONDITION: We assume neither &lhs nor &rhs == this
// if it did, use preMult or postMult instead
@@ -232,31 +339,33 @@ void Matrix::postMult( const Matrix& other )
#undef SET_ROW
#undef INNER_PRODUCT
bool Matrix::invert( const Matrix& other )
bool Matrix::invert( const Matrix& invm )
{
if (&other==this)
{
Matrix tm(other);
if (&invm==this) {
Matrix tm(invm);
return invert(tm);
}
// if ( other._mat[0][3] == 0.0
// && other._mat[1][3] == 0.0
// && other._mat[2][3] == 0.0
// && other._mat[3][3] == 1.0 )
// {
// return invertAffine( other );
// }
/*
if ( invm._mat[0][3] == 0.0
&& invm._mat[1][3] == 0.0
&& invm._mat[2][3] == 0.0
&& invm._mat[3][3] == 1.0 )
{
return inverseAffine( invm,*this );
}
*/
return inverse(invm,*this);
// code lifted from VR Juggler.
// not cleanly added, but seems to work. RO.
const float* a = reinterpret_cast<const float*>(other._mat);
const float* a = reinterpret_cast<const float*>(invm._mat);
float* b = reinterpret_cast<float*>(_mat);
int n = 4;
int i, j, k;
int r[ 4], c[ 4], row[ 4], col[ 4];
float m[ 4][ 4*2], pivot, maxother, tmpother, fac;
float m[ 4][ 4*2], pivot, max_m, tmp_m, fac;
/* Initialization */
for ( i = 0; i < n; i ++ )
@@ -265,7 +374,7 @@ bool Matrix::invert( const Matrix& other )
row[ i] = col[ i] = 0;
}
/* Set working Matrix */
/* Set working matrix */
for ( i = 0; i < n; i++ )
{
for ( j = 0; j < n; j++ )
@@ -279,16 +388,16 @@ bool Matrix::invert( const Matrix& other )
for ( k = 0; k < n; k++ )
{
/* Choosing the pivot */
for ( i = 0, maxother = 0; i < n; i++ )
for ( i = 0, max_m = 0; i < n; i++ )
{
if ( row[ i] ) continue;
for ( j = 0; j < n; j++ )
{
if ( col[ j] ) continue;
tmpother = fabs( m[ i][j]);
if ( tmpother > maxother)
tmp_m = fabs( m[ i][j]);
if ( tmp_m > max_m)
{
maxother = tmpother;
max_m = tmp_m;
r[ k] = i;
c[ k] = j;
}
@@ -299,9 +408,9 @@ bool Matrix::invert( const Matrix& other )
if ( fabs( pivot) <= 1e-20)
{
notify(WARN) << "Warning: pivot = "<< pivot <<" in Matrix::invert(), cannot compute inverse."<<std::endl;
notify(WARN) << "input matrix to Matrix::invert() was = "<< other <<std::endl;
makeIdentity();
notify(INFO) << "Warning: pivot = "<< pivot <<" in Matrix::invert(), cannot compute inverse."<<std::endl;
notify(INFO) << "input matrix to Matrix::invert() was = "<< invm <<std::endl;
//makeIdentity();
return false;
}
@@ -330,7 +439,7 @@ bool Matrix::invert( const Matrix& other )
}
}
/* Assign invers to a Matrix */
/* Assign invers to a matrix */
for ( i = 0; i < n; i++ )
for ( j = 0; j < n; j++ )
row[ i] = ( c[ j] == i ) ? r[j] : row[ i];
@@ -340,6 +449,8 @@ bool Matrix::invert( const Matrix& other )
b[ i * n + j] = m[ row[ i]][j + n];
return true; // It worked
}
const double PRECISION_LIMIT = 1.0e-15;
@@ -412,3 +523,4 @@ bool Matrix::invertAffine( const Matrix& other )
return true;
}