-
Notifications
You must be signed in to change notification settings - Fork 54
/
zernike.py
305 lines (257 loc) · 12.2 KB
/
zernike.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
#!/usr/bin/python
"""
Compute the Zernike moments of a collection of points.
Authors:
- Arthur Mikhno, 2013, Columbia University (original MATLAB code)
- Brian Rossa, 2013, Tank Think Labs, LLC (port to Python)
- Arno Klein, 2013 (arno@mindboggle.info) http://binarybottle.com
Copyright 2013, Mindboggle team (http://mindboggle.info), Apache v2.0 License
"""
def zernike_moments(points, faces, order=10, scale_input=True,
decimate_fraction=0, decimate_smooth=0, verbose=False):
"""
Compute the Zernike moments of a surface patch of points and faces.
Optionally decimate the input mesh.
Note::
Decimation sometimes leads to an error of "Segmentation fault: 11"
(Twins-2-1 left label 14 gives such an error only when decimated.)
Parameters
----------
points : list of lists of 3 floats
x,y,z coordinates for each vertex
faces : list of lists of 3 integers
each list contains indices to vertices that form a triangle on a mesh
order : integer
order of the moments being calculated
scale_input : bool
translate and scale each object so it is bounded by a unit sphere?
(this is the expected input to zernike_moments())
decimate_fraction : float
fraction of mesh faces to remove for decimation (0 for no decimation)
decimate_smooth : integer
number of smoothing steps for decimation
verbose : bool
print statements?
Returns
-------
descriptors : list of floats
Zernike descriptors
Examples
--------
>>> # Example 1: simple cube (decimation results in a Segmentation Fault):
>>> import numpy as np
>>> from mindboggle.shapes.zernike.zernike import zernike_moments
>>> points = [[0,0,0], [1,0,0], [0,0,1], [0,1,1],
... [1,0,1], [0,1,0], [1,1,1], [1,1,0]]
>>> faces = [[0,2,4], [0,1,4], [2,3,4], [3,4,5], [3,5,6], [0,1,7]]
>>> order = 3
>>> scale_input = True
>>> decimate_fraction = 0
>>> decimate_smooth = 0
>>> verbose = False
>>> descriptors = zernike_moments(points, faces, order, scale_input,
... decimate_fraction, decimate_smooth, verbose)
>>> [np.float("{0:.{1}f}".format(x, 5)) for x in descriptors]
[0.09189, 0.09357, 0.04309, 0.06466, 0.0382, 0.04138]
Example 2: Twins-2-1 left postcentral pial surface -- NO decimation:
(zernike_moments took 142 seconds for order = 3 with no decimation)
>>> from mindboggle.shapes.zernike.zernike import zernike_moments
>>> from mindboggle.mio.vtks import read_vtk
>>> from mindboggle.guts.mesh import keep_faces
>>> from mindboggle.mio.fetch_data import prep_tests
>>> urls, fetch_data = prep_tests()
>>> label_file = fetch_data(urls['left_freesurfer_labels'], '', '.vtk')
>>> points, f1,f2, faces, labels, f3,f4,f5 = read_vtk(label_file)
>>> I22 = [i for i,x in enumerate(labels) if x==1022] # postcentral
>>> faces = keep_faces(faces, I22)
>>> order = 3
>>> scale_input = True
>>> decimate_fraction = 0
>>> decimate_smooth = 0
>>> verbose = False
>>> descriptors = zernike_moments(points, faces, order, scale_input,
... decimate_fraction, decimate_smooth, verbose)
>>> [np.float("{0:.{1}f}".format(x, 5)) for x in descriptors]
[0.00471, 0.0084, 0.00295, 0.00762, 0.0014, 0.00076]
Example 3: left postcentral + pars triangularis pial surfaces:
>>> from mindboggle.mio.vtks import read_vtk, write_vtk
>>> points, f1,f2, faces, labels, f3,f4,f5 = read_vtk(label_file)
>>> I20 = [i for i,x in enumerate(labels) if x==1020] # pars triangularis
>>> I22 = [i for i,x in enumerate(labels) if x==1022] # postcentral
>>> I22.extend(I20)
>>> faces = keep_faces(faces, I22)
>>> order = 3
>>> scale_input = True
>>> decimate_fraction = 0
>>> decimate_smooth = 0
>>> verbose = False
>>> descriptors = zernike_moments(points, faces, order, scale_input,
... decimate_fraction, decimate_smooth, verbose)
>>> [np.float("{0:.{1}f}".format(x, 5)) for x in descriptors]
[0.00586, 0.00973, 0.00322, 0.00818, 0.0013, 0.00131]
View both segments (skip test):
>>> from mindboggle.mio.plots import plot_surfaces # doctest: +SKIP
>>> from mindboggle.mio.vtks import rewrite_scalars # doctest: +SKIP
>>> scalars = -1 * np.ones(np.shape(labels)) # doctest: +SKIP
>>> scalars[I22] = 1 # doctest: +SKIP
>>> rewrite_scalars(label_file, 'test_two_labels.vtk', scalars,
... 'two_labels', scalars) # doctest: +SKIP
>>> plot_surfaces(vtk_file) # doctest: +SKIP
"""
import numpy as np
from mindboggle.guts.mesh import reindex_faces_0to1
from mindboggle.guts.mesh import decimate
from mindboggle.shapes.zernike.pipelines import DefaultPipeline as Pipeline
# Convert 0-indices (Python) to 1-indices (Matlab) for all face indices:
index1 = False # already done elsewhere in the code
if index1:
faces = reindex_faces_0to1(faces)
# Convert lists to numpy arrays:
if isinstance(points, list):
points = np.array(points)
if isinstance(faces, list):
faces = np.array(faces)
# ------------------------------------------------------------------------
# Translate all points so that they are centered at their mean,
# and scale them so that they are bounded by a unit sphere:
# ------------------------------------------------------------------------
if scale_input:
center = np.mean(points, axis=0)
points = points - center
maxd = np.max(np.sqrt(np.sum(points**2, axis=1)))
points /= maxd
# ------------------------------------------------------------------------
# Decimate surface:
# ------------------------------------------------------------------------
if 0 < decimate_fraction < 1:
points, faces, u1,u2 = decimate(points, faces,
decimate_fraction, decimate_smooth, [], save_vtk=False)
# Convert lists to numpy arrays:
points = np.array(points)
faces = np.array(faces)
# ------------------------------------------------------------------------
# Multiprocessor pipeline:
# ------------------------------------------------------------------------
pl = Pipeline()
# ------------------------------------------------------------------------
# Geometric moments:
# ------------------------------------------------------------------------
G = pl.geometric_moments_exact(points, faces, order)
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
Z = pl.zernike(G, order)
# ------------------------------------------------------------------------
# Extract Zernike descriptors:
# ------------------------------------------------------------------------
descriptors = pl.feature_extraction(Z, order).tolist()
if verbose:
print("Zernike moments: {0}".format(descriptors))
return descriptors
def zernike_moments_per_label(vtk_file, order=10, exclude_labels=[-1],
scale_input=True, decimate_fraction=0,
decimate_smooth=25, verbose=False):
"""
Compute the Zernike moments per labeled region in a file.
Optionally decimate the input mesh.
Parameters
----------
vtk_file : string
name of VTK surface mesh file containing index scalars (labels)
order : integer
number of moments to compute
exclude_labels : list of integers
labels to be excluded
scale_input : bool
translate and scale each object so it is bounded by a unit sphere?
(this is the expected input to zernike_moments())
decimate_fraction : float
fraction of mesh faces to remove for decimation (1 for no decimation)
decimate_smooth : integer
number of smoothing steps for decimation
verbose : bool
print statements?
Returns
-------
descriptors_lists : list of lists of floats
Zernike descriptors per label
label_list : list of integers
list of unique labels for which moments are computed
Examples
--------
>>> # Zernike moments per label of a FreeSurfer-labeled left cortex.
>>> # Uncomment "if label==22:" below to run example
>>> # for left postcentral (22) pial surface:
>>> import numpy as np
>>> from mindboggle.shapes.zernike.zernike import zernike_moments_per_label
>>> from mindboggle.mio.fetch_data import prep_tests
>>> urls, fetch_data = prep_tests()
>>> vtk_file = fetch_data(urls['left_freesurfer_labels'], '', '.vtk')
>>> order = 3
>>> exclude_labels = [-1]
>>> scale_input = True
>>> verbose = False
>>> descriptors_lists, label_list = zernike_moments_per_label(vtk_file,
... order, exclude_labels, scale_input, verbose)
>>> label_list[0:10]
[999, 1001, 1002, 1003, 1005, 1006, 1007, 1008, 1009, 1010]
>>> [np.float("{0:.{1}f}".format(x, 5)) for x in descriptors_lists[0]]
[0.00587, 0.01143, 0.0031, 0.00881, 0.00107, 0.00041]
>>> [np.float("{0:.{1}f}".format(x, 5)) for x in descriptors_lists[1]]
[4e-05, 9e-05, 3e-05, 9e-05, 2e-05, 1e-05]
>>> [np.float("{0:.{1}f}".format(x, 5)) for x in descriptors_lists[2]]
[0.00144, 0.00232, 0.00128, 0.00304, 0.00084, 0.00051]
>>> [np.float("{0:.{1}f}".format(x, 5)) for x in descriptors_lists[3]]
[0.00393, 0.006, 0.00371, 0.00852, 0.00251, 0.00153]
>>> [np.float("{0:.{1}f}".format(x, 5)) for x in descriptors_lists[4]]
[0.00043, 0.0003, 0.00095, 0.00051, 0.00115, 0.00116]
"""
import numpy as np
from mindboggle.mio.vtks import read_vtk
from mindboggle.guts.mesh import keep_faces
from mindboggle.shapes.zernike.zernike import zernike_moments
min_points_faces = 4
# ------------------------------------------------------------------------
# Read VTK surface mesh file:
# ------------------------------------------------------------------------
points, indices, lines, faces, labels, scalar_names, npoints, \
input_vtk = read_vtk(vtk_file)
# ------------------------------------------------------------------------
# Loop through labeled regions:
# ------------------------------------------------------------------------
ulabels = [x for x in np.unique(labels) if x not in exclude_labels]
label_list = []
descriptors_lists = []
for label in ulabels:
#if label == 1022: # 22:
# print("DEBUG: COMPUTE FOR ONLY ONE LABEL")
# --------------------------------------------------------------------
# Determine the indices per label:
# --------------------------------------------------------------------
Ilabel = [i for i,x in enumerate(labels) if x == label]
if verbose:
print(' {0} vertices for label {1}'.format(len(Ilabel), label))
if len(Ilabel) > min_points_faces:
# ----------------------------------------------------------------
# Remove background faces:
# ----------------------------------------------------------------
pick_faces = keep_faces(faces, Ilabel)
if len(pick_faces) > min_points_faces:
# ------------------------------------------------------------
# Compute Zernike moments for the label:
# ------------------------------------------------------------
descriptors = zernike_moments(points, pick_faces,
order, scale_input,
decimate_fraction,
decimate_smooth, verbose)
# ------------------------------------------------------------
# Append to a list of lists of spectra:
# ------------------------------------------------------------
descriptors_lists.append(descriptors)
label_list.append(label)
return descriptors_lists, label_list
# ============================================================================
# Doctests
# ============================================================================
if __name__ == "__main__":
import doctest
doctest.testmod(verbose=True) # py.test --doctest-modules