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,"
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;
}

// 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 );
}
}

// scale the adjoint matrix to get the inverse
for (int i = 0; i < 16; i++)
{
outElements[i] /= det;
}
return true;
}
``````
Categories: MathVTK

Article Rating
Subscribe
Notify of