DGtal 2.0.0
Loading...
Searching...
No Matches
fullConvexityShortestPaths3D.cpp
Go to the documentation of this file.
1
60#include <iostream>
61#include <queue>
62#include "DGtal/base/Common.h"
63#include "DGtal/io/viewers/PolyscopeViewer.h"
64#include "DGtal/io/Color.h"
65#include "DGtal/io/colormaps/SimpleDistanceColorMap.h"
66#include "DGtal/shapes/Shapes.h"
67#include "DGtal/helpers/StdDefs.h"
68#include "DGtal/helpers/Shortcuts.h"
69#include "DGtal/images/ImageContainerBySTLVector.h"
71#include "DGtal/geometry/volumes/TangencyComputer.h"
73#include "ConfigExamples.h"
74
76
77using namespace std;
78using namespace DGtal;
87
88// Called when an user clicks on a surfel.
89int reaction( void* viewer, DGtal::int32_t name, void* data )
90{
91 ((void) viewer);
92
93 DGtal::int32_t* selected = (DGtal::int32_t*) data;
94 *selected = name;
95 std::cout << "Selected surfel=" << *selected << std::endl;
96 return 0;
97}
98
99int main( int argc, char** argv )
100{
101 trace.info() << "Usage: " << argv[ 0 ] << " <input.vol> <m> <M> <opt>" << std::endl;
102 trace.info() << "\tComputes shortest paths to a source point" << std::endl;
103 trace.info() << "\t- input.vol: choose your favorite shape" << std::endl;
104 trace.info() << "\t- m [==0], M [==255]: used to threshold input vol image" << std::endl;
105 trace.info() << "\t- opt >= sqrt(3): secure shortest paths, 0: fast" << std::endl;
106 string inputFilename = examplesPath + "samples/Al.100.vol";
107 std::string fn= argc > 1 ? argv[ 1 ] : inputFilename; //< vol filename
108 int m = argc > 2 ? atoi( argv[ 2 ] ) : 0; //< low for thresholding
109 int M = argc > 3 ? atoi( argv[ 3 ] ) : 255; //< up for thresholding
110 double opt = argc > 4 ? atof( argv[ 4 ] ) : sqrt(3.0); //< exact (sqrt(3)) or inexact (0) computations
111
112 PolyscopeViewer<> viewer;
113
114 // Set up shortcuts parameters.
115 auto params = SH3::defaultParameters();
116 params( "thresholdMin", m )( "thresholdMax", M );
117 params( "surfaceComponents" , "All" );
118
119 // Domain creation from two bounding points.
120 trace.info() << "Building set or importing vol ... ";
121 KSpace K;
122 auto bimage = SH3::makeBinaryImage( fn, params );
123 K = SH3::getKSpace( bimage );
124 trace.info() << " [Done]" << std::endl;
125
126 // Compute surface
127 auto surface = SH3::makeDigitalSurface( bimage, K, params );
128
129 // Compute interior boundary points
130 // They are less immediate interior points than surfels.
131 std::vector< Point > points;
132 std::map< SCell, int > surfel2idx;
133 std::map< Point, int > point2idx;
134 int idx = 0;
135 for ( auto s : (*surface) )
136 {
137 // get inside point on the border of the shape.
138 Dimension k = K.sOrthDir( s );
139 auto voxel = K.sIncident( s, k, K.sDirect( s, k ) );
140 Point p = K.sCoords( voxel );
141 auto it = point2idx.find( p );
142 if ( it == point2idx.end() )
143 {
144 points.push_back( p );
145 surfel2idx[ s ] = idx;
146 point2idx [ p ] = idx++;
147 }
148 else
149 surfel2idx[ s ] = it->second;
150 }
151 trace.info() << "Shape has " << points.size() << " interior boundary points"
152 << std::endl;
153
154 // Select a starting point.
155 DGtal::int32_t selected_surfels[ 2 ] = { 0, 0 };
156 typedef PolyscopeViewer<> MViewer3D;
157 auto surfels = SH3::getSurfelRange ( surface );
158 for ( int i = 0; i < 2; i++ )
159 {
160 MViewer3D viewerCore( K );
161 Color colSurfel( 200, 200, 255, 255 );
162 Color colStart( 255, 0, 0, 255 );
163 viewerCore.drawColor( colSurfel );
164 for ( auto && s : surfels ) viewerCore << s;
165
166 viewerCore.show();
167 }
168
169 // Get selected surfel/point
170 const auto s0 = surfels[ selected_surfels[ 0 ] ];
171 Dimension k0 = K.sOrthDir( s0 );
172 auto voxel0 = K.sIncident( s0, k0, K.sDirect( s0, k0 ) );
173 Point p0 = K.sCoords( voxel0 );
174 auto start0 = point2idx[ p0 ];
175 std::cout << "Start0 index is " << start0 << std::endl;
176 const auto s1 = surfels[ selected_surfels[ 1 ] ];
177 Dimension k1 = K.sOrthDir( s1 );
178 auto voxel1 = K.sIncident( s1, k1, K.sDirect( s1, k1 ) );
179 Point p1 = K.sCoords( voxel1 );
180 auto start1 = point2idx[ p1 ];
181 std::cout << "Start1 index is " << start1 << std::endl;
182
183 // (I) Extracts shortest paths to a target
184
188 TC.init( points.cbegin(), points.cend() );
189 auto SP = TC.makeShortestPaths( opt );
190 SP.init( start0 ); //< set source
191 double last_distance = 0.0;
192 while ( ! SP.finished() )
193 {
194 last_distance = std::get<2>( SP.current() );
195 SP.expand();
196 }
197 std::cout << "Max distance is " << last_distance << std::endl;
199
200 {
201 const int nb_repetitions = 10;
202 const double period = last_distance / nb_repetitions;
203 SimpleDistanceColorMap< double > cmap( 0.0, period, false );
204 MViewer3D viewerCore;
205 Color colSurfel( 200, 200, 255, 128 );
206 Color colStart( 255, 0, 0, 128 );
207
208 for ( auto i = 0; i < points.size(); ++i )
209 {
210 const double d_s = SP.distance( i );
211 Color c_s = cmap( fmod( d_s, period ) );
212 viewerCore.drawColor( c_s );
213 viewerCore.drawBall( RealPoint( points[ i ][ 0 ],
214 points[ i ][ 1 ],
215 points[ i ][ 2 ] ) );
216 }
217
218 // JOL: Left if you wish to display it as a surface, and to display paths.
219
220 auto surfels = SH3::getSurfelRange ( surface );
221 for ( auto && s : surfels )
222 {
223 const double d_s = SP.distance( surfel2idx[ s ] );
224 Color c_s = cmap( fmod( d_s, period ) );
225 viewerCore.drawColor( c_s );
226 viewerCore << s;
227 }
228 viewerCore.drawColor( colStart );
229 viewerCore.drawColor( Color::Green );
230 for ( Index i = 0; i < SP.size(); i++ ) {
231 Point p1 = SP.point( i );
232 Point p2 = SP.point( SP.ancestor( i ) );
233 viewerCore.drawLine( p1, p2 );
234 }
235 viewerCore.show();
236 }
237
238 // (II) Extracts a shortest path between two points.
239
241 auto SP0 = TC.makeShortestPaths( opt );
242 auto SP1 = TC.makeShortestPaths( opt );
243 SP0.init( start0 ); //< source point
244 SP1.init( start1 ); //< target point
245 std::vector< Index > Q; //< the output shortest path
246 while ( ! SP0.finished() && ! SP1.finished() )
247 {
248 auto n0_ = SP0.current(); ((void) n0_);
249 auto n1_ = SP1.current();
250 // auto p0_ = std::get<0>( n0_ );
251 auto p1_ = std::get<0>( n1_ );
252 SP0.expand();
253 SP1.expand();
254 if ( SP0.isVisited( p1_ ) )
255 {
256 auto c0 = SP0.pathToSource( p1_ );
257 auto c1 = SP1.pathToSource( p1_ );
258 std::copy(c0.rbegin(), c0.rend(), std::back_inserter(Q));
259 Q.pop_back();
260 std::copy(c1.begin(), c1.end(), std::back_inserter(Q));
261 break;
262 }
263 }
264 // Q is empty if there is no path.
266
267 // Display computed distances and shortest path
268 {
269 const int nb_repetitions = 10;
270 const double period = last_distance / nb_repetitions;
271 SimpleDistanceColorMap< double > cmap( 0.0, period, false );
272 MViewer3D viewerCore;
273 Color colSurfel( 200, 200, 255, 128 );
274 Color colStart( 255, 0, 0, 128 );
275 for ( auto i = 0; i < points.size(); ++i )
276 {
277 const double d_s0 = SP0.isVisited( i ) ? SP0.distance( i ) : SP0.infinity();
278 const double d_s1 = SP1.isVisited( i ) ? SP1.distance( i ) : SP1.infinity();
279 const double d_s = std::min( d_s0, d_s1 );
280 Color c_s = ( d_s != SP0.infinity() )
281 ? cmap( fmod( d_s, period ) )
282 : Color::Black;
283 viewerCore.drawColor( c_s );
284 viewerCore.drawBall( RealPoint( points[ i ][ 0 ],
285 points[ i ][ 1 ],
286 points[ i ][ 2 ] ) );
287 }
288
289 viewerCore.drawColor( Color::Green );
290 for ( auto i = 1; i < Q.size(); i++ )
291 {
292 Point p1_ = TC.point( Q[ i-1 ] );
293 Point p2_ = TC.point( Q[ i ] );
294 viewerCore.drawLine( p1_, p2_ );
295 }
296 viewerCore.show();
297 }
298
299 // (III) Extracts multiple shortest paths between many sources and two targets.
300
301 std::vector< Index > sources;
302 std::vector< Index > dests;
303 for ( int i = 0; i < 20; i++ )
304 sources.push_back( rand() % TC.size() );
305 dests.push_back( start0 );
306 dests.push_back( start1 );
307 auto paths = TC.shortestPaths( sources, dests, opt );
308
309 // Display them.
310 {
311 MViewer3D viewerCore;
312 Color colSurfel( 200, 200, 255, 128 );
313 Color colStart( 255, 0, 0, 128 );
314 for ( auto i = 0; i < points.size(); ++i )
315 {
316 viewerCore.drawColor( Color( 150, 150, 150, 255 ) );
317 viewerCore.drawBall( RealPoint( points[ i ][ 0 ],
318 points[ i ][ 1 ],
319 points[ i ][ 2 ] ) );
320 }
321 viewerCore.drawColor( Color::Green );
322 for ( auto path : paths )
323 {
324 for ( auto i = 1; i < path.size(); i++ )
325 {
326 Point p1_ = TC.point( path[ i-1 ] );
327 Point p2_ = TC.point( path[ i ] );
328 viewerCore.drawLine( p1_, p2_ );
329 }
330 trace.info() << "length=" << TC.length( path ) << std::endl;
331 }
332 viewerCore.show();
333 }
334
335 return 0;
336}
337// //
339
Structure representing an RGB triple with alpha component.
Definition Color.h:77
static const Color Black
Definition Color.h:422
static const Color Green
Definition Color.h:426
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Dimension sOrthDir(const SCell &s) const
Given a signed surfel [s], returns its orthogonal direction (ie, the coordinate where the surfel is c...
Point sCoords(const SCell &c) const
Return its digital coordinates.
bool sDirect(const SCell &p, Dimension k) const
Return 'true' if the direct orientation of [p] along [k] is in the positive coordinate direction.
SCell sIncident(const SCell &c, Dimension k, bool up) const
Return the forward or backward signed cell incident to [c] along axis [k], depending on [up].
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition Shortcuts.h:102
Aim: simple blue to red colormap for distance information for instance.
PointVector< dim, Integer > Point
Points in DGtal::SpaceND.
Definition SpaceND.h:110
PointVector< dim, double > RealPoint
Definition SpaceND.h:117
PointVector< dim, Integer > Vector
Vectors in DGtal::SpaceND.
Definition SpaceND.h:113
Aim: A class that computes tangency to a given digital set. It provides services to compute all the c...
std::ostream & info()
CountedPtr< SH3::DigitalSurface > surface
Space::Point Point
Z3i::KSpace KSpace
int reaction(void *viewer, DGtal::int32_t name, void *data)
Z3i::SCell SCell
Space::Vector Vector
Shortcuts< KSpace > SH3
Z3i::Space Space
Space::RealPoint RealPoint
Z3i::Domain Domain
SMesh::Index Index
DGtal is the top-level namespace which contains all DGtal functions and types.
std::int32_t int32_t
signed 32-bit integer.
Definition BasicTypes.h:71
DGtal::uint32_t Dimension
Definition Common.h:119
Trace trace
STL namespace.
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