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