DGtal 1.3.0
Loading...
Searching...
No Matches
trackImplicitPolynomialSurfaceToOFF.cpp
Go to the documentation of this file.
1
57#include <iostream>
58#include "DGtal/helpers/StdDefs.h"
59#include "DGtal/topology/helpers/Surfaces.h"
60#include "DGtal/topology/DigitalSurface.h"
61#include "DGtal/topology/SetOfSurfels.h"
62#include "DGtal/math/MPolynomial.h"
63#include "DGtal/shapes/GaussDigitizer.h"
64#include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
65#include "DGtal/shapes/implicit/ImplicitFunctionDiff1LinearCellEmbedder.h"
66#include "DGtal/io/readers/MPolynomialReader.h"
68
70
71using namespace std;
72using namespace DGtal;
73using namespace Z3i;
74
76
77
78void usage( int /*argc*/, char** argv )
79{
80 std::cerr << "Usage: " << argv[ 0 ] << " <Polynomial> <Px> <Py> <Pz> <Qx> <Qy> <Qz> <step>" << std::endl;
81 std::cerr << "\t - displays the boundary of a shape defined implicitly by a 3-polynomial <Polynomial>." << std::endl;
82 std::cerr << "\t - P and Q defines the bounding box." << std::endl;
83 std::cerr << "\t - step is the grid step." << std::endl;
84 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;
85 std::cerr << "\t - See http://www.freigeist.cc/gallery.html" << std::endl;
86}
87
88int main( int argc, char** argv )
89{
90 if ( argc < 9 )
91 {
92 usage( argc, argv );
93 return 1;
94 }
95 double p1[ 3 ];
96 double p2[ 3 ];
97 for ( unsigned int i = 0; i < 3; ++i )
98 {
99 p1[ i ] = atof( argv[ 2+i ] );
100 p2[ i ] = atof( argv[ 5+i ] );
101 }
102 double step = atof( argv[ 8 ] );
103
105 trace.beginBlock( "Making polynomial surface." );
107 typedef RealPoint::Coordinate Ring;
113
114 // See http://www.freigeist.cc/gallery.html
115 Polynomial3 P;
116 Polynomial3Reader reader;
117 std::string poly_str = argv[ 1 ];
118 std::string::const_iterator iter
119 = reader.read( P, poly_str.begin(), poly_str.end() );
120 if ( iter != poly_str.end() )
121 {
122 std::cerr << "ERROR: I read only <"
123 << poly_str.substr( 0, iter - poly_str.begin() )
124 << ">, and I built P=" << P << std::endl;
125 return 1;
126 }
127 // Durchblick x3y+xz3+y3z+z3+5z = 0
128
129 // MPolynomial<3, double> P = mmonomial<double>( 3, 1, 0 )
130 // + mmonomial<double>( 1, 0, 3 )
131 // + mmonomial<double>( 0, 3, 1 )
132 // + mmonomial<double>( 0, 0, 3 )
133 // + 5 * mmonomial<double>( 0, 0, 1 );
134 // Crixxi (y2+z2-1)2 +(x2+y2-1)3 = 0
135 // developed = y4 +2y2z2+z4-2z2 -y2 + x6+3x4y2+3x2y4+y6-3x4-6x2y2-3y4+3x2
136 // MPolynomial<3, double> P = mmonomial<double>(0,4,0)
137 // + 2 * mmonomial<double>(0,2,2)
138 // + mmonomial<double>(0,2,0)
139 // + mmonomial<double>(0,0,4)
140 // - 2 * mmonomial<double>(0,0,2)
141 // + mmonomial<double>(6,0,0)
142 // + 3 * mmonomial<double>(4,2,0)
143 // + 3 * mmonomial<double>(2,4,0)
144 // + mmonomial<double>(0,6,0)
145 // - 3 * mmonomial<double>(4,0,0)
146 // - 6 * mmonomial<double>(2,2,0)
147 // - 3 * mmonomial<double>(0,4,0)
148 // + 3 * mmonomial<double>(2,0,0);
149
150 trace.info() << "P( X_0, X_1, X_2 ) = " << P << std::endl;
151 ImplicitShape ishape( P );
152 DigitalShape dshape;
153 dshape.attach( ishape );
154 dshape.init( RealPoint( p1 ), RealPoint( p2 ), step );
155 Domain domain = dshape.getDomain();
156 trace.endBlock();
158
160 // Construct the Khalimsky space from the image domain
161 KSpace K;
162 // NB: it is \b necessary to work with a \b closed cellular space
163 // since umbrellas use separators and pivots, which must exist for
164 // arbitrary surfels.
165 bool space_ok = K.init( domain.lowerBound(),
166 domain.upperBound(), true // necessary
167 );
168 if (!space_ok)
169 {
170 trace.error() << "Error in the Khamisky space construction."<<std::endl;
171 return 2;
172 }
174
176 typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
177 MySurfelAdjacency surfAdj( true ); // interior in all directions.
179
181 trace.beginBlock( "Extracting boundary by tracking the space. " );
182 typedef KSpace::Surfel Surfel;
184 typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
186
187 // The surfels are tracked from an initial surfel, which is found by
188 // try/error.
189 MySetOfSurfels theSetOfSurfels( K, surfAdj );
190 Surfel bel = Surfaces<KSpace>::findABel( K, dshape, 100000 );
191 Surfaces<KSpace>::trackBoundary( theSetOfSurfels.surfelSet(),
192 K, surfAdj,
193 dshape, bel );
194 trace.endBlock();
196
198 trace.beginBlock( "Making OFF surface <marching-cube.off>. " );
199 MyDigitalSurface digSurf( theSetOfSurfels );
200 trace.info() << "Digital surface has " << digSurf.size() << " surfels."
201 << std::endl;
202 // The cell embedder is used to place vertices closer to the set
203 // P(x,y,z)=0
204 typedef
208 CellEmbedder;
209 CellEmbedder cellEmbedder;
210 cellEmbedder.init( K, ishape, dshape.pointEmbedder() );
211 ofstream out( "marching-cube.off" );
212 if ( out.good() )
213 digSurf.exportEmbeddedSurfaceAs3DOFF( out, cellEmbedder );
214 out.close();
215 trace.endBlock();
217 trace.beginBlock( "Making NOFF surface <marching-cube.noff>. " );
218 ofstream out2( "marching-cube.noff" );
219 if ( out2.good() )
220 digSurf.exportEmbeddedSurfaceAs3DNOFF( out2, cellEmbedder );
221 out2.close();
222 trace.endBlock();
223
224 return 0;
225}
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
void exportEmbeddedSurfaceAs3DOFF(std::ostream &out, const CellEmbedder &cembedder) const
void exportEmbeddedSurfaceAs3DNOFF(std::ostream &out, const CellEmbedder &cembedder) const
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
void attach(ConstAlias< EuclideanShape > shape)
const PointEmbedder & pointEmbedder() const
void init(const RealPoint &xLow, const RealPoint &xUp, typename RealVector::Component gridStep)
Domain getDomain() const
const Point & lowerBound() const
const Point & upperBound() const
Aim: a cellular embedder for implicit functions, (default constructible, copy constructible,...
Aim: model of CEuclideanOrientedShape concepts to create a shape from a polynomial.
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
std::set< SCell > SurfelSet
Preferred type for defining a set of surfels (always signed cells).
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
Aim: This class converts a string polynomial expression in a multivariate polynomial.
Iterator read(Polynomial &p, Iterator begin, Iterator end)
Aim: Represents a multivariate polynomial, i.e. an element of , where K is some ring or field.
Definition: MPolynomial.h:955
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
Component Coordinate
Type for Point elements.
Definition: PointVector.h:617
void init(typename RealVector::Component gridStep)
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels....
Definition: SetOfSurfels.h:74
static void trackBoundary(SCellSet &surface, const KSpace &K, const SurfelAdjacency< KSpace::dimension > &surfel_adj, const PointPredicate &pp, const SCell &start_surfel)
static SCell findABel(const KSpace &K, const PointPredicate &pp, unsigned int nbtries=1000)
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
void beginBlock(const std::string &keyword="")
std::ostream & error()
std::ostream & info()
double endBlock()
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
MyDigitalSurface::SurfelSet SurfelSet
KhalimskySpaceND< 3, Integer > KSpace
Definition: StdDefs.h:146
Space::RealPoint RealPoint
Definition: StdDefs.h:170
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
STL namespace.
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
int main()
Definition: testBits.cpp:56
KSpace K
Domain domain