DGtal  0.9.3beta
trackImplicitPolynomialSurfaceToOFF.cpp
1 
30 #include <iostream>
33 #include "DGtal/helpers/StdDefs.h"
34 #include "DGtal/topology/helpers/Surfaces.h"
35 #include "DGtal/topology/DigitalSurface.h"
36 #include "DGtal/topology/SetOfSurfels.h"
37 #include "DGtal/math/MPolynomial.h"
38 #include "DGtal/shapes/GaussDigitizer.h"
39 #include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
40 #include "DGtal/shapes/implicit/ImplicitFunctionDiff1LinearCellEmbedder.h"
41 #include "DGtal/io/readers/MPolynomialReader.h"
43 
45 
46 using namespace std;
47 using namespace DGtal;
48 using namespace Z3i;
49 
51 
52 
53 void usage( int /*argc*/, char** argv )
54 {
55  std::cerr << "Usage: " << argv[ 0 ] << " <Polynomial> <Px> <Py> <Pz> <Qx> <Qy> <Qz> <step>" << std::endl;
56  std::cerr << "\t - displays the boundary of a shape defined implicitly by a 3-polynomial <Polynomial>." << std::endl;
57  std::cerr << "\t - P and Q defines the bounding box." << std::endl;
58  std::cerr << "\t - step is the grid step." << std::endl;
59  std::cerr << "\t - You may try x^3y+xz^3+y^3z+z^3+5z or (y^2+z^2-1)^2 +(x^2+y^2-1)^3 " << std::endl;
60  std::cerr << "\t - See http://www.freigeist.cc/gallery.html" << std::endl;
61 }
62 
63 int main( int argc, char** argv )
64 {
65  if ( argc < 9 )
66  {
67  usage( argc, argv );
68  return 1;
69  }
70  double p1[ 3 ];
71  double p2[ 3 ];
72  for ( unsigned int i = 0; i < 3; ++i )
73  {
74  p1[ i ] = atof( argv[ 2+i ] );
75  p2[ i ] = atof( argv[ 5+i ] );
76  }
77  double step = atof( argv[ 8 ] );
78 
80  trace.beginBlock( "Making polynomial surface." );
82  typedef RealPoint::Coordinate Ring;
88 
89  // See http://www.freigeist.cc/gallery.html
90  Polynomial3 P;
91  Polynomial3Reader reader;
92  std::string poly_str = argv[ 1 ];
93  std::string::const_iterator iter
94  = reader.read( P, poly_str.begin(), poly_str.end() );
95  if ( iter != poly_str.end() )
96  {
97  std::cerr << "ERROR: I read only <"
98  << poly_str.substr( 0, iter - poly_str.begin() )
99  << ">, and I built P=" << P << std::endl;
100  return 1;
101  }
102  // Durchblick x3y+xz3+y3z+z3+5z = 0
103 
104  // MPolynomial<3, double> P = mmonomial<double>( 3, 1, 0 )
105  // + mmonomial<double>( 1, 0, 3 )
106  // + mmonomial<double>( 0, 3, 1 )
107  // + mmonomial<double>( 0, 0, 3 )
108  // + 5 * mmonomial<double>( 0, 0, 1 );
109  // Crixxi (y2+z2-1)2 +(x2+y2-1)3 = 0
110  // developed = y4 +2y2z2+z4-2z2 -y2 + x6+3x4y2+3x2y4+y6-3x4-6x2y2-3y4+3x2
111  // MPolynomial<3, double> P = mmonomial<double>(0,4,0)
112  // + 2 * mmonomial<double>(0,2,2)
113  // + mmonomial<double>(0,2,0)
114  // + mmonomial<double>(0,0,4)
115  // - 2 * mmonomial<double>(0,0,2)
116  // + mmonomial<double>(6,0,0)
117  // + 3 * mmonomial<double>(4,2,0)
118  // + 3 * mmonomial<double>(2,4,0)
119  // + mmonomial<double>(0,6,0)
120  // - 3 * mmonomial<double>(4,0,0)
121  // - 6 * mmonomial<double>(2,2,0)
122  // - 3 * mmonomial<double>(0,4,0)
123  // + 3 * mmonomial<double>(2,0,0);
124 
125  trace.info() << "P( X_0, X_1, X_2 ) = " << P << std::endl;
126  ImplicitShape ishape( P );
127  DigitalShape dshape;
128  dshape.attach( ishape );
129  dshape.init( RealPoint( p1 ), RealPoint( p2 ), step );
130  Domain domain = dshape.getDomain();
131  trace.endBlock();
133 
135  // Construct the Khalimsky space from the image domain
136  KSpace K;
137  // NB: it is \b necessary to work with a \b closed cellular space
138  // since umbrellas use separators and pivots, which must exist for
139  // arbitrary surfels.
140  bool space_ok = K.init( domain.lowerBound(),
141  domain.upperBound(), true // necessary
142  );
143  if (!space_ok)
144  {
145  trace.error() << "Error in the Khamisky space construction."<<std::endl;
146  return 2;
147  }
149 
151  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
152  MySurfelAdjacency surfAdj( true ); // interior in all directions.
154 
156  trace.beginBlock( "Extracting boundary by tracking the space. " );
157  typedef KSpace::Surfel Surfel;
158  typedef KSpace::SurfelSet SurfelSet;
159  typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
161 
162  // The surfels are tracked from an initial surfel, which is found by
163  // try/error.
164  MySetOfSurfels theSetOfSurfels( K, surfAdj );
165  Surfel bel = Surfaces<KSpace>::findABel( K, dshape, 100000 );
166  Surfaces<KSpace>::trackBoundary( theSetOfSurfels.surfelSet(),
167  K, surfAdj,
168  dshape, bel );
169  trace.endBlock();
171 
173  trace.beginBlock( "Making OFF surface <marching-cube.off>. " );
174  MyDigitalSurface digSurf( theSetOfSurfels );
175  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
176  << std::endl;
177  // The cell embedder is used to place vertices closer to the set
178  // P(x,y,z)=0
179  typedef
181  ImplicitShape,
182  DigitalEmbedder >
183  CellEmbedder;
184  CellEmbedder cellEmbedder;
185  cellEmbedder.init( K, ishape, dshape.pointEmbedder() );
186  ofstream out( "marching-cube.off" );
187  if ( out.good() )
188  digSurf.exportEmbeddedSurfaceAs3DOFF( out, cellEmbedder );
189  out.close();
190  trace.endBlock();
192  trace.beginBlock( "Making NOFF surface <marching-cube.noff>. " );
193  ofstream out2( "marching-cube.noff" );
194  if ( out2.good() )
195  digSurf.exportEmbeddedSurfaceAs3DNOFF( out2, cellEmbedder );
196  out2.close();
197  trace.endBlock();
198 
199  return 0;
200 }
void beginBlock(const std::string &keyword="")
Iterator read(Polynomial &p, Iterator begin, Iterator end)
Aim: a cellular embedder for implicit functions, (default constructible, copy constructible, assignable). Model of CCellEmbedder and CWithGradientMap.
Trace trace
Definition: Common.h:137
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
Definition: Surfaces.h:78
KhalimskySpaceND< 2, Integer > KSpace
Definition: StdDefs.h:77
Aim: This class converts a string polynomial expression in a multivariate polynomial.
STL namespace.
double endBlock()
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:141
void attach(const EuclideanShape &shape)
void init(const RealPoint &xLow, const RealPoint &xUp, typename RealVector::Component gridStep)
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
Domain getDomain() const
Aim: model of CEuclideanOrientedShape concepts to create a shape from a polynomial.
const Point & upperBound() const
DGtal is the top-level namespace which contains all DGtal functions and types.
std::ostream & info()
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels...
Definition: SetOfSurfels.h:73
Space::RealPoint RealPoint
Definition: StdDefs.h:97
const Point & lowerBound() const
std::ostream & error()