DGtal  1.3.beta
fullConvexitySphereGeodesics.cpp
Go to the documentation of this file.
1 
90 #include <iostream>
92 #include <queue>
93 #include "DGtal/base/Common.h"
94 #include "DGtal/shapes/Shapes.h"
95 #include "DGtal/shapes/SurfaceMesh.h"
96 #include "DGtal/helpers/StdDefs.h"
97 #include "DGtal/helpers/Shortcuts.h"
98 #include "DGtal/images/ImageContainerBySTLVector.h"
99 #include "DGtal/io/writers/SurfaceMeshWriter.h"
100 #include "DGtal/geometry/volumes/TangencyComputer.h"
101 
103 
104 using namespace std;
105 using namespace DGtal;
119 
120 void saveToObj( const std::string& output,
121  const SMesh& surfmesh,
122  const std::vector< double >& vvalues,
123  int nb_isolines_per_unit = 10,
124  const double thickness = 0.1 )
125 {
126  std::string cobjname = output;
127  std::string cisoname = output + "-iso";
128  auto quantify = [] ( double v, double m, double nb )
129  { return round( v/m*nb )*m/nb; };
130  trace.beginBlock( "Output Corrected HD OBJ files" );
131  const auto fvalues = surfmesh.computeFaceValuesFromVertexValues( vvalues );
132  double maxValue = * ( std::max_element( vvalues.cbegin(), vvalues.cend() ) );
133  double minValue = * ( std::min_element( vvalues.cbegin(), vvalues.cend() ) );
134  double maxDist = maxValue - minValue;
135  trace.info() << "Max distance is " << maxDist << std::endl;
136  auto cmap = SH3::getColorMap( 0.0, maxDist );
137  std::vector< Color > fcolors( surfmesh.nbFaces() );
138  for ( Index f = 0; f < fvalues.size(); ++f )
139  fcolors[ f ] = cmap( quantify( fvalues[ f ] - minValue, maxDist, 50 ) );
140  SMeshWriter::writeOBJ( cobjname, surfmesh, fcolors );
141  double unit = pow( 10.0, floor( log( maxDist ) / log( 10.0 ) ) - 1.0 );
142  const int N = 10 * nb_isolines_per_unit;
143  std::vector< double > isolines( N );
144  std::vector< Color > isocolors( N );
145  for ( int i = 0; i < N; i++ )
146  {
147  isolines [ i ] = (double) i * 10.0 * unit / (double) nb_isolines_per_unit
148  + minValue;
149  isocolors[ i ] = ( i % nb_isolines_per_unit == 0 )
150  ? Color::Red : Color::Black;
151  }
152  SMeshWriter::writeIsoLinesOBJ( cisoname, surfmesh, fvalues, vvalues,
153  isolines, thickness, isocolors );
154  trace.endBlock();
155 }
156 
157 
158 int main( int argc, char** argv )
159 {
160  trace.info() << "Usage: " << argv[ 0 ] << " <h> <opt>" << std::endl;
161  trace.info() << "\tComputes shortest paths to a source point on a sphere digitized with gridstep <h>." << std::endl;
162  trace.info() << "\t- h [==1.0]: digitization gridstep" << std::endl;
163  trace.info() << "\t- opt [==sqrt(3)]: >= sqrt(3): secure shortest paths, 0: fast" << std::endl;
164  double h = argc > 1 ? atof( argv[ 1 ] ) : 0.0625; //< exact (sqrt(3)) or inexact (0) computations
165  double opt = argc > 2 ? atof( argv[ 2 ] ) : sqrt(3.0); //< exact (sqrt(3)) or inexact (0) computations
166 
167  // Domain creation from two bounding points.
168  trace.beginBlock( "Building sphere1 shape ... " );
169  auto params = SH3::defaultParameters();
170  params( "polynomial", "sphere1" )( "gridstep", h );
171  params( "minAABB", -2)( "maxAABB", 2)( "offset", 1.0 )( "closed", 1 );
172  auto implicit_shape = SH3::makeImplicitShape3D ( params );
173  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
174  auto K = SH3::getKSpace( params );
175  auto binary_image = SH3::makeBinaryImage(digitized_shape,
177  params );
178  trace.endBlock();
179 
180  trace.beginBlock( "Build mesh from primal surface" );
181  // Compute surface
182  auto surface = SH3::makeDigitalSurface( binary_image, K, params );
183  // Build a mesh
184  SMesh smesh;
185  auto embedder = SH3::getCellEmbedder( K );
186  std::vector< Point > lattice_points;
187  SH3::RealPoints vertices;
188  std::vector< Vertices > faces;
189  SH3::Cell2Index c2i;
190  auto pointels = SH3::getPointelRange( c2i, surface );
191  for ( auto p : pointels ) lattice_points.push_back( K.uCoords( p ) );
192  trace.info() << "#surfels =" << surface->size() << std::endl;
193  trace.info() << "#pointels=" << pointels.size() << std::endl;
194  vertices = SH3::RealPoints( pointels.size() );
195  std::transform( pointels.cbegin(), pointels.cend(), vertices.begin(),
196  [&] (const SH3::Cell& c) { return h * embedder( c ); } );
197  // Build faces
198  for ( auto&& surfel : *surface )
199  {
200  const auto primal_surfel_vtcs = SH3::getPointelRange( K, surfel );
201  std::vector< Index > face;
202  for ( auto&& primal_vtx : primal_surfel_vtcs )
203  face.push_back( c2i[ primal_vtx ] );
204  faces.push_back( face );
205  }
206  smesh.init( vertices.cbegin(), vertices.cend(),
207  faces.cbegin(), faces.cend() );
208  trace.info() << smesh << std::endl;
209  trace.endBlock();
210 
211  // Find lowest and uppest point.
212  const Index nb = lattice_points.size();
213  Index lowest = 0;
214  Index uppest = 0;
215  for ( Index i = 1; i < nb; i++ )
216  {
217  if ( lattice_points[ i ] < lattice_points[ lowest ] ) lowest = i;
218  if ( lattice_points[ uppest ] < lattice_points[ i ] ) uppest = i;
219  }
220 
221  // Extracts shortest paths to a target
223  trace.beginBlock( "Compute geodesics" );
225  TC.init( lattice_points.cbegin(), lattice_points.cend() );
226  auto SP = TC.makeShortestPaths( opt );
227  SP.init( lowest ); //< set source
228  double last_distance = 0.0;
229  Index last = 0;
230  while ( ! SP.finished() )
231  {
232  last = std::get<0>( SP.current() );
233  last_distance = std::get<2>( SP.current() );
234  SP.expand();
235  }
236  double time = trace.endBlock();
237  std::cout << "Max distance is " << last_distance*h << std::endl;
238  std::cout << "Comput. time is " << time << std::endl;
239  std::cout << "Last index is " << last << std::endl;
240  std::cout << "Uppest index is " << uppest << std::endl;
241 
242  // Export surface for display
243  std::vector<double> distances = SP.distances();
244  for ( Index i = 0; i < distances.size(); i++ )
245  distances[ i ] *= h;
246  saveToObj( "sphere1-geodesics", smesh, distances, 10, 0.1 );
247  return 0;
248 }
249 // //
251 
DGtal::TangencyComputer::Index
std::size_t Index
Definition: TangencyComputer.h:83
DGtal::Shortcuts::Cell2Index
std::map< Cell, IdxVertex > Cell2Index
Definition: Shortcuts.h:189
DGtal::KhalimskySpaceND::upperBound
const Point & upperBound() const
Return the upper bound for digital points in this space.
DGtal::HyperRectDomain< Space >
DGtal::Trace::endBlock
double endBlock()
DGtal::Shortcuts::getCellEmbedder
static CanonicCellEmbedder< KSpace > getCellEmbedder(const KSpace &K)
Definition: Shortcuts.h:438
DGtal::Shortcuts::Cell
LightDigitalSurface::Cell Cell
Definition: Shortcuts.h:162
Index
SMesh::Index Index
Definition: fullConvexitySphereGeodesics.cpp:117
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::TangencyComputer
Aim: A class that computes tangency to a given digital set. It provides services to compute all the c...
Definition: TangencyComputer.h:72
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::SurfaceMesh
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:91
DGtal::SurfaceMeshWriter
Aim: An helper class for writing mesh file formats (Waverfront OBJ at this point) and creating a Surf...
Definition: SurfaceMeshWriter.h:64
DGtal::Shortcuts::makeDigitalSurface
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
main
int main(int argc, char **argv)
Definition: fullConvexitySphereGeodesics.cpp:158
DGtal::KhalimskySpaceND::lowerBound
const Point & lowerBound() const
Return the lower bound for digital points in this space.
KSpace
Z3i::KSpace KSpace
Definition: fullConvexitySphereGeodesics.cpp:107
DGtal::Shortcuts
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:104
DGtal::SignedKhalimskyCell
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Definition: KhalimskySpaceND.h:208
DGtal::TangencyComputer::makeShortestPaths
ShortestPaths makeShortestPaths(double secure=sqrt(KSpace::dimension)) const
DGtal::SpaceND
Definition: SpaceND.h:95
DGtal::SurfaceMesh::Vertices
std::vector< Vertex > Vertices
The type that defines a list/range of vertices (e.g. to define faces)
Definition: SurfaceMesh.h:112
DGtal::Trace::info
std::ostream & info()
Vertices
SMesh::Vertices Vertices
Definition: fullConvexitySphereGeodesics.cpp:118
Space
Z3i::Space Space
Definition: fullConvexitySphereGeodesics.cpp:106
DGtal::Shortcuts::defaultParameters
static Parameters defaultParameters()
Definition: Shortcuts.h:203
Vector
Space::Vector Vector
Definition: fullConvexitySphereGeodesics.cpp:114
DGtal::Shortcuts::makeImplicitShape3D
static CountedPtr< ImplicitShape3D > makeImplicitShape3D(const Parameters &params=parametersImplicitShape3D())
Definition: Shortcuts.h:282
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
SCell
Z3i::SCell SCell
Definition: fullConvexitySphereGeodesics.cpp:109
SMeshWriter
SurfaceMeshWriter< RealPoint, RealVector > SMeshWriter
Definition: fullConvexitySphereGeodesics.cpp:116
saveToObj
void saveToObj(const std::string &output, const SMesh &surfmesh, const std::vector< double > &vvalues, int nb_isolines_per_unit=10, const double thickness=0.1)
Definition: fullConvexitySphereGeodesics.cpp:120
DGtal::KhalimskySpaceND::uCoords
Point uCoords(const Cell &c) const
Return its digital coordinates.
SH3
Shortcuts< KSpace > SH3
Definition: fullConvexitySphereGeodesics.cpp:110
DGtal::Shortcuts::getColorMap
static ColorMap getColorMap(Scalar min, Scalar max, const Parameters &params=parametersUtilities())
Definition: Shortcuts.h:2668
RealPoint
Space::RealPoint RealPoint
Definition: fullConvexitySphereGeodesics.cpp:112
DGtal::Shortcuts::makeDigitizedImplicitShape3D
static CountedPtr< DigitizedImplicitShape3D > makeDigitizedImplicitShape3D(CountedPtr< ImplicitShape3D > shape, Parameters params=parametersDigitizedImplicitShape3D())
Definition: Shortcuts.h:523
DGtal::TangencyComputer::init
void init(PointIterator itB, PointIterator itE)
DGtal::Shortcuts::RealPoints
std::vector< RealPoint > RealPoints
Definition: Shortcuts.h:180
Domain
Z3i::Domain Domain
Definition: fullConvexitySphereGeodesics.cpp:108
Point
Space::Point Point
Definition: fullConvexitySphereGeodesics.cpp:111
DGtal::SurfaceMesh::computeFaceValuesFromVertexValues
std::vector< AnyRing > computeFaceValuesFromVertexValues(const std::vector< AnyRing > &vvalues) const
DGtal::Shortcuts::getKSpace
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
RealVector
Space::RealVector RealVector
Definition: fullConvexitySphereGeodesics.cpp:113
DGtal::Shortcuts::getPointelRange
static PointelRange getPointelRange(Cell2Index &c2i, CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer > > surface)
Definition: Shortcuts.h:1469
DGtal::PointVector< dim, Integer >
SMesh
SurfaceMesh< RealPoint, RealVector > SMesh
Definition: fullConvexitySphereGeodesics.cpp:115
DGtal::Shortcuts::makeBinaryImage
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
DGtal::SurfaceMesh::init
bool init(RealPointIterator itPos, RealPointIterator itPosEnd, VerticesIterator itVertices, VerticesIterator itVerticesEnd)
DGtal::SurfaceMesh::nbFaces
Size nbFaces() const
Definition: SurfaceMesh.h:288
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:64
DGtal::SurfaceMesh::Index
std::size_t Index
The type used for numbering vertices and faces.
Definition: SurfaceMesh.h:105