------------------------------------------------------------------------------
| 45 | |
| 46 | //------------------------------------------------------------------------------ |
| 47 | vtkSmartPointer<vtkMatrix3x3> PolarDecomposition(vtkMatrix3x3* M) |
| 48 | { |
| 49 | // Compute the SVD of the input |
| 50 | double arrayM[3][3]; |
| 51 | for (int i = 0; i < 3; ++i) |
| 52 | { |
| 53 | for (int j = 0; j < 3; ++j) |
| 54 | { |
| 55 | arrayM[i][j] = M->GetElement(i, j); |
| 56 | } |
| 57 | } |
| 58 | double U[3][3], Vt[3][3], S[3]; |
| 59 | vtkMath::SingularValueDecomposition3x3(arrayM, U, S, Vt); |
| 60 | |
| 61 | // Estimate the scale factor by averaging the singular values |
| 62 | double scale = (S[0] + S[1] + S[2]) / 3.0; |
| 63 | |
| 64 | // Determine the polar matrix |
| 65 | vtkNew<vtkMatrix3x3> polar; |
| 66 | vtkNew<vtkMatrix3x3> matrixU; |
| 67 | vtkNew<vtkMatrix3x3> matrixVt; |
| 68 | matrixU->SetData(U[0]); |
| 69 | matrixVt->SetData(Vt[0]); |
| 70 | vtkMatrix3x3::Multiply3x3(matrixU, matrixVt, polar); |
| 71 | for (int i = 0; i < 3; ++i) |
| 72 | { |
| 73 | for (int j = 0; j < 3; ++j) |
| 74 | { |
| 75 | polar->SetElement(i, j, scale * polar->GetElement(i, j)); |
| 76 | } |
| 77 | } |
| 78 | |
| 79 | return polar; |
| 80 | } |
| 81 | |
| 82 | } // anonymous namespace |
| 83 |
no test coverage detected