Use Singular Value Decomposition In VTK

Published by StephenWei on

Affine transformation can be decomposed to translation, rotation, and scale.
As we know, the transformation can be represented by matrix.

    \[M_1 = \begin{pmatrix} X_0 & Y_0 & Z_0 & P_0 \\ X_1 & Y_1 & Z_1 & P_1 \\ X_2 & Y_2 & Z_2 & P_2 \\ 0 & 0 & 0 & 1 \end{pmatrix}\]

Related post: https://www.weiy.city/2021/11/the-releationship-between-local-transform-and-pose-transform/

It can be computed by translate matrix and rotate & scale matrix.

    \[\begin{pmatrix} X_0 & Y_0 & Z_0 & P_0 \\ X_1 & Y_1 & Z_1 & P_1 \\ X_2 & Y_2 & Z_2 & P_2 \\ 0 & 0 & 0 & 1 \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 & P_0 \\ 0 & 1 & 0 & P_1 \\ 0 & 0 & 1 & P_2 \\ 0 & 0 & 0 & 1 \end{pmatrix} \times \begin{pmatrix} X_0 & Y_0 & Z_0 & 0 \\ X_1 & Y_1 & Z_1 & 0 \\ X_2 & Y_2 & Z_2 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}\]

Example: I try to do decomposition for a matrix.

    \[\begin{pmatrix} 2 & -2 & 0 & 1 \\ 0 & 0 & -1 & 1 \\ 2 & 2 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \times \begin{pmatrix} 2 & -2 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 2 & 2 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}\]

The matrix comes from the three basic linear transforms, scale, rotate and translate.

#include <vtkPointData.h>
#include <vtkSmartPointer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkPolyData.h>
#include <vtkProperty.h>
#include <vtkSphereSource.h>
#include <vtkPolyDataMapper.h>
#include <vtkCommand.h>
#include <vtkSliderWidget.h>
#include <vtkSliderRepresentation.h>
#include <vtkTransform.h>
#include <vtkSliderRepresentation3D.h>
#include <vtkWidgetEventTranslator.h>
#include <vtkWidgetEvent.h>
#include <vtkCubeSource.h>
#include <vtkTransformFilter.h>
#include <vtkAxesActor.h>
#include "./tool.h"

#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();

int main(int, char *[])
{
    vtkSPtrNew( source, vtkCubeSource );
    source->Update();

    vtkSPtrNew( scaleM, vtkTransform );
    scaleM->Scale( 2, 2, 1 );
    scaleM->Update();

    double elements[16] = { 1, -1, 0, 0,
                            0, 0, -1, 0,
                            1, 1, 0, 0,
                            0, 0, 0, 1 };
    vtkSPtrNew( rotateM, vtkTransform );
    rotateM->SetMatrix( elements );
    rotateM->Update();

    vtkSPtrNew( moveM, vtkTransform );
    moveM->Translate( 1, 1, 0 );
    moveM->Update();

    vtkSPtrNew( finalM, vtkTransform );
    finalM->Concatenate( moveM );
    finalM->Concatenate( rotateM );
    finalM->Concatenate( scaleM );
    finalM->Update();
    finalM->PrintSelf( std::cout, vtkIndent() );

    vtkSPtrNew( axes, vtkAxesActor );
    axes->SetTotalLength( 1, 1, 1 ); // change length of three axis
    axes->SetUserTransform( finalM );

    double *scaleFactor = finalM->GetScale();
    cout << "scaleFactor: " << scaleFactor[0] << ", " << scaleFactor[1] << ", " << scaleFactor[2] << endl;

    vtkSPtrNew( renderer, vtkRenderer );
    renderer->AddActor( axes );
    renderer->SetBackground( 0, 0, 0 );

    vtkSPtrNew( renderWindow, vtkRenderWindow );
    renderWindow->AddRenderer( renderer );

    vtkSPtrNew( renderWindowInteractor, vtkRenderWindowInteractor );
    renderWindowInteractor->SetRenderWindow( renderWindow );

    renderer->ResetCamera();
    renderWindow->Render();
    renderWindowInteractor->Start();
    return EXIT_SUCCESS;
}

Define variable M =

    \[\begin{pmatrix} 2 & -2 & 0 \\ 0 & 0 & -1 \\ 2 & 2 & 0 \end{pmatrix}\]

eigenvalue decomposition

We have to review the eigenvalue decomposition.

If there is a constant value \lambda and non-zero column vector \vec {x} has the relationship AX = \lambda X, we can get the new representation A = W \Sigma W^{-1}.
The eigenvalues are going to fill the diagonal, and the eigenvectors are going to form W.

    \[AX = \lambda X\]

    \[\lambda E X - A X = (\lambda E - A)X =0\]

When the eigenvalues polynomial is equal to 0, it’s called the eigenvalues equation of A.

    \[|\lambda E - A| = 0\]

For example, A =

    \[\begin{pmatrix} 1 & -2 \\ 1 & 4 \end{pmatrix}\]

    \[|\lambda E - A| = \begin{vmatrix} \lambda - 1 & 2 \\ -1 & \lambda - 4 \end{vmatrix} = (\lambda - 1)(\lambda - 4) +2 = 0\]

Solve it, we get: \lambda_1 = 2, \lambda_2 = 3.

When \lambda = 2, we have:

    \[(\lambda E - A)X = \begin{pmatrix} 1 & 2 \\ -1 & -2 \end{pmatrix} \times \begin{pmatrix} x \\ y \end{pmatrix} = 0\]

Set X=1, we get Y = -0.5. So an eigenvector is (1, -0.5).
In the case \lambda = 3, we get the eigenvector(1, -1).

So we have :

    \[A = \begin{pmatrix} 1 & -2 \\ 1 & 4 \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ -0.5 & -1 \end{pmatrix} \times \begin{pmatrix} 2 & 0 \\ 0 & 3 \end{pmatrix} \times \begin{pmatrix} 2 & 2 \\ -1 & -2 \end{pmatrix} = W \Sigma W^{-1}\]

singular value decomposition

A more popular decomposition way for matrix is singular value decomposition, SVD. It says, A = U \times W \times V^T. It can be proved by eigenvalue decomposition. SVD does not require the matrix to be decomposed as a square matrix. Let’s say our matrix A is an m x n matrix.

We can use SVD algorithm to decomposite the example’s matrix M.

// decomposition.
    double factors[3];
    auto matrix = finalM->GetMatrix();
    double U[3][3], VT[3][3];
    for( int i = 0; i < 3; ++i )
    {
        U[0][i] = matrix->GetElement( 0, i );
        U[1][i] = matrix->GetElement( 1, i );
        U[2][i] = matrix->GetElement( 2, i );
    }
    vtkMath::SingularValueDecomposition3x3( U, U, factors, VT );
    vtkSPtrNew( matrix_U, vtkMatrix3x3 );
    matrix_U->DeepCopy( (const double *)U );
    matrix_U->PrintSelf( std::cout, vtkIndent() );
    /*
    0.707107    -0.707107   0
    0   0   -1
    0.707107    0.707107    0
*/

    vtkSPtrNew( matrix_W, vtkMatrix3x3 );
    matrix_W->Identity();
    for( int i = 0; i < 3; ++i )
        matrix_W->SetElement( i, i, factors[i] );
    matrix_W->PrintSelf( std::cout, vtkIndent() );
    /*
    2.82843 0   0
    0   2.82843 0
    0   0   1
*/

    vtkSPtrNew( matrix_VT, vtkMatrix3x3 );
    matrix_VT->DeepCopy( (const double *)VT );
    matrix_VT->PrintSelf( std::cout, vtkIndent() );
    /*
    1   0   0
    0   1   0
    0   0   1
*/

So we get:

M= \begin{pmatrix} 2 & -2 & 0 \\ 0 & 0 & -1 \\ 2 & 2 & 0 \end{pmatrix} = \begin{pmatrix} 0.707107 & -0.707107 & 0 \\ 0 & 0 & -1 \\ 0.707107 & 0.707107 & 0 \end{pmatrix} \times \begin{pmatrix} 2.82843 & 0 & 0 \\ 0 & 2.82843 & 0 \\ 0 & 0 & 1 \end{pmatrix} \times \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}

The interface finalM->GetScale(); returns singular values of original matrix.

Categories: MathVTK
Tags:
0 0 votes
Article Rating
Subscribe
Notify of
guest

0 Comments
Inline Feedbacks
View all comments

Related Posts

VTK

Size Of VTK Render Window Will Not Match The External Canvas

I found an interesting issue about vtkSDL2OpenGLRenderWindow::SetSize when changed canvas for VTK WASM 3D window.The related post on the website is Change Canvas For VTK WASM Project. The interface SetSize can work if it was used before m_RenderWindowInteractor->SetRenderWindow( newWin );.It means Read more…

VTK

VTK WASM Error – Uncaught TypeError: Cannot read properties of null (reading ‘bindFramebuffer’) at _glBindFramebuffer

There are some errors and warnings after I changed canvas for VTK wasm. The related post: https://www.weiy.city/2023/11/change-canvas-for-vtk-wasm-project/ We need to free up graphics resources with the old window before we render again. worker.cpp:

VTK

Pick Operation Of vtkPropPicker On vtkSDL2OpenGLRenderWindow Is Not Working

Environment:VTK Branch_v9.1.0emcc (Emscripten gcc/clang-like replacement + linker emulating GNU ld) 3.1.39 I tested two interfaces for pick operation of vtkPropPicker on web page. The result value is awalys 0. Invoking code: There are warning on Read more…

Content Summary
 : Input your strings, the tool can get a brief summary of the content for you.

Got it!
X
0
Would love your thoughts, please comment.x
()
x