DGtal  1.4.beta
geometry/volumes/fullConvexitySphereGeodesics.cpp

This example shows how to use tangency to compute exact or approximate geodesic shortest paths on 3D digital objects, here a unit sphere digitized at your chosen gridstep.

Approximated shortest paths to speed-up computation

For instance, you may call it with a digitization gridstep 0.0625 and parameter 1.8 (greater than sqrt(3), so guarantees exact shortest paths).

./examples/geometry/volumes/fullConvexitySphereGeodesics 0.0625 1.8


and outputs

Usage: ./examples/geometry/volumes/fullConvexitySphereGeodesics <h> <opt>
Computes shortest paths to a source point on a sphere digitized with gridstep <h>.
- h [==1.0]: digitization gridstep
- opt [==sqrt(3)]: >= sqrt(3): secure shortest paths, 0: fast
New Block [Building sphere1 shape ... ]
EndBlock [Building sphere1 shape ... ] (13.2603 ms)
New Block [Build mesh from primal surface]
#surfels =4782
#pointels=4784
[SurfaceMesh (OK) #V=4784 #VN=0 #E=9564 #F=4782 #FN=0 E[IF]=4 E[IV]=3.99833 E[IFE]=2]
EndBlock [Build mesh from primal surface] (111.024 ms)
New Block [Compute geodesics]
EndBlock [Compute geodesics] (11958.8 ms)
Max distance is 3.02389
Comput. time is 11958.8
Last index is   3698
Uppest index is 3698
New Block [Output Corrected HD OBJ files]
Max distance is 3.02389
EndBlock [Output Corrected HD OBJ files] (494.499 ms)


It computes all shortest paths to the lowest point of the digitized sphere, outputs the maximal distance (so less than pi), and checks that the furthest point is indeed antipodal to the source point. Two OBJ files named "sphere1-geodesics.obj "and "sphere1-geodesics-iso.obj" are also outputed, and you may use them to render the result like below.

 Exact geodesic distances on unit sphere digitized at h=0.0625 (max d=3.02389) Approximate geodesic distances on unit sphere digitized at h=0.01 (max d=3.1269)
#include <iostream>
#include <queue>
#include "DGtal/base/Common.h"
#include "DGtal/shapes/Shapes.h"
#include "DGtal/shapes/SurfaceMesh.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/helpers/Shortcuts.h"
#include "DGtal/images/ImageContainerBySTLVector.h"
#include "DGtal/io/writers/SurfaceMeshWriter.h"
#include "DGtal/geometry/volumes/TangencyComputer.h"
using namespace std;
using namespace DGtal;
typedef Z3i::Space Space;
typedef Z3i::SCell SCell;
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 )
{
std::string cobjname = output;
std::string cisoname = output + "-iso";
auto quantify = [] ( double v, double m, double nb )
{ return round( v/m*nb )*m/nb; };
trace.beginBlock( "Output Corrected HD OBJ files" );
const auto fvalues = surfmesh.computeFaceValuesFromVertexValues( vvalues );
double maxValue = * ( std::max_element( vvalues.cbegin(), vvalues.cend() ) );
double minValue = * ( std::min_element( vvalues.cbegin(), vvalues.cend() ) );
double maxDist = maxValue - minValue;
trace.info() << "Max distance is " << maxDist << std::endl;
auto cmap = SH3::getColorMap( 0.0, maxDist );
std::vector< Color > fcolors( surfmesh.nbFaces() );
for ( Index f = 0; f < fvalues.size(); ++f )
fcolors[ f ] = cmap( quantify( fvalues[ f ] - minValue, maxDist, 50 ) );
SMeshWriter::writeOBJ( cobjname, surfmesh, fcolors );
double unit = pow( 10.0, floor( log( maxDist ) / log( 10.0 ) ) - 1.0 );
const int N = 10 * nb_isolines_per_unit;
std::vector< double > isolines( N );
std::vector< Color > isocolors( N );
for ( int i = 0; i < N; i++ )
{
isolines [ i ] = (double) i * 10.0 * unit / (double) nb_isolines_per_unit
+ minValue;
isocolors[ i ] = ( i % nb_isolines_per_unit == 0 )
? Color::Red : Color::Black;
}
SMeshWriter::writeIsoLinesOBJ( cisoname, surfmesh, fvalues, vvalues,
isolines, thickness, isocolors );
}
int main( int argc, char** argv )
{
trace.info() << "Usage: " << argv[ 0 ] << " <h> <opt>" << std::endl;
trace.info() << "\tComputes shortest paths to a source point on a sphere digitized with gridstep <h>." << std::endl;
trace.info() << "\t- h [==1.0]: digitization gridstep" << std::endl;
trace.info() << "\t- opt [==sqrt(3)]: >= sqrt(3): secure shortest paths, 0: fast" << std::endl;
double h = argc > 1 ? atof( argv[ 1 ] ) : 0.0625; //< exact (sqrt(3)) or inexact (0) computations
double opt = argc > 2 ? atof( argv[ 2 ] ) : sqrt(3.0); //< exact (sqrt(3)) or inexact (0) computations
// Domain creation from two bounding points.
trace.beginBlock( "Building sphere1 shape ... " );
auto params = SH3::defaultParameters();
params( "polynomial", "sphere1" )( "gridstep", h );
params( "minAABB", -2)( "maxAABB", 2)( "offset", 1.0 )( "closed", 1 );
auto implicit_shape = SH3::makeImplicitShape3D ( params );
auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
auto K = SH3::getKSpace( params );
auto binary_image = SH3::makeBinaryImage(digitized_shape,
params );
trace.beginBlock( "Build mesh from primal surface" );
// Compute surface
auto surface = SH3::makeDigitalSurface( binary_image, K, params );
// Build a mesh
SMesh smesh;
auto embedder = SH3::getCellEmbedder( K );
std::vector< Point > lattice_points;
std::vector< Vertices > faces;
auto pointels = SH3::getPointelRange( c2i, surface );
for ( auto p : pointels ) lattice_points.push_back( K.uCoords( p ) );
trace.info() << "#surfels =" << surface->size() << std::endl;
trace.info() << "#pointels=" << pointels.size() << std::endl;
vertices = SH3::RealPoints( pointels.size() );
std::transform( pointels.cbegin(), pointels.cend(), vertices.begin(),
[&] (const SH3::Cell& c) { return h * embedder( c ); } );
// Build faces
for ( auto&& surfel : *surface )
{
const auto primal_surfel_vtcs = SH3::getPointelRange( K, surfel );
std::vector< Index > face;
for ( auto&& primal_vtx : primal_surfel_vtcs )
face.push_back( c2i[ primal_vtx ] );
faces.push_back( face );
}
smesh.init( vertices.cbegin(), vertices.cend(),
faces.cbegin(), faces.cend() );
trace.info() << smesh << std::endl;
// Find lowest and uppest point.
const Index nb = lattice_points.size();
Index lowest = 0;
Index uppest = 0;
for ( Index i = 1; i < nb; i++ )
{
if ( lattice_points[ i ] < lattice_points[ lowest ] ) lowest = i;
if ( lattice_points[ uppest ] < lattice_points[ i ] ) uppest = i;
}
// Extracts shortest paths to a target
trace.beginBlock( "Compute geodesics" );
TC.init( lattice_points.cbegin(), lattice_points.cend() );
auto SP = TC.makeShortestPaths( opt );
SP.init( lowest ); //< set source
double last_distance = 0.0;
Index last = 0;
while ( ! SP.finished() )
{
last = std::get<0>( SP.current() );
last_distance = std::get<2>( SP.current() );
SP.expand();
}
double time = trace.endBlock();
std::cout << "Max distance is " << last_distance*h << std::endl;
std::cout << "Comput. time is " << time << std::endl;
std::cout << "Last index is " << last << std::endl;
std::cout << "Uppest index is " << uppest << std::endl;
// Export surface for display
std::vector<double> distances = SP.distances();
for ( Index i = 0; i < distances.size(); i++ )
distances[ i ] *= h;
saveToObj( "sphere1-geodesics", smesh, distances, 10, 0.1 );
return 0;
}
// //
DGtal::TangencyComputer::Index
std::size_t Index
Definition: TangencyComputer.h:84
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
SH3
Shortcuts< KSpace > SH3
Definition: testArithmeticalDSSComputerOnSurfels.cpp:49
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:73
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="")
DGtal::KhalimskySpaceND::lowerBound
const Point & lowerBound() const
Return the lower bound for digital points in this space.
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
boost::vertices
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator > vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
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
KSpace
Z3i::KSpace KSpace
Definition: testArithmeticalDSSComputerOnSurfels.cpp:48
DGtal::Trace::info
std::ostream & info()
Vertices
SMesh::Vertices Vertices
Definition: fullConvexitySphereGeodesics.cpp:118
DGtal::Shortcuts::defaultParameters
static Parameters defaultParameters()
Definition: Shortcuts.h:203
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.
SMeshWriter
SurfaceMeshWriter< RealPoint, RealVector > SMeshWriter
Definition: fullConvexitySphereGeodesics.cpp:116
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
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.
DGtal::Shortcuts::getColorMap
static ColorMap getColorMap(Scalar min, Scalar max, const Parameters &params=parametersUtilities())
Definition: Shortcuts.h:2668
DGtal::Shortcuts::makeDigitizedImplicitShape3D
static CountedPtr< DigitizedImplicitShape3D > makeDigitizedImplicitShape3D(CountedPtr< ImplicitShape3D > shape, Parameters params=parametersDigitizedImplicitShape3D())
Definition: Shortcuts.h:523
Vector
FreemanChain< int >::Vector Vector
Definition: testCombinDSS.cpp:60
RealVector
Space::RealVector RealVector
Definition: digitalPolyhedronBuilder3D.cpp:82
DGtal::Shortcuts::RealPoints
std::vector< RealPoint > RealPoints
Definition: Shortcuts.h:180
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
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
DGtal::Shortcuts::getPointelRange
static PointelRange getPointelRange(Cell2Index &c2i, CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer > > surface)
Definition: Shortcuts.h:1469
DGtal::PointVector< dim, Integer >
Space
SpaceND< 2 > Space
Definition: testSimpleRandomAccessRangeFromPoint.cpp:42
SMesh
SurfaceMesh< RealPoint, RealVector > SMesh
Definition: fullConvexitySphereGeodesics.cpp:115
DGtal::Shortcuts::makeBinaryImage
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
SCell
Z3i::SCell SCell
Definition: digitalPolyhedronBuilder3D.cpp:79
Point
MyPointD Point
Definition: testClone2.cpp:383
RealPoint
Z2i::RealPoint RealPoint
Definition: testAstroid2D.cpp:46
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