-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathsupplementary_code.py
217 lines (163 loc) · 5.83 KB
/
supplementary_code.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
import matplotlib.pyplot as plt
import numpy as np
from ipywidgets import interact
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from skimage import exposure, io, measure
def show_plane(axis, plane, cmap="gray", title=None):
"""Shows a specific plane within 3D data.
"""
axis.imshow(plane, cmap=cmap)
axis.set_xticks([])
axis.set_yticks([])
if title:
axis.set_title(title)
return None
def slice_in_3d(axis, shape, plane):
"""Draws a cube in a 3D plot.
Parameters
----------
axis : matplotlib.Axes
A matplotlib axis to be drawn.
shape : tuple or array (1, 3)
Shape of the input data.
plane : int
Number of the plane to be drawn.
Notes
-----
Originally from:
https://stackoverflow.com/questions/44881885/python-draw-parallelepiped
"""
Z = np.array([[0, 0, 0],
[1, 0, 0],
[1, 1, 0],
[0, 1, 0],
[0, 0, 1],
[1, 0, 1],
[1, 1, 1],
[0, 1, 1]])
Z = Z * shape
r = [-1, 1]
X, Y = np.meshgrid(r, r)
# plotting vertices
axis.scatter3D(Z[:, 0], Z[:, 1], Z[:, 2])
# list of sides' polygons of figure
verts = [[Z[0], Z[1], Z[2], Z[3]],
[Z[4], Z[5], Z[6], Z[7]],
[Z[0], Z[1], Z[5], Z[4]],
[Z[2], Z[3], Z[7], Z[6]],
[Z[1], Z[2], Z[6], Z[5]],
[Z[4], Z[7], Z[3], Z[0]],
[Z[2], Z[3], Z[7], Z[6]]]
# plotting sides
axis.add_collection3d(
Poly3DCollection(verts,
facecolors=(0, 1, 1, 0.25),
linewidths=1,
edgecolors='darkblue')
)
verts = np.array([[[0, 0, 0],
[0, 0, 1],
[0, 1, 1],
[0, 1, 0]]])
verts = verts * shape
verts += [plane, 0, 0]
axis.add_collection3d(
Poly3DCollection(verts,
facecolors='magenta',
linewidths=1,
edgecolors='black')
)
axis.set_xlabel('plane')
axis.set_ylabel('col')
axis.set_zlabel('row')
# auto-scale plot axes
scaling = np.array([getattr(axis, 'get_{}lim'.format(dim))() for dim in 'xyz'])
axis.auto_scale_xyz(*[[np.min(scaling), np.max(scaling)]] * 3)
return None
def slice_explorer(data, cmap='gray'):
"""Allows to explore 2D slices in 3D data.
Parameters
----------
data : array (M, N, P)
3D interest image.
cmap : str (optional)
A string referring to one of matplotlib's colormaps.
"""
data_len = len(data)
@interact(plane=(0, data_len-1), continuous_update=False)
def display_slice(plane=data_len/2):
fig, axis = plt.subplots(figsize=(20, 7))
axis_3d = fig.add_subplot(133, projection='3d')
show_plane(axis, data[plane], title='Plane {}'.format(plane), cmap=cmap)
slice_in_3d(axis=axis_3d, shape=data.shape, plane=plane)
plt.show()
return display_slice
def display(data, cmap="gray", step=2):
_, axes = plt.subplots(nrows=5, ncols=6, figsize=(16, 14))
# getting data min and max to plot limits
vmin, vmax = data.min(), data.max()
for axis, image in zip(axes.flatten(), data[::step]):
axis.imshow(image, cmap=cmap, vmin=vmin, vmax=vmax)
axis.set_xticks([])
axis.set_yticks([])
return None
def plot_hist(axis, data, title=None):
"""Helper function for plotting histograms.
"""
axis.hist(data.ravel(), bins=256)
axis.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
if title:
axis.set_title(title)
return None
def plot_3d_surface(data, labels, region=3, spacing=(1.0, 1.0, 1.0)):
"""Generates a 3D surface plot for the specified region.
Parameters
----------
data : array (M, N, P)
3D interest image.
labels : array (M, N, P)
Labels corresponding to data, obtained by measure.label.
region : int, optional
The region of interest to be plotted.
spacing : array (1, 3)
Spacing information, set distances between pixels.
Notes
-----
The volume is visualized using the mesh vertexes and faces.
"""
properties = measure.regionprops(labels, intensity_image=data)
# skimage.measure.marching_cubes expects ordering (row, col, plane).
# We need to transpose the data:
volume = (labels == properties[region].label).transpose(1, 2, 0)
verts_px, faces_px, _, _ = measure.marching_cubes_lewiner(volume, level=0, spacing=(1.0, 1.0, 1.0))
surface_area_pixels = measure.mesh_surface_area(verts_px, faces_px)
verts_actual, faces_actual, _, _ = measure.marching_cubes_lewiner(volume, level=0, spacing=tuple(spacing))
surface_area_actual = measure.mesh_surface_area(verts_actual, faces_actual)
print('Surface area\n')
print(' * Total pixels: {:0.2f}'.format(surface_area_pixels))
print(' * Actual: {:0.2f}'.format(surface_area_actual))
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
mesh = Poly3DCollection(verts_px[faces_px])
mesh.set_edgecolor('black')
ax.add_collection3d(mesh)
ax.set_xlabel('col')
ax.set_ylabel('row')
ax.set_zlabel('plane')
min_pln, min_row, min_col, max_pln, max_row, max_col = properties[region].bbox
ax.set_xlim(min_row, max_row)
ax.set_ylim(min_col, max_col)
ax.set_zlim(min_pln, max_pln)
plt.tight_layout()
plt.show()
return None
def results_from_part_1():
data = io.imread("data/cells.tif")
vmin, vmax = np.percentile(data, q=(0.5, 99.5))
rescaled = exposure.rescale_intensity(
data,
in_range=(vmin, vmax),
out_range=np.float32
)
equalized = exposure.equalize_hist(data)
return data, rescaled, equalized