1 // This is mul/vil3d/algo/vil3d_find_blobs.cxx
2 #include <vector>
3 #include <iostream>
4 #include <algorithm>
5 #include "vil3d_find_blobs.h"
6 //:
7 // \file
8 // \brief Identify and enumerate all disjoint blobs in a binary image.
9 // \author Ian Scott, Kevin de Souza
10
11 #include <cassert>
12 #ifdef _MSC_VER
13 # include "vcl_msvc_warnings.h"
14 #endif
15
16 #include <vil3d/vil3d_image_view.h>
17 #include <vil3d/vil3d_math.h>
18
19 namespace
20 {
21 //: Manages the relabelling structure.
22 // Got the idea from http://en.wikipedia.org/wiki/Disjoint-set_data_structure Mar 2011.
23 // although no code was copied.
24 class disjoint_sets
25 {
26 typedef unsigned LABEL;
27 typedef unsigned LEN;
28 struct node
29 {
30 LABEL parent;
31 LEN rank;
32 };
33 std::vector<node> store_;
34 public:
disjoint_sets()35 disjoint_sets(): store_(1u)
36 { // renumber 0->0
37 store_.front().parent=0;
38 store_.front().rank=0;
39 }
40
41 //: Get the root label for label v;
root(LABEL v)42 LABEL root(LABEL v)
43 {
44 node & n=store_[v];
45 if (n.parent == v)
46 return v;
47 else
48 {
49 n.parent=root(n.parent); // relabel to speed up later searches
50 return n.parent;
51 }
52 }
53
54 //: Merge two sets containing labels left and right.
merge_labels(LABEL left_label,LABEL right_label)55 void merge_labels(LABEL left_label, LABEL right_label)
56 {
57 LABEL left_root = root(left_label);
58 LABEL right_root = root(right_label);
59 if (left_root == right_root) return; // already merged.
60 node& left_root_node = store_[left_root];
61 node& right_root_node = store_[right_root];
62 if (left_root_node.rank > right_root_node.rank) // Find the larger tree.
63 { // add the right tree to the left
64 right_root_node.parent = left_root_node.parent;
65 }
66 else
67 { // add the left tree to the right
68 left_root_node.parent = right_root_node.parent;
69 if (left_root_node.rank == right_root_node.rank)
70 right_root_node.rank++;
71 }
72 }
73
74 //: Create a new label;
new_label()75 LABEL new_label()
76 {
77 node n = {(LABEL)store_.size(), 0};
78 store_.push_back(n);
79 return n.parent;
80 }
81
size()82 LEN size()
83 {
84 return store_.size();
85 }
86 };
87 }
88
89 // Identify and enumerate all disjoint blobs in a binary image.
vil3d_find_blobs(const vil3d_image_view<bool> & src,vil3d_find_blob_connectivity conn,vil3d_image_view<unsigned> & dst)90 void vil3d_find_blobs(const vil3d_image_view<bool>& src,
91 vil3d_find_blob_connectivity conn,
92 vil3d_image_view<unsigned>& dst)
93 {
94 unsigned ni=src.ni();
95 unsigned nj=src.nj();
96 unsigned nk=src.nk();
97 dst.set_size(ni, nj, nk);
98 dst.fill(0);
99
100 disjoint_sets merge_list;
101 std::vector<unsigned> neighbouring_labels;
102
103 unsigned n_neighbours;
104 switch (conn)
105 {
106 case vil3d_find_blob_connectivity_6_conn:
107 n_neighbours=3;
108 break;
109 case vil3d_find_blob_connectivity_26_conn:
110 n_neighbours=13;
111 break;
112 default:
113 n_neighbours=0; assert(!"unknown connectivity");
114 }
115
116 // The 3-prev-(6)-neighbourhood are the first three entries.
117 // The 13-prev-(26)-neighbourhood are those three plus the rest.
118 int neighbourhood_ii[] = { -1, 0, 0, -1, 0, +1, -1, +1, -1, 0, +1, -1, +1};
119 int neighbourhood_jj[] = { 0, -1, 0, -1, -1, -1, 0, 0, +1, +1, +1, -1, -1};
120 int neighbourhood_kk[] = { 0, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0};
121
122
123 for (unsigned k=0; k<nk; ++k)
124 for (unsigned j=0; j<nj; ++j)
125 for (unsigned i=0; i<ni; ++i)
126 {
127 if (!src(i,j,k)) continue;
128 neighbouring_labels.clear();
129 for (unsigned l=0; l<n_neighbours; ++l)
130 {
131 unsigned ii = i + neighbourhood_ii[l];
132 if (ii >= ni) continue; // rely on wraparound to find -ne overruns.
133 unsigned jj = j + neighbourhood_jj[l];
134 if (jj >= nj) continue;
135 unsigned kk = k + neighbourhood_kk[l];
136 if (kk >= nk) continue;
137 unsigned d = dst(ii,jj,kk);
138 if (d==0) continue;
139 neighbouring_labels.push_back(d);
140 }
141 if (neighbouring_labels.empty())
142 {
143 unsigned new_label = merge_list.new_label();
144 dst(i,j,k) = new_label;
145 }
146 else
147 {
148 // See how many unique labels neighbouring labels we have
149 std::sort(neighbouring_labels.begin(), neighbouring_labels.end());
150 auto end = std::unique(neighbouring_labels.begin(), neighbouring_labels.end());
151 // don't bother erasing unique suffix, just keeping the end iterator
152 // will be enough.
153 auto it=neighbouring_labels.begin();
154 unsigned label = *it++;
155 dst(i,j,k) = label;
156
157 // If we have neighbours with multiple labels.
158 // then record merges of the previously disjointly labelled blocks.
159 // If there was only a single unique label, the following for loop
160 // will not execute.
161 for (; it!=end; ++it)
162 merge_list.merge_labels(*it, label);
163 }
164 }
165
166 unsigned n_merge=merge_list.size();
167 std::vector<unsigned> renumbering(n_merge, 0u);
168 // Convert the merge lsit into a simple renumbering array,
169 // and change to root of each disjoint set to its lowest member.
170 // The reinstates label order to the original raster order.
171 for (unsigned l=1; l<n_merge; ++l)
172 {
173 if (renumbering[l]!=0) continue;
174 unsigned root_label = merge_list.root(l);
175 unsigned root_label_renumber = renumbering[root_label];
176 if (root_label_renumber==0)
177 {
178 renumbering[root_label]=l;
179 renumbering[l]=l;
180 }
181 else
182 renumbering[l]=renumbering[root_label];
183 }
184
185 // Now due to the renumbering, the set of labels may not compactly occupy
186 // the number line. So renumber the renumbering array.
187 std::vector<unsigned> labels(renumbering.begin(), renumbering.end());
188 std::sort(labels.begin(), labels.end());
189 labels.erase(std::unique(labels.begin(), labels.end()), labels.end());
190 const auto dodgy = static_cast<unsigned>(-1);
191 std::vector<unsigned> renum_renum(renumbering.size(), dodgy);
192 renum_renum[0]=0;
193 for (unsigned l=0, n=labels.size(); l<n; ++l)
194 renum_renum[labels[l]] = l;
195
196 for (unsigned int & it : renumbering)
197 it=renum_renum[it];
198
199 // Check than no DODGY values got into the renumbering.
200 assert(std::find(renumbering.begin(), renumbering.end(), dodgy)
201 == renumbering.end() );
202
203 // Renumber the labels, to merge connected blobs, with a compact set of labels.
204
205 for (unsigned k=0; k<nk; ++k)
206 for (unsigned j=0; j<nj; ++j)
207 for (unsigned i=0; i<ni; ++i)
208 dst(i,j,k) = renumbering[dst(i,j,k)];
209 }
210
211
212
213 //: Convert a label image into a list of chorded regions.
214 // A blob label value of n will be returned in dest_regions[n-1].
vil3d_blob_labels_to_regions(const vil3d_image_view<unsigned> & src_label,std::vector<vil3d_region> & regions)215 void vil3d_blob_labels_to_regions(const vil3d_image_view<unsigned>& src_label,
216 std::vector<vil3d_region>& regions)
217 {
218 regions.clear();
219 unsigned ni=src_label.ni();
220 unsigned nj=src_label.nj();
221 unsigned nk=src_label.nk();
222
223 unsigned min_v,max_v;
224 vil3d_math_value_range(src_label,min_v,max_v);
225 if (max_v==0) return; // No blobs.
226
227 regions.resize(max_v);
228
229 for (unsigned k=0; k<nk; ++k)
230 for (unsigned j=0; j<nj; ++j)
231 for (unsigned i=0; i<ni;) // don't auto increment i, since the loop body does it.
232 {
233 unsigned label = src_label(i,j,k);
234 if (!label)
235 { // if not a label - go to next pixel.
236 ++i;
237 continue;
238 }
239 unsigned i_start=i;
240 // Find end of chord.
241 while (++i < ni && src_label(i,j,k)==label);
242 regions[label-1].push_back(vil3d_chord(i_start, i-1, j,k));
243 }
244 }
245