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