Wed, 05 Aug 2009 17:10:56 +0100
add README
1 /*
2 #
3 # File : skeleton.h
4 # ( C++ header file - CImg plug-in )
5 #
6 # Description : CImg plugin that implements the computation of the Hamilton-Jacobi skeletons
7 # using Siddiqi algorithm with the correction proposed by Torsello,
8 # as described in :
9 #
10 # [SBTZ02] K. Siddiqi, S. Bouix, A. Tannenbaum and S.W. Zucker. Hamilton-Jacobi Skeletons
11 # International Journal of Computer Vision, 48(3):215-231, 2002
12 #
13 # [TH03] A. Torsello and E. R. Hancock. Curvature Correction of the Hamilton-Jacobi Skeleton
14 # IEEE Computer Vision and Pattern Recognition, 2003
15 #
16 # [BST05] S. Bouix, K. Siddiqi and A. Tannenbaum. Flux driven automatic centerline
17 # extraction. Medical Image Analysis, 9:209-221, 2005
18 #
19 # IMPORTANT WARNING : You must include STL's <queue> before plugin inclusion to make it working !
20 #
21 # Copyright : Francois-Xavier Dupe
22 # ( http://www.greyc.ensicaen.fr/~fdupe/ )
23 #
24 # This software is governed by the CeCILL license under French law and
25 # abiding by the rules of distribution of free software. You can use,
26 # modify and/or redistribute the software under the terms of the CeCILL
27 # license as circulated by CEA, CNRS and INRIA at the following URL
28 # "http://www.cecill.info".
29 #
30 # As a counterpart to the access to the source code and rights to copy,
31 # modify and redistribute granted by the license, users are provided only
32 # with a limited warranty and the software's author, the holder of the
33 # economic rights, and the successive licensors have only limited
34 # liability.
35 #
36 # In this respect, the user's attention is drawn to the risks associated
37 # with loading, using, modifying and/or developing or reproducing the
38 # software by the user in light of its specific status of free software,
39 # that may mean that it is complicated to manipulate, and that also
40 # therefore means that it is reserved for developers and experienced
41 # professionals having in-depth computer knowledge. Users are therefore
42 # encouraged to load and test the software's suitability as regards their
43 # requirements in conditions enabling the security of their systems and/or
44 # data to be ensured and, more generally, to use and operate it in the
45 # same conditions as regards security.
46 #
47 # The fact that you are presently reading this means that you have had
48 # knowledge of the CeCILL license and that you accept its terms.
49 #
50 */
51 #ifndef cimg_plugin_skeleton
52 #define cimg_plugin_skeleton
54 /**
55 * Compute the flux of the gradient
56 * @param grad the gradient of the distance function
57 * @param sY the sampling size in Y
58 * @param sZ the sampling size in Z
59 * @return the flux
60 */
61 CImg<floatT> get_flux ( const CImgList<floatT> & grad,
62 float sY = 1.0f, float sZ = 1.0f ) const
63 {
64 int stop = 0; // Stop flag
65 float f = 0; // The current flux
66 int count = 0; // Counter
67 CImg<floatT> flux ((*this).dimx(), (*this).dimy(), (*this).dimz());
68 flux.fill(0);
70 cimg_forXYZ((*this),x,y,z)
71 {
72 // If the point is the backgroung continue
73 if ( (*this)(x,y,z) == 0 )
74 continue;
75 // Look at the neigthboorhound and compute the flux
76 stop = 0;
77 f = 0;
78 count = 0;
80 for ( int k = -1; k <= 1; ++k )
81 for ( int l = -1; l <= 1; ++l )
82 for ( int m = -1; m <= 1; ++m )
83 {
84 if ( stop == 1 )
85 continue;
87 // Protection
88 if (( x+k < 0 ) || ( x+k >= (*this).dimx() ) || ( y+l < 0 ) || ( y+l >= (*this).dimy() ) ||
89 ( z+m < 0 ) || ( z+m >= (*this).dimz() ) || ( k==0 && l==0 && m==0 ))
90 continue;
92 ++count;
94 //Test if the point is in the interior
95 if ( (*this)(x+k,y+l,z+m) == 0 )
96 {
97 stop = 1;
98 continue;
99 }
101 // Compute the flux
102 f += ( grad(0,x+k,y+l,z+m)*k + grad(1,x+k,y+l,z+m)*l/sY + grad(2,x+k,y+l,z+m)*m/sZ ) / cimg_std::sqrt((float)(k*k+l*l+m*m));
103 }
105 //Update
106 if ( stop == 1 || count == 0 )
107 flux(x,y,z) = 0;
108 else
109 flux(x,y,z) = f / count;
111 }
113 return flux;
114 }
116 /**
117 * Definition of a point with his flux value
118 */
119 struct _PointFlux
120 {
121 int pos [3];
122 float flux;
123 float dist;
124 };
126 /**
127 * Class for the priority queue
128 */
129 class _compare_point
130 {
131 /**
132 * Create medial curves
133 */
134 bool curve;
136 public:
137 _compare_point ( bool curve = false )
138 {
139 this->curve = curve;
140 }
142 bool operator() ( const _PointFlux & p1, const _PointFlux & p2 ) const
143 {
144 if ( curve )
145 {
146 if ( p1.dist > p2.dist )
147 return true;
148 else if ( p1.dist == p2.dist && p1.flux < p2.flux )
149 return true;
150 }
151 else
152 {
153 if ( p1.flux < p2.flux )
154 return true;
155 else if ( p1.flux == p2.flux && p1.dist > p2.dist )
156 return true;
157 }
159 return false;
160 }
161 };
163 /**
164 * Compute the log-density using the algorithm from Torsello
165 * @param dist the distance map
166 * @param grad the gradient of the distance map, e.g. the flux
167 * @param flux the divergence map
168 * @param delta the threshold for the division
169 * @return the logdensity \rho
170 */
171 CImg<floatT> get_logdensity ( const CImg<floatT> & dist,
172 const CImgList<floatT> & grad,
173 const CImg<floatT> & flux, float delta = 0.1 ) const
174 {
175 std::priority_queue< _PointFlux, std::vector<_PointFlux>, _compare_point > pqueue(true);
176 CImg<floatT> logdensity ((*this).dimx(), (*this).dimy(), (*this).dimz(), 1, 0);
178 // 1 - Put all the pixel inside the priority queue
179 cimg_forXYZ(dist,x,y,z)
180 if ( dist(x,y,z) != 0 )
181 {
182 _PointFlux p;
183 p.pos[0] = x;
184 p.pos[1] = y;
185 p.pos[2] = z;
186 p.flux = 0;
187 p.dist = dist(x,y,z);
188 pqueue.push(p);
189 }
191 // 2 - Compute the logdensity
192 while ( ! pqueue.empty() )
193 {
194 _PointFlux p = pqueue.top();
195 pqueue.pop();
197 float Fx = grad(0,p.pos[0],p.pos[1],p.pos[2]);
198 float Fy = grad(1,p.pos[0],p.pos[1],p.pos[2]);
199 float Fz = grad(2,p.pos[0],p.pos[1],p.pos[2]);
201 logdensity(p.pos[0],p.pos[1],p.pos[2]) = logdensity.linear_atXYZ(p.pos[0]-Fx,p.pos[1]-Fy,p.pos[2]-Fz)
202 - 0.5f * (flux(p.pos[0],p.pos[1],p.pos[2])+flux.linear_atXYZ(p.pos[0]-Fx,p.pos[1]-Fy,p.pos[2]-Fz));
204 float tmp = 1.0f - (1.0f-fabs(Fx)) * (1.0f-fabs(Fy)) * (1.0f-fabs(Fz));
205 if ( tmp > delta )
206 logdensity(p.pos[0],p.pos[1],p.pos[2]) /= tmp;
207 else if ( delta < 1 )
208 logdensity(p.pos[0],p.pos[1],p.pos[2]) = 0;
209 }
211 return logdensity;
212 }
214 /**
215 * Computed the corrected divergence map using Torsello formula and idea
216 * @param logdensity the log density map
217 * @param grad the gradient of the distance map
218 * @param flux the flux using siddiqi formula
219 * @param delta the discrete step
220 * @return the corrected divergence map
221 */
222 CImg<floatT> get_corrected_flux ( const CImg<floatT> & logdensity,
223 const CImgList<floatT> & grad,
224 const CImg<floatT> & flux,
225 float delta = 1.0 ) const
226 {
227 CImg<floatT> corr_map ((*this).dimx(), (*this).dimy(), (*this).dimz(), 1, 0);
229 cimg_forXYZ(corr_map,x,y,z)
230 {
231 float Fx = grad(0,x,y,z);
232 float Fy = grad(1,x,y,z);
233 float Fz = grad(2,x,y,z);
235 corr_map(x,y,z) = (logdensity(x,y,z) - logdensity.linear_atXYZ(x-Fx,y-Fy,z-Fz)) * expf(logdensity(x,y,z) - 0.5f * delta) +
236 0.5f * ( flux.linear_atXYZ(x-Fx,y-Fy,z-Fz)*expf(logdensity.linear_atXYZ(x-Fx,y-Fy,z-Fz)) + flux(x,y,z)*expf(logdensity(x,y,z)));
237 }
239 return corr_map;
240 }
242 /**
243 * Test if a point is simple using Euler number for 2D case
244 * or using Malandain criterion for 3D case
245 * @param img the image
246 * @param x the x coordinate
247 * @param y the y coordinate
248 * @param z the z coordinate
249 * @return true if simple
250 */
251 bool _isSimple ( const CImg<T> & img, int x, int y, int z ) const
252 {
253 if ( img.dimz() == 1 ) // 2D case
254 {
255 int V = 0; // Number of vertices
256 int E = 0; // Number of edges
258 for ( int k = -1; k <= 1; ++k )
259 for ( int l = -1; l <= 1; ++l )
260 {
261 // Protection
262 if ( x+k < 0 || x+k >= img.dimx() || y+l < 0 || y+l >= img.dimy() )
263 continue;
265 // Count the number of vertices
266 if ( img(x+k,y+l) != 0 && ! ( k == 0 && l == 0 ))
267 {
268 ++V;
270 // Count the number of edges
271 for ( int k1 = -1; k1 <= 1; ++k1 )
272 for ( int l1 = -1; l1 <= 1; ++l1 )
273 {
274 // Protection
275 if ( x+k+k1 < 0 || x+k+k1 >= img.dimx() || y+l+l1 < 0 || y+l+l1 >= img.dimy() )
276 continue;
278 if ( !(k1 == 0 && l1 == 0) && img(x+k+k1,y+l+l1) != 0 && k+k1 > -2 && l+l1 > -2
279 && k+k1 < 2 && l+l1 < 2 && !( k+k1 == 0 && l+l1 == 0 ))
280 ++E;
281 }
282 }
283 }
285 // Remove the corner if exists
286 if ( x-1 >= 0 && y-1 >= 0 && img(x-1,y-1) != 0 && img(x,y-1) != 0 && img(x-1,y) != 0 )
287 E -= 2;
289 if ( x-1 >= 0 && y+1 < img.dimy() && img(x-1,y+1) != 0 && img(x,y+1) != 0 && img(x-1,y) != 0 )
290 E -= 2;
292 if ( x+1 < img.dimx() && y-1 >= 0 && img(x+1,y-1) != 0 && img(x,y-1) != 0 && img(x+1,y) != 0 )
293 E -= 2;
295 if ( x+1 < img.dimx() && y+1 < img.dimy() && img(x+1,y+1) != 0 && img(x,y+1) != 0 && img(x+1,y) != 0 )
296 E -= 2;
298 // Final return true if it is a tree (eg euler number equal to 1)
299 if (( V - E/2 ) == 1 )
300 return true;
301 }
302 else // 3D case
303 {
304 int C_asterix = 0;
305 int C_bar = 0;
306 CImg<intT> visit ( 3, 3, 3, 1, 0 ); // Visitor table
307 int count = 0;
309 visit(1,1,1) = -1;
311 // Compute C^*
313 // Seeking for a component
314 for ( int k = -1; k <= 1; ++k )
315 for ( int l = -1; l <= 1; ++l )
316 for ( int m = -1; m <= 1; ++m )
317 {
318 int label = 0;
320 // Protection
321 if ( x+k < 0 || x+k >= img.dimx() ||
322 y+l < 0 || y+l >= img.dimy() ||
323 z+m < 0 || z+m >= img.dimz() ||
324 ( k==0 && l==0 && m==0 ))
325 continue;
327 if ( visit(1+k,1+l,1+m) == 0 && img(x+k,y+l,z+m) != 0 )
328 {
329 // Look after the neightbor
330 for ( int k1 = -1; k1 <= 1; ++k1 )
331 for ( int l1 = -1; l1 <= 1; ++l1 )
332 for ( int m1 = -1; m1 <= 1; ++m1 )
333 {
334 // Protection
335 if ( x+k+k1 < 0 || x+k+k1 >= img.dimx() ||
336 y+l+l1 < 0 || y+l+l1 >= img.dimy() ||
337 z+m+m1 < 0 || z+m+m1 >= img.dimz() ||
338 k+k1 > 1 || k+k1 < -1 ||
339 l+l1 > 1 || l+l1 < -1 ||
340 m+m1 > 1 || m+m1 < -1 )
341 continue;
343 // Search for a already knew component
344 if ( visit(1+k+k1,1+l+l1,1+m+m1) > 0 &&
345 img(x+k+k1,y+l+l1,z+m+m1) != 0 )
346 {
347 if ( label == 0 )
348 label = visit(1+k+k1,1+l+l1,1+m+m1);
349 else if ( label != visit(1+k+k1,1+l+l1,1+m+m1) )
350 {
351 // Meld component
352 --C_asterix;
354 int C = visit(1+k+k1,1+l+l1,1+m+m1);
355 cimg_forXYZ(visit,a,b,c)
356 if ( visit(a,b,c) == C )
357 visit(a,b,c) = label;
358 }
359 }
360 }
362 // Label the point
363 if ( label == 0 )
364 {
365 // Find a new component
366 ++C_asterix;
367 ++count;
368 visit(1+k,1+l,1+m) = count;
369 }
370 else
371 {
372 visit(1+k,1+l,1+m) = label;
373 }
374 }
375 }
377 if ( C_asterix != 1 )
378 return false;
380 // Compute \bar{C}
382 // Reinit visit
383 visit.fill(0);
384 visit(1,1,1) = -1;
386 // Seeking for a component
388 // Look at X-axis
389 { for ( int k = -1; k <= 1; ++k )
390 {
391 if ( x+k < 0 || x+k >= img.dimx() )
392 continue;
394 if ( img(x+k,y,z) == 0 && visit(1+k,1,1) == 0 )
395 {
396 ++C_bar;
397 ++count;
398 visit(1+k,1,1) = count;
400 // Follow component
401 for ( int l = -1; l <= 1; ++l )
402 {
403 if ( y+l < img.dimy() && y+l >= 0 && img(x+k,y+l,z) == 0 && visit(1+k,1+l,1) == 0 )
404 visit(1+k,1+l,1) = count;
405 if ( z+l < img.dimz() && z+l >= 0 && img(x+k,y,z+l) == 0 && visit(1+k,1,1+l) == 0 )
406 visit(1+k,1,1+l) = count;
407 }
408 }
409 }
410 }
412 // Look at Y-axis
413 { for ( int k = -1; k <= 1; ++k )
414 {
415 if ( y+k < 0 || y+k >= img.dimy() )
416 continue;
418 if ( img(x,y+k,z) == 0 && visit(1,1+k,1) == 0 )
419 {
420 int label = 0;
421 ++C_bar;
422 ++count;
423 visit(1,1+k,1) = count;
424 label = count;
426 // Follow component
427 for ( int l = -1; l <= 1; ++l )
428 {
429 if ( l == 0 )
430 continue;
432 if ( x+l < img.dimx() && x+l >= 0 && img(x+l,y+k,z) == 0 )
433 {
434 if ( visit(1+l,1+k,1) != 0 )
435 {
436 if ( label != visit(1+l,1+k,1) )
437 {
438 // Meld component
439 --C_bar;
441 int C = visit(1+l,1+k,1);
442 cimg_forXYZ(visit,a,b,c)
443 if ( visit(a,b,c) == C )
444 visit(a,b,c) = label;
445 }
446 }
447 else
448 visit(1+l,1+k,1) = label;
449 }
451 if ( z+l < img.dimz() && z+l >= 0 && img(x,y+k,z+l) == 0 )
452 {
453 if ( visit(1,1+k,1+l) != 0 )
454 {
455 if ( label != visit(1,1+k,1+l) )
456 {
457 // Meld component
458 --C_bar;
460 int C = visit(1,1+k,1+l);
461 cimg_forXYZ(visit,a,b,c)
462 if ( visit(a,b,c) == C )
463 visit(a,b,c) = label;
464 }
465 }
466 else
467 visit(1,1+k,1+l) = label;
468 }
469 }
470 }
471 }
472 }
474 // Look at Z-axis
475 { for ( int k = -1; k <= 1; ++k )
476 {
477 if ( z+k < 0 || z+k >= img.dimz() )
478 continue;
480 if ( img(x,y,z+k) == 0 && visit(1,1,1+k) == 0 )
481 {
482 int label = 0;
483 ++C_bar;
484 ++count;
485 visit(1,1,1+k) = count;
486 label = count;
488 // Follow component
489 for ( int l = -1; l <= 1; ++l )
490 {
491 if ( l == 0 )
492 continue;
494 if ( x+l < img.dimx() && x+l >= 0 && img(x+l,y,z+k) == 0 )
495 {
496 if ( visit(1+l,1,1+k) != 0 )
497 {
498 if ( label != visit(1+l,1,1+k) )
499 {
500 // Meld component
501 --C_bar;
503 int C = visit(1+l,1,1+k);
504 cimg_forXYZ(visit,a,b,c)
505 if ( visit(a,b,c) == C )
506 visit(a,b,c) = label;
507 }
508 }
509 else
510 visit(1+l,1,1+k) = label;
511 }
513 if ( y+l < img.dimy() && y+l >= 0 && img(x,y+l,z+k) == 0 )
514 {
515 if ( visit(1,1+l,1+k) != 0 )
516 {
517 if ( label != visit(1,1+l,1+k) )
518 {
519 // Meld component
520 --C_bar;
522 int C = visit(1,1+l,1+k);
523 cimg_forXYZ(visit,a,b,c)
524 if ( visit(a,b,c) == C )
525 visit(a,b,c) = label;
526 }
527 }
528 else
529 visit(1,1+l,1+k) = label;
530 }
531 }
532 }
533 }
534 }
536 if ( C_bar == 1 )
537 return true;
538 }
540 return false;
541 }
543 /**
544 * Test if a point is a end point
545 * @param img the image
546 * @param label the table of labels
547 * @param curve set it to true for having medial curve
548 * @param x the x coordinate
549 * @param y the y coordinate
550 * @param z the z coordinate
551 * @return true if simple
552 */
553 bool _isEndPoint ( const CImg<T> & img, const CImg<T> & label,
554 bool curve, int x, int y, int z ) const
555 {
556 if ( label(x,y,z) == 1 )
557 return true;
559 if (( ! curve ) && ( img.dimz() != 1 )) // 3D case with medial surface
560 {
561 // Use Pudney specification with the 9 plans
562 int plan9 [9][8][3] = { { {-1,0,-1}, {0,0,-1}, {1,0,-1}, {-1,0,0}, {1,0,0}, {-1,0,1}, {0,0,1}, {1,0,1} }, // Plan 1
563 { {-1,1,0}, {0,1,0}, {1,1,0}, {-1,0,0}, {1,0,0}, {-1,-1,0}, {0,-1,0}, {1,-1,0} }, // Plan 2
564 { {0,-1,-1}, {0,0,-1}, {0,1,-1}, {0,-1,0}, {0,1,0}, {0,-1,1}, {0,0,1}, {0,1,1} }, // Plan 3
565 { {1,1,1}, {0,1,0}, {-1,1,-1}, {1,0,1}, {-1,0,-1}, {-1,-1,-1}, {0,-1,0}, {1,-1,1} }, // Plan 4
566 { {-1,1,1}, {0,1,0}, {1,1,-1}, {-1,0,1}, {1,0,-1}, {-1,-1,1}, {0,-1,0}, {1,-1,-1} }, // Plan 5
567 { {-1,1,1}, {0,1,1}, {1,1,1}, {-1,0,0}, {1,0,0}, {-1,-1,-1}, {0,-1,-1}, {1,-1,-1} }, // Plan 6
568 { {-1,1,-1}, {0,1,-1}, {1,1,-1}, {-1,0,0}, {1,0,0}, {-1,-1,1}, {0,-1,1}, {1,-1,1} }, // Plan 7
569 { {-1,1,-1}, {-1,1,0}, {-1,1,1}, {0,0,-1}, {0,0,1}, {1,-1,-1}, {1,-1,0}, {1,-1,1} }, // Plan 8
570 { {1,1,-1}, {1,1,0}, {1,1,1}, {0,0,-1}, {0,0,1}, {-1,-1,-1}, {-1,-1,0}, {-1,-1,1} } // Plan 9
571 };
573 // Count the number of neighbors on each plan
574 for ( int k = 0; k < 9; ++k )
575 {
576 int count = 0;
578 for ( int l = 0; l < 8; ++l )
579 {
580 if ( x+plan9[k][l][0] < 0 || x+plan9[k][l][0] >= img.dimx() ||
581 y+plan9[k][l][1] < 0 || y+plan9[k][l][1] >= img.dimy() ||
582 z+plan9[k][l][2] < 0 || z+plan9[k][l][2] >= img.dimz() )
583 continue;
585 if ( img(x+plan9[k][l][0],y+plan9[k][l][1],z+plan9[k][l][2]) != 0 )
586 ++count;
587 }
589 if ( count < 2 )
590 return true;
591 }
592 }
593 else // 2D or 3D case with medial curve
594 {
595 int isb = 0;
597 for ( int k = -1; k <= 1; ++k )
598 for ( int l = -1; l <= 1; ++l )
599 for ( int m = -1; m <= 1; ++m )
600 {
601 // Protection
602 if ( x+k < 0 || x+k >= img.dimx() || y+l < 0 || y+l >= img.dimy() ||
603 z+m < 0 || z+m >= img.dimz() )
604 continue;
606 if ( img(x+k,y+l,z+m) != 0 )
607 ++isb;
608 }
610 if ( isb == 2 ) // The pixel with one neighbor
611 return true;
612 }
614 // Else it's not...
615 return false;
616 }
618 /**
619 * Compute the skeleton of the shape using Hamilton-Jacobi scheme
620 * @param flux the flux of the distance gradient
621 * @param dist the euclidean distance of the object
622 * @param curve create or not medial curve
623 * @param thres the threshold on the flux
624 * @return the skeleton
625 */
626 CImg<T> get_skeleton ( const CImg<floatT> & flux,
627 const CImg<floatT> & dist, bool curve ,float thres ) const
628 {
629 CImg<T> skeleton ( *this ); // The skeleton
630 CImg<T> label ( (*this).dimx(), (*this).dimy(), (*this).dimz(), 1, 0 ); // Save label
631 CImg<T> count ( (*this).dimx(), (*this).dimy(), (*this).dimz(), 1, 0 ); // A counter for the queue
632 std::priority_queue< _PointFlux, std::vector<_PointFlux>, _compare_point > pqueue(curve);
633 int isb = 0;
635 // 1 - Init get the bound points
636 cimg_forXYZ(*this,x,y,z)
637 {
638 if ( skeleton(x,y,z) == 0 )
639 continue;
641 // Test bound condition
642 isb = 0;
643 for ( int k = -1; k <= 1; ++k )
644 for ( int l = -1; l <= 1; ++l )
645 for ( int m = -1; m <= 1; ++m )
646 {
647 // Protection
648 if ( x+k < 0 || x+k >= (*this).dimx() || y+l < 0 || y+l >= (*this).dimy() ||
649 z+m < 0 || z+m >= (*this).dimz() )
650 continue;
652 if ( skeleton(x+k,y+l,z+m) == 0 )
653 isb = 1;
654 }
656 if ( isb == 1 && _isSimple(skeleton,x,y,z) )
657 {
658 _PointFlux p;
659 p.pos[0] = x;
660 p.pos[1] = y;
661 p.pos[2] = z;
662 p.flux = flux(x,y,z);
663 p.dist = dist(x,y,z);
664 pqueue.push(p);
665 count(x,y,z) = 1;
666 }
668 }
670 // 2 - Compute the skeleton
671 while ( ! pqueue.empty() )
672 {
673 _PointFlux p = pqueue.top(); // Get the point with the max flux
674 pqueue.pop(); // Remove the point from the queue
675 count(p.pos[0],p.pos[1],p.pos[2]) = 0; // Reinit counter
677 // Test if the point is simple
678 if ( _isSimple(skeleton,p.pos[0],p.pos[1],p.pos[2]) )
679 {
680 if (( ! _isEndPoint(skeleton,label,curve,p.pos[0],p.pos[1],p.pos[2]) ) || p.flux > thres )
681 {
682 skeleton(p.pos[0],p.pos[1],p.pos[2]) = 0; // Remove the point
684 for ( int k = -1; k <= 1; ++k )
685 for ( int l = -1; l <= 1; ++l )
686 for ( int m = -1; m <= 1; ++m )
687 {
688 // Protection
689 if ( p.pos[0]+k < 0 || p.pos[0]+k >= (*this).dimx() ||
690 p.pos[1]+l < 0 || p.pos[1]+l >= (*this).dimy() ||
691 p.pos[2]+m < 0 || p.pos[2]+m >= (*this).dimz() )
692 continue;
694 if ( skeleton(p.pos[0]+k,p.pos[1]+l,p.pos[2]+m) != 0 &&
695 count(p.pos[0]+k,p.pos[1]+l,p.pos[2]+m) < 1 &&
696 _isSimple(skeleton,p.pos[0]+k,p.pos[1]+l,p.pos[2]+m) )
697 {
698 _PointFlux p1;
699 p1.pos[0] = p.pos[0]+k;
700 p1.pos[1] = p.pos[1]+l;
701 p1.pos[2] = p.pos[2]+m;
702 p1.flux = flux(p.pos[0]+k,p.pos[1]+l,p.pos[2]+m);
703 p1.dist = dist(p.pos[0]+k,p.pos[1]+l,p.pos[2]+m);
704 pqueue.push(p1);
705 count(p.pos[0]+k,p.pos[1]+l,p.pos[2]+m) = 1;
706 }
707 }
708 }
709 else
710 {
711 label(p.pos[0],p.pos[1],p.pos[2]) = 1; // Mark the point as skeletal
712 }
713 }
714 }
716 return skeleton;
717 }
719 /**
720 * In place version of get_skeleton
721 */
722 CImg<T> skeleton ( const CImg<floatT> & flux,
723 const CImg<floatT> & dist, bool curve ,float thres )
724 {
725 return get_skeleton(flux,dist,curve,thres).transfer_to(*this);
726 }
728 #endif /* cimg_skeleton_plugin */