DGtal  1.0.0
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  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 
140 double rescale( double x )
141 {
142  return ( x - 1.0 ) * 0.5 + 0.5;
143 }
144 
145 template <typename Viewer,
146  typename Vector>
147 void viewPolygons
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 
173 template <typename Vector>
174 unsigned 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 
182 template <typename Vector>
183 unsigned 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 }
189 template <typename Vector>
190 Vector lower( const Vector & z, unsigned int k )
191 {
192  Vector z2( z );
193  --z2[ k ];
194  return z2;
195 }
196 template <typename Vector>
197 Vector upper( const Vector & z, unsigned int k )
198 {
199  Vector z2( z );
200  ++z2[ k ];
201  return z2;
202 }
203 
204 template <typename Vector>
205 unsigned 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.
217 template <typename Vector>
218 bool 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 
252 template <typename Vector>
253 bool 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 }
269 template <typename Vector>
270 bool 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 }
288 template <typename Vector>
289 bool 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 
308 template <typename Vector>
309 bool 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 
336 template <typename Vector>
337 void 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 
365 template <typename Vector>
366 void localDualVolume( std::vector<Vector> & points,
367  std::map< Vector, bool > & f,
368  const Vector & z )
369 {
370  points.clear();
371  Z3i::Domain domain( z, z + Vector::diagonal(1) );
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 
379 template <typename Vector>
380 struct 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 
392 int 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 }
virtual void setFillColor(DGtal::Color aColor)
const ConstIterator & end() const
virtual DGtal::Color getFillColor()
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)
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
Domain domain
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:165
PointVector< 3, double > RealPoint
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 init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
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)
std::set< Cell > CellSet
Preferred type for defining a set of Cell(s).
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
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
Cell uTranslation(const Cell &p, const Vector &vec) const
Add the vector [vec] to [p].
Vector upper(const Vector &z, unsigned int k)
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
unsigned int nbLighted(std::map< Vector, bool > &f, const Vector &z)