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

VTK – How to make top renderer background color visible?

I added two renderers in renderWindow, and set a blue background for the top renderer. But the color is not visible, here is source code. We can find the description for the usage of the Read more…

VTK

Change ColorLevel Of vtkImageMapper To Adjust The Exposure

If you have used vtkImageActor or vtkImageMapper + vtkActor2D to display an image in VTK renderer window, you can change the value of ColorLevel of vtkImageMapper to adjust the exposure to get a better effect. Read more…

VTK

Make vtkImageActor Object Always Face Camera And Be Under vtkFollower Object

Here is a simple example about make vtkImageActor always face current active camera and keep it under a vtkFollower object. It’s impossible to send image data to vtkFollower directly, so we have to add the Read more…

Tex To PDF
 : convert the Latex file which suffix is tex to a PDF file

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