DGtal 1.3.0
Loading...
Searching...
No Matches
testPolynomial.cpp
1
30
31
32#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"
42#include "DGtal/topology/SCellsFunctors.h"
43#include "DGtal/topology/helpers/BoundaryPredicate.h"
44#include "DGtal/topology/SetOfSurfels.h"
45#include "DGtal/io/viewers/Viewer3D.h"
46#include "DGtal/io/colormaps/GradientColorMap.h"
47#include <boost/math/special_functions/fpclassify.hpp>
48
50
51using namespace std;
52using namespace DGtal;
53using namespace Z3i;
54
56
57
58void usage( int /*argc*/, char** argv )
59{
60 std::cerr << "Usage: " << argv[ 0 ] << " <Polynomial> <Px> <Py> <Pz> <Qx> <Qy> <Qz> <step>" << std::endl;
61 std::cerr << "\t - displays the boundary of a shape defined implicitly by a 3-polynomial <Polynomial>." << std::endl;
62 std::cerr << "\t - P and Q defines the bounding box." << std::endl;
63 std::cerr << "\t - step is the grid step." << std::endl;
64 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;
65 std::cerr << "\t - See http://www.freigeist.cc/gallery.html" << std::endl;
66}
67
69 typedef RealPointT::Coordinate Ring;
75
76
77
78int main( int argc, char** argv )
79{
80 if ( argc < 9 )
81 {
82 usage( argc, argv );
83 return 1;
84 }
85 double p1[ 3 ];
86 double p2[ 3 ];
87 for ( unsigned int i = 0; i < 3; ++i )
88 {
89 p1[ i ] = atof( argv[ 2 + i ] );
90 p2[ i ] = atof( argv[ 5 + i ] );
91 }
92 double step = atof( argv[ 8 ] );
93
94
96 Polynomial3Reader reader;
97 std::string poly_str = argv[ 1 ];
98 std::string::const_iterator iter
99 = reader.read( P, poly_str.begin(), poly_str.end() );
100 if ( iter != poly_str.end() )
101 {
102 std::cerr << "ERROR: I read only <"
103 << poly_str.substr( 0, iter - poly_str.begin() )
104 << ">, and I built P=" << P << std::endl;
105 return 1;
106 }
107
108
109 ImplicitShape ishape( P );
110 DigitalShape dshape;
111 dshape.attach( ishape );
112 dshape.init( RealPointT( p1 ), RealPointT( p2 ), step );
113 Domain domain = dshape.getDomain();
114
115
116 KSpace K;
117
118 bool space_ok = K.init( domain.lowerBound(),
119 domain.upperBound(), true
120 );
121 if ( !space_ok )
122 {
123 trace.error() << "Error in the Khamisky space construction." << std::endl;
124 return 2;
125 }
126
127 typedef SurfelAdjacency< KSpace::dimension > MySurfelAdjacency;
128 MySurfelAdjacency surfAdj( true ); // interior in all directions.
129
130
131 typedef KSpace::Surfel Surfel;
133 typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
134
135 MySetOfSurfels theSetOfSurfels( K, surfAdj );
136 Surfel bel = Surfaces< KSpace >::findABel( K, dshape, 100000 );
137 Surfaces< KSpace >::trackBoundary( theSetOfSurfels.surfelSet(),
138 K, surfAdj,
139 dshape, bel );
140
141
142
143
144 QApplication application( argc, argv );
145 Viewer3D<> viewer;
146 viewer.show();
147 viewer << SetMode3D( domain.className(), "BoundingBox" ) << domain;
148
149
150
151
152 //-----------------------------------------------------------------------
153 // Looking for the min and max values
154
155 double minCurv = 1;
156 double maxCurv = 0;
158 for ( std::set< SCell >::iterator it = theSetOfSurfels.begin(), it_end = theSetOfSurfels.end();
159 it != it_end; ++it)
160 {
161
162 RealPointT A = midpoint( *it ) * step;
163 A = ishape.nearestPoint (A,0.01,200,0.1*step);
164 double a = ishape.meanCurvature( A );
165// double a=ishape.gaussianCurvature(A);
166 if ( !boost::math::isnan( a ))
167 {
168
169 if ( a > maxCurv )
170 {
171 maxCurv = a;
172 }
173 else
174 if ( a < minCurv )
175 {
176 minCurv = a;
177 }
178 }
179 }
180 trace.info() << " Min = " << minCurv << std::endl;
181 trace.info() << " Max = " << maxCurv << std::endl;
182
183
184 //-----------------------------------------------------------------------
185 //Specifing a color map
186
187 GradientColorMap< double > cmap_grad( minCurv, maxCurv );
188 cmap_grad.addColor( Color( 50, 50, 255 ) );
189 cmap_grad.addColor( Color( 255, 0, 0 ) );
190 cmap_grad.addColor( Color( 255, 255, 10 ) );
191
192 //------------------------------------------------------------------------------------
193 //drawing
194 unsigned int nbSurfels = 0;
195
196 for ( std::set<SCell>::iterator it = theSetOfSurfels.begin(),
197 it_end = theSetOfSurfels.end();
198 it != it_end; ++it, ++nbSurfels )
199 {
200
201
202 RealPointT A = midpoint( *it ) * step;
203 A = ishape.nearestPoint (A,0.01,200,0.1*step);
204// double a=ishape.gaussianCurvature(A);
205 double a = ishape.meanCurvature( A );
206 if ( boost::math::isnan( a ))
207 {
208 a = 0;
209 }
210
211 viewer << CustomColors3D( Color::Black, cmap_grad( a ));
212 viewer << *it;
213 }
214
215 viewer << Viewer3D<>::updateDisplay;
216
217 return application.exec();
218}
Structure representing an RGB triple with alpha component.
Definition: Color.h:68
static const Color Black
Definition: Color.h:413
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
void attach(ConstAlias< EuclideanShape > shape)
void init(const RealPoint &xLow, const RealPoint &xUp, typename RealVector::Component gridStep)
Domain getDomain() const
Aim: This class template may be used to (linearly) convert scalar values in a given range into a colo...
const Point & lowerBound() const
const Point & upperBound() const
std::string className() const
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
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)
std::ostream & error()
std::ostream & info()
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
MyDigitalSurface::SurfelSet SurfelSet
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
STL namespace.
Modifier class in a Display3D stream. Useful to choose your own mode for a given class....
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