DGtal  0.9.2
testPolynomial.cpp
1 
29 
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 
51 using namespace std;
52 using namespace DGtal;
53 using namespace Z3i;
54 
56 
57 
58 void 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 
78 int 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 
95  Polynomial3 P;
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;
132  typedef KSpace::SurfelSet SurfelSet;
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;
157  CanonicSCellEmbedder< KSpace > midpoint( K );
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 }
Iterator read(Polynomial &p, Iterator begin, Iterator end)
Aim: This class template may be used to (linearly) convert scalar values in a given range into a colo...
Trace trace
Definition: Common.h:130
Component Coordinate
Type for Point elements.
Definition: PointVector.h:158
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
Definition: Surfaces.h:78
std::set< SCell > SurfelSet
Preferred type for defining a set of surfels (always signed cells).
Aim: This class converts a string polynomial expression in a multivariate polynomial.
STL namespace.
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:141
void attach(const EuclideanShape &shape)
std::string className() const
void init(const RealPoint &xLow, const RealPoint &xUp, typename RealVector::Component gridStep)
bool init(const Point &lower, const Point &upper, bool isClosed)
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()
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value...
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels...
Definition: SetOfSurfels.h:73
Structure representing an RGB triple with alpha component.
Definition: Color.h:66
const Point & lowerBound() const
std::ostream & error()
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex...
Modifier class in a Display3D stream. Useful to choose your own mode for a given class. Realizes the concept CDrawableWithDisplay3D.