1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: TestDelaunay2DMeshes.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15
16 // Test meshes obtained with vtkDelaunay2D.
17
18 #include "vtkCellArray.h"
19 #include "vtkDelaunay2D.h"
20 #include "vtkNew.h"
21 #include "vtkPoints.h"
22 #include "vtkPolyData.h"
23 #include "vtkPolyDataReader.h"
24 #include "vtkPolyDataWriter.h"
25 #include "vtkSmartPointer.h"
26 #include "vtkTestUtilities.h"
27 #include "vtkTransform.h"
28 #include "vtkTriangle.h"
29 #include "vtkXMLPolyDataReader.h"
30
31 #define VTK_FAILURE 1
32
CompareMeshes(vtkPolyData * p1,vtkPolyData * p2)33 bool CompareMeshes(vtkPolyData* p1, vtkPolyData* p2)
34 {
35 vtkIdType nbPoints1 = p1->GetNumberOfPoints();
36 vtkIdType nbPoints2 = p2->GetNumberOfPoints();
37 vtkIdType nbCells1 = p1->GetNumberOfCells();
38 vtkIdType nbCells2 = p2->GetNumberOfCells();
39 if (nbPoints1 != nbPoints2 || nbCells1 != nbCells2)
40 {
41 return false;
42 }
43
44 vtkCellArray* polys1 = p1->GetPolys();
45 vtkCellArray* polys2 = p2->GetPolys();
46 polys1->InitTraversal();
47 polys2->InitTraversal();
48 vtkIdType npts1;
49 vtkIdType npts2;
50 const vtkIdType* pts1;
51 const vtkIdType* pts2;
52 while (polys1->GetNextCell(npts1, pts1) && polys2->GetNextCell(npts2, pts2))
53 {
54 if (npts1 != npts2)
55 {
56 return false;
57 }
58 for (vtkIdType i = 0; i < npts1; i++)
59 {
60 if (pts1[i] != pts2[i])
61 {
62 return false;
63 }
64 }
65 }
66
67 return true;
68 }
69
DumpMesh(vtkPolyData * mesh)70 void DumpMesh(vtkPolyData* mesh)
71 {
72 vtkNew<vtkPolyDataWriter> writer;
73 writer->SetInputData(mesh);
74 writer->WriteToOutputStringOn();
75 writer->Write();
76 std::cerr << writer->GetOutputString() << std::endl;
77 }
78
TriangulationTest(const std::string & filePath)79 bool TriangulationTest(const std::string& filePath)
80 {
81 vtkNew<vtkPolyDataReader> inputReader;
82 inputReader->SetFileName((filePath + "-Input.vtk").c_str());
83 inputReader->Update();
84
85 vtkNew<vtkDelaunay2D> delaunay2D;
86 delaunay2D->SetInputConnection(inputReader->GetOutputPort());
87 delaunay2D->SetSourceConnection(inputReader->GetOutputPort());
88 delaunay2D->Update();
89
90 vtkPolyData* obtainedMesh = delaunay2D->GetOutput();
91
92 vtkNew<vtkPolyDataReader> outputReader;
93 outputReader->SetFileName((filePath + "-Output.vtk").c_str());
94 outputReader->Update();
95
96 vtkPolyData* validMesh = outputReader->GetOutput();
97
98 if (!CompareMeshes(validMesh, obtainedMesh))
99 {
100 std::cerr << "Obtained mesh is different from expected! "
101 "Its VTK file follows:"
102 << std::endl;
103 DumpMesh(obtainedMesh);
104 return false;
105 }
106
107 return true;
108 }
109
GetTransform(vtkTransform * transform,vtkPoints * points)110 void GetTransform(vtkTransform* transform, vtkPoints* points)
111 {
112 double zaxis[3] = { 0., 0., 1. };
113 double pt0[3], pt1[3], pt2[3], normal[3];
114 points->GetPoint(0, pt0);
115 points->GetPoint(1, pt1);
116 points->GetPoint(2, pt2);
117 vtkTriangle::ComputeNormal(pt0, pt1, pt2, normal);
118
119 double rotationAxis[3], center[3], rotationAngle;
120 double dotZAxis = vtkMath::Dot(normal, zaxis);
121 if (fabs(1.0 - dotZAxis) < 1e-6)
122 {
123 // Aligned with z-axis
124 rotationAxis[0] = 1.0;
125 rotationAxis[1] = 0.0;
126 rotationAxis[2] = 0.0;
127 rotationAngle = 0.0;
128 }
129 else if (fabs(1.0 + dotZAxis) < 1e-6)
130 {
131 // Co-linear with z-axis, but reversed sense.
132 // Aligned with z-axis
133 rotationAxis[0] = 1.0;
134 rotationAxis[1] = 0.0;
135 rotationAxis[2] = 0.0;
136 rotationAngle = 180.0;
137 }
138 else
139 {
140 // The general case
141 vtkMath::Cross(normal, zaxis, rotationAxis);
142 vtkMath::Normalize(rotationAxis);
143 rotationAngle = vtkMath::DegreesFromRadians(acos(vtkMath::Dot(zaxis, normal)));
144 }
145
146 transform->PreMultiply();
147 transform->Identity();
148 transform->RotateWXYZ(rotationAngle, rotationAxis[0], rotationAxis[1], rotationAxis[2]);
149
150 vtkTriangle::TriangleCenter(pt0, pt1, pt2, center);
151 transform->Translate(-center[0], -center[1], -center[2]);
152 }
153
TessellationTestWithTransform(const std::string & dataPath)154 bool TessellationTestWithTransform(const std::string& dataPath)
155 {
156 std::string transformFilePath = dataPath + "-Transform.vtp";
157 std::string boundaryFilePath = dataPath + "-Input.vtp";
158
159 vtkNew<vtkXMLPolyDataReader> reader;
160 reader->SetFileName(transformFilePath.c_str());
161 reader->Update();
162
163 vtkNew<vtkTransform> transform;
164 vtkPoints* points = reader->GetOutput()->GetPoints();
165 GetTransform(transform, points);
166
167 reader->SetFileName(boundaryFilePath.c_str());
168 reader->Update();
169 vtkPolyData* boundaryPoly = reader->GetOutput();
170
171 vtkNew<vtkDelaunay2D> del2D;
172 del2D->SetInputData(boundaryPoly);
173 del2D->SetSourceData(boundaryPoly);
174 del2D->SetTolerance(0.0);
175 del2D->SetAlpha(0.0);
176 del2D->SetOffset(0);
177 del2D->SetProjectionPlaneMode(VTK_SET_TRANSFORM_PLANE);
178 del2D->SetTransform(transform);
179 del2D->BoundingTriangulationOff();
180 del2D->Update();
181
182 vtkPolyData* outPoly = del2D->GetOutput();
183
184 if (outPoly->GetNumberOfCells() != boundaryPoly->GetNumberOfPoints() - 2)
185 {
186 std::cerr << "Bad triangulation for " << dataPath << "!" << std::endl;
187 std::cerr << "Output has " << outPoly->GetNumberOfCells() << " cells instead of "
188 << boundaryPoly->GetNumberOfPoints() - 2 << std::endl;
189 return false;
190 }
191 return true;
192 }
193
TestDelaunay2DMeshes(int argc,char * argv[])194 int TestDelaunay2DMeshes(int argc, char* argv[])
195 {
196
197 char* data_dir = vtkTestUtilities::GetDataRoot(argc, argv);
198 if (!data_dir)
199 {
200 cerr << "Could not determine data directory." << endl;
201 return VTK_FAILURE;
202 }
203
204 std::string dataPath = std::string(data_dir) + "/Data/Delaunay/";
205 delete[] data_dir;
206
207 bool result = true;
208 result &= TriangulationTest(dataPath + "DomainWithHole");
209 result &= TessellationTestWithTransform(dataPath + "Test1");
210 result &= TessellationTestWithTransform(dataPath + "Test2");
211 result &= TessellationTestWithTransform(dataPath + "Test3");
212 result &= TessellationTestWithTransform(dataPath + "Test4");
213 result &= TessellationTestWithTransform(dataPath + "Test5");
214
215 return result ? EXIT_SUCCESS : EXIT_FAILURE;
216 }
217