DGtal  1.4.beta
testSurfaceMesh.cpp
Go to the documentation of this file.
1 
31 #include <iostream>
32 #include <sstream>
33 #include <algorithm>
34 #include "DGtal/base/Common.h"
35 #include "ConfigTest.h"
36 #include "DGtalCatch.h"
37 #include "DGtal/helpers/StdDefs.h"
38 #include "DGtal/kernel/PointVector.h"
39 #include "DGtal/graph/CUndirectedSimpleGraph.h"
40 #include "DGtal/graph/BreadthFirstVisitor.h"
41 #include "DGtal/shapes/SurfaceMesh.h"
42 #include "DGtal/shapes/SurfaceMeshHelper.h"
43 #include "DGtal/io/readers/SurfaceMeshReader.h"
44 #include "DGtal/io/writers/SurfaceMeshWriter.h"
46 
47 using namespace std;
48 using namespace DGtal;
49 
51 // Functions for testing class SurfaceMesh.
53 
54 
55 
58 {
61  typedef SurfaceMesh< RealPoint, RealVector > PolygonMesh;
63  std::vector< RealPoint > positions;
64  std::vector< Vertices > faces;
65  positions.push_back( RealPoint( 0, 0, 0 ) );
66  positions.push_back( RealPoint( 1, 0, 0 ) );
67  positions.push_back( RealPoint( 0, 1, 0 ) );
68  positions.push_back( RealPoint( 1, 1, 0 ) );
69  positions.push_back( RealPoint( 0, 0, 1 ) );
70  positions.push_back( RealPoint( 1, 0, 1 ) );
71  positions.push_back( RealPoint( 0, 1, 1 ) );
72  positions.push_back( RealPoint( 1, 1, 1 ) );
73  positions.push_back( RealPoint( 1, 0, 2 ) );
74  positions.push_back( RealPoint( 0, 0, 2 ) );
75  faces.push_back( { 1, 0, 2, 3 } );
76  faces.push_back( { 0, 1, 5, 4 } );
77  faces.push_back( { 1, 3, 7, 5 } );
78  faces.push_back( { 3, 2, 6, 7 } );
79  faces.push_back( { 2, 0, 4, 6 } );
80  faces.push_back( { 4, 5, 8, 9 } );
81  return PolygonMesh( positions.cbegin(), positions.cend(),
82  faces.cbegin(), faces.cend() );
83 }
84 
85 
88 {
91  typedef SurfaceMesh< RealPoint, RealVector > PolygonMesh;
93  std::vector< RealPoint > positions;
94  std::vector< Vertices > faces;
95  positions.push_back( RealPoint( 0, 0, 1 ) );
96  positions.push_back( RealPoint( 0, -1, 0 ) );
97  positions.push_back( RealPoint( 1, 0, 0 ) );
98  positions.push_back( RealPoint( 0, 1, 0 ) );
99  positions.push_back( RealPoint( 0, 0, 0 ) );
100  faces.push_back( { 0, 4, 1 } );
101  faces.push_back( { 0, 4, 2 } );
102  faces.push_back( { 0, 4, 3 } );
103  return PolygonMesh( positions.cbegin(), positions.cend(),
104  faces.cbegin(), faces.cend() );
105 }
106 
107 SCENARIO( "SurfaceMesh< RealPoint3 > concept check tests", "[surfmesh][concepts]" )
108 {
111  typedef SurfaceMesh< RealPoint, RealVector > PolygonMesh;
112  BOOST_CONCEPT_ASSERT(( concepts::CUndirectedSimpleGraph< PolygonMesh > ));
113 }
114 
115 SCENARIO( "SurfaceMesh< RealPoint3 > build tests", "[surfmesh][build]" )
116 {
119  typedef SurfaceMesh< RealPoint, RealVector > PolygonMesh;
121  typedef PolygonMesh::Edge Edge;
122  typedef PolygonMesh::Vertex Vertex;
123  GIVEN( "A box with an open side" ) {
124  PolygonMesh polymesh = makeBox();
125  THEN( "The mesh has 10 vertices, v0 has 3 neighbors, v1 has 3 neighbors, etc" ) {
126  REQUIRE( polymesh.size() == 10 );
127  REQUIRE( polymesh.degree( 0 ) == 3 );
128  REQUIRE( polymesh.degree( 1 ) == 3 );
129  REQUIRE( polymesh.degree( 2 ) == 3 );
130  REQUIRE( polymesh.degree( 3 ) == 3 );
131  REQUIRE( polymesh.degree( 4 ) == 4 );
132  REQUIRE( polymesh.degree( 5 ) == 4 );
133  REQUIRE( polymesh.degree( 6 ) == 3 );
134  REQUIRE( polymesh.degree( 7 ) == 3 );
135  REQUIRE( polymesh.degree( 8 ) == 2 );
136  REQUIRE( polymesh.degree( 9 ) == 2 );
137  }
138  THEN( "Euler number is 1 as is the Euler number of a disk." )
139  {
140  REQUIRE( polymesh.nbVertices() == 10 );
141  REQUIRE( polymesh.nbEdges() == 15 );
142  REQUIRE( polymesh.nbFaces() == 6 );
143  REQUIRE( polymesh.Euler() == 1 );
144  }
145  THEN( "Checking distances." )
146  {
147  REQUIRE( polymesh.distance(0,0) == Approx(0.0) );
148  REQUIRE( polymesh.distance(0,7) == Approx(std::sqrt(3)));
149  }
150  THEN( "Breadth-first visiting the mesh from vertex 0, visit {0}, then {1,2,4}, then {3,5,6,9}, then {7,8}." )
151  {
152  BreadthFirstVisitor< PolygonMesh > visitor( polymesh, 0 );
153  std::vector<unsigned long> vertices;
154  std::vector<unsigned long> distances;
155  while ( ! visitor.finished() )
156  {
157  vertices.push_back( visitor.current().first );
158  distances.push_back( visitor.current().second );
159  visitor.expand();
160  }
161  REQUIRE( vertices.size() == 10 );
162  REQUIRE( distances.size() == 10 );
163  int expected_vertices[] = { 0, 1, 2, 4, 3, 5, 6, 9, 7, 8 };
164  int expected_distance[] = { 0, 1, 1, 1, 2, 2, 2, 2, 3, 3 };
165  auto itb = vertices.begin();
166  std::sort( itb+1, itb+4 );
167  std::sort( itb+4, itb+8 );
168  std::sort( itb+8, itb+10 );
169  bool vertices_ok
170  = std::equal( vertices.begin(), vertices.end(), expected_vertices );
171  REQUIRE( vertices_ok );
172  bool distances_ok
173  = std::equal( distances.begin(), distances.end(), expected_distance );
174  REQUIRE( distances_ok );
175  }
176  THEN( "The mesh has 6 boundary edges and 9 manifold inner consistent edges, the boundary is a 1d manifold" ) {
177  auto mani_bdry = polymesh.computeManifoldBoundaryEdges();
178  auto mani_inner = polymesh.computeManifoldInnerEdges();
179  auto mani_inner_c = polymesh.computeManifoldInnerConsistentEdges();
180  auto mani_inner_u = polymesh.computeManifoldInnerUnconsistentEdges();
181  auto non_mani = polymesh.computeNonManifoldEdges();
182  CAPTURE( polymesh );
183  REQUIRE( mani_bdry.size() == 6 );
184  REQUIRE( mani_inner.size() == 9 );
185  REQUIRE( mani_inner_c.size() == 9 );
186  REQUIRE( mani_inner_u.size() == 0 );
187  REQUIRE( non_mani.size() == 0 );
188  }
189  THEN( "The face along (1,3) is a quadrangle (1,3,7,5)" ) {
190  Edge e13 = polymesh.makeEdge( 1, 3 );
191  auto lfs = polymesh.edgeLeftFaces( e13 );
192  Vertices T = polymesh.incidentVertices( lfs[ 0 ] );
193  REQUIRE( T.size() == 4 );
194  std::sort( T.begin(), T.end() );
195  REQUIRE( T[ 0 ] == 1 );
196  REQUIRE( T[ 1 ] == 3 );
197  REQUIRE( T[ 2 ] == 5 );
198  REQUIRE( T[ 3 ] == 7 );
199  }
200  THEN( "The face along (3,1) is a quadrangle (3,1,0,2)" ) {
201  Edge e13 = polymesh.makeEdge( 1, 3 );
202  auto rfs = polymesh.edgeRightFaces( e13 );
203  Vertices T = polymesh.incidentVertices( rfs[ 0 ] );
204  REQUIRE( T.size() == 4 );
205  std::sort( T.begin(), T.end() );
206  REQUIRE( T[ 0 ] == 0 );
207  REQUIRE( T[ 1 ] == 1 );
208  REQUIRE( T[ 2 ] == 2 );
209  REQUIRE( T[ 3 ] == 3 );
210  }
211  THEN( "The lower part of the mesh has the barycenter (0.5, 0.5, 0.5) " ) {
212  auto positions = polymesh.positions();
213  RealPoint b;
214  for ( Vertex v = 0; v < 8; ++v )
215  b += positions[ v ];
216  b /= 8;
217  REQUIRE( b[ 0 ] == 0.5 );
218  REQUIRE( b[ 1 ] == 0.5 );
219  REQUIRE( b[ 2 ] == 0.5 );
220  }
221  THEN( "We can iterate over the vertices" ) {
222  auto positions = polymesh.positions();
223  RealPoint exp_positions[] = { { 0,0,0 }, { 1,0,0 }, { 0,1,0 }, { 1,1,0 },
224  { 0,0,1 }, { 1,0,1 }, { 0,1,1 }, { 1,1,1 },
225  { 1,0,2 }, { 0,0,2 } };
226  for ( auto it = polymesh.begin(), itE = polymesh.end(); it != itE; ++it ) {
227  REQUIRE( positions[ *it ] == exp_positions[ *it ] );
228  }
229  }
230  }
231 }
232 
233 
234 SCENARIO( "SurfaceMesh< RealPoint3 > mesh helper tests", "[surfmesh][helper]" )
235 {
238  typedef SurfaceMeshHelper< RealPoint, RealVector > PolygonMeshHelper;
239  typedef PolygonMeshHelper::NormalsType NormalsType;
240  GIVEN( "A sphere of radius 10" ) {
241  auto polymesh = PolygonMeshHelper::makeSphere( 3.0, RealPoint::zero,
242  10, 10, NormalsType::NO_NORMALS );
243  THEN( "The mesh has Euler characteristic 2" ) {
244  REQUIRE( polymesh.Euler() == 2 );
245  }
246  THEN( "It is a consistent manifold without boundary" ) {
247  auto mani_bdry = polymesh.computeManifoldBoundaryEdges();
248  auto mani_inner = polymesh.computeManifoldInnerEdges();
249  auto mani_inner_c = polymesh.computeManifoldInnerConsistentEdges();
250  auto mani_inner_u = polymesh.computeManifoldInnerUnconsistentEdges();
251  auto non_mani = polymesh.computeNonManifoldEdges();
252  CAPTURE( polymesh );
253  REQUIRE( mani_bdry.size() == 0 );
254  REQUIRE( mani_inner.size() == mani_inner_c.size() );
255  REQUIRE( mani_inner_u.size() == 0 );
256  REQUIRE( non_mani.size() == 0 );
257  }
258  }
259  GIVEN( "A torus with radii 3 and 1" ) {
260  auto polymesh = PolygonMeshHelper::makeTorus( 3.0, 1.0, RealPoint::zero,
261  10, 10, 0, NormalsType::NO_NORMALS );
262  THEN( "The mesh has Euler characteristic 0" ) {
263  REQUIRE( polymesh.Euler() == 0 );
264  }
265  THEN( "It is a consistent manifold without boundary" ) {
266  auto mani_bdry = polymesh.computeManifoldBoundaryEdges();
267  auto mani_inner = polymesh.computeManifoldInnerEdges();
268  auto mani_inner_c = polymesh.computeManifoldInnerConsistentEdges();
269  auto mani_inner_u = polymesh.computeManifoldInnerUnconsistentEdges();
270  auto non_mani = polymesh.computeNonManifoldEdges();
271  CAPTURE( polymesh );
272  REQUIRE( mani_bdry.size() == 0 );
273  REQUIRE( mani_inner.size() == mani_inner_c.size() );
274  REQUIRE( mani_inner_u.size() == 0 );
275  REQUIRE( non_mani.size() == 0 );
276  }
277  }
278  GIVEN( "A lantern with radii 3" ) {
279  auto polymesh = PolygonMeshHelper::makeLantern( 3.0, 3.0, RealPoint::zero,
280  10, 10, NormalsType::NO_NORMALS );
281  THEN( "The mesh has Euler characteristic 0" ) {
282  REQUIRE( polymesh.Euler() == 0 );
283  }
284  THEN( "It is a consistent manifold with boundary" ) {
285  auto mani_bdry = polymesh.computeManifoldBoundaryEdges();
286  auto mani_inner = polymesh.computeManifoldInnerEdges();
287  auto mani_inner_c = polymesh.computeManifoldInnerConsistentEdges();
288  auto mani_inner_u = polymesh.computeManifoldInnerUnconsistentEdges();
289  auto non_mani = polymesh.computeNonManifoldEdges();
290  CAPTURE( polymesh );
291  REQUIRE( mani_bdry.size() == 20 );
292  REQUIRE( mani_inner.size() == mani_inner_c.size() );
293  REQUIRE( mani_inner_u.size() == 0 );
294  REQUIRE( non_mani.size() == 0 );
295  }
296  }
297 }
298 
299 SCENARIO( "SurfaceMesh< RealPoint3 > reader/writer tests", "[surfmesh][io]" )
300 {
303  typedef SurfaceMesh< RealPoint, RealVector > PolygonMesh;
304  typedef SurfaceMeshHelper< RealPoint, RealVector > PolygonMeshHelper;
305  typedef SurfaceMeshReader< RealPoint, RealVector > PolygonMeshReader;
306  typedef SurfaceMeshWriter< RealPoint, RealVector > PolygonMeshWriter;
307  typedef PolygonMeshHelper::NormalsType NormalsType;
308  auto polymesh = PolygonMeshHelper::makeSphere( 3.0, RealPoint::zero,
309  10, 10, NormalsType::VERTEX_NORMALS );
310  WHEN( "Writing the mesh as an OBJ file and reading into another mesh" ) {
311  PolygonMesh readmesh;
312  std::ostringstream output;
313  bool okw = PolygonMeshWriter::writeOBJ( output, polymesh );
314  std::string file = output.str();
315  std::istringstream input( file );
316  bool okr = PolygonMeshReader::readOBJ ( input, readmesh );
317  THEN( "The read mesh is the same as the original one" ) {
318  CAPTURE( file );
319  CAPTURE( polymesh );
320  CAPTURE( readmesh );
321  REQUIRE( okw );
322  REQUIRE( okr );
323  REQUIRE( polymesh.Euler() == readmesh.Euler() );
324  REQUIRE( polymesh.nbVertices() == readmesh.nbVertices() );
325  REQUIRE( polymesh.nbEdges() == readmesh.nbEdges() );
326  REQUIRE( polymesh.nbFaces() == readmesh.nbFaces() );
327  REQUIRE( polymesh.neighborVertices( 0 ).size()
328  == readmesh.neighborVertices( 0 ).size() );
329  REQUIRE( polymesh.neighborVertices( 20 ).size()
330  == readmesh.neighborVertices( 20 ).size() );
331  REQUIRE( polymesh.vertexNormals().size() == readmesh.vertexNormals().size() );
332  }
333  }
334 }
335 
336 SCENARIO( "SurfaceMesh< RealPoint3 > boundary tests", "[surfmesh][boundary]" )
337 {
338  auto polymesh = makeNonManifoldBoundary();
339  auto polymesh2 = makeBox();
340  WHEN( "Checking the topolopgy of the mesh boundary" ) {
341  auto chains = polymesh2.computeManifoldBoundaryChains();
342  THEN( "The box as a manifold boundary" ) {
343  CAPTURE(chains);
344  REQUIRE( polymesh2.isBoundariesManifold() == true);
345  REQUIRE( polymesh2.isBoundariesManifold(false) == true);
346  REQUIRE( chains.size() == 1);
347  REQUIRE( chains[0].size() == 6);
348  }
349  THEN( "The extra mesh does not have a manifold boundary" ) {
350  REQUIRE( polymesh.isBoundariesManifold() == false);
351  }
352  }
353 }
Aim: This class is useful to perform a breadth-first exploration of a graph given a starting point or...
const Node & current() const
Space::RealVector RealVector
SMesh::Vertices Vertices
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: An helper class for building classical meshes.
Aim: An helper class for reading mesh files (Wavefront OBJ at this point) and creating a SurfaceMesh.
Aim: An helper class for writing mesh file formats (Waverfront OBJ at this point) and creating a Surf...
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:92
Aim: Represents the concept of local graph: each vertex has neighboring vertices, but we do not neces...
CAPTURE(thicknessHV)
GIVEN("A cubical complex with random 3-cells")
HalfEdgeDataStructure::Edge Edge
REQUIRE(domain.isInside(aPoint))
SCENARIO("SurfaceMesh< RealPoint3 > concept check tests", "[surfmesh][concepts]")
SurfaceMesh< PointVector< 3, double >, PointVector< 3, double > > makeBox()
SurfaceMesh< PointVector< 3, double >, PointVector< 3, double > > makeNonManifoldBoundary()
PointVector< 3, double > RealPoint
TriMesh::Vertex Vertex