# 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]]

# 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]]]

# 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:
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))

## 2D Images node_array/node_index structure

1 for img_name in imgs_list:
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))

## 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
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