DGtal  1.3.beta
topology/trackImplicitPolynomialSurfaceToOFF.cpp

Implicit polynomial surface defined on the command-line by the user (as "(x^2+y^2+2z^2-1)*(z^2x-0.1)"), then extracted using digital surface tracking and converted into the corresponding combinatorial surface.

See also
Input and output for multivariate polynomials
Application to export surface in OFF format
* $ ./examples/topology/trackImplicitPolynomialSurfaceToOFF "(x^2+y^2+2*z^2-1)*(z^2x-0.1)" -2 -2 -2 2 2 2 0.02
* # Digital surface has 112826 surfels.
* # output in marching-cube.off (in 1809ms)
* # You may display it with your favorite OFF displayer (like geomview, etc).
* $ ctmviewer marching-cube.off
* 
Implicit polynomial surface (x^2+y^2+2*z^2-1)*(z^2x-0.1) between [-2,-2,-2] and [2,2,2], step 0.02.
#include <iostream>
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/topology/helpers/Surfaces.h"
#include "DGtal/topology/DigitalSurface.h"
#include "DGtal/topology/SetOfSurfels.h"
#include "DGtal/math/MPolynomial.h"
#include "DGtal/shapes/GaussDigitizer.h"
#include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
#include "DGtal/shapes/implicit/ImplicitFunctionDiff1LinearCellEmbedder.h"
#include "DGtal/io/readers/MPolynomialReader.h"
using namespace std;
using namespace DGtal;
using namespace Z3i;
void usage( int /*argc*/, char** argv )
{
std::cerr << "Usage: " << argv[ 0 ] << " <Polynomial> <Px> <Py> <Pz> <Qx> <Qy> <Qz> <step>" << std::endl;
std::cerr << "\t - displays the boundary of a shape defined implicitly by a 3-polynomial <Polynomial>." << std::endl;
std::cerr << "\t - P and Q defines the bounding box." << std::endl;
std::cerr << "\t - step is the grid step." << std::endl;
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;
std::cerr << "\t - See http://www.freigeist.cc/gallery.html" << std::endl;
}
int main( int argc, char** argv )
{
if ( argc < 9 )
{
usage( argc, argv );
return 1;
}
double p1[ 3 ];
double p2[ 3 ];
for ( unsigned int i = 0; i < 3; ++i )
{
p1[ i ] = atof( argv[ 2+i ] );
p2[ i ] = atof( argv[ 5+i ] );
}
double step = atof( argv[ 8 ] );
trace.beginBlock( "Making polynomial surface." );
typedef RealPoint::Coordinate Ring;
// See http://www.freigeist.cc/gallery.html
std::string poly_str = argv[ 1 ];
std::string::const_iterator iter
= reader.read( P, poly_str.begin(), poly_str.end() );
if ( iter != poly_str.end() )
{
std::cerr << "ERROR: I read only <"
<< poly_str.substr( 0, iter - poly_str.begin() )
<< ">, and I built P=" << P << std::endl;
return 1;
}
// Durchblick x3y+xz3+y3z+z3+5z = 0
// MPolynomial<3, double> P = mmonomial<double>( 3, 1, 0 )
// + mmonomial<double>( 1, 0, 3 )
// + mmonomial<double>( 0, 3, 1 )
// + mmonomial<double>( 0, 0, 3 )
// + 5 * mmonomial<double>( 0, 0, 1 );
// Crixxi (y2+z2-1)2 +(x2+y2-1)3 = 0
// developed = y4 +2y2z2+z4-2z2 -y2 + x6+3x4y2+3x2y4+y6-3x4-6x2y2-3y4+3x2
// MPolynomial<3, double> P = mmonomial<double>(0,4,0)
// + 2 * mmonomial<double>(0,2,2)
// + mmonomial<double>(0,2,0)
// + mmonomial<double>(0,0,4)
// - 2 * mmonomial<double>(0,0,2)
// + mmonomial<double>(6,0,0)
// + 3 * mmonomial<double>(4,2,0)
// + 3 * mmonomial<double>(2,4,0)
// + mmonomial<double>(0,6,0)
// - 3 * mmonomial<double>(4,0,0)
// - 6 * mmonomial<double>(2,2,0)
// - 3 * mmonomial<double>(0,4,0)
// + 3 * mmonomial<double>(2,0,0);
trace.info() << "P( X_0, X_1, X_2 ) = " << P << std::endl;
ImplicitShape ishape( P );
DigitalShape dshape;
dshape.attach( ishape );
dshape.init( RealPoint( p1 ), RealPoint( p2 ), step );
Domain domain = dshape.getDomain();
// Construct the Khalimsky space from the image domain
// NB: it is \b necessary to work with a \b closed cellular space
// since umbrellas use separators and pivots, which must exist for
// arbitrary surfels.
bool space_ok = K.init( domain.lowerBound(),
domain.upperBound(), true // necessary
);
if (!space_ok)
{
trace.error() << "Error in the Khamisky space construction."<<std::endl;
return 2;
}
typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
MySurfelAdjacency surfAdj( true ); // interior in all directions.
trace.beginBlock( "Extracting boundary by tracking the space. " );
typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
// The surfels are tracked from an initial surfel, which is found by
// try/error.
MySetOfSurfels theSetOfSurfels( K, surfAdj );
Surfel bel = Surfaces<KSpace>::findABel( K, dshape, 100000 );
Surfaces<KSpace>::trackBoundary( theSetOfSurfels.surfelSet(),
K, surfAdj,
dshape, bel );
trace.beginBlock( "Making OFF surface <marching-cube.off>. " );
MyDigitalSurface digSurf( theSetOfSurfels );
trace.info() << "Digital surface has " << digSurf.size() << " surfels."
<< std::endl;
// The cell embedder is used to place vertices closer to the set
// P(x,y,z)=0
typedef
CellEmbedder;
CellEmbedder cellEmbedder;
cellEmbedder.init( K, ishape, dshape.pointEmbedder() );
ofstream out( "marching-cube.off" );
if ( out.good() )
digSurf.exportEmbeddedSurfaceAs3DOFF( out, cellEmbedder );
out.close();
trace.beginBlock( "Making NOFF surface <marching-cube.noff>. " );
ofstream out2( "marching-cube.noff" );
if ( out2.good() )
digSurf.exportEmbeddedSurfaceAs3DNOFF( out2, cellEmbedder );
out2.close();
return 0;
}
DGtal::GaussDigitizer
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
Definition: GaussDigitizer.h:79
DGtal::Surfaces
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
Definition: Surfaces.h:78
DGtal::GaussDigitizer::getDomain
Domain getDomain() const
DGtal::HyperRectDomain< Space >
DGtal::Trace::endBlock
double endBlock()
DGtal::KhalimskySpaceND::SurfelSet
std::set< SCell > SurfelSet
Preferred type for defining a set of surfels (always signed cells).
Definition: KhalimskySpaceND.h:450
DGtal::PointVector::Coordinate
Component Coordinate
Type for Point elements.
Definition: PointVector.h:617
DGtal::DigitalSurface
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Definition: DigitalSurface.h:139
DGtal::SurfelAdjacency< KSpace::dimension >
DGtal::Trace::error
std::ostream & error()
DGtal::HyperRectDomain::upperBound
const Point & upperBound() const
DGtal::KhalimskySpaceND::init
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
DGtal::RegularPointEmbedder::init
void init(typename RealVector::Component gridStep)
DGtal::GaussDigitizer::pointEmbedder
const PointEmbedder & pointEmbedder() const
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::GaussDigitizer::attach
void attach(ConstAlias< EuclideanShape > shape)
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::SignedKhalimskyCell
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Definition: KhalimskySpaceND.h:208
DGtal::RegularPointEmbedder< Space >
DGtal::MPolynomialReader::read
Iterator read(Polynomial &p, Iterator begin, Iterator end)
Definition: MPolynomialReader.h:284
DGtal::GaussDigitizer::init
void init(const RealPoint &xLow, const RealPoint &xUp, typename RealVector::Component gridStep)
KSpace
Z3i::KSpace KSpace
Definition: testArithmeticalDSSComputerOnSurfels.cpp:48
DGtal::Trace::info
std::ostream & info()
Surfel
KSpace::SCell Surfel
Definition: testArithmeticalDSSComputerOnSurfels.cpp:50
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::ImplicitPolynomial3Shape
Aim: model of CEuclideanOrientedShape concepts to create a shape from a polynomial.
Definition: ImplicitPolynomial3Shape.h:67
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
DGtal::MPolynomialReader
Aim: This class converts a string polynomial expression in a multivariate polynomial.
Definition: MPolynomialReader.h:258
usage
void usage(int, char **argv)
Definition: approximation.cpp:64
DGtal::ImplicitFunctionDiff1LinearCellEmbedder
Aim: a cellular embedder for implicit functions, (default constructible, copy constructible,...
Definition: ImplicitFunctionDiff1LinearCellEmbedder.h:79
domain
Domain domain
Definition: testProjection.cpp:88
DGtal::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:165
DGtal::SetOfSurfels
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels....
Definition: SetOfSurfels.h:73
DGtal::MPolynomial< 3, Ring >
RealPoint
Z2i::RealPoint RealPoint
Definition: testAstroid2D.cpp:46
DGtal::HyperRectDomain::lowerBound
const Point & lowerBound() const
MyDigitalSurface
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
Definition: greedy-plane-segmentation-ex2.cpp:92
SurfelSet
MyDigitalSurface::SurfelSet SurfelSet
Definition: greedy-plane-segmentation-ex2.cpp:95
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:64