I found a bug when we do matrix inverse calculation by vtkTransform object. The version 8.2.0 and 9.1.0 of VTK has the same potential risk in vtkMatrix4x4::Invert(const double inElements[16], double outElements[16])
.
Let’s see the source code, it judges whether the determinant value is zero by ==
. The equal comparison operation for float data type in C++ is not reliable. It can be dangerous if the zero value cheats the judgment expression. So my program go in an infinite loop after it.
//------------------------------------------------------------------------------
// Matrix Inversion (adapted from Richard Carling in "Graphics Gems,"
// Academic Press, 1990).
void vtkMatrix4x4::Invert(const double inElements[16], double outElements[16])
{
// inverse( original_matrix, inverse_matrix )
// calculate the inverse of a 4x4 matrix
//
// -1
// A = ___1__ adjoint A
// det A
//
// calculate the 4x4 determinent
// if the determinent is zero,
// then the inverse matrix is not unique.
double det = vtkMatrix4x4::Determinant(inElements);
if (det == 0.0)
{
return;
}
// calculate the adjoint matrix
vtkMatrix4x4::Adjoint(inElements, outElements);
// scale the adjoint matrix to get the inverse
for (int i = 0; i < 16; i++)
{
outElements[i] /= det;
}
}
Rewrite it to fix the potential bug.
bool Invert(const vtkSmartPointer<vtkMatrix4x4> inMatrix, double outElements[])
{
// calculate the inverse of a 4x4 matrix
//
// -1
// A = ___1__ adjoint A
// det A
//
double det = inMatrix->Determinant();
if( fabs( det ) < 1e-5 )
{
LOG( ERR, "jump: zero det!" );
return false;
}
double inElements[16];
for( int i = 0; i < 4; ++i )
{
for( int j = 0; j < 4; ++j )
{
inElements[i*4+j] = inMatrix->GetElement( i, j );
}
}
// calculate the adjoint matrix
vtkMatrix4x4::Adjoint(inElements, outElements);
// scale the adjoint matrix to get the inverse
for (int i = 0; i < 16; i++)
{
outElements[i] /= det;
}
return true;
}