DGtal  0.9.3
viewDualSurface.cpp
Go to the documentation of this file.
1 
36 #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 
47 using namespace std;
48 using namespace DGtal;
49 using namespace Z3i;
50 
51 
52 
53 template <typename Vector>
54 Vector 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 
61 template <typename Vector>
62 struct 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.
79 template <typename Vector>
80 void naiveConvexHull
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 ) ) )
105  ++nbBadPos;
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  unsigned int 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  // std::cerr << "----------- " << std::endl;
139  // for ( unsigned int i = 0; i < indices.size(); ++i )
140  // {
141  // std::cerr << "( ";
142  // for ( unsigned int j = 0; j < indices[ i ].size(); ++j )
143  // std::cerr << indices[ i ][ j ] << " ";
144  // std::cerr << ")" << std::endl;
145  // }
146 }
147 
148 double rescale( double x )
149 {
150  return ( x - 1.0 ) * 0.5 + 0.5;
151 }
152 
153 template <typename Viewer,
154  typename Vector>
155 void viewPolygons
156 ( Viewer & viewer,
157  const DGtal::Color & color,
158  const std::vector< std::vector< unsigned int > > & indices,
159  const std::vector<Vector> & points )
160 {
161  typedef typename Viewer::RealPoint RealPoint;
162  std::vector<RealPoint> pts3d;
163  DGtal::Color fillColorSave = viewer.getFillColor();
164  for ( unsigned int f = 0; f < indices.size(); ++f )
165  {
166  pts3d.clear();
167  RealPoint P;
168  for ( unsigned int v = 0; v < indices[ f ].size(); ++v )
169  {
170  unsigned int i = indices[ f ][ v ];
171  P[0] = rescale( points[ i ][ 0 ] );
172  P[1] = rescale( points[ i ][ 1 ] );
173  P[2] = rescale( points[ i ][ 2 ] );
174  pts3d.push_back( P );
175  }
176  viewer.setFillColor(color);
177  viewer.addPolygon( pts3d );
178  }
179 }
180 
181 template <typename Vector>
182 unsigned int dim( const Vector & z )
183 {
184  unsigned int d = 0;
185  for ( unsigned int i = 0; i < Vector::dimension; ++i )
186  if ( ( z[ i ] % 2 ) == 1 ) ++d;
187  return d;
188 }
189 
190 template <typename Vector>
191 unsigned int openDim( const Vector & z )
192 {
193  for ( unsigned int i = 0; i < Vector::dimension; ++i )
194  if ( ( z[ i ] % 2 ) == 1 ) return i;
195  return Vector::dimension;
196 }
197 template <typename Vector>
198 Vector lower( const Vector & z, unsigned int k )
199 {
200  Vector z2( z );
201  --z2[ k ];
202  return z2;
203 }
204 template <typename Vector>
205 Vector upper( const Vector & z, unsigned int k )
206 {
207  Vector z2( z );
208  ++z2[ k ];
209  return z2;
210 }
211 
212 template <typename Vector>
213 unsigned int nbLighted( std::map< Vector, bool > & f,
214  const Vector & z )
215 { // z of dim >=2
216  unsigned int d = dim( z );
217  if ( d == 0 ) return f[ z ] ? 1 : 0;
218  unsigned int i = openDim( z );
219  return nbLighted( f, lower( z, i ) )
220  + nbLighted( f, upper( z, i ) );
221 }
222 
223 
224 // Most similar to convex hull... but not exactly, e.g. cfg 31.
225 template <typename Vector>
226 bool lightBetween( std::map< Vector, bool > & f,
227  const Vector & z )
228 {
229  unsigned int d = dim( z );
230  if ( d == 0 ) return f[ z ];
231  else if ( d == 1 )
232  {
233  unsigned int i = openDim( z );
234  return f[ lower( z, i ) ] || f[ upper( z, i ) ];
235  }
236  else
237  {
238  Vector v1, v2;
239  for ( unsigned int i = 0; i < Vector::dimension; ++i )
240  {
241  v1[ i ] = ( ( z[ i ] % 2 ) == 1 ) ? z[ i ] - 1 : z[ i ];
242  v2[ i ] = ( ( z[ i ] % 2 ) == 1 ) ? z[ i ] + 1 : z[ i ];
243  }
244  Domain domain( v1, v2 );
245  for ( Domain::ConstIterator it = domain.begin(), itE = domain.end();
246  it != itE; ++it )
247  {
248  if ( *it == z ) break;
249  Point zp = z*2 - *it;
250  // std::cerr << *it << " <--> " << zp << std::endl;
251  if ( lightBetween( f, *it ) && lightBetween( f, zp ) )
252  return true;
253  }
254  return false;
255  }
256 
257 }
258 
259 
260 template <typename Vector>
261 bool lightMax( std::map< Vector, bool > & f,
262  const Vector & z )
263 {
264  unsigned int d = dim( z );
265  if ( d == 0 ) return f[ z ];
266  else if ( d == 1 )
267  {
268  unsigned int i = openDim( z );
269  return f[ lower( z, i ) ] || f[ upper( z, i ) ];
270  }
271  else // if ( d > 1 )
272  {
273  unsigned int n = nbLighted( f, z );
274  return n >= 2;
275  }
276 }
277 template <typename Vector>
278 bool lightMinMax( std::map< Vector, bool > & f,
279  const Vector & z )
280 {
281  unsigned int d = dim( z );
282  if ( d == 0 ) return f[ z ];
283  else
284  {
285  Vector tmp( z );
286  bool val = true;
287  for ( unsigned i = 0; i < d; ++i )
288  {
289  unsigned int k = openDim( tmp );
290  tmp = lower( tmp, k );
291  val = val && ( lightMinMax( f, lower( z, k ) ) || lightMinMax( f, upper( z, k ) ) );
292  }
293  return val;
294  }
295 }
296 template <typename Vector>
297 bool lightMaxMin( std::map< Vector, bool > & f,
298  const Vector & z )
299 {
300  unsigned int d = dim( z );
301  if ( d == 0 ) return f[ z ];
302  else
303  {
304  Vector tmp( z );
305  bool val = false;
306  for ( unsigned i = 0; i < d; ++i )
307  {
308  unsigned int k = openDim( tmp );
309  tmp = lower( tmp, k );
310  val = val || ( lightMaxMin( f, lower( z, k ) ) && lightMaxMin( f, upper( z, k ) ) );
311  }
312  return val;
313  }
314 }
315 
316 template <typename Vector>
317 bool lightEpsilon( std::map< Vector, bool > & f,
318  const Vector & z,
319  unsigned int epsilon )
320 {
321  unsigned int d = dim( z );
322  if ( d == 0 ) return f[ z ];
323  else
324  {
325  Vector tmp( z );
326  bool eps_d = ( ( epsilon >> (d-1)) & 1 ) != 0;
327  bool val = eps_d ? true : false;
328  for ( unsigned i = 0; i < d; ++i )
329  {
330  unsigned int k = openDim( tmp );
331  tmp = lower( tmp, k );
332  if ( eps_d )
333  val = val && ( lightEpsilon( f, lower( z, k ), epsilon )
334  || lightEpsilon( f, upper( z, k ), epsilon ) );
335  else
336  val = val || ( lightEpsilon( f, lower( z, k ), epsilon )
337  && lightEpsilon( f, upper( z, k ), epsilon ) );
338  }
339  return val;
340  }
341 }
342 
343 
344 template <typename Vector>
345 void fillCfg( std::map< Vector, bool > & f,
346  const Vector & z,
347  unsigned int cfg )
348 {
349  unsigned int d = dim( z );
350  if ( d == 0 )
351  {
352  f[ z ] = (cfg == 1);
353  //std::cerr << "f[" << z << "] = " << f[ z ] << std::endl;
354  }
355  else
356  {
357  unsigned n = 1 << ( d - 1 );
358  unsigned int cfgLow = 0;
359  unsigned int cfgUp = 0;
360  for ( unsigned int j = 0; j < n; ++j )
361  {
362  cfgLow += ( cfg & 1 ) << j;
363  cfg >>= 1;
364  cfgUp += ( cfg & 1 ) << j;
365  cfg >>= 1;
366  }
367  unsigned int i = openDim( z );
368  fillCfg( f, lower( z, i ), cfgLow );
369  fillCfg( f, upper( z, i ), cfgUp );
370  }
371 }
372 
373 template <typename Vector>
374 void localDualVolume( std::vector<Vector> & points,
375  std::map< Vector, bool > & f,
376  const Vector & z )
377 {
378  points.clear();
379  Z3i::Domain domain( z, z + Vector::diagonal(1) );
380  for ( Z3i::Domain::ConstIterator it = domain.begin(), itE = domain.end();
381  it != itE; ++it )
382  {
383  if ( f[ *it ] ) points.push_back( *it );
384  }
385 }
386 
387 template <typename Vector>
388 struct ConfigPointPredicate
389 {
390  std::map< Vector, bool > & fct;
391  Vector base;
392  ConfigPointPredicate( std::map< Vector, bool > & aFct, Vector aBase )
393  : fct( aFct ), base( aBase ) {}
394  bool operator()( const Vector & p ) const
395  {
396  return fct[ p * 2 + base];
397  }
398 };
399 
400 int main( int argc, char** argv )
401 {
402  typedef KSpace::CellSet CellSet;
403  QApplication application(argc,argv);
404 
405  KSpace KS;
406 
408  viewer.show();
409  DGtal::Color fillColor( 200, 200, 220, 255 );
410  DGtal::Color surfelColor( 255, 0, 0, 150 );
411  DGtal::Color voxelColor( 150, 150, 0, 150 );
412 
413  std::vector<Vector> pts;
414 
415  unsigned int cfg = argc > 1 ? atoi( argv[1] ) : 0;
416  unsigned int cfg2 = argc > 2 ? atoi( argv[2] ) : 255;
417  std::map< Vector, bool > f;
418  for ( unsigned int y = 0; (y < 16) && (cfg <= cfg2); ++y )
419  for ( unsigned int x = 0; (x < 16) && (cfg <= cfg2); ++x, ++cfg )
420  {
421  Vector offset( x*6, y*6, 0 );
422  fillCfg( f, offset + Vector( 1, 1, 1 ), cfg );
423  Domain domain( offset + Vector( 0, 0, 0), offset + Vector( 2, 2, 2 ) );
424  KSpace K;
425  K.init( Vector( 0, 0, 0), Vector( 2, 2, 2 ), true );
426  ConfigPointPredicate<Vector> cpp( f, offset );
427  CellSet aBoundary;
428  Surfaces<KSpace>::uMakeBoundary( aBoundary, K, cpp, Vector( 0, 0, 0), Vector( 1, 1, 1 ) );
429  for ( CellSet::const_iterator it = aBoundary.begin(), itE = aBoundary.end();
430  it != itE; ++it )
431  {
432  viewer << CustomColors3D( surfelColor, surfelColor );
433  viewer << KS.uTranslation( *it, offset/2 );
434  }
435  for ( Domain::ConstIterator it = domain.begin(), itE = domain.end();
436  it != itE; ++it )
437  {
438  // lightEpsilon( f, *it, 5 ); // {1,-1,1}=5 // interesting
439  f[ *it ] = lightBetween( f, *it );
440  }
441  viewer << CustomColors3D( DGtal::Color( 255, 0, 0, 255 ), fillColor );
442  std::vector< std::vector< unsigned int > > indices;
443  Domain domain2( offset + Vector( 0, 0, 0), offset + Vector( 1, 1, 1 ) );
444 
445  for ( Domain::ConstIterator it = domain.begin(), itE = domain.end();
446  it != itE; ++it )
447  {
448  localDualVolume( pts, f, *it );
449  indices.clear();
450  naiveConvexHull( indices, pts, false ); // right_handed
451  viewPolygons( viewer, fillColor, indices, pts );
452  }
453  }
454  viewer << Viewer3D<>::updateDisplay;
455 
456  return application.exec();
457 }
virtual void setFillColor(DGtal::Color aColor)
const ConstIterator & end() const
virtual DGtal::Color getFillColor()
const Domain domain(Point(1, 2), Point(6, 5))
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
Definition: Surfaces.h:78
Vector lower(const Vector &z, unsigned int k)
STL namespace.
void naiveConvexHull(std::vector< std::vector< unsigned int > > &indices, const std::vector< Vector > &points, bool left_handed)
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:141
Vector wedge(const Vector &v1, const Vector &v2)
double rescale(double x)
int main(int argc, char **argv)
unsigned int openDim(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)
const ConstIterator & begin() const
void viewPolygons(Viewer &viewer, const DGtal::Color &color, const std::vector< std::vector< unsigned int > > &indices, const std::vector< Vector > &points)
bool lightBetween(std::map< Vector, bool > &f, const Vector &z)
DGtal is the top-level namespace which contains all DGtal functions and types.
MyPointD Point
Definition: testClone2.cpp:383
bool lightMaxMin(std::map< Vector, bool > &f, const Vector &z)
unsigned int dim(const Vector &z)
FreemanChain< int >::Vector Vector
bool lightMinMax(std::map< Vector, bool > &f, const Vector &z)
void addPolygon(const std::vector< RealPoint > &vertices)
Structure representing an RGB triple with alpha component.
Definition: Color.h:66
bool lightMax(std::map< Vector, bool > &f, const Vector &z)
void fillCfg(std::map< Vector, bool > &f, const Vector &z, unsigned int cfg)
KSpace K
Vector upper(const Vector &z, unsigned int k)
unsigned int nbLighted(std::map< Vector, bool > &f, const Vector &z)