Skip to content

Commit 20da539

Browse files
g-gisbertg-gisbert
committed
Add simple transform filter based on 4 input points
1 parent d0693b8 commit 20da539

File tree

3 files changed

+333
-0
lines changed

3 files changed

+333
-0
lines changed

‎Filters/General/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,7 @@ set(classes
8585
vtkPassArrays
8686
vtkPassSelectedArrays
8787
vtkPointConnectivityFilter
88+
vtkPointsMatchingTransformFilter
8889
vtkPolyDataStreamer
8990
vtkPolyDataToReebGraphFilter
9091
vtkProbePolyhedron

‎Filters/General/vtkPointsMatchingTransformFilter.cxx‎

Lines changed: 231 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,231 @@
1+
// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
2+
// SPDX-License-Identifier: BSD-3-Clause
3+
#include "vtkPointsMatchingTransformFilter.h"
4+
5+
#include "vtkInformation.h"
6+
#include "vtkInformationVector.h"
7+
#include "vtkMatrix3x3.h"
8+
#include "vtkMatrix4x4.h"
9+
#include "vtkPointSet.h"
10+
#include "vtkSmartPointer.h"
11+
#include "vtkSphericalPointIterator.h"
12+
#include "vtkTransform.h"
13+
#include "vtkTransformFilter.h"
14+
15+
#include <array>
16+
17+
VTK_ABI_NAMESPACE_BEGIN
18+
vtkStandardNewMacro(vtkPointsMatchingTransformFilter);
19+
vtkCxxSetSmartPointerMacro(vtkPointsMatchingTransformFilter, SourceMatrix, vtkMatrix4x4);
20+
vtkCxxSetSmartPointerMacro(vtkPointsMatchingTransformFilter, TargetMatrix, vtkMatrix4x4);
21+
22+
namespace
23+
{
24+
//------------------------------------------------------------------------------
25+
vtkSmartPointer<vtkMatrix3x3> ExtractRotationFromMatrix4x4(vtkMatrix4x4* M)
26+
{
27+
vtkNew<vtkMatrix3x3> rotation;
28+
double data[9];
29+
for (int i = 0; i < 9; ++i)
30+
{
31+
data[i] = M->GetElement(i / 3, i % 3);
32+
}
33+
rotation->SetData(data);
34+
return rotation;
35+
}
36+
37+
//------------------------------------------------------------------------------
38+
void SetRotationInMatrix4x4(vtkMatrix3x3* source, vtkMatrix4x4* dest)
39+
{
40+
for (int i = 0; i < 9; ++i)
41+
{
42+
dest->SetElement(i / 3, i % 3, source->GetElement(i / 3, i % 3));
43+
}
44+
}
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+
84+
//------------------------------------------------------------------------------
85+
vtkPointsMatchingTransformFilter::vtkPointsMatchingTransformFilter()
86+
{
87+
this->SourceMatrix = vtkSmartPointer<vtkMatrix4x4>::New();
88+
this->TargetMatrix = vtkSmartPointer<vtkMatrix4x4>::New();
89+
constexpr std::array<double, 16> init = { 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1 };
90+
this->SourceMatrix->SetData(init.data());
91+
this->TargetMatrix->SetData(init.data());
92+
}
93+
94+
//------------------------------------------------------------------------------
95+
int vtkPointsMatchingTransformFilter::RequestData(vtkInformation* vtkNotUsed(request),
96+
vtkInformationVector** inputVector, vtkInformationVector* outputVector)
97+
{
98+
vtkSmartPointer<vtkPointSet> input = vtkPointSet::GetData(inputVector[0]);
99+
vtkSmartPointer<vtkPointSet> output = vtkPointSet::GetData(outputVector);
100+
if (!input)
101+
{
102+
vtkErrorMacro(<< "Invalid or missing input");
103+
return 0;
104+
}
105+
106+
// Build transform matrix
107+
vtkNew<vtkMatrix4x4> srcInv;
108+
srcInv->DeepCopy(this->SourceMatrix);
109+
if (std::abs(srcInv->Determinant()) < 10e-4)
110+
{
111+
vtkWarningMacro("Source matrix is not invertible. Source points are likely coplanar.");
112+
output->ShallowCopy(input);
113+
return 1;
114+
}
115+
srcInv->Invert();
116+
vtkNew<vtkMatrix4x4> transformMatrix;
117+
vtkMatrix4x4::Multiply4x4(this->TargetMatrix, srcInv, transformMatrix);
118+
119+
if (this->RigidTransform)
120+
{
121+
vtkSmartPointer<vtkMatrix3x3> rotation = ::ExtractRotationFromMatrix4x4(transformMatrix);
122+
rotation = ::PolarDecomposition(rotation);
123+
::SetRotationInMatrix4x4(rotation, transformMatrix);
124+
}
125+
126+
// Apply the transform
127+
vtkNew<vtkTransform> transform;
128+
transform->SetMatrix(transformMatrix);
129+
vtkNew<vtkTransformFilter> transformFilter;
130+
transformFilter->SetInputData(input);
131+
transformFilter->SetTransform(transform);
132+
transformFilter->Update();
133+
134+
output->ShallowCopy(transformFilter->GetOutput());
135+
return 1;
136+
}
137+
138+
//------------------------------------------------------------------------------
139+
void vtkPointsMatchingTransformFilter::SetSourcePoint1(double x, double y, double z)
140+
{
141+
this->SetSourcePoint(0, x, y, z);
142+
}
143+
144+
//------------------------------------------------------------------------------
145+
void vtkPointsMatchingTransformFilter::SetSourcePoint2(double x, double y, double z)
146+
{
147+
this->SetSourcePoint(1, x, y, z);
148+
}
149+
150+
//------------------------------------------------------------------------------
151+
void vtkPointsMatchingTransformFilter::SetSourcePoint3(double x, double y, double z)
152+
{
153+
this->SetSourcePoint(2, x, y, z);
154+
}
155+
156+
//------------------------------------------------------------------------------
157+
void vtkPointsMatchingTransformFilter::SetSourcePoint4(double x, double y, double z)
158+
{
159+
this->SetSourcePoint(3, x, y, z);
160+
}
161+
162+
//------------------------------------------------------------------------------
163+
void vtkPointsMatchingTransformFilter::SetTargetPoint1(double x, double y, double z)
164+
{
165+
this->SetTargetPoint(0, x, y, z);
166+
}
167+
168+
//------------------------------------------------------------------------------
169+
void vtkPointsMatchingTransformFilter::SetTargetPoint2(double x, double y, double z)
170+
{
171+
this->SetTargetPoint(1, x, y, z);
172+
}
173+
174+
//------------------------------------------------------------------------------
175+
void vtkPointsMatchingTransformFilter::SetTargetPoint3(double x, double y, double z)
176+
{
177+
this->SetTargetPoint(2, x, y, z);
178+
}
179+
180+
//------------------------------------------------------------------------------
181+
void vtkPointsMatchingTransformFilter::SetTargetPoint4(double x, double y, double z)
182+
{
183+
this->SetTargetPoint(3, x, y, z);
184+
}
185+
186+
//------------------------------------------------------------------------------
187+
void vtkPointsMatchingTransformFilter::SetSourcePoint(int index, double x, double y, double z)
188+
{
189+
this->SourceMatrix->SetElement(0, index, x);
190+
this->SourceMatrix->SetElement(1, index, y);
191+
this->SourceMatrix->SetElement(2, index, z);
192+
}
193+
194+
//------------------------------------------------------------------------------
195+
void vtkPointsMatchingTransformFilter::SetTargetPoint(int index, double x, double y, double z)
196+
{
197+
this->TargetMatrix->SetElement(0, index, x);
198+
this->TargetMatrix->SetElement(1, index, y);
199+
this->TargetMatrix->SetElement(2, index, z);
200+
}
201+
202+
//------------------------------------------------------------------------------
203+
vtkMTimeType vtkPointsMatchingTransformFilter::GetMTime()
204+
{
205+
vtkMTimeType mTime = this->Superclass::GetMTime();
206+
vtkMTimeType time;
207+
208+
if (this->SourceMatrix != nullptr)
209+
{
210+
time = this->SourceMatrix->GetMTime();
211+
mTime = std::max(time, mTime);
212+
}
213+
if (this->TargetMatrix != nullptr)
214+
{
215+
time = this->TargetMatrix->GetMTime();
216+
mTime = std::max(time, mTime);
217+
}
218+
219+
return mTime;
220+
}
221+
222+
//------------------------------------------------------------------------------
223+
void vtkPointsMatchingTransformFilter::PrintSelf(ostream& os, vtkIndent indent)
224+
{
225+
this->Superclass::PrintSelf(os, indent);
226+
227+
os << indent << "Source Matrix: " << this->SourceMatrix << "\n";
228+
os << indent << "Target Matrix: " << this->SourceMatrix << "\n";
229+
os << indent << "Rigid Transform: " << this->RigidTransform << "\n";
230+
}
231+
VTK_ABI_NAMESPACE_END

‎Filters/General/vtkPointsMatchingTransformFilter.h‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
2+
// SPDX-License-Identifier: BSD-3-Clause
3+
/**
4+
* @class vtkPointsMatchingTransformFilter
5+
* @brief affine transform a vtkPointSet based on 4 pairs of points
6+
*
7+
* Allow the user to specify 4 pairs of points in space, to compute
8+
* the affine transform between the two defined coordinates systems and
9+
* to perform the transform on the input.
10+
*/
11+
12+
#ifndef vtkPointsMatchingTransformFilter_h
13+
#define vtkPointsMatchingTransformFilter_h
14+
15+
#include "vtkFiltersGeneralModule.h" // For export macro
16+
#include "vtkPointSetAlgorithm.h"
17+
18+
#include "vtkMatrix4x4.h" // For vtkMatrix4x4
19+
20+
VTK_ABI_NAMESPACE_BEGIN
21+
class vtkMatrix3x3;
22+
23+
class VTKFILTERSGENERAL_EXPORT vtkPointsMatchingTransformFilter : public vtkPointSetAlgorithm
24+
{
25+
public:
26+
///@{
27+
/**
28+
* Standard methods for instantiation, obtaining type information, and
29+
* printing.
30+
*/
31+
static vtkPointsMatchingTransformFilter* New();
32+
vtkTypeMacro(vtkPointsMatchingTransformFilter, vtkPointSetAlgorithm);
33+
void PrintSelf(ostream& os, vtkIndent indent) override;
34+
///@}
35+
36+
///@{
37+
/**
38+
* Sets one of the four point pairs that define the transformation.
39+
* These methods directly update the corresponding column in the Source
40+
* or Target matrix used to compute the transform.
41+
*/
42+
void SetSourcePoint1(double x, double y, double z);
43+
void SetSourcePoint2(double x, double y, double z);
44+
void SetSourcePoint3(double x, double y, double z);
45+
void SetSourcePoint4(double x, double y, double z);
46+
void SetTargetPoint1(double x, double y, double z);
47+
void SetTargetPoint2(double x, double y, double z);
48+
void SetTargetPoint3(double x, double y, double z);
49+
void SetTargetPoint4(double x, double y, double z);
50+
void SetSourcePoint(int index, double x, double y, double z);
51+
void SetTargetPoint(int index, double x, double y, double z);
52+
///@}
53+
54+
///@{
55+
/**
56+
* Set/Get the source matrix.
57+
*/
58+
void SetSourceMatrix(vtkMatrix4x4*);
59+
vtkGetObjectMacro(SourceMatrix, vtkMatrix4x4);
60+
///@}
61+
62+
///@{
63+
/**
64+
* Set/Get the target matrix.
65+
*/
66+
void SetTargetMatrix(vtkMatrix4x4*);
67+
vtkGetObjectMacro(TargetMatrix, vtkMatrix4x4);
68+
///@}
69+
70+
///@{
71+
/**
72+
* Set/Get the RigidTransform option. If true, approximate the transform with the closest rigid +
73+
* scale transform. Default is false.
74+
*/
75+
vtkSetMacro(RigidTransform, bool);
76+
vtkGetMacro(RigidTransform, bool);
77+
vtkBooleanMacro(RigidTransform, bool);
78+
///@}
79+
80+
/**
81+
* Return the mtime also considering the source and target matrices.
82+
*/
83+
vtkMTimeType GetMTime() override;
84+
85+
protected:
86+
vtkPointsMatchingTransformFilter();
87+
~vtkPointsMatchingTransformFilter() override = default;
88+
89+
int RequestData(vtkInformation*, vtkInformationVector**, vtkInformationVector*) override;
90+
91+
private:
92+
vtkPointsMatchingTransformFilter(const vtkPointsMatchingTransformFilter&) = delete;
93+
void operator=(const vtkPointsMatchingTransformFilter&) = delete;
94+
95+
vtkSmartPointer<vtkMatrix4x4> SourceMatrix;
96+
vtkSmartPointer<vtkMatrix4x4> TargetMatrix;
97+
bool RigidTransform = false;
98+
};
99+
100+
VTK_ABI_NAMESPACE_END
101+
#endif

0 commit comments

Comments
 (0)