DGtal 1.3.0
Searching...
No Matches
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;
typedef Space::Point Point;
typedef Space::RealPoint RealPoint;
typedef Space::RealVector RealVector;
typedef Space::Vector Vector;
typedef SMesh::Vertices Vertices;
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 );
trace.endBlock();
}
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,
SH3::Domain(K.lowerBound(),K.upperBound()),
params );
trace.endBlock();
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;
SH3::RealPoints vertices;
std::vector< Vertices > faces;
SH3::Cell2Index c2i;
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;
}
// //
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:105
Aim: A class that computes tangency to a given digital set. It provides services to compute all the c...
void beginBlock(const std::string &keyword="")
std::ostream & info()
double endBlock()
Z3i::SCell SCell
SurfaceMesh< RealPoint, RealVector > SMesh
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)
SurfaceMeshWriter< RealPoint, RealVector > SMeshWriter
SMesh::Vertices Vertices
SMesh::Index Index
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator > vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
STL namespace.
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Aim: An helper class for writing mesh file formats (Waverfront OBJ at this point) and creating a Surf...
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:92
std::size_t Index
The type used for numbering vertices and faces.
Definition: SurfaceMesh.h:105
Shortcuts< KSpace > SH3
Z2i::RealPoint RealPoint
int main()
Definition: testBits.cpp:56
MyPointD Point
Definition: testClone2.cpp:383
FreemanChain< int >::Vector Vector
KSpace K
HyperRectDomain< Space > Domain