DGtal 1.3.0
Searching...
No Matches
viewDualSurface.cpp
Go to the documentation of this file.
1
37#include <iostream>
38#include <algorithm>
39
40#include "DGtal/base/Common.h"
41#include "DGtal/helpers/StdDefs.h"
42#include "DGtal/topology/helpers/Surfaces.h"
43#include "ConfigExamples.h"
44#include "DGtal/io/viewers/Viewer3D.h"
45
46
47using namespace std;
48using namespace DGtal;
49using namespace Z3i;
50
51
52
53template <typename Vector>
54Vector wedge( const Vector & v1, const Vector & v2 )
55{
56 return Vector( v1[ 1 ] * v2[ 2 ] - v1[ 2 ] * v2[ 1 ],
57 v1[ 2 ] * v2[ 0 ] - v1[ 0 ] * v2[ 2 ],
58 v1[ 0 ] * v2[ 1 ] - v1[ 1 ] * v2[ 0 ] );
59}
60
61template <typename Vector>
62struct LessThanOnFace
63{
64 Vector N; // expected normal
65 Vector P; // origin or first point
66 const std::vector< Vector > & pts;
67 inline LessThanOnFace( const Vector & aN, const Vector & aP,
68 const std::vector< Vector > & aPts )
69 : N( aN ), P( aP ), pts( aPts )
70 {}
71 inline bool operator()( unsigned int i1, unsigned int i2 ) const
72 {
73 return N.dot( wedge( pts[ i1 ] - P, pts[ i2 ] - P ) ) > 0;
74 }
75};
76
77// Very naive convex hull algorithm. O(n^4) complexity ! But very
78// short. Takes care also of polygonal faces.
79template <typename Vector>
81( std::vector< std::vector< unsigned int > > & indices,
82 const std::vector<Vector> & points, bool left_handed )
83{
84 typedef typename Vector::Component Scalar;
85 // Checks all triplets of points.
86 std::vector< unsigned int > aFace;
87 for ( unsigned int i1 = 0; i1 < points.size(); ++i1 )
88 for ( unsigned int i2 = 0; i2 < points.size(); ++i2 )
89 if ( i1 != i2 )
90 for ( unsigned int i3 = i1 > i2 ? i1+1 : i2+1; i3 < points.size(); ++i3 )
91 {
92 Vector P12 = points[ i2 ] - points[ i1 ];
93 Vector P13 = points[ i3 ] - points[ i1 ];
94 Vector N = wedge( P12, P13 );
95 if ( N == Vector::zero ) continue;
96
97 unsigned int nbBadPos = 0;
98 for ( unsigned int i4 = 0; i4 < points.size(); ++i4 )
99 {
100 Vector P14 = points[ i4 ] - points[ i1 ];
101 Scalar c = N.dot( P14 );
102 if ( c == 0 ) aFace.push_back( i4 );
103 else if ( ( left_handed && ( c > 0 ) )
104 || ( ! left_handed && ( c < 0 ) ) )
106 }
107 if ( nbBadPos == 0 )
108 {
109 LessThanOnFace<Vector> LTOF( N, points[ aFace[ 0 ] ], points );
110 std::sort( ++aFace.begin(), aFace.end(), LTOF );
111 indices.push_back( aFace );
112 }
113 aFace.clear();
114 }
115 // purge faces.
116 for ( unsigned int i = 0; i < indices.size(); ++i )
117 {
118 auto s = indices[ i ].size();
119 for ( unsigned int j = i+1; j < indices.size(); )
120 {
121 if ( indices[ j ].size() == s )
122 {
123 bool equal = true;
124 for ( unsigned int k = 0; equal && ( k < s ); ++k )
125 if ( indices[ i ][ k ] != indices[ j ][ k ] )
126 equal = false;
127 if ( equal )
128 {
129 std::swap( indices[ j ], indices.back() );
130 indices.pop_back();
131 }
132 else
133 ++j;
134 }
135 else ++j;
136 }
137 }
138}
139
140double rescale( double x )
141{
142 return ( x - 1.0 ) * 0.5 + 0.5;
143}
144
145template <typename Viewer,
146 typename Vector>
148( Viewer & viewer,
149 const DGtal::Color & color,
150 const std::vector< std::vector< unsigned int > > & indices,
151 const std::vector<Vector> & points )
152{
153 typedef typename Viewer::RealPoint RealPoint;
154 std::vector<RealPoint> pts3d;
155 DGtal::Color fillColorSave = viewer.getFillColor();
156 for ( unsigned int f = 0; f < indices.size(); ++f )
157 {
158 pts3d.clear();
159 RealPoint P;
160 for ( unsigned int v = 0; v < indices[ f ].size(); ++v )
161 {
162 unsigned int i = indices[ f ][ v ];
163 P[0] = rescale( points[ i ][ 0 ] );
164 P[1] = rescale( points[ i ][ 1 ] );
165 P[2] = rescale( points[ i ][ 2 ] );
166 pts3d.push_back( P );
167 }
168 viewer.setFillColor(color);
169 viewer.addPolygon( pts3d );
170 }
171}
172
173template <typename Vector>
174unsigned int dim( const Vector & z )
175{
176 unsigned int d = 0;
177 for ( unsigned int i = 0; i < Vector::dimension; ++i )
178 if ( ( z[ i ] % 2 ) == 1 ) ++d;
179 return d;
180}
181
182template <typename Vector>
183unsigned int openDim( const Vector & z )
184{
185 for ( unsigned int i = 0; i < Vector::dimension; ++i )
186 if ( ( z[ i ] % 2 ) == 1 ) return i;
187 return Vector::dimension;
188}
189template <typename Vector>
190Vector lower( const Vector & z, unsigned int k )
191{
192 Vector z2( z );
193 --z2[ k ];
194 return z2;
195}
196template <typename Vector>
197Vector upper( const Vector & z, unsigned int k )
198{
199 Vector z2( z );
200 ++z2[ k ];
201 return z2;
202}
203
204template <typename Vector>
205unsigned int nbLighted( std::map< Vector, bool > & f,
206 const Vector & z )
207{ // z of dim >=2
208 unsigned int d = dim( z );
209 if ( d == 0 ) return f[ z ] ? 1 : 0;
210 unsigned int i = openDim( z );
211 return nbLighted( f, lower( z, i ) )
212 + nbLighted( f, upper( z, i ) );
213}
214
215
216// Most similar to convex hull... but not exactly, e.g. cfg 31.
217template <typename Vector>
218bool lightBetween( std::map< Vector, bool > & f,
219 const Vector & z )
220{
221 unsigned int d = dim( z );
222 if ( d == 0 ) return f[ z ];
223 else if ( d == 1 )
224 {
225 unsigned int i = openDim( z );
226 return f[ lower( z, i ) ] || f[ upper( z, i ) ];
227 }
228 else
229 {
230 Vector v1, v2;
231 for ( unsigned int i = 0; i < Vector::dimension; ++i )
232 {
233 v1[ i ] = ( ( z[ i ] % 2 ) == 1 ) ? z[ i ] - 1 : z[ i ];
234 v2[ i ] = ( ( z[ i ] % 2 ) == 1 ) ? z[ i ] + 1 : z[ i ];
235 }
236 Domain domain( v1, v2 );
237 for ( Domain::ConstIterator it = domain.begin(), itE = domain.end();
238 it != itE; ++it )
239 {
240 if ( *it == z ) break;
241 Point zp = z*2 - *it;
242 // std::cerr << *it << " <--> " << zp << std::endl;
243 if ( lightBetween( f, *it ) && lightBetween( f, zp ) )
244 return true;
245 }
246 return false;
247 }
248
249}
250
251
252template <typename Vector>
253bool lightMax( std::map< Vector, bool > & f,
254 const Vector & z )
255{
256 unsigned int d = dim( z );
257 if ( d == 0 ) return f[ z ];
258 else if ( d == 1 )
259 {
260 unsigned int i = openDim( z );
261 return f[ lower( z, i ) ] || f[ upper( z, i ) ];
262 }
263 else // if ( d > 1 )
264 {
265 unsigned int n = nbLighted( f, z );
266 return n >= 2;
267 }
268}
269template <typename Vector>
270bool lightMinMax( std::map< Vector, bool > & f,
271 const Vector & z )
272{
273 unsigned int d = dim( z );
274 if ( d == 0 ) return f[ z ];
275 else
276 {
277 Vector tmp( z );
278 bool val = true;
279 for ( unsigned i = 0; i < d; ++i )
280 {
281 unsigned int k = openDim( tmp );
282 tmp = lower( tmp, k );
283 val = val && ( lightMinMax( f, lower( z, k ) ) || lightMinMax( f, upper( z, k ) ) );
284 }
285 return val;
286 }
287}
288template <typename Vector>
289bool lightMaxMin( std::map< Vector, bool > & f,
290 const Vector & z )
291{
292 unsigned int d = dim( z );
293 if ( d == 0 ) return f[ z ];
294 else
295 {
296 Vector tmp( z );
297 bool val = false;
298 for ( unsigned i = 0; i < d; ++i )
299 {
300 unsigned int k = openDim( tmp );
301 tmp = lower( tmp, k );
302 val = val || ( lightMaxMin( f, lower( z, k ) ) && lightMaxMin( f, upper( z, k ) ) );
303 }
304 return val;
305 }
306}
307
308template <typename Vector>
309bool lightEpsilon( std::map< Vector, bool > & f,
310 const Vector & z,
311 unsigned int epsilon )
312{
313 unsigned int d = dim( z );
314 if ( d == 0 ) return f[ z ];
315 else
316 {
317 Vector tmp( z );
318 bool eps_d = ( ( epsilon >> (d-1)) & 1 ) != 0;
319 bool val = eps_d ? true : false;
320 for ( unsigned i = 0; i < d; ++i )
321 {
322 unsigned int k = openDim( tmp );
323 tmp = lower( tmp, k );
324 if ( eps_d )
325 val = val && ( lightEpsilon( f, lower( z, k ), epsilon )
326 || lightEpsilon( f, upper( z, k ), epsilon ) );
327 else
328 val = val || ( lightEpsilon( f, lower( z, k ), epsilon )
329 && lightEpsilon( f, upper( z, k ), epsilon ) );
330 }
331 return val;
332 }
333}
334
335
336template <typename Vector>
337void fillCfg( std::map< Vector, bool > & f,
338 const Vector & z,
339 unsigned int cfg )
340{
341 unsigned int d = dim( z );
342 if ( d == 0 )
343 {
344 f[ z ] = (cfg == 1);
345 //std::cerr << "f[" << z << "] = " << f[ z ] << std::endl;
346 }
347 else
348 {
349 unsigned n = 1 << ( d - 1 );
350 unsigned int cfgLow = 0;
351 unsigned int cfgUp = 0;
352 for ( unsigned int j = 0; j < n; ++j )
353 {
354 cfgLow += ( cfg & 1 ) << j;
355 cfg >>= 1;
356 cfgUp += ( cfg & 1 ) << j;
357 cfg >>= 1;
358 }
359 unsigned int i = openDim( z );
360 fillCfg( f, lower( z, i ), cfgLow );
361 fillCfg( f, upper( z, i ), cfgUp );
362 }
363}
364
365template <typename Vector>
366void localDualVolume( std::vector<Vector> & points,
367 std::map< Vector, bool > & f,
368 const Vector & z )
369{
370 points.clear();
372 for ( Z3i::Domain::ConstIterator it = domain.begin(), itE = domain.end();
373 it != itE; ++it )
374 {
375 if ( f[ *it ] ) points.push_back( *it );
376 }
377}
378
379template <typename Vector>
380struct ConfigPointPredicate
381{
382 std::map< Vector, bool > & fct;
383 Vector base;
384 ConfigPointPredicate( std::map< Vector, bool > & aFct, Vector aBase )
385 : fct( aFct ), base( aBase ) {}
386 bool operator()( const Vector & p ) const
387 {
388 return fct[ p * 2 + base];
389 }
390};
391
392int main( int argc, char** argv )
393{
394 typedef KSpace::CellSet CellSet;
395 QApplication application(argc,argv);
396
397 KSpace KS;
398
400 viewer.show();
401 DGtal::Color fillColor( 200, 200, 220, 255 );
402 DGtal::Color surfelColor( 255, 0, 0, 150 );
403 DGtal::Color voxelColor( 150, 150, 0, 150 );
404
405 std::vector<Vector> pts;
406
407 unsigned int cfg = argc > 1 ? atoi( argv[1] ) : 0;
408 unsigned int cfg2 = argc > 2 ? atoi( argv[2] ) : 255;
409 std::map< Vector, bool > f;
410 for ( unsigned int y = 0; (y < 16) && (cfg <= cfg2); ++y )
411 for ( unsigned int x = 0; (x < 16) && (cfg <= cfg2); ++x, ++cfg )
412 {
413 Vector offset( x*6, y*6, 0 );
414 fillCfg( f, offset + Vector( 1, 1, 1 ), cfg );
415 Domain domain( offset + Vector( 0, 0, 0), offset + Vector( 2, 2, 2 ) );
416 KSpace K;
417 K.init( Vector( 0, 0, 0), Vector( 2, 2, 2 ), true );
418 ConfigPointPredicate<Vector> cpp( f, offset );
419 CellSet aBoundary;
420 Surfaces<KSpace>::uMakeBoundary( aBoundary, K, cpp, Vector( 0, 0, 0), Vector( 1, 1, 1 ) );
421 for ( CellSet::const_iterator it = aBoundary.begin(), itE = aBoundary.end();
422 it != itE; ++it )
423 {
424 viewer << CustomColors3D( surfelColor, surfelColor );
425 viewer << KS.uTranslation( *it, offset/2 );
426 }
427 for ( Domain::ConstIterator it = domain.begin(), itE = domain.end();
428 it != itE; ++it )
429 {
430 // lightEpsilon( f, *it, 5 ); // {1,-1,1}=5 // interesting
431 f[ *it ] = lightBetween( f, *it );
432 }
433 viewer << CustomColors3D( DGtal::Color( 255, 0, 0, 255 ), fillColor );
434 std::vector< std::vector< unsigned int > > indices;
435 Domain domain2( offset + Vector( 0, 0, 0), offset + Vector( 1, 1, 1 ) );
436
437 for ( Domain::ConstIterator it = domain.begin(), itE = domain.end();
438 it != itE; ++it )
439 {
440 localDualVolume( pts, f, *it );
441 indices.clear();
442 naiveConvexHull( indices, pts, false ); // right_handed
443 viewPolygons( viewer, fillColor, indices, pts );
444 }
445 }
446 viewer << Viewer3D<>::updateDisplay;
447
448 return application.exec();
449}
Structure representing an RGB triple with alpha component.
Definition: Color.h:68
virtual DGtal::Color getFillColor()
virtual void setFillColor(DGtal::Color aColor)
void addPolygon(const std::vector< RealPoint > &vertices)
Iterator for HyperRectDomain.
const ConstIterator & begin() const
const ConstIterator & end() const
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
Cell uTranslation(const Cell &p, const Vector &vec) const
Add the vector [vec] to [p].
std::set< Cell > CellSet
Preferred type for defining a set of Cell(s).
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
Integer Component
Type for Vector elements.
Definition: PointVector.h:614
static Self diagonal(Component val=1)
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:1595
static const Dimension dimension
Copy of the static dimension of the Point/Vector.
Definition: PointVector.h:626
static void uMakeBoundary(CellSet &aBoundary, const KSpace &aKSpace, const PointPredicate &pp, const Point &aLowerBound, const Point &aUpperBound)
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
Space::RealPoint RealPoint
Definition: StdDefs.h:170
Space::Vector Vector
Definition: StdDefs.h:169
DGtal is the top-level namespace which contains all DGtal functions and types.
STL namespace.
int main()
Definition: testBits.cpp:56
KSpace K
Domain domain
void naiveConvexHull(std::vector< std::vector< unsigned int > > &indices, const std::vector< Vector > &points, bool left_handed)
void viewPolygons(Viewer &viewer, const DGtal::Color &color, const std::vector< std::vector< unsigned int > > &indices, const std::vector< Vector > &points)
bool lightMaxMin(std::map< Vector, bool > &f, const Vector &z)
Vector lower(const Vector &z, unsigned int k)
bool lightBetween(std::map< Vector, bool > &f, const Vector &z)
unsigned int openDim(const Vector &z)
void fillCfg(std::map< Vector, bool > &f, const Vector &z, unsigned int cfg)
double rescale(double x)
Vector wedge(const Vector &v1, const Vector &v2)
unsigned int nbLighted(std::map< Vector, bool > &f, const Vector &z)
Vector upper(const Vector &z, unsigned int k)
bool lightMax(std::map< Vector, bool > &f, const Vector &z)
bool lightEpsilon(std::map< Vector, bool > &f, const Vector &z, unsigned int epsilon)
void localDualVolume(std::vector< Vector > &points, std::map< Vector, bool > &f, const Vector &z)
bool lightMinMax(std::map< Vector, bool > &f, const Vector &z)