Max-tree implementation C++

MainPage

Implementation based on the union-find with level compression algorithm, which is described in: A Comparative Review of Component Tree Computation Algorithms

OK
OK [C/C++ extension is up-to-date]

Python Implementation

 1 #from max_tree_c_01 import counting_sort_c, canonicalize_c, union_find2d_c, \
 2 #                                 union_find3d_c, computeNodeArray2d_c,\
 3 #                                 computeNodeArray3d_c
 4 
 5 from se2off import se2off
 6 import numpy as np
 7 
 8 
 9 def build_max_tree(f, Bc, option = 0):
10 
11     ndim = f.ndim
12     off = se2off(Bc) # Array of offsets
13 
14     parent = np.empty(f.size, dtype = np.int32)
15     parent[:] = -1
16 
17     zpar = np.empty(f.size, dtype = np.int32)
18 
19     ftype = f.dtype
20     if (ftype == np.uint8):
21         fu16 = f.astype(np.uint16)
22     else:
23         fu16 = f
24 
25     flat_img = fu16.ravel()
26     MAX = flat_img.max()
27     S_rev = counting_sort_c(int(MAX),flat_img)
28 
29     if ndim == 2:
30         H,W = f.shape
31         union_find2d_c(H,W,off,parent,zpar,S_rev,flat_img)
32     elif ndim == 3:
33         L,M,N = f.shape
34         union_find3d_c(L,M,N,off,parent,zpar,S_rev,flat_img)
35     else:
36         print "Invalid option"
37         return
38 
39     canonicalize_c(flat_img,parent,S_rev)
40 
41     if option == 0:
42         return parent,S_rev
43     else:
44         node_index = np.empty(f.shape, dtype = np.int32)
45         node_index[:] = -1
46         if ndim == 2:
47             node_array = computeNodeArray2d_c(parent,flat_img,S_rev,node_index)
48         else:
49             node_array = computeNodeArray3d_c(parent,flat_img,S_rev,node_index)
50 
51         node_array = node_array.T.copy()
52         return parent,S_rev,node_array,node_index

Simple 2D Example

 1 from max_tree_c_01 import build_max_tree
 2 import numpy as np
 3 import time
 4 from generate_graph import generateGraph
 5 
 6 img = np.array([[9,10, 3, 3, 3, 5, 8],
 7                 [2, 0, 4, 0, 0, 0, 7],
 8                 [8, 0, 5, 0, 6, 0, 9],
 9                 [9, 0, 6, 0, 8, 4, 2]],dtype = np.uint8)
10 
11 print "2D test image:\n", img
12 
13 Bc = np.array([[1,1,1],[1,1,1],[1,1,1]], dtype = bool)
14 print "\n2D structuring Element Bc:\n", Bc.astype(int)
15 
16 parent,S_rev,node_array,node_index = build_max_tree(img,Bc, option = 1)
17 print '\nNode_array:\n', node_array
18 print '\nNode_index:\n', node_index
19 graph = generateGraph(parent,img)
20 mmgraphviz(graph,"Pixel oriented max-tree")
2D test image:
[[ 9 10  3  3  3  5  8]
 [ 2  0  4  0  0  0  7]
 [ 8  0  5  0  6  0  9]
 [ 9  0  6  0  8  4  2]]

2D structuring Element Bc:
[[1 1 1]
 [1 1 1]
 [1 1 1]]

Node_array:
[[ 0  0  1  1  3  3  4  4  5  5  8  6  7 11 11  2  9]
 [ 1  2  1  2  2  2  1  1  1  1  0  2  0  0  0  0  0]
 [ 0  2  8  3  4  4  5  6  5  9  6  7  8  9  8  9 10]
 [28 19  2 15  7  5  4  2  2  2  1  3  1  1  1  1  1]
 [ 8  7 14  2 26  9  5 18 16  0 23 13 25 20  6 21  1]
 [42 26  5 17 11  6  3  5  5  0  3  3  3  2  0  3  0]
 [ 0  0  2  0  0  0  0  2  2  0  3  0  3  2  0  3  0]
 [ 3  3  3  3  3  3  2  3  3  0  3  2  3  2  0  3  0]
 [84 58  0 52 36  7 23  8  4  1  2 18  4  6  6  0  1]
 [ 0  0  0  0  4  0  5  4  2  0  2  6  4  6  6  0  1]
 [ 6  6  0  6  6  2  6  4  2  1  2  6  4  6  6  0  1]]

Node_index:
[[ 9 16  3  3  3  6 14]
 [ 1  0  5  0  0  0 11]
 [ 2  0  8  0  7  0 13]
 [15  0 10  0 12  4  1]]
/media/_xsb/iamxt/max_tree_c_01/GRVIZ12946_001.png

Pixel oriented max-tree

Simple 3D Example

 1 img2 = np.zeros((5,5,5), dtype = np.uint8)
 2 img2[2,1:-1,1:-1] = 1
 3 img2[2,2,2] = 2
 4 img2[1,1:-1,1:-1] = 1
 5 img2[3,1:-1,1:-1] = 1
 6 
 7 print "3D test image:\n", img2
 8 
 9 Bc2 = np.ones((3,3,3), dtype = bool)
10 print "3D structuring element:\n", Bc2.astype(int)
11 
12 parent2,S_rev2,node_array2,node_index2 = build_max_tree(img2,Bc2, option = 1)
13 print "\nNode array:\n", node_array2
14 print "\nNode index:\n", node_index2
15 mmgraphviz(generateGraph(parent2,img2.ravel()), 'Max-tree of the 3D image')
3D test image:
[[[0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]]

 [[0 0 0 0 0]
  [0 1 1 1 0]
  [0 1 1 1 0]
  [0 1 1 1 0]
  [0 0 0 0 0]]

 [[0 0 0 0 0]
  [0 1 1 1 0]
  [0 1 2 1 0]
  [0 1 1 1 0]
  [0 0 0 0 0]]

 [[0 0 0 0 0]
  [0 1 1 1 0]
  [0 1 1 1 0]
  [0 1 1 1 0]
  [0 0 0 0 0]]

 [[0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]]]
3D structuring element:
[[[1 1 1]
  [1 1 1]
  [1 1 1]]

 [[1 1 1]
  [1 1 1]
  [1 1 1]]

 [[1 1 1]
  [1 1 1]
  [1 1 1]]]

Node array:
[[  0   0   1]
 [  1   1   0]
 [  0   1   2]
 [125  27   1]
 [  0  31  62]
 [250  54   2]
 [  0   1   2]
 [  4   3   2]
 [250  54   2]
 [  0   1   2]
 [  4   3   2]
 [250  54   2]
 [  0   1   2]
 [  4   3   2]]

Node index:
[[[0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]]

 [[0 0 0 0 0]
  [0 1 1 1 0]
  [0 1 1 1 0]
  [0 1 1 1 0]
  [0 0 0 0 0]]

 [[0 0 0 0 0]
  [0 1 1 1 0]
  [0 1 2 1 0]
  [0 1 1 1 0]
  [0 0 0 0 0]]

 [[0 0 0 0 0]
  [0 1 1 1 0]
  [0 1 1 1 0]
  [0 1 1 1 0]
  [0 0 0 0 0]]

 [[0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]]]
/media/_xsb/iamxt/max_tree_c_01/GRVIZ12946_002.png

Max-tree of the 3D image

Processing Time

2D Images usual structure

 1 import time
 2 
 3 imgs_list = [ 'ia636/woodlog.tif', 'ia636/lenina.pgm', 'ia636/cameraman.pgm', \
 4               'MVBook/bamboo.png', 'itajaiacu/itajaiacu01.jpg']
 5 
 6 for img_name in imgs_list:
 7   img = adreadgray(img_name)
 8   par,S_rev = 0,0
 9   t = time.time()
10   par,S_rev = build_max_tree(img,Bc)
11   t = time.time() - t
12   name = img_name.split('/')[-1]
13   H,W = img.shape
14   adshow(img,"%s, shape = %d x %d, processing time = %fs" %(name,H,W,t))

woodlog.tif, shape = 256 x 256, processing time = 0.009978s

lenina.pgm, shape = 256 x 256, processing time = 0.008347s

cameraman.pgm, shape = 256 x 256, processing time = 0.008955s

bamboo.png, shape = 512 x 512, processing time = 0.034291s

itajaiacu01.jpg, shape = 1200 x 1600, processing time = 0.306124s

2D Images node_array/node_index structure

1 for img_name in imgs_list:
2   img = adreadgray(img_name)
3   par,S_rev,node_array,node_index = 0,0,0,0
4   t = time.time()
5   parent,S_rev,node_array,node_index = build_max_tree(img,Bc, option = 1)
6   t = time.time() - t
7   name = img_name.split('/')[-1]
8   H,W = img.shape
9   adshow(img,"%s, shape = %d x %d, processing time = %fs" %(name,H,W,t))

woodlog.tif, shape = 256 x 256, processing time = 0.014886s

lenina.pgm, shape = 256 x 256, processing time = 0.010476s

cameraman.pgm, shape = 256 x 256, processing time = 0.011146s

bamboo.png, shape = 512 x 512, processing time = 0.044598s

itajaiacu01.jpg, shape = 1200 x 1600, processing time = 0.529856s

3D Image usual structure

 1 from toolbox.fattach import find_attachment_files
 2 import dicom
 3 import os
 4 import ia636 as ia
 5 import homedanilo.library_dpereira84_dicom as lib
 6 import nibabel as nib
 7 
 8 
 9 
10 filename = '/awmedia3/q/orig/LPBA40/native_space/S21/S21.native.mri.img'
11 
12 img1 = nib.load(filename)
13 data = img1.get_data()[:,:,:,0]
14 data = data.astype(np.uint16)
15 print data.shape
16 print "Image size:", data.size
17 Bc = np.ones((3,3,3), dtype = bool)
18 a = np.array([[0,0,0],[0,1,0],[0,0,0]], dtype = bool)
19 b = np.array([[0,1,0],[1,1,1],[0,1,0]], dtype = bool)
20 Bc2 = np.array([a,b,a], dtype = bool)
21 
22 parent,S = 0,0
23 t = time.time()
24 parent,S_rev = build_max_tree(data,Bc2)
25 t = time.time() - t
26 print "Processing time connectivity C6:", t
27 
28 parent,S = 0,0
29 t = time.time()
30 parent,S_rev = build_max_tree(data,Bc)
31 t = time.time() - t
32 print "Processing time connectivity C26:", t
(256, 124, 256)
Image size: 8126464
Processing time connectivity C6: 1.79462909698
Processing time connectivity C26: 3.56044101715

3D Image node_array/node_index structure

 1 print "Image size:", data.size
 2 parent,S,node_array,node_index = 0,0,0,0
 3 t = time.time()
 4 parent,S_rev,node_array,node_index = build_max_tree(data,Bc2, option =1)
 5 t = time.time() - t
 6 print "Processing time connectivity C6:", t
 7 
 8 parent,S,node_array,node_index = 0,0,0,0
 9 t = time.time()
10 parent,S_rev,node_array,node_index = build_max_tree(data,Bc, option = 1)
11 t = time.time() - t
12 print "Processing time connectivity C26:", t
Image size: 8126464
Processing time connectivity C6: 2.48807787895
Processing time connectivity C26: 3.94429898262