-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcolor_separate_functions.py
185 lines (138 loc) · 6.76 KB
/
color_separate_functions.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
"""
This module contains functions for color separation, nuclei segmentation, and analysis of TMA core images.
Functions:
- `analyze_TMA_core`: Performs color separation, nuclei segmentation, and intensity analysis on a TMA core image.
- `color_separate_only`: Performs color separation without nuclei segmentation, useful for previewing channels.
- `color_separate`: Separates H, E, D stains and optionally combines H and D for fluorescence effect.
- `segment_nuclei`: Segments nuclei in the H component using StarDist.
- `filter_objects`: Filters segmented nuclei based on size.
- `brown_intensity_from_nuclei`: Extracts brown intensity from nuclei locations.
"""
import numpy as np
from skimage import img_as_ubyte
from skimage.color import rgb2hed, hed2rgb
from skimage.exposure import rescale_intensity
from stardist.models import StarDist2D
from stardist.plot import render_label
from csbdeep.utils import normalize
from skimage import measure
import pandas as pd
def analyze_TMA_core(image, min_size=30, max_size=1000):
"""
Analyzes a TMA core image by performing color separation, nuclei segmentation, and filtering based on size.
Args:
- image: Input TMA core image.
- min_size: Minimum size of nuclei to consider.
- max_size: Maximum size of nuclei to consider.
Returns:
- Tuple of image, H component, brown component, filtered segmentation image, means, and standard deviations.
"""
print("Now color separating ...")
H, _, brown, _ = color_separate(image)
print("Now segmenting nuclei using StarDist ...")
nuclei_segm = segment_nuclei(H)
print("Now filtering nuclei based on user-defined min and max ...")
filtered_segm_image = filter_objects(nuclei_segm, image, min_size, max_size)
print("Now extracting mean and std dev. brown intensity values from Nuclei ...")
means, stdevs = brown_intensity_from_nuclei(filtered_segm_image, brown)
return image, H, brown, filtered_segm_image, means, stdevs
def color_separate_only(image):
"""
Performs color separation without segmentation, useful for previewing separated channels.
Args:
- image: Input RGB image.
Returns:
- Tuple of input image, H component, and brown component.
"""
print("Now color separating ...")
H, _, brown, _ = color_separate(image)
return image, H, brown
def color_separate(image_rgb):
"""
Performs color separation on an RGB image.
Args:
- image_rgb: Input RGB image.
Returns:
- Tuple of H, E, D components, and a combined image for visualization.
"""
# Convert the RGB image to HED using the prebuilt skimage method
image_hed = rgb2hed(image_rgb)
# Create an RGB image for each of the separated stains
# Convert them to ubyte for easy saving to drive as an image
null = np.zeros_like(image_hed[:, :, 0])
image_h = img_as_ubyte(hed2rgb(np.stack((image_hed[:, :, 0], null, null), axis=-1)))
image_e = img_as_ubyte(hed2rgb(np.stack((null, image_hed[:, :, 1], null), axis=-1)))
image_d = img_as_ubyte(hed2rgb(np.stack((null, null, image_hed[:, :, 2]), axis=-1)))
# Optional fun exercise of combining H and DAB stains into a single image with fluorescence look
h = rescale_intensity(image_hed[:, :, 0], out_range=(0, 1),
in_range=(0, np.percentile(image_hed[:, :, 0], 99)))
d = rescale_intensity(image_hed[:, :, 2], out_range=(0, 1),
in_range=(0, np.percentile(image_hed[:, :, 2], 99)))
# Cast the two channels into an RGB image, as the blue and green channels
# Convert to ubyte for easy saving as image to local drive
zdh = img_as_ubyte(np.dstack((null, d, h))) # DAB in green and H in Blue
return image_h, image_e, image_d, zdh
def segment_nuclei(H_img):
"""
Segments nuclei using StarDist.
Args:
- H_img: Input H component image.
Returns:
- Nuclei segmentation labels.
"""
# Define a pretrained model to segment nuclei in fluorescence images (download from pretrained)
model = StarDist2D.from_pretrained('2D_versatile_fluo')
H_inverted = np.invert(H_img)
H_gray = H_inverted[:, :, 0]
H_labels, _ = model.predict_instances(normalize(H_gray), nms_thresh=0.3, prob_thresh=0.6)
return H_labels
def filter_objects(H_labels, intensity_image, min_size=30, max_size=1000):
"""
Filters segmented objects based on size.
Args:
- H_labels: Nuclei segmentation labels.
- intensity_image: Original intensity image.
- min_size: Minimum size of nuclei to consider.
- max_size: Maximum size of nuclei to consider.
Returns:
- Filtered segmentation labels.
"""
props = measure.regionprops_table(H_labels, intensity_image,
properties=['label',
'area', 'equivalent_diameter',
'mean_intensity', 'solidity', 'centroid'])
df = pd.DataFrame(props)
# Filter objects by size
df_filtered = df[df.area > min_size]
df_filtered = df_filtered[df_filtered.area < max_size]
# Filter objects from the labeled image
useful_labels = df_filtered.label.values # Labels of objects that passed our filter criteria
filtered_H_labels = np.zeros_like(H_labels) # Array same size as labeled image but with 0s, to be filled later
for i in range(H_labels.shape[0]):
for j in range(H_labels.shape[1]):
if (H_labels[i, j] in useful_labels) == True:
filtered_H_labels[i, j] = H_labels[i, j]
else:
filtered_H_labels[i, j] = 0
return filtered_H_labels
def brown_intensity_from_nuclei(filtered_H_labels, brown_image):
"""
Extracts brown intensity from nuclei locations.
Args:
- filtered_H_labels: Filtered nuclei segmentation labels.
- brown_image: Brown intensity image.
Returns:
- Tuple of means and standard deviations of brown intensity.
"""
DAB_intensity_from_H_filtered = measure.regionprops_table(filtered_H_labels, brown_image,
properties=['label', 'mean_intensity'])
df_DAB_intensity = pd.DataFrame(DAB_intensity_from_H_filtered)
mean_R = df_DAB_intensity["mean_intensity-0"].mean()
mean_G = df_DAB_intensity["mean_intensity-1"].mean()
mean_B = df_DAB_intensity["mean_intensity-2"].mean()
std_R = df_DAB_intensity["mean_intensity-0"].std()
std_G = df_DAB_intensity["mean_intensity-1"].std()
std_B = df_DAB_intensity["mean_intensity-2"].std()
means = [mean_R, mean_G, mean_B]
stdevs = [std_R, std_G, std_B]
return means, stdevs