DGtal 1.3.0
Loading...
Searching...
No Matches
standardDigitalPolyhedronBuilder3D.cpp
Go to the documentation of this file.
1
30namespace DGtal {
61} // namespace DGtal {
62
64#include <iostream>
65#include <queue>
66#include "DGtal/base/Common.h"
67#include "DGtal/helpers/StdDefs.h"
68#include "DGtal/io/viewers/Viewer3D.h"
69#include "DGtal/shapes/Shapes.h"
70#include "DGtal/shapes/SurfaceMesh.h"
71#include "DGtal/io/readers/SurfaceMeshReader.h"
72#include "DGtal/geometry/volumes/DigitalConvexity.h"
73#include "ConfigExamples.h"
74
76
77using namespace std;
78using namespace DGtal;
88typedef std::vector<Point> PointRange;
89
90// Convenient class to represent different types of arithmetic planes as a predicate.
91template < bool Naive, bool Symmetric >
92struct MedianPlane {
93 Vector N;
94 Integer mu;
95 Integer omega;
96 MedianPlane() = default;
97 MedianPlane( const MedianPlane& other ) = default;
98 MedianPlane( MedianPlane&& other ) = default;
99 MedianPlane& operator=( const MedianPlane& other ) = default;
100 MedianPlane& operator=( MedianPlane&& other ) = default;
101 MedianPlane( Point p, Point q, Point r )
102 : N ( ( q - p ).crossProduct( r - p ) )
103 {
104 mu = N.dot( p );
105 omega = Naive ? N.norm( N.L_infty ) : N.norm( N.L_1 );
106 if ( Symmetric && ( ( omega & 1 ) == 0 ) ) omega += 1;
107 mu -= omega / 2;
108 }
109 bool operator()( const Point& p ) const
110 {
111 auto r = N.dot( p );
112 return ( mu <= r ) && ( r < mu+omega );
113 }
114};
115
116// Choose your plane !
117// typedef MedianPlane< true, false > Plane; //< Naive, thinnest possible
118// typedef MedianPlane< true, true > Plane; //< Naive, Symmetric
119// typedef MedianPlane< false, false > Plane; //< Standard
120typedef MedianPlane< false, true > Plane; //< Standard, Symmetric, thickest here
121
122int main( int argc, char** argv )
123{
124 trace.info() << "Usage: " << argv[ 0 ] << " <input.obj> <h> <view>" << std::endl;
125 trace.info() << "\tComputes a digital polyhedron from an OBJ file" << std::endl;
126 trace.info() << "\t- input.obj: choose your favorite mesh" << std::endl;
127 trace.info() << "\t- h [==1]: the digitization gridstep" << std::endl;
128 trace.info() << "\t- view [==31]: display vertices(1), common edges(2), positive side f edges(4), negative side f edges (8), faces(16)" << std::endl;
129 string filename = examplesPath + "samples/lion.obj";
130 std::string fn = argc > 1 ? argv[ 1 ] : filename; //< vol filename
131 double h = argc > 2 ? atof( argv[ 2 ] ) : 1.0;
132 int view = argc > 3 ? atoi( argv[ 3 ] ) : 31;
133 // Read OBJ file
134 std::ifstream input( fn.c_str() );
136 bool ok = SurfaceMeshReader< RealPoint, RealVector >::readOBJ( input, surfmesh );
137 if ( ! ok )
138 {
139 trace.error() << "Unable to read obj file : " << fn << std::endl;
140 return 1;
141 }
142
143 QApplication application(argc,argv);
144 typedef Viewer3D<Space,KSpace> MViewer;
145 MViewer viewer;
146 viewer.setWindowTitle("standardDigitalPolyhedronBuilder3D");
147 viewer.show();
148
149 Point lo(-500,-500,-500);
150 Point up(500,500,500);
151 DigitalConvexity< KSpace > dconv( lo, up );
153
154 auto vertices = std::vector<Point>( surfmesh.nbVertices() );
155 for ( auto v : surfmesh )
156 {
157 RealPoint p = (1.0 / h) * surfmesh.position( v );
158 Point q ( (Integer) round( p[ 0 ] ),
159 (Integer) round( p[ 1 ] ),
160 (Integer) round( p[ 2 ] ) );
161 vertices[ v ] = q;
162 }
163 std::set< Point > faces_set, pos_edges_set, neg_edges_set;
164 auto faceVertices = surfmesh.allIncidentVertices();
165
166 trace.beginBlock( "Checking face planarity" );
167 std::vector< Plane > face_planes;
168 face_planes.resize( surfmesh.nbFaces() );
169 bool planarity = true;
170 for ( int f = 0; f < surfmesh.nbFaces() && planarity; ++f )
171 {
172 PointRange X;
173 for ( auto v : faceVertices[ f ] )
174 X.push_back( vertices[ v ] );
175 face_planes[ f ] = Plane( X[ 0 ], X[ 1 ], X[ 2 ] );
176 for ( int v = 3; v < X.size(); v++ )
177 if ( ! face_planes[ f ]( X[ v ] ) )
178 {
179 trace.error() << "Face " << f << " is not planar." << std::endl;
180 planarity = false; break;
181 }
182 }
183 trace.endBlock();
184 if ( ! planarity ) return 1;
185 trace.beginBlock( "Computing polyhedron" );
186 for ( int f = 0; f < surfmesh.nbFaces(); ++f )
187 {
188 PointRange X;
189 for ( auto v : faceVertices[ f ] )
190 X.push_back( vertices[ v ] );
191 auto F = dconv.relativeEnvelope( X, face_planes[ f ], Algorithm::DIRECT );
192 faces_set.insert( F.cbegin(), F.cend() );
193 for ( int i = 0; i < X.size(); i++ )
194 {
195 PointRange Y { X[ i ], X[ (i+1)%X.size() ] };
196 if ( Y[ 1 ] < Y[ 0 ] ) std::swap( Y[ 0 ], Y[ 1 ] );
197 int idx1 = faceVertices[ f ][ i ];
198 int idx2 = faceVertices[ f ][ (i+1)%X.size() ];
199 // Variant (1): edges of both sides have many points in common
200 // auto A = dconv.relativeEnvelope( Y, face_planes[ f ], Algorithm::DIRECT );
201 // Variant (2): edges of both sides have much less points in common
202 auto A = dconv.relativeEnvelope( Y, F, Algorithm::DIRECT );
203 bool pos = idx1 < idx2;
204 (pos ? pos_edges_set : neg_edges_set).insert( A.cbegin(), A.cend() );
205 }
206 }
207 trace.endBlock();
208 std::vector< Point > face_points, common_edge_points, arc_points, final_arc_points ;
209 std::vector< Point > pos_edge_points, neg_edge_points, both_edge_points;
210 std::vector< Point > vertex_points = vertices;
211 std::sort( vertex_points.begin(), vertex_points.end() );
212 std::set_symmetric_difference( pos_edges_set.cbegin(), pos_edges_set.cend(),
213 neg_edges_set.cbegin(), neg_edges_set.cend(),
214 std::back_inserter( arc_points ) );
215 std::set_intersection( pos_edges_set.cbegin(), pos_edges_set.cend(),
216 neg_edges_set.cbegin(), neg_edges_set.cend(),
217 std::back_inserter( common_edge_points ) );
218 std::set_union( pos_edges_set.cbegin(), pos_edges_set.cend(),
219 neg_edges_set.cbegin(), neg_edges_set.cend(),
220 std::back_inserter( both_edge_points ) );
221 std::set_difference( faces_set.cbegin(), faces_set.cend(),
222 both_edge_points.cbegin(), both_edge_points.cend(),
223 std::back_inserter( face_points ) );
224 std::set_difference( pos_edges_set.cbegin(), pos_edges_set.cend(),
225 common_edge_points.cbegin(), common_edge_points.cend(),
226 std::back_inserter( pos_edge_points ) );
227 std::set_difference( neg_edges_set.cbegin(), neg_edges_set.cend(),
228 common_edge_points.cbegin(), common_edge_points.cend(),
229 std::back_inserter( neg_edge_points ) );
230 std::set_difference( common_edge_points.cbegin(), common_edge_points.cend(),
231 vertex_points.cbegin(), vertex_points.cend(),
232 std::back_inserter( final_arc_points ) );
233 auto total = vertex_points.size() + pos_edge_points.size()
234 + neg_edge_points.size()
235 + final_arc_points.size() + face_points.size();
236 trace.info() << "#vertex points=" << vertex_points.size() << std::endl;
237 trace.info() << "#pos edge points=" << pos_edge_points.size() << std::endl;
238 trace.info() << "#neg edge points=" << neg_edge_points.size() << std::endl;
239 trace.info() << "#arc points=" << final_arc_points.size() << std::endl;
240 trace.info() << "#face points=" << face_points.size() << std::endl;
241 trace.info() << "#total points=" << total << std::endl;
242
243 // display everything
244 Color colors[] =
246 Color::Magenta, Color( 200, 200, 200 ) };
247 if ( view & 0x1 )
248 {
249 viewer.setLineColor( colors[ 0 ] );
250 viewer.setFillColor( colors[ 0 ] );
251 for ( auto p : vertices ) viewer << p;
252 }
253 if ( view & 0x2 )
254 {
255 viewer.setLineColor( colors[ 3 ] );
256 viewer.setFillColor( colors[ 3 ] );
257 for ( auto p : final_arc_points ) viewer << p;
258 }
259 if ( view & 0x4 )
260 {
261 viewer.setLineColor( colors[ 1 ] );
262 viewer.setFillColor( colors[ 1 ] );
263 for ( auto p : pos_edge_points ) viewer << p;
264 }
265 if ( view & 0x8 )
266 {
267 viewer.setLineColor( colors[ 2 ] );
268 viewer.setFillColor( colors[ 2 ] );
269 for ( auto p : neg_edge_points ) viewer << p;
270 }
271 if ( view & 0x10 )
272 {
273 viewer.setLineColor( colors[ 4 ] );
274 viewer.setFillColor( colors[ 4 ] );
275 for ( auto p : face_points ) viewer << p;
276 }
277 viewer << MViewer::updateDisplay;
278 return application.exec();
279
280}
281// //
283
Structure representing an RGB triple with alpha component.
Definition: Color.h:68
static const Color Red
Definition: Color.h:416
static const Color Blue
Definition: Color.h:419
static const Color Black
Definition: Color.h:413
static const Color Magenta
Definition: Color.h:421
PointRange relativeEnvelope(const PointRange &Z, const PointRange &Y, EnvelopeAlgorithm algo=EnvelopeAlgorithm::DIRECT) const
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
void beginBlock(const std::string &keyword="")
std::ostream & error()
std::ostream & info()
double endBlock()
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
std::vector< Point > PointRange
DGtal::int32_t Integer
Definition: StdDefs.h:143
DGtal is the top-level namespace which contains all DGtal functions and types.
auto crossProduct(PointVector< 3, LeftEuclideanRing, LeftContainer > const &lhs, PointVector< 3, RightEuclideanRing, RightContainer > const &rhs) -> decltype(DGtal::constructFromArithmeticConversion(lhs, rhs))
Cross product of two 3D Points/Vectors.
Trace trace
Definition: Common.h:154
STL namespace.
Z3i::Integer Integer
std::vector< Point > PointRange
Space::Vector Vector
MedianPlane< false, true > Plane
Space::RealPoint RealPoint
Space::RealVector RealVector
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
static bool readOBJ(std::istream &input, SurfaceMesh &smesh)
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:92
Size nbFaces() const
Definition: SurfaceMesh.h:288
const std::vector< Vertices > & allIncidentVertices() const
Definition: SurfaceMesh.h:362
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:637
Size nbVertices() const
Definition: SurfaceMesh.h:280
int main()
Definition: testBits.cpp:56
MyPointD Point
Definition: testClone2.cpp:383
FreemanChain< int >::Vector Vector
void insert(VContainer1 &c1, LContainer2 &c2, unsigned int idx, double v)