In Viewer1D, use more structured data as example.
[stefan/lulu.git] / lulu / ccomp.pyx
1 # -*- python -*-
2 #cython: cdivision=True
3
4 import numpy as np
5 cimport numpy as np
6
7 """
8 See also:
9
10   Christophe Fiorio and Jens Gustedt,
11   "Two linear time Union-Find strategies for image processing",
12   Theoretical Computer Science 154 (1996), pp. 165-181.
13
14   Kensheng Wu, Ekow Otoo and Arie Shoshani,
15   "Optimizing connected component labeling algorithms",
16   Paper LBNL-56864, 2005,
17   Lawrence Berkeley National Laboratory
18   (University of California),
19   http://repositories.cdlib.org/lbnl/LBNL-56864.
20
21 """
22
23 # Tree operations implemented by an array as described in Wu et al.
24
25 DTYPE = np.int
26 ctypedef np.int_t DTYPE_t
27
28 cdef DTYPE_t find_root(np.int_t *work, np.int_t n):
29     """Find the root of node n.
30
31     """
32     cdef np.int_t root = n
33     while (work[root] < root):
34         root = work[root]
35     return root
36
37 cdef set_root(np.int_t *work, np.int_t n, np.int_t root):
38     """
39     Set all nodes on a path to point to new_root.
40
41     """
42     cdef np.int_t j
43     while (work[n] < n):
44         j = work[n]
45         work[n] = root
46         n = j
47
48     work[n] = root
49
50
51 cdef join_trees(np.int_t *work, np.int_t n, np.int_t m):
52     """Join two trees containing nodes n and m.
53
54     """
55     cdef np.int_t root = find_root(work, n)
56     cdef np.int_t root_m
57
58     if (n != m):
59         root_m = find_root(work, m)
60
61         if (root > root_m):
62             root = root_m
63
64         set_root(work, n, root)
65         set_root(work, m, root)
66
67 # Connected components search as described in Fiorio et al.
68
69 def label(np.ndarray[DTYPE_t, ndim=2] input):
70     """Label connected regions of an integer array.
71
72     """
73     cdef np.int_t rows = input.shape[0]
74     cdef np.int_t cols = input.shape[1]
75
76     cdef np.ndarray[DTYPE_t, ndim=2] data = input.copy()
77     cdef np.ndarray[DTYPE_t, ndim=2] work
78
79     work = np.arange(data.size, dtype=DTYPE).reshape((rows, cols))
80
81     cdef np.int_t *work_p = <np.int_t*>work.data
82     cdef np.int_t *data_p = <np.int_t*>data.data
83
84     cdef np.int_t i, j
85
86     # Initialize the first row
87     for j in range(1, cols):
88         if data[0, j] == data[0, j-1]:
89             join_trees(work_p, j, j-1)
90
91     for i in range(1, rows):
92         # Handle the first column
93         if data[i, 0] == data[i-1, 0]:
94             join_trees(work_p, i*cols, (i-1)*cols)
95
96         if data[i, 0] == data[i-1, 1]:
97             join_trees(work_p, i*cols, (i-1)*cols + 1)
98
99         for j in range(1, cols):
100             if data[i, j] == data[i-1, j-1]:
101                 join_trees(work_p, i*cols + j, (i-1)*cols + j - 1)
102
103             if data[i, j] == data[i-1, j]:
104                 join_trees(work_p, i*cols + j, (i-1)*cols + j)
105
106             if j < cols - 1:
107                 if data[i, j] == data[i - 1, j + 1]:
108                     join_trees(work_p, i*cols + j, (i-1)*cols + j + 1)
109
110             if data[i, j] == data[i, j-1]:
111                 join_trees(work_p, i*cols + j, i*cols + j - 1)
112
113     # Label output
114
115     cdef np.int_t ctr = 0
116     for i in range(rows):
117         for j in range(cols):
118             if (i*cols + j) == work[i, j]:
119                 data[i, j] = ctr
120                 ctr = ctr + 1
121             else:
122                 data[i, j] = data_p[work[i, j]]
123
124     return data