1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18
19 #include "itkMesh.h"
20 #include "itkTriangleCell.h"
21
22 #include <iostream>
23
24
itkTriangleCellTest(int,char * [])25 int itkTriangleCellTest(int, char* [] )
26 {
27
28
29 /**
30 * Define a mesh type that stores a PixelType of "int". Use the defaults for
31 * the other template parameters.
32 */
33 using MeshType = itk::Mesh<int>;
34 using CellTraits = MeshType::CellTraits;
35
36 /**
37 * Define a few cell types which uses a PixelType of "int". Again,
38 * use the defaults for the other parameters. Note that a cell's template
39 * parameters must match those of the mesh into which it is inserted.
40 */
41 using CellInterfaceType = itk::CellInterface< int, CellTraits >;
42 using TriangleCellType = itk::TriangleCell<CellInterfaceType>;
43
44 class TriangleHelper : public TriangleCellType
45 {
46 using Superclass = TriangleCellType;
47 using CoordRepType = Superclass::CoordRepType;
48 using PointsContainer = Superclass::PointsContainer;
49 using InterpolationWeightType = Superclass::InterpolationWeightType;
50
51 public:
52 bool EvaluatePosition(CoordRepType* inputPoint,
53 PointsContainer* points,
54 CoordRepType* closestPoint,
55 CoordRepType pcoord [],
56 double * distance,
57 InterpolationWeightType* weights) override
58 {
59 return this->Superclass::EvaluatePosition( inputPoint,
60 points, closestPoint, pcoord, distance, weights );
61 }
62
63 };
64
65
66 /**
67 * Typedef the generic cell type for the mesh. It is an abstract class,
68 * so we can only use information from it, like get its pointer type.
69 */
70 using CellType = MeshType::CellType;
71 using CellAutoPointer = CellType::CellAutoPointer;
72
73 /**
74 * The type of point stored in the mesh. Because mesh was instantiated
75 * with defaults (itkDefaultStaticMeshTraits), the point dimension is 3 and
76 * the coordinate representation is float.
77 */
78 using PointType = MeshType::PointType;
79
80
81 /**
82 * Create the mesh through its object factory.
83 */
84 MeshType::Pointer mesh = MeshType::New();
85 mesh->DebugOn();
86
87 constexpr unsigned int numberOfPoints = 4;
88 /**
89 * Define the 3d geometric positions for 4 points in a square.
90 */
91 MeshType::CoordRepType testPointCoords[numberOfPoints][3]
92 = { {0,0,0}, {10,0,0}, {10,10,0}, {0,10,0} };
93
94 /**
95 * Add our test points to the mesh.
96 * mesh->SetPoint(pointId, point)
97 * Note that the constructor for Point is public, and takes an array
98 * of coordinates for the point.
99 */
100 for(unsigned int i=0; i < numberOfPoints; ++i)
101 {
102 mesh->SetPoint(i, PointType( testPointCoords[i] ) );
103 }
104
105 /**
106 * Specify the method used for allocating cells
107 */
108 mesh->SetCellsAllocationMethod( MeshType::CellsAllocatedDynamicallyCellByCell );
109
110 /**
111 * Create the test cell. Note that testCell is a generic auto
112 * pointer to a cell; in this example it ends up pointing to
113 * different types of cells.
114 */
115 CellAutoPointer testCell;
116 auto * newcell = new TriangleHelper;
117 testCell.TakeOwnership( newcell ); // polymorphism
118
119 /**
120 * List the points that the polygon will use from the mesh.
121 */
122 MeshType::PointIdentifier polygon1Points1[3] = {2,0,1};
123
124 /**
125 * Assign the points to the tetrahedron through their identifiers.
126 */
127 testCell->SetPointIds(polygon1Points1);
128
129 /**
130 * Add the test cell to the mesh.
131 * mesh->SetCell(cellId, cell)
132 */
133 mesh->SetCell(0, testCell ); // Transfer ownership to the mesh
134 std::cout << "TriangleCell pointer = " << (void const *)testCell.GetPointer() << std::endl;
135 std::cout << "TriangleCell Owner = " << testCell.IsOwner() << std::endl;
136
137 {
138 std::cout << "Test MakeCopy" << std::endl;
139
140 CellAutoPointer anotherCell;
141
142 testCell->MakeCopy( anotherCell );
143
144 if( anotherCell->GetNumberOfPoints() != testCell->GetNumberOfPoints() )
145 {
146 std::cerr << "Make Copy failed !" << std::endl;
147 return EXIT_FAILURE;
148 }
149 }
150
151 //
152 // Exercise the EvaluatePosition() method of the TriangleCell
153 //
154 TriangleCellType::CoordRepType inputPoint[3];
155 TriangleCellType::PointsContainer * points = mesh->GetPoints();
156 TriangleCellType::CoordRepType closestPoint[3];
157 TriangleCellType::CoordRepType pcoords[3];
158 double distance;
159 TriangleCellType::InterpolationWeightType weights[3];
160
161 const double tolerance = 1e-5;
162
163 bool isInside;
164
165 // Test 1:
166 inputPoint[0] = 5.0;
167 inputPoint[1] = 3.0;
168 inputPoint[2] = 0.0;
169
170 std::cout << "Calling EvaluatePosition for ";
171 std::cout << inputPoint[0] << ", ";
172 std::cout << inputPoint[1] << ", ";
173 std::cout << inputPoint[2] << std::endl;
174
175 isInside = testCell->EvaluatePosition(inputPoint,
176 points, closestPoint, pcoords , &distance, weights);
177
178 if( !isInside )
179 {
180 std::cerr << "Error: point should be reported as being inside" << std::endl;
181 return EXIT_FAILURE;
182 }
183
184 if( ( itk::Math::abs( pcoords[0] - 0.3 ) > tolerance ) ||
185 ( itk::Math::abs( pcoords[1] - 0.5 ) > tolerance ) ||
186 ( itk::Math::abs( pcoords[2] - 0.2 ) > tolerance ) )
187 {
188 std::cerr << "Error: pcoords computed incorrectly" << std::endl;
189 std::cerr << "pcoords[0] = " << pcoords[0] << std::endl;
190 std::cerr << "pcoords[1] = " << pcoords[1] << std::endl;
191 std::cerr << "pcoords[2] = " << pcoords[2] << std::endl;
192 return EXIT_FAILURE;
193 }
194
195 if( ( itk::Math::abs( weights[0] - 0.3 ) > tolerance ) ||
196 ( itk::Math::abs( weights[1] - 0.5 ) > tolerance ) ||
197 ( itk::Math::abs( weights[2] - 0.2 ) > tolerance ) )
198 {
199 std::cerr << "Error: weights computed incorrectly" << std::endl;
200 std::cerr << "weights[0] = " << weights[0] << std::endl;
201 std::cerr << "weights[1] = " << weights[1] << std::endl;
202 std::cerr << "weights[2] = " << weights[2] << std::endl;
203 return EXIT_FAILURE;
204 }
205
206 if ((itk::Math::abs(closestPoint[0] - 5.0) > tolerance) ||
207 (itk::Math::abs(closestPoint[1] - 3.0) > tolerance) ||
208 (itk::Math::abs(closestPoint[2] - 0.0) > tolerance))
209 {
210 std::cerr << "Error: closestPoint computed incorrectly" << std::endl;
211 std::cerr << "closestPoint[0] = " << closestPoint[0] << std::endl;
212 std::cerr << "closestPoint[1] = " << closestPoint[1] << std::endl;
213 std::cerr << "closestPoint[2] = " << closestPoint[2] << std::endl;
214 return EXIT_FAILURE;
215 }
216
217
218 // Test 2:
219 inputPoint[0] = 15.0;
220 inputPoint[1] = 5.0;
221 inputPoint[2] = 0.0;
222
223 std::cout << "Calling EvaluatePosition for ";
224 std::cout << inputPoint[0] << ", ";
225 std::cout << inputPoint[1] << ", ";
226 std::cout << inputPoint[2] << std::endl;
227
228 isInside = testCell->EvaluatePosition(inputPoint,
229 points, closestPoint, pcoords , &distance, weights);
230
231 if( isInside )
232 {
233 std::cerr << "Error: point should be reported as being outside" << std::endl;
234 return EXIT_FAILURE;
235 }
236
237 if( ( itk::Math::abs( pcoords[0] - 0.5 ) > tolerance ) ||
238 ( itk::Math::abs( pcoords[1] + 0.5 ) > tolerance ) ||
239 ( itk::Math::abs( pcoords[2] - 1.0 ) > tolerance ) )
240 {
241 std::cerr << "Error: pcoords computed incorrectly" << std::endl;
242 std::cerr << "pcoords[0] = " << pcoords[0] << std::endl;
243 std::cerr << "pcoords[1] = " << pcoords[1] << std::endl;
244 std::cerr << "pcoords[2] = " << pcoords[2] << std::endl;
245 return EXIT_FAILURE;
246 }
247
248 if ((itk::Math::abs(closestPoint[0] - 10.0) > tolerance) ||
249 (itk::Math::abs(closestPoint[1] - 5.0) > tolerance) ||
250 (itk::Math::abs(closestPoint[2] - 0.0) > tolerance))
251 {
252 std::cerr << "Error: closestPoint computed incorrectly" << std::endl;
253 std::cerr << "closestPoint[0] = " << closestPoint[0] << std::endl;
254 std::cerr << "closestPoint[1] = " << closestPoint[1] << std::endl;
255 std::cerr << "closestPoint[2] = " << closestPoint[2] << std::endl;
256 return EXIT_FAILURE;
257 }
258
259 //
260 // NOTE: Outside points don't get their weights computed.
261 //
262
263
264 // Test 3:
265 inputPoint[0] = 0.0;
266 inputPoint[1] = 10.0;
267 inputPoint[2] = 0.0;
268
269 std::cout << "Calling EvaluatePosition for ";
270 std::cout << inputPoint[0] << ", ";
271 std::cout << inputPoint[1] << ", ";
272 std::cout << inputPoint[2] << std::endl;
273
274 isInside = testCell->EvaluatePosition(inputPoint,
275 points, closestPoint, pcoords , &distance, weights);
276
277 if( isInside )
278 {
279 std::cerr << "Error: point should be reported as being outside" << std::endl;
280 return EXIT_FAILURE;
281 }
282
283 if( ( itk::Math::abs( pcoords[0] - 1.0 ) > tolerance ) ||
284 ( itk::Math::abs( pcoords[1] - 1.0 ) > tolerance ) ||
285 ( itk::Math::abs( pcoords[2] + 1.0 ) > tolerance ) )
286 {
287 std::cerr << "Error: pcoords computed incorrectly" << std::endl;
288 std::cerr << "pcoords[0] = " << pcoords[0] << std::endl;
289 std::cerr << "pcoords[1] = " << pcoords[1] << std::endl;
290 std::cerr << "pcoords[2] = " << pcoords[2] << std::endl;
291 return EXIT_FAILURE;
292 }
293
294 //note should be half way up the line (0,0,0) -> (10,10,0)
295 if ((itk::Math::abs(closestPoint[0] - 5.0) > tolerance) ||
296 (itk::Math::abs(closestPoint[1] - 5.0) > tolerance) ||
297 (itk::Math::abs(closestPoint[2] - 0.0) > tolerance))
298 {
299 std::cerr << "Error: closestPoint computed incorrectly" << std::endl;
300 std::cerr << "closestPoint[0] = " << closestPoint[0] << std::endl;
301 std::cerr << "closestPoint[1] = " << closestPoint[1] << std::endl;
302 std::cerr << "closestPoint[2] = " << closestPoint[2] << std::endl;
303 return EXIT_FAILURE;
304 }
305
306 //
307 // NOTE: Outside points don't get their weights computed.
308 //
309
310
311 // Test 4:
312 inputPoint[0] = 5.0;
313 inputPoint[1] = -5.0;
314 inputPoint[2] = 0.0;
315
316 std::cout << "Calling EvaluatePosition for ";
317 std::cout << inputPoint[0] << ", ";
318 std::cout << inputPoint[1] << ", ";
319 std::cout << inputPoint[2] << std::endl;
320
321 isInside = testCell->EvaluatePosition(inputPoint,
322 points, closestPoint, pcoords , &distance, weights);
323
324 if( isInside )
325 {
326 std::cerr << "Error: point should be reported as being outside" << std::endl;
327 return EXIT_FAILURE;
328 }
329
330 if( ( itk::Math::abs( pcoords[0] + 0.5 ) > tolerance ) ||
331 ( itk::Math::abs( pcoords[1] - 0.5 ) > tolerance ) ||
332 ( itk::Math::abs( pcoords[2] - 1.0 ) > tolerance ) )
333 {
334 std::cerr << "Error: pcoords computed incorrectly" << std::endl;
335 std::cerr << "pcoords[0] = " << pcoords[0] << std::endl;
336 std::cerr << "pcoords[1] = " << pcoords[1] << std::endl;
337 std::cerr << "pcoords[2] = " << pcoords[2] << std::endl;
338 return EXIT_FAILURE;
339 }
340
341 if ((itk::Math::abs(closestPoint[0] - 5.0) > tolerance) ||
342 (itk::Math::abs(closestPoint[1] - 0.0) > tolerance) ||
343 (itk::Math::abs(closestPoint[2] - 0.0) > tolerance))
344 {
345 std::cerr << "Error: closestPoint computed incorrectly" << std::endl;
346 std::cerr << "closestPoint[0] = " << closestPoint[0] << std::endl;
347 std::cerr << "closestPoint[1] = " << closestPoint[1] << std::endl;
348 std::cerr << "closestPoint[2] = " << closestPoint[2] << std::endl;
349 return EXIT_FAILURE;
350 }
351
352 //
353 // NOTE: Outside points don't get their weights computed.
354 //
355
356
357 // Test 5:
358 inputPoint[0] = -5.0;
359 inputPoint[1] = -3.0;
360 inputPoint[2] = 0.0;
361
362 std::cout << "Calling EvaluatePosition for ";
363 std::cout << inputPoint[0] << ", ";
364 std::cout << inputPoint[1] << ", ";
365 std::cout << inputPoint[2] << std::endl;
366
367 isInside = testCell->EvaluatePosition(inputPoint,
368 points, closestPoint, pcoords , &distance, weights);
369
370 if( isInside )
371 {
372 std::cerr << "Error: point should be reported as being outside" << std::endl;
373 return EXIT_FAILURE;
374 }
375
376 if( ( itk::Math::abs( pcoords[0] + 0.3 ) > tolerance ) ||
377 ( itk::Math::abs( pcoords[1] - 1.5 ) > tolerance ) ||
378 ( itk::Math::abs( pcoords[2] + 0.2 ) > tolerance ) )
379 {
380 std::cerr << "Error: pcoords computed incorrectly" << std::endl;
381 std::cerr << "pcoords[0] = " << pcoords[0] << std::endl;
382 std::cerr << "pcoords[1] = " << pcoords[1] << std::endl;
383 std::cerr << "pcoords[2] = " << pcoords[2] << std::endl;
384 return EXIT_FAILURE;
385 }
386
387 if ((itk::Math::abs(closestPoint[0] - 0.0) > tolerance) ||
388 (itk::Math::abs(closestPoint[1] - 0.0) > tolerance) ||
389 (itk::Math::abs(closestPoint[2] - 0.0) > tolerance))
390 {
391 std::cerr << "Error: closestPoint computed incorrectly" << std::endl;
392 std::cerr << "closestPoint[0] = " << closestPoint[0] << std::endl;
393 std::cerr << "closestPoint[1] = " << closestPoint[1] << std::endl;
394 std::cerr << "closestPoint[2] = " << closestPoint[2] << std::endl;
395 return EXIT_FAILURE;
396 }
397
398 //
399 // NOTE: Outside points don't get their weights computed.
400 //
401
402 return EXIT_SUCCESS;
403 }
404