1.1 diff -r 5edfbd3e7a46 -r 1204ebf9340d PTdecode/CImg-1.3.0/plugins/skeleton.h 1.2 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.3 +++ b/PTdecode/CImg-1.3.0/plugins/skeleton.h Mon Aug 03 14:09:20 2009 +0100 1.4 @@ -0,0 +1,728 @@ 1.5 +/* 1.6 + # 1.7 + # File : skeleton.h 1.8 + # ( C++ header file - CImg plug-in ) 1.9 + # 1.10 + # Description : CImg plugin that implements the computation of the Hamilton-Jacobi skeletons 1.11 + # using Siddiqi algorithm with the correction proposed by Torsello, 1.12 + # as described in : 1.13 + # 1.14 + # [SBTZ02] K. Siddiqi, S. Bouix, A. Tannenbaum and S.W. Zucker. Hamilton-Jacobi Skeletons 1.15 + # International Journal of Computer Vision, 48(3):215-231, 2002 1.16 + # 1.17 + # [TH03] A. Torsello and E. R. Hancock. Curvature Correction of the Hamilton-Jacobi Skeleton 1.18 + # IEEE Computer Vision and Pattern Recognition, 2003 1.19 + # 1.20 + # [BST05] S. Bouix, K. Siddiqi and A. Tannenbaum. Flux driven automatic centerline 1.21 + # extraction. Medical Image Analysis, 9:209-221, 2005 1.22 + # 1.23 + # IMPORTANT WARNING : You must include STL's <queue> before plugin inclusion to make it working ! 1.24 + # 1.25 + # Copyright : Francois-Xavier Dupe 1.26 + # ( http://www.greyc.ensicaen.fr/~fdupe/ ) 1.27 + # 1.28 + # This software is governed by the CeCILL license under French law and 1.29 + # abiding by the rules of distribution of free software. You can use, 1.30 + # modify and/or redistribute the software under the terms of the CeCILL 1.31 + # license as circulated by CEA, CNRS and INRIA at the following URL 1.32 + # "http://www.cecill.info". 1.33 + # 1.34 + # As a counterpart to the access to the source code and rights to copy, 1.35 + # modify and redistribute granted by the license, users are provided only 1.36 + # with a limited warranty and the software's author, the holder of the 1.37 + # economic rights, and the successive licensors have only limited 1.38 + # liability. 1.39 + # 1.40 + # In this respect, the user's attention is drawn to the risks associated 1.41 + # with loading, using, modifying and/or developing or reproducing the 1.42 + # software by the user in light of its specific status of free software, 1.43 + # that may mean that it is complicated to manipulate, and that also 1.44 + # therefore means that it is reserved for developers and experienced 1.45 + # professionals having in-depth computer knowledge. Users are therefore 1.46 + # encouraged to load and test the software's suitability as regards their 1.47 + # requirements in conditions enabling the security of their systems and/or 1.48 + # data to be ensured and, more generally, to use and operate it in the 1.49 + # same conditions as regards security. 1.50 + # 1.51 + # The fact that you are presently reading this means that you have had 1.52 + # knowledge of the CeCILL license and that you accept its terms. 1.53 + # 1.54 +*/ 1.55 +#ifndef cimg_plugin_skeleton 1.56 +#define cimg_plugin_skeleton 1.57 + 1.58 +/** 1.59 + * Compute the flux of the gradient 1.60 + * @param grad the gradient of the distance function 1.61 + * @param sY the sampling size in Y 1.62 + * @param sZ the sampling size in Z 1.63 + * @return the flux 1.64 + */ 1.65 +CImg<floatT> get_flux ( const CImgList<floatT> & grad, 1.66 + float sY = 1.0f, float sZ = 1.0f ) const 1.67 +{ 1.68 + int stop = 0; // Stop flag 1.69 + float f = 0; // The current flux 1.70 + int count = 0; // Counter 1.71 + CImg<floatT> flux ((*this).dimx(), (*this).dimy(), (*this).dimz()); 1.72 + flux.fill(0); 1.73 + 1.74 + cimg_forXYZ((*this),x,y,z) 1.75 + { 1.76 + // If the point is the backgroung continue 1.77 + if ( (*this)(x,y,z) == 0 ) 1.78 + continue; 1.79 + // Look at the neigthboorhound and compute the flux 1.80 + stop = 0; 1.81 + f = 0; 1.82 + count = 0; 1.83 + 1.84 + for ( int k = -1; k <= 1; ++k ) 1.85 + for ( int l = -1; l <= 1; ++l ) 1.86 + for ( int m = -1; m <= 1; ++m ) 1.87 + { 1.88 + if ( stop == 1 ) 1.89 + continue; 1.90 + 1.91 + // Protection 1.92 + if (( x+k < 0 ) || ( x+k >= (*this).dimx() ) || ( y+l < 0 ) || ( y+l >= (*this).dimy() ) || 1.93 + ( z+m < 0 ) || ( z+m >= (*this).dimz() ) || ( k==0 && l==0 && m==0 )) 1.94 + continue; 1.95 + 1.96 + ++count; 1.97 + 1.98 + //Test if the point is in the interior 1.99 + if ( (*this)(x+k,y+l,z+m) == 0 ) 1.100 + { 1.101 + stop = 1; 1.102 + continue; 1.103 + } 1.104 + 1.105 + // Compute the flux 1.106 + 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)); 1.107 + } 1.108 + 1.109 + //Update 1.110 + if ( stop == 1 || count == 0 ) 1.111 + flux(x,y,z) = 0; 1.112 + else 1.113 + flux(x,y,z) = f / count; 1.114 + 1.115 + } 1.116 + 1.117 + return flux; 1.118 +} 1.119 + 1.120 +/** 1.121 + * Definition of a point with his flux value 1.122 + */ 1.123 +struct _PointFlux 1.124 +{ 1.125 + int pos [3]; 1.126 + float flux; 1.127 + float dist; 1.128 +}; 1.129 + 1.130 +/** 1.131 + * Class for the priority queue 1.132 + */ 1.133 +class _compare_point 1.134 +{ 1.135 + /** 1.136 + * Create medial curves 1.137 + */ 1.138 + bool curve; 1.139 + 1.140 + public: 1.141 + _compare_point ( bool curve = false ) 1.142 + { 1.143 + this->curve = curve; 1.144 + } 1.145 + 1.146 + bool operator() ( const _PointFlux & p1, const _PointFlux & p2 ) const 1.147 + { 1.148 + if ( curve ) 1.149 + { 1.150 + if ( p1.dist > p2.dist ) 1.151 + return true; 1.152 + else if ( p1.dist == p2.dist && p1.flux < p2.flux ) 1.153 + return true; 1.154 + } 1.155 + else 1.156 + { 1.157 + if ( p1.flux < p2.flux ) 1.158 + return true; 1.159 + else if ( p1.flux == p2.flux && p1.dist > p2.dist ) 1.160 + return true; 1.161 + } 1.162 + 1.163 + return false; 1.164 + } 1.165 +}; 1.166 + 1.167 +/** 1.168 + * Compute the log-density using the algorithm from Torsello 1.169 + * @param dist the distance map 1.170 + * @param grad the gradient of the distance map, e.g. the flux 1.171 + * @param flux the divergence map 1.172 + * @param delta the threshold for the division 1.173 + * @return the logdensity \rho 1.174 + */ 1.175 +CImg<floatT> get_logdensity ( const CImg<floatT> & dist, 1.176 + const CImgList<floatT> & grad, 1.177 + const CImg<floatT> & flux, float delta = 0.1 ) const 1.178 +{ 1.179 + std::priority_queue< _PointFlux, std::vector<_PointFlux>, _compare_point > pqueue(true); 1.180 + CImg<floatT> logdensity ((*this).dimx(), (*this).dimy(), (*this).dimz(), 1, 0); 1.181 + 1.182 + // 1 - Put all the pixel inside the priority queue 1.183 + cimg_forXYZ(dist,x,y,z) 1.184 + if ( dist(x,y,z) != 0 ) 1.185 + { 1.186 + _PointFlux p; 1.187 + p.pos[0] = x; 1.188 + p.pos[1] = y; 1.189 + p.pos[2] = z; 1.190 + p.flux = 0; 1.191 + p.dist = dist(x,y,z); 1.192 + pqueue.push(p); 1.193 + } 1.194 + 1.195 + // 2 - Compute the logdensity 1.196 + while ( ! pqueue.empty() ) 1.197 + { 1.198 + _PointFlux p = pqueue.top(); 1.199 + pqueue.pop(); 1.200 + 1.201 + float Fx = grad(0,p.pos[0],p.pos[1],p.pos[2]); 1.202 + float Fy = grad(1,p.pos[0],p.pos[1],p.pos[2]); 1.203 + float Fz = grad(2,p.pos[0],p.pos[1],p.pos[2]); 1.204 + 1.205 + 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) 1.206 + - 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)); 1.207 + 1.208 + float tmp = 1.0f - (1.0f-fabs(Fx)) * (1.0f-fabs(Fy)) * (1.0f-fabs(Fz)); 1.209 + if ( tmp > delta ) 1.210 + logdensity(p.pos[0],p.pos[1],p.pos[2]) /= tmp; 1.211 + else if ( delta < 1 ) 1.212 + logdensity(p.pos[0],p.pos[1],p.pos[2]) = 0; 1.213 + } 1.214 + 1.215 + return logdensity; 1.216 +} 1.217 + 1.218 +/** 1.219 + * Computed the corrected divergence map using Torsello formula and idea 1.220 + * @param logdensity the log density map 1.221 + * @param grad the gradient of the distance map 1.222 + * @param flux the flux using siddiqi formula 1.223 + * @param delta the discrete step 1.224 + * @return the corrected divergence map 1.225 + */ 1.226 +CImg<floatT> get_corrected_flux ( const CImg<floatT> & logdensity, 1.227 + const CImgList<floatT> & grad, 1.228 + const CImg<floatT> & flux, 1.229 + float delta = 1.0 ) const 1.230 +{ 1.231 + CImg<floatT> corr_map ((*this).dimx(), (*this).dimy(), (*this).dimz(), 1, 0); 1.232 + 1.233 + cimg_forXYZ(corr_map,x,y,z) 1.234 + { 1.235 + float Fx = grad(0,x,y,z); 1.236 + float Fy = grad(1,x,y,z); 1.237 + float Fz = grad(2,x,y,z); 1.238 + 1.239 + 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) + 1.240 + 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))); 1.241 + } 1.242 + 1.243 + return corr_map; 1.244 +} 1.245 + 1.246 +/** 1.247 + * Test if a point is simple using Euler number for 2D case 1.248 + * or using Malandain criterion for 3D case 1.249 + * @param img the image 1.250 + * @param x the x coordinate 1.251 + * @param y the y coordinate 1.252 + * @param z the z coordinate 1.253 + * @return true if simple 1.254 + */ 1.255 +bool _isSimple ( const CImg<T> & img, int x, int y, int z ) const 1.256 +{ 1.257 + if ( img.dimz() == 1 ) // 2D case 1.258 + { 1.259 + int V = 0; // Number of vertices 1.260 + int E = 0; // Number of edges 1.261 + 1.262 + for ( int k = -1; k <= 1; ++k ) 1.263 + for ( int l = -1; l <= 1; ++l ) 1.264 + { 1.265 + // Protection 1.266 + if ( x+k < 0 || x+k >= img.dimx() || y+l < 0 || y+l >= img.dimy() ) 1.267 + continue; 1.268 + 1.269 + // Count the number of vertices 1.270 + if ( img(x+k,y+l) != 0 && ! ( k == 0 && l == 0 )) 1.271 + { 1.272 + ++V; 1.273 + 1.274 + // Count the number of edges 1.275 + for ( int k1 = -1; k1 <= 1; ++k1 ) 1.276 + for ( int l1 = -1; l1 <= 1; ++l1 ) 1.277 + { 1.278 + // Protection 1.279 + if ( x+k+k1 < 0 || x+k+k1 >= img.dimx() || y+l+l1 < 0 || y+l+l1 >= img.dimy() ) 1.280 + continue; 1.281 + 1.282 + if ( !(k1 == 0 && l1 == 0) && img(x+k+k1,y+l+l1) != 0 && k+k1 > -2 && l+l1 > -2 1.283 + && k+k1 < 2 && l+l1 < 2 && !( k+k1 == 0 && l+l1 == 0 )) 1.284 + ++E; 1.285 + } 1.286 + } 1.287 + } 1.288 + 1.289 + // Remove the corner if exists 1.290 + if ( x-1 >= 0 && y-1 >= 0 && img(x-1,y-1) != 0 && img(x,y-1) != 0 && img(x-1,y) != 0 ) 1.291 + E -= 2; 1.292 + 1.293 + if ( x-1 >= 0 && y+1 < img.dimy() && img(x-1,y+1) != 0 && img(x,y+1) != 0 && img(x-1,y) != 0 ) 1.294 + E -= 2; 1.295 + 1.296 + if ( x+1 < img.dimx() && y-1 >= 0 && img(x+1,y-1) != 0 && img(x,y-1) != 0 && img(x+1,y) != 0 ) 1.297 + E -= 2; 1.298 + 1.299 + 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 ) 1.300 + E -= 2; 1.301 + 1.302 + // Final return true if it is a tree (eg euler number equal to 1) 1.303 + if (( V - E/2 ) == 1 ) 1.304 + return true; 1.305 + } 1.306 + else // 3D case 1.307 + { 1.308 + int C_asterix = 0; 1.309 + int C_bar = 0; 1.310 + CImg<intT> visit ( 3, 3, 3, 1, 0 ); // Visitor table 1.311 + int count = 0; 1.312 + 1.313 + visit(1,1,1) = -1; 1.314 + 1.315 + // Compute C^* 1.316 + 1.317 + // Seeking for a component 1.318 + for ( int k = -1; k <= 1; ++k ) 1.319 + for ( int l = -1; l <= 1; ++l ) 1.320 + for ( int m = -1; m <= 1; ++m ) 1.321 + { 1.322 + int label = 0; 1.323 + 1.324 + // Protection 1.325 + if ( x+k < 0 || x+k >= img.dimx() || 1.326 + y+l < 0 || y+l >= img.dimy() || 1.327 + z+m < 0 || z+m >= img.dimz() || 1.328 + ( k==0 && l==0 && m==0 )) 1.329 + continue; 1.330 + 1.331 + if ( visit(1+k,1+l,1+m) == 0 && img(x+k,y+l,z+m) != 0 ) 1.332 + { 1.333 + // Look after the neightbor 1.334 + for ( int k1 = -1; k1 <= 1; ++k1 ) 1.335 + for ( int l1 = -1; l1 <= 1; ++l1 ) 1.336 + for ( int m1 = -1; m1 <= 1; ++m1 ) 1.337 + { 1.338 + // Protection 1.339 + if ( x+k+k1 < 0 || x+k+k1 >= img.dimx() || 1.340 + y+l+l1 < 0 || y+l+l1 >= img.dimy() || 1.341 + z+m+m1 < 0 || z+m+m1 >= img.dimz() || 1.342 + k+k1 > 1 || k+k1 < -1 || 1.343 + l+l1 > 1 || l+l1 < -1 || 1.344 + m+m1 > 1 || m+m1 < -1 ) 1.345 + continue; 1.346 + 1.347 + // Search for a already knew component 1.348 + if ( visit(1+k+k1,1+l+l1,1+m+m1) > 0 && 1.349 + img(x+k+k1,y+l+l1,z+m+m1) != 0 ) 1.350 + { 1.351 + if ( label == 0 ) 1.352 + label = visit(1+k+k1,1+l+l1,1+m+m1); 1.353 + else if ( label != visit(1+k+k1,1+l+l1,1+m+m1) ) 1.354 + { 1.355 + // Meld component 1.356 + --C_asterix; 1.357 + 1.358 + int C = visit(1+k+k1,1+l+l1,1+m+m1); 1.359 + cimg_forXYZ(visit,a,b,c) 1.360 + if ( visit(a,b,c) == C ) 1.361 + visit(a,b,c) = label; 1.362 + } 1.363 + } 1.364 + } 1.365 + 1.366 + // Label the point 1.367 + if ( label == 0 ) 1.368 + { 1.369 + // Find a new component 1.370 + ++C_asterix; 1.371 + ++count; 1.372 + visit(1+k,1+l,1+m) = count; 1.373 + } 1.374 + else 1.375 + { 1.376 + visit(1+k,1+l,1+m) = label; 1.377 + } 1.378 + } 1.379 + } 1.380 + 1.381 + if ( C_asterix != 1 ) 1.382 + return false; 1.383 + 1.384 + // Compute \bar{C} 1.385 + 1.386 + // Reinit visit 1.387 + visit.fill(0); 1.388 + visit(1,1,1) = -1; 1.389 + 1.390 + // Seeking for a component 1.391 + 1.392 + // Look at X-axis 1.393 + { for ( int k = -1; k <= 1; ++k ) 1.394 + { 1.395 + if ( x+k < 0 || x+k >= img.dimx() ) 1.396 + continue; 1.397 + 1.398 + if ( img(x+k,y,z) == 0 && visit(1+k,1,1) == 0 ) 1.399 + { 1.400 + ++C_bar; 1.401 + ++count; 1.402 + visit(1+k,1,1) = count; 1.403 + 1.404 + // Follow component 1.405 + for ( int l = -1; l <= 1; ++l ) 1.406 + { 1.407 + if ( y+l < img.dimy() && y+l >= 0 && img(x+k,y+l,z) == 0 && visit(1+k,1+l,1) == 0 ) 1.408 + visit(1+k,1+l,1) = count; 1.409 + if ( z+l < img.dimz() && z+l >= 0 && img(x+k,y,z+l) == 0 && visit(1+k,1,1+l) == 0 ) 1.410 + visit(1+k,1,1+l) = count; 1.411 + } 1.412 + } 1.413 + } 1.414 + } 1.415 + 1.416 + // Look at Y-axis 1.417 + { for ( int k = -1; k <= 1; ++k ) 1.418 + { 1.419 + if ( y+k < 0 || y+k >= img.dimy() ) 1.420 + continue; 1.421 + 1.422 + if ( img(x,y+k,z) == 0 && visit(1,1+k,1) == 0 ) 1.423 + { 1.424 + int label = 0; 1.425 + ++C_bar; 1.426 + ++count; 1.427 + visit(1,1+k,1) = count; 1.428 + label = count; 1.429 + 1.430 + // Follow component 1.431 + for ( int l = -1; l <= 1; ++l ) 1.432 + { 1.433 + if ( l == 0 ) 1.434 + continue; 1.435 + 1.436 + if ( x+l < img.dimx() && x+l >= 0 && img(x+l,y+k,z) == 0 ) 1.437 + { 1.438 + if ( visit(1+l,1+k,1) != 0 ) 1.439 + { 1.440 + if ( label != visit(1+l,1+k,1) ) 1.441 + { 1.442 + // Meld component 1.443 + --C_bar; 1.444 + 1.445 + int C = visit(1+l,1+k,1); 1.446 + cimg_forXYZ(visit,a,b,c) 1.447 + if ( visit(a,b,c) == C ) 1.448 + visit(a,b,c) = label; 1.449 + } 1.450 + } 1.451 + else 1.452 + visit(1+l,1+k,1) = label; 1.453 + } 1.454 + 1.455 + if ( z+l < img.dimz() && z+l >= 0 && img(x,y+k,z+l) == 0 ) 1.456 + { 1.457 + if ( visit(1,1+k,1+l) != 0 ) 1.458 + { 1.459 + if ( label != visit(1,1+k,1+l) ) 1.460 + { 1.461 + // Meld component 1.462 + --C_bar; 1.463 + 1.464 + int C = visit(1,1+k,1+l); 1.465 + cimg_forXYZ(visit,a,b,c) 1.466 + if ( visit(a,b,c) == C ) 1.467 + visit(a,b,c) = label; 1.468 + } 1.469 + } 1.470 + else 1.471 + visit(1,1+k,1+l) = label; 1.472 + } 1.473 + } 1.474 + } 1.475 + } 1.476 + } 1.477 + 1.478 + // Look at Z-axis 1.479 + { for ( int k = -1; k <= 1; ++k ) 1.480 + { 1.481 + if ( z+k < 0 || z+k >= img.dimz() ) 1.482 + continue; 1.483 + 1.484 + if ( img(x,y,z+k) == 0 && visit(1,1,1+k) == 0 ) 1.485 + { 1.486 + int label = 0; 1.487 + ++C_bar; 1.488 + ++count; 1.489 + visit(1,1,1+k) = count; 1.490 + label = count; 1.491 + 1.492 + // Follow component 1.493 + for ( int l = -1; l <= 1; ++l ) 1.494 + { 1.495 + if ( l == 0 ) 1.496 + continue; 1.497 + 1.498 + if ( x+l < img.dimx() && x+l >= 0 && img(x+l,y,z+k) == 0 ) 1.499 + { 1.500 + if ( visit(1+l,1,1+k) != 0 ) 1.501 + { 1.502 + if ( label != visit(1+l,1,1+k) ) 1.503 + { 1.504 + // Meld component 1.505 + --C_bar; 1.506 + 1.507 + int C = visit(1+l,1,1+k); 1.508 + cimg_forXYZ(visit,a,b,c) 1.509 + if ( visit(a,b,c) == C ) 1.510 + visit(a,b,c) = label; 1.511 + } 1.512 + } 1.513 + else 1.514 + visit(1+l,1,1+k) = label; 1.515 + } 1.516 + 1.517 + if ( y+l < img.dimy() && y+l >= 0 && img(x,y+l,z+k) == 0 ) 1.518 + { 1.519 + if ( visit(1,1+l,1+k) != 0 ) 1.520 + { 1.521 + if ( label != visit(1,1+l,1+k) ) 1.522 + { 1.523 + // Meld component 1.524 + --C_bar; 1.525 + 1.526 + int C = visit(1,1+l,1+k); 1.527 + cimg_forXYZ(visit,a,b,c) 1.528 + if ( visit(a,b,c) == C ) 1.529 + visit(a,b,c) = label; 1.530 + } 1.531 + } 1.532 + else 1.533 + visit(1,1+l,1+k) = label; 1.534 + } 1.535 + } 1.536 + } 1.537 + } 1.538 + } 1.539 + 1.540 + if ( C_bar == 1 ) 1.541 + return true; 1.542 + } 1.543 + 1.544 + return false; 1.545 +} 1.546 + 1.547 +/** 1.548 + * Test if a point is a end point 1.549 + * @param img the image 1.550 + * @param label the table of labels 1.551 + * @param curve set it to true for having medial curve 1.552 + * @param x the x coordinate 1.553 + * @param y the y coordinate 1.554 + * @param z the z coordinate 1.555 + * @return true if simple 1.556 + */ 1.557 +bool _isEndPoint ( const CImg<T> & img, const CImg<T> & label, 1.558 + bool curve, int x, int y, int z ) const 1.559 +{ 1.560 + if ( label(x,y,z) == 1 ) 1.561 + return true; 1.562 + 1.563 + if (( ! curve ) && ( img.dimz() != 1 )) // 3D case with medial surface 1.564 + { 1.565 + // Use Pudney specification with the 9 plans 1.566 + 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 1.567 + { {-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 1.568 + { {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 1.569 + { {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 1.570 + { {-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 1.571 + { {-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 1.572 + { {-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 1.573 + { {-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 1.574 + { {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 1.575 + }; 1.576 + 1.577 + // Count the number of neighbors on each plan 1.578 + for ( int k = 0; k < 9; ++k ) 1.579 + { 1.580 + int count = 0; 1.581 + 1.582 + for ( int l = 0; l < 8; ++l ) 1.583 + { 1.584 + if ( x+plan9[k][l][0] < 0 || x+plan9[k][l][0] >= img.dimx() || 1.585 + y+plan9[k][l][1] < 0 || y+plan9[k][l][1] >= img.dimy() || 1.586 + z+plan9[k][l][2] < 0 || z+plan9[k][l][2] >= img.dimz() ) 1.587 + continue; 1.588 + 1.589 + if ( img(x+plan9[k][l][0],y+plan9[k][l][1],z+plan9[k][l][2]) != 0 ) 1.590 + ++count; 1.591 + } 1.592 + 1.593 + if ( count < 2 ) 1.594 + return true; 1.595 + } 1.596 + } 1.597 + else // 2D or 3D case with medial curve 1.598 + { 1.599 + int isb = 0; 1.600 + 1.601 + for ( int k = -1; k <= 1; ++k ) 1.602 + for ( int l = -1; l <= 1; ++l ) 1.603 + for ( int m = -1; m <= 1; ++m ) 1.604 + { 1.605 + // Protection 1.606 + if ( x+k < 0 || x+k >= img.dimx() || y+l < 0 || y+l >= img.dimy() || 1.607 + z+m < 0 || z+m >= img.dimz() ) 1.608 + continue; 1.609 + 1.610 + if ( img(x+k,y+l,z+m) != 0 ) 1.611 + ++isb; 1.612 + } 1.613 + 1.614 + if ( isb == 2 ) // The pixel with one neighbor 1.615 + return true; 1.616 + } 1.617 + 1.618 + // Else it's not... 1.619 + return false; 1.620 +} 1.621 + 1.622 +/** 1.623 + * Compute the skeleton of the shape using Hamilton-Jacobi scheme 1.624 + * @param flux the flux of the distance gradient 1.625 + * @param dist the euclidean distance of the object 1.626 + * @param curve create or not medial curve 1.627 + * @param thres the threshold on the flux 1.628 + * @return the skeleton 1.629 + */ 1.630 +CImg<T> get_skeleton ( const CImg<floatT> & flux, 1.631 + const CImg<floatT> & dist, bool curve ,float thres ) const 1.632 +{ 1.633 + CImg<T> skeleton ( *this ); // The skeleton 1.634 + CImg<T> label ( (*this).dimx(), (*this).dimy(), (*this).dimz(), 1, 0 ); // Save label 1.635 + CImg<T> count ( (*this).dimx(), (*this).dimy(), (*this).dimz(), 1, 0 ); // A counter for the queue 1.636 + std::priority_queue< _PointFlux, std::vector<_PointFlux>, _compare_point > pqueue(curve); 1.637 + int isb = 0; 1.638 + 1.639 + // 1 - Init get the bound points 1.640 + cimg_forXYZ(*this,x,y,z) 1.641 + { 1.642 + if ( skeleton(x,y,z) == 0 ) 1.643 + continue; 1.644 + 1.645 + // Test bound condition 1.646 + isb = 0; 1.647 + for ( int k = -1; k <= 1; ++k ) 1.648 + for ( int l = -1; l <= 1; ++l ) 1.649 + for ( int m = -1; m <= 1; ++m ) 1.650 + { 1.651 + // Protection 1.652 + if ( x+k < 0 || x+k >= (*this).dimx() || y+l < 0 || y+l >= (*this).dimy() || 1.653 + z+m < 0 || z+m >= (*this).dimz() ) 1.654 + continue; 1.655 + 1.656 + if ( skeleton(x+k,y+l,z+m) == 0 ) 1.657 + isb = 1; 1.658 + } 1.659 + 1.660 + if ( isb == 1 && _isSimple(skeleton,x,y,z) ) 1.661 + { 1.662 + _PointFlux p; 1.663 + p.pos[0] = x; 1.664 + p.pos[1] = y; 1.665 + p.pos[2] = z; 1.666 + p.flux = flux(x,y,z); 1.667 + p.dist = dist(x,y,z); 1.668 + pqueue.push(p); 1.669 + count(x,y,z) = 1; 1.670 + } 1.671 + 1.672 + } 1.673 + 1.674 + // 2 - Compute the skeleton 1.675 + while ( ! pqueue.empty() ) 1.676 + { 1.677 + _PointFlux p = pqueue.top(); // Get the point with the max flux 1.678 + pqueue.pop(); // Remove the point from the queue 1.679 + count(p.pos[0],p.pos[1],p.pos[2]) = 0; // Reinit counter 1.680 + 1.681 + // Test if the point is simple 1.682 + if ( _isSimple(skeleton,p.pos[0],p.pos[1],p.pos[2]) ) 1.683 + { 1.684 + if (( ! _isEndPoint(skeleton,label,curve,p.pos[0],p.pos[1],p.pos[2]) ) || p.flux > thres ) 1.685 + { 1.686 + skeleton(p.pos[0],p.pos[1],p.pos[2]) = 0; // Remove the point 1.687 + 1.688 + for ( int k = -1; k <= 1; ++k ) 1.689 + for ( int l = -1; l <= 1; ++l ) 1.690 + for ( int m = -1; m <= 1; ++m ) 1.691 + { 1.692 + // Protection 1.693 + if ( p.pos[0]+k < 0 || p.pos[0]+k >= (*this).dimx() || 1.694 + p.pos[1]+l < 0 || p.pos[1]+l >= (*this).dimy() || 1.695 + p.pos[2]+m < 0 || p.pos[2]+m >= (*this).dimz() ) 1.696 + continue; 1.697 + 1.698 + if ( skeleton(p.pos[0]+k,p.pos[1]+l,p.pos[2]+m) != 0 && 1.699 + count(p.pos[0]+k,p.pos[1]+l,p.pos[2]+m) < 1 && 1.700 + _isSimple(skeleton,p.pos[0]+k,p.pos[1]+l,p.pos[2]+m) ) 1.701 + { 1.702 + _PointFlux p1; 1.703 + p1.pos[0] = p.pos[0]+k; 1.704 + p1.pos[1] = p.pos[1]+l; 1.705 + p1.pos[2] = p.pos[2]+m; 1.706 + p1.flux = flux(p.pos[0]+k,p.pos[1]+l,p.pos[2]+m); 1.707 + p1.dist = dist(p.pos[0]+k,p.pos[1]+l,p.pos[2]+m); 1.708 + pqueue.push(p1); 1.709 + count(p.pos[0]+k,p.pos[1]+l,p.pos[2]+m) = 1; 1.710 + } 1.711 + } 1.712 + } 1.713 + else 1.714 + { 1.715 + label(p.pos[0],p.pos[1],p.pos[2]) = 1; // Mark the point as skeletal 1.716 + } 1.717 + } 1.718 + } 1.719 + 1.720 + return skeleton; 1.721 +} 1.722 + 1.723 +/** 1.724 + * In place version of get_skeleton 1.725 + */ 1.726 +CImg<T> skeleton ( const CImg<floatT> & flux, 1.727 + const CImg<floatT> & dist, bool curve ,float thres ) 1.728 +{ 1.729 + return get_skeleton(flux,dist,curve,thres).transfer_to(*this); 1.730 +} 1.731 + 1.732 +#endif /* cimg_skeleton_plugin */