-
Notifications
You must be signed in to change notification settings - Fork 260
/
conncomp.py
381 lines (320 loc) · 16.6 KB
/
conncomp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
"""Class / utilities for connected components."""
############################################################
# Program is part of MintPy #
# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi #
# Author: Zhang Yunjun, 2018 #
############################################################
# Recommend import:
# from mintpy.objects.conncomp import connectComponent
import itertools
import time
import numpy as np
from scipy.sparse import csgraph as csg
from scipy.spatial import cKDTree
from skimage import measure, morphology as morph, segmentation as seg
from mintpy.objects.ramp import deramp
######################################## utilities functions ##############################
def label_conn_comp(mask, min_area=2.5e3, erosion_size=5, print_msg=False):
"""Label / clean up the conn comp (mask)
Parameters: mask - 2D np.ndarray of bool/int
min_area - float, minimum region/area size
erosion_size - int (odd number), size of erosion structure
set to 0 to turn it off.
Returns: label_img - 2d np.ndarray of int, labeled array where all
connected regions are assigned the same value
num_label - int, number of labeled regions
"""
# label
label_img, num_label = measure.label(mask, connectivity=1, return_num=True)
## remove small regions
min_area = min(min_area, label_img.size * 3e-3)
if print_msg:
print(f'remove regions with area < {int(min_area)}')
mask = morph.remove_small_objects(label_img, min_size=min_area, connectivity=1)
label_img[mask == 0] = 0
# update label
label_img, num_label = measure.label(label_img, connectivity=1, return_num=True) # re-label
## remove regions that would disappear after erosion
# to ensure the consistency between label_img and label_bound
if erosion_size > 0:
erosion_structure = np.ones((erosion_size, erosion_size))
label_erosion_img = morph.erosion(label_img, erosion_structure).astype(np.uint8)
erosion_regions = measure.regionprops(label_erosion_img)
if len(erosion_regions) < num_label:
if print_msg:
print('regions lost during morphological erosion operation:')
label_erosion = [reg.label for reg in erosion_regions]
for orig_reg in measure.regionprops(label_img):
if orig_reg.label not in label_erosion:
label_img[label_img == orig_reg.label] = 0
if print_msg:
print('label: {}, area: {}, bbox: {}'.format(orig_reg.label,
orig_reg.area,
orig_reg.bbox))
# update label
label_img, num_label = measure.label(label_img, connectivity=1, return_num=True)
return label_img, num_label
def label_boundary(label_img, num_label, erosion_size=5, print_msg=False):
"""Label the boundary of the labeled array
Parameters: label_img - 2d np.ndarray of int, labeled array where all connected regions are assigned the same value
num_label - int, number of labeled regions
Returns: label_img - 2d np.ndarray of int, labeled array where all connected regions are assigned the same value
num_label - int, number of labeled regions
label_bound - 2d np.ndarrary of bool, where True represent a boundary pixel.
"""
if erosion_size > 0:
# remove regions that would disappear after erosion
# to ensure the consistency between label_img and label_bound
erosion_structure = np.ones((erosion_size, erosion_size))
label_erosion_img = morph.erosion(label_img, erosion_structure).astype(np.uint8)
erosion_regions = measure.regionprops(label_erosion_img)
if len(erosion_regions) < num_label:
if print_msg:
print('regions lost during morphological erosion operation:')
label_erosion = [reg.label for reg in erosion_regions]
for orig_reg in measure.regionprops(label_img):
if orig_reg.label not in label_erosion:
label_img[label_img == orig_reg.label] = 0
if print_msg:
print('label: {}, area: {}, bbox: {}'.format(orig_reg.label,
orig_reg.area,
orig_reg.bbox))
# update label
label_img, num_label = measure.label(label_img, connectivity=1, return_num=True) # re-label
# get label boundaries to facilitate bridge finding
label_bound = seg.find_boundaries(label_erosion_img, mode='thick').astype(np.uint8)
label_bound *= label_erosion_img
return label_img, num_label, label_bound
######################################## beginning of connectComponent class ##############################
class connectComponent:
""" Object for bridging connected components.
Example:
unw_file = 'filt_fine.unw'
# prepare connectComponent object
atr = readfile.read_attribute(unw_file)
conncomp = readfile.read(unw_file+'.conncomp')[0]
cc = connectComponent(conncomp=conncomp, metadata=atr)
cc.label()
cc.find_mst_bridge()
# run bridging
unw = readfile.read(unw_file)[0]
bdg_unw = cc.unwrap_conn_comp(unw, ramp_type='linear')
# write output file
writefile.write(bdg_unw, 'bdg_'+unw_file, atr)
"""
def __init__(self, conncomp, metadata):
"""Parameters: conncomp : 2D np.ndarray in np.bool_ format
metadata : dict, attributes
"""
if type(conncomp).__module__ != np.__name__:
raise ValueError(f'Input conncomp is not np.ndarray: {type(conncomp).__module__}')
self.conncomp = conncomp
self.metadata = metadata
if 'REF_Y' in metadata.keys():
self.refY = int(self.metadata['REF_Y'])
self.refX = int(self.metadata['REF_X'])
else:
self.refY = None
self.refX = None
self.length, self.width = self.conncomp.shape
def label(self, min_area=2.5e3, erosion_size=5, print_msg=False):
""" Label the connected components
Returns: self.labelImg - 2D np.ndarray in int64 to mask areas to be corrected
self.labelBound - 2D np.ndarray in uint8 for label boundaries to find bridges
"""
# label conn comp
(self.labelImg,
self.numLabel) = label_conn_comp(self.conncomp,
min_area=min_area,
print_msg=print_msg)
# label conn comp boundaries
(self.labelImg,
self.numLabel,
self.labelBound) = label_boundary(self.labelImg,
self.numLabel,
erosion_size=erosion_size,
print_msg=print_msg)
# reference label (ref_y/x or the largest one)
if self.refY is not None:
self.labelRef = self.labelImg[self.refY, self.refX]
if self.labelRef == 0:
raise ValueError('input reference point is NOT included in the connectComponent.')
else:
regions = measure.regionprops(self.labelImg)
idx = np.argmax([region.area for region in regions])
self.labelRef = regions[idx].label
return
def get_all_bridge(self):
""" Search all possible connections among labeled regions
Returns: connDict : dict of connection, i.e.:
{'1_2': {'1': array([1232, 345]),
'2': array([868, 239]),
'distance': 379.1200337623956},
'1_3': {'1': array([1232, 345]),
'3': array([1089, 191]),
'distance': 210.1547049199708},
'1_4': {'1': array([1204, 1143]),
'4': array([1217, 1157]),
'distance': 19.1049731745428},
'1_5': {'1': array([1263, 557]),
'5': array([1270, 565]),
'distance': 10.63014581273465},
'2_3': {'2': array([868, 239]),
'3': array([891, 249]),
'distance': 25.079872407968907},
'2_4': {'2': array([868, 239]),
'4': array([1273, 1103]),
'distance': 954.2122405419037},
'2_5': {'2': array([868, 239]),
'5': array([1269, 566]),
'distance': 517.4263232577175},
'3_4': {'3': array([996, 275]),
'4': array([1319, 1085]),
'distance': 872.0258023705492},
'3_5': {'3': array([1015, 264]),
'5': array([1289, 545]),
'distance': 392.4754769409167},
'4_5': {'4': array([1319, 1085]),
'5': array([1305, 670]),
'distance': 415.2360774306587}
}
distMat : 2D np.array in size of (nLabel, nLabel), i.e.:
array([[ 0. , 379.12003 , 210.15471 , 19.104973, 10.630146],
[379.12003 , 0. , 25.079872, 954.2122 , 517.42633 ],
[210.15471 , 25.079872, 0. , 872.0258 , 392.47546 ],
[ 19.104973, 954.2122 , 872.0258 , 0. , 415.23608 ],
[ 10.630146, 517.42633 , 392.47546 , 415.23608 , 0. ]],
dtype=float32)
"""
regions = measure.regionprops(self.labelBound)
trees = []
for i in range(self.numLabel):
trees.append(cKDTree(regions[i].coords))
self.connDict = dict()
self.distMat = np.zeros((self.numLabel, self.numLabel), dtype=np.float32)
for i, j in itertools.combinations(range(self.numLabel), 2):
# find shortest bridge
dist, idx = trees[i].query(regions[j].coords)
idx_min = np.argmin(dist)
yxj = regions[j].coords[idx_min,:]
yxi = regions[i].coords[idx[idx_min],:]
dist_min = dist[idx_min]
# save
n0, n1 = str(i+1), str(j+1)
conn = dict()
conn[n0] = yxi
conn[n1] = yxj
conn['distance'] = dist_min
self.connDict[f'{n0}_{n1}'] = conn
self.distMat[i,j] = self.distMat[j,i] = dist_min
return self.connDict, self.distMat
def find_mst_bridge(self):
"""Search for bridges to connect all labeled areas using the minimum spanning tree algorithm
Returns: bridges : list of dict, i.e.:
[{'label0': 1, 'label1': 3, 'x0': 345, 'x1': 191, 'y0': 1232, 'y1': 1089},
{'label0': 1, 'label1': 4, 'x0': 1143, 'x1': 1157, 'y0': 1204, 'y1': 1217},
{'label0': 1, 'label1': 5, 'x0': 557, 'x1': 565, 'y0': 1263, 'y1': 1270},
{'label0': 3, 'label1': 2, 'x0': 249, 'x1': 239, 'y0': 891, 'y1': 868}]
"""
if not hasattr(self, 'distMat'):
self.get_all_bridge()
# MST bridges with breadth_first_order
distMatMst = csg.minimum_spanning_tree(self.distMat)
succs, preds = csg.breadth_first_order(distMatMst, i_start=self.labelRef-1, directed=False)
# save to self.bridges
self.bridges = []
for i in range(1, succs.size):
n0 = preds[succs[i]] + 1
n1 = succs[i] + 1
# read conn
if n0 > n1:
nn = [str(n1), str(n0)]
else:
nn = [str(n0), str(n1)]
conn = self.connDict[f'{nn[0]}_{nn[1]}']
y0, x0 = conn[str(n0)]
y1, x1 = conn[str(n1)]
# save bdg
bridge = dict()
bridge['x0'] = x0
bridge['y0'] = y0
bridge['x1'] = x1
bridge['y1'] = y1
bridge['label0'] = n0
bridge['label1'] = n1
bridge['distance'] = ((x1 - x0)**2 + (y1 - y0)**2)**0.5
self.bridges.append(bridge)
self.num_bridge = len(self.bridges)
return self.bridges
def get_bridge_endpoint_aoi_mask(self, bridge, radius=50):
# get AOI mask
x0, y0 = bridge['x0'], bridge['y0']
x1, y1 = bridge['x1'], bridge['y1']
x00 = max(0, x0 - radius); x01 = min(self.width, x0 + radius)
y00 = max(0, y0 - radius); y01 = min(self.length, y0 + radius)
x10 = max(0, x1 - radius); x11 = min(self.width, x1 + radius)
y10 = max(0, y1 - radius); y11 = min(self.length, y1 + radius)
aoi_mask0 = np.zeros(self.labelImg.shape, dtype=np.bool_)
aoi_mask1 = np.zeros(self.labelImg.shape, dtype=np.bool_)
aoi_mask0[y00:y01, x00:x01] = True
aoi_mask1[y10:y11, x10:x11] = True
return aoi_mask0, aoi_mask1
def unwrap_conn_comp(self, unw, radius=50, ramp_type=None, print_msg=False):
start_time = time.time()
radius = int(min(radius, min(self.conncomp.shape)*0.05))
unw = np.array(unw, dtype=np.float32)
if self.refY is not None:
unw[unw != 0.] -= unw[self.refY, self.refX]
if ramp_type is not None:
if print_msg:
print(f'estimate a {ramp_type} ramp')
ramp_mask = (self.labelImg == self.labelRef)
unw, ramp = deramp(unw, ramp_mask, ramp_type, metadata=self.metadata)
for bridge in self.bridges:
# prepare masks
aoi_mask0, aoi_mask1 = self.get_bridge_endpoint_aoi_mask(bridge, radius=radius)
label_mask0 = self.labelImg == bridge['label0']
label_mask1 = self.labelImg == bridge['label1']
# get phase difference
value0 = np.nanmedian(unw[aoi_mask0 * label_mask0])
value1 = np.nanmedian(unw[aoi_mask1 * label_mask1])
diff_value = value1 - value0
# estimate integer number of phase jump
num_jump = (np.abs(diff_value) + np.pi) // (2.*np.pi)
if diff_value > 0:
num_jump *= -1
# add phase jump
unw[label_mask1] += 2.* np.pi * num_jump
if print_msg:
print(('phase diff {}_{}: {:04.1f} rad --> '
'num of jump: {}').format(bridge['label1'],
bridge['label0'],
diff_value,
num_jump))
# add ramp back
if ramp_type is not None:
unw += ramp
if print_msg:
print(f'time used: {time.time()-start_time:.2f} secs.')
return unw
def plot_bridge(self, ax, cmap='jet', radius=50):
# background label
ax.imshow(self.labelImg, cmap=cmap, interpolation='nearest')
for bridge in self.bridges:
# bridges
ax.plot([bridge['x0'], bridge['x1']],
[bridge['y0'], bridge['y1']], 'w-', lw=1)
# endpoint window
if radius > 0:
aoi_mask0, aoi_mask1 = self.get_bridge_endpoint_aoi_mask(bridge, radius=radius)
label_mask0 = self.labelImg == bridge['label0']
label_mask1 = self.labelImg == bridge['label1']
# Note by Emre, Mar 2024: overlay bridge regions directly using plot() function,
# to save memory while calling this func in a loop by avoiding creating separate
# mask arrays (https://github.com/insarlab/MintPy/pull/1155)
ax.plot(np.nonzero(aoi_mask0*label_mask0)[1], np.nonzero(aoi_mask0*label_mask0)[0], 'gray', alpha=0.3)
ax.plot(np.nonzero(aoi_mask1*label_mask1)[1], np.nonzero(aoi_mask1*label_mask1)[0], 'gray', alpha=0.3)
# reference pixel
ax.plot(self.refX, self.refY, 'ks', ms=2)
return ax
######################################## end of connectComponent class ####################################