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