-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathsegment_3D.py
230 lines (177 loc) · 8.32 KB
/
segment_3D.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
import operator
import shutil
from pathlib import Path
import numpy as np
import pandas as pd
import skimage.filters
from cellpose import models
from roifile import ImagejRoi
from scipy import ndimage
from skimage.feature import peak_local_max
from skimage.measure import find_contours, regionprops
from skimage.morphology import watershed, remove_small_objects
def segment_watershed(img, sigma: float = 2.0, threshold: float = 25.0, min_size: float = 45.0, do_3d: bool = True):
""" filter and get labeled 3D mask with a watershed segmentation
Parameters
---------
img : array
image should be of shape Z*Y*X (if 3D) or Y*X (if 2D)
sigma : float (optional, default = 2.0)
standard deviation for gaussian filtering
threshold : float (optional, default = 25.0)
value of pixel intensity to create a binary image
min_size : float (optional, default = 45.0)
the smallest allowable object size
do_3d : bool (default =True)
specify if the segmentation is done in 3D or in 2D
Return
---------
single 3D array with labeled objects
0 = no masks // 1,2,... = mask labels
"""
if do_3d:
img_filter = skimage.filters.gaussian(img, sigma, preserve_range=True, multichannel=False) # gaussian filter
masks = []
for z in range(0, img_filter.shape[0]):
mask = img_filter[z, :, :] > threshold
mask = ndimage.binary_fill_holes(mask)
mask = remove_small_objects(mask.astype(bool), min_size)
masks.append(mask)
mask3d = np.stack(masks)
# watershed segmentation
distance = ndimage.distance_transform_edt(mask3d)
local_maxi = peak_local_max(distance, indices=False ,footprint=np.ones((15,15,15)), labels=mask3d)
markers = ndimage.label(local_maxi, structure=np.ones((3, 3, 3)))[0]
labels = watershed(-distance, markers, mask=mask3d)
else:
img_filter = skimage.filters.gaussian(img, sigma, preserve_range=True, multichannel=False)
mask = img_filter > threshold
mask = ndimage.binary_fill_holes(mask)
mask = remove_small_objects(mask.astype(bool), min_size)
# watershed segmentation
distance = ndimage.distance_transform_edt(mask)
local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((15,15)), labels=mask)
markers = ndimage.label(local_maxi, structure=np.ones((3, 3)))[0]
labels = watershed(-distance, markers, mask=mask)
return labels
def segment_cellpose(img, mode: str = 'cyto', do_3d: bool = True):
""" filter and get labeled 2D or 3D mask with cellpose
Parameter
---------
img : array
image should be of shape Z*Y*X (if 3D) or Y*X (if 2D)
mode : str (default = 'cyto')
Run either 'cyto' segmentation or 'nuclei' segmentation
do_3d : bool (default =True)
specify if the segmentation is done in 3D or in 2D
Return
---------
single 3D array with labeled objects
0 = no masks // 1,2,... = mask labels
"""
model = models.Cellpose(gpu=False, model_type=mode, torch=True)
channels = [0, 0]
masks, flows, styles, diams = model.eval(img, diameter=None, channels=channels, do_3D=do_3d, z_axis=0)
return masks
def ROIs_archive(img_label, rois_path: str, filename: str = 'ROIset', do_3d: bool = True):
""" save ROIs archive, ROIs are recovered from their most representative Z (plane)
saved to filename + '.zip'
Parameter
---------
img_label : array
2D or 3D array with labeled objects, image should be of shape Z*Y*X (if 3D) or Y*X (if 2D)
rois_path : str
where the archive will be saved
filename : str (optional, default = 'ROIset')
name of archive
do_3d : bool (default = True)
specify if img_label is 3D or 2D array
"""
""" find areas of each label for every Z if 3D array"""
if do_3d:
dico = {}
for z in range(0, img_label.shape[0]):
props = regionprops(img_label[z, :, :].astype(int))
for obj in props:
if obj.label not in dico:
dico[obj.label] = {}
dico[obj.label][z] = obj.area
dico = sorted(dico.items(), key=lambda t: t[0])
""" find the most representative Z for each label """
typicalZ = []
for label in dico:
maxi = max(label[1].items(), key=operator.itemgetter(1))[0]
typicalZ.append(maxi)
labels_mask = np.swapaxes(img_label, 1, 2)
""" save filename.zip at rois_path """
directory = rois_path
path = Path(directory)
pathROIset = path / filename
pathROIset.mkdir(parents=True, exist_ok=False)
for i, value in enumerate(typicalZ, np.unique(labels_mask.astype(int))[1]):
contour = find_contours(labels_mask[value, :, :] == i, level=0)
roi = ImagejRoi.frompoints(contour[0])
roi.tofile(str(path) + '/' + f'roi_{i:04d}.roi')
filepath = str(path / f'roi_{i:04d}.roi')
shutil.move(filepath, pathROIset)
shutil.make_archive(str(path) + '/' + filename, 'zip', str(path) + '/' + filename)
shutil.rmtree(str(path) + '/' + filename)
else:
labels_mask = np.swapaxes(img_label, 0, 1)
directory = rois_path
path = Path(directory)
pathROIset = path / filename
pathROIset.mkdir(parents=True, exist_ok=False)
for i in np.unique(labels_mask)[1:]:
contour = find_contours(labels_mask == i, level=0)
roi = ImagejRoi.frompoints(contour[0])
roi.tofile(str(path) + '/' + f'roi_{i:04d}.roi')
filepath = str(path / f'roi_{i:04d}.roi')
shutil.move(filepath, pathROIset)
shutil.make_archive(str(path) + '/' + filename, 'zip', str(path) + '/' + filename)
shutil.rmtree(str(path) + '/' + filename)
def result_file(mask_label, channel0, channel1, result_path: str = "", filename: str = 'result',
do_3d: bool = True):
""" save ROIs archive, ROIs are recovered from their most representative Z (plane)
saved to filename + '.csv'
Parameters
---------
mask_label : array
2D or 3D array with labeled objects, image should be of shape Z*Y*X (if 3D) or Y*X (if 2D)
channel0 : array
first channel to measure area and mean intensity
2D or 3D array, image should be of shape Z*Y*X (if 3D) or Y*X (if 2D)
channel1 : array
second channel to measure area and mean intensity
2D or 3D array, image should be of shape Z*Y*X (if 3D) or Y*X (if 2D)
result_path : str (optional, default = 'result')
filename : str (optional, default = 'ROIset')
name of file
do_3d : bool (default = True)
specify if img_label is 3D or 2D array
"""
mask_label = mask_label.astype(int)
if do_3d:
# Generate .csv file with intensities of 2 channels
region0 = skimage.measure.regionprops_table(mask_label, intensity_image=channel0,
properties=['label', 'mean_intensity', 'area'])
region1 = skimage.measure.regionprops_table(mask_label, intensity_image=channel1,
properties=['label', 'mean_intensity', 'area'])
region1['mean_intensity1'] = region1.pop('mean_intensity')
region1['area1'] = region1.pop('area')
region0.update(region1)
df = pd.DataFrame(region0)
df.to_csv(result_path + '/' + filename + '.csv', float_format='%.4f',
header=['ID', 'Int_TYPE1', 'Area_TYPE1', 'Int_TYPE2_N', 'Area_TYPE1'], index=False, sep='\t')
else:
#Generate .csv file with intensities of 2 channels
region0 = skimage.measure.regionprops_table(mask_label, intensity_image=channel0,
properties=['label', 'mean_intensity', 'area'])
region1 = skimage.measure.regionprops_table(mask_label, intensity_image=channel1,
properties=['label', 'mean_intensity', 'area'])
region1['mean_intensity1'] = region1.pop('mean_intensity')
region1['area1'] = region1.pop('area')
region0.update(region1)
df = pd.DataFrame(region0)
df.to_csv(result_path + '/' + filename + '.csv', float_format='%.4f',
header=['ID', 'Int_TYPE1', 'Area_TYPE1', 'Int_TYPE2_N', 'Area_TYPE1'], index=False, sep='\t')