-
Notifications
You must be signed in to change notification settings - Fork 14
/
surface_eigenmodes.py
244 lines (197 loc) · 10.2 KB
/
surface_eigenmodes.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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Calculate the eigenmodes of a cortical surface
@author: James Pang and Kevin Aquino, Monash University, 2022
"""
# Import all the libraries
from lapy import Solver, TriaIO
import numpy as np
import nibabel as nib
import brainspace.mesh as mesh
import os
from argparse import ArgumentParser
def calc_eig(tria, num_modes):
"""Calculate the eigenvalues and eigenmodes of a surface.
Parameters
----------
tria : lapy compatible object
Loaded vtk object corresponding to a surface triangular mesh
num_modes : int
Number of eigenmodes to be calculated
Returns
------
evals : array (num_modes x 1)
Eigenvalues
emodes : array (number of surface points x num_modes)
Eigenmodes
"""
fem = Solver(tria)
evals, emodes = fem.eigs(k=num_modes)
return evals, emodes
def create_temp_surface(surface_input, surface_output_filename):
"""Write surface to a new vtk file.
Parameters
----------
surface_input : brainspace compatible object
Loaded vtk object corresponding to a surface triangular mesh
surface_output_filename : str
Filename of surface to be saved
"""
f = open(surface_output_filename, 'w')
f.write('# vtk DataFile Version 2.0\n')
f.write(surface_output_filename + '\n')
f.write('ASCII\n')
f.write('DATASET POLYDATA\n')
f.write('POINTS ' + str(np.shape(surface_input.Points)[0]) + ' float\n')
for i in range(np.shape(surface_input.Points)[0]):
f.write(' '.join(map(str, np.array(surface_input.Points[i, :]))))
f.write('\n')
f.write('\n')
f.write('POLYGONS ' + str(np.shape(surface_input.polys2D)[0]) + ' ' + str(4* np.shape(surface_input.polys2D)[0]) + '\n')
for i in range(np.shape(surface_input.polys2D)[0]):
f.write(' '.join(map(str, np.append(3, np.array(surface_input.polys2D[i, :])))))
f.write('\n')
f.close()
def get_indices(surface_original, surface_new):
"""Extract indices of vertices of the two surfaces that match.
Parameters
----------
surface_original : brainspace compatible object
Loaded vtk object corresponding to a surface triangular mesh
surface_new : brainspace compatible object
Loaded vtk object corresponding to a surface triangular mesh
Returns
------
indices : array
indices of vertices
"""
indices = np.zeros([np.shape(surface_new.Points)[0],1])
for i in range(np.shape(surface_new.Points)[0]):
indices[i] = np.where(np.all(np.equal(surface_new.Points[i,:],surface_original.Points), axis=1))[0][0]
indices = indices.astype(int)
return indices
def calc_surface_eigenmodes(surface_input_filename, mask_input_filename, output_eval_filename, output_emode_filename, save_cut, num_modes):
"""Main function to calculate the eigenmodes of a cortical surface with application of a mask (e.g., to remove the medial wall).
Parameters
----------
surface_input_filename : str
Filename of input surface
mask_input_filename : str
Filename of mask to be applied on the surface (e.g., cortex without medial wall, values = 1 for mask and 0 elsewhere)
output_eval_filename : str
Filename of text file where the output eigenvalues will be stored
output_emode_filename : str
Filename of text file where the output eigenmodes will be stored
save_cut : boolean
Boolean to decide if the new surface with mask applied will be saved to a new surface file
num_modes : int
Number of eigenmodes to be calculated
"""
# load surface (as a brainspace object)
surface_orig = mesh.mesh_io.read_surface(surface_input_filename)
# load mask
# can be any ROI (even whole cortex)
mask_input_file_main, mask_input_file_ext = os.path.splitext(mask_input_filename)
if mask_input_file_ext == '.txt':
mask = np.loadtxt(mask_input_filename)
elif mask_input_file_ext == '.gii':
mask = nib.load(mask_input_filename).darrays[0].data
# create temporary suface based on mask
surface_cut = mesh.mesh_operations.mask_points(surface_orig, mask)
if save_cut == 1:
# old method: save vtk of surface_cut and open via lapy TriaIO
# The writing phase of this process is very slow especially for large surfaces
temp_cut_filename='temp_cut.vtk'
create_temp_surface(surface_cut, temp_cut_filename)
# load surface (as a lapy object)
tria = TriaIO.import_vtk(temp_cut_filename)
else:
# new method: replace v and t of surface_orig with v and t of surface_cut
# faster version without the need to write the vtk file
# load surface (as a lapy object)
tria = TriaIO.import_vtk(surface_input_filename)
tria.v = surface_cut.Points
tria.t = np.reshape(surface_cut.Polygons, [surface_cut.n_cells, 4])[:,1:4]
# calculate eigenvalues and eigenmodes
evals, emodes = calc_eig(tria, num_modes)
# get indices of vertices of surface_orig that match surface_cut
indices = get_indices(surface_orig, surface_cut)
# reshape emodes to match vertices of original surface
emodes_reshaped = np.zeros([surface_orig.n_points,np.shape(emodes)[1]])
for mode in range(np.shape(emodes)[1]):
emodes_reshaped[indices,mode] = np.expand_dims(emodes[:,mode], axis=1);
# save results to text files
np.savetxt(output_eval_filename, evals)
np.savetxt(output_emode_filename, emodes_reshaped)
if save_cut == 0:
if os.path.exists('temp_cut.vtk'):
os.remove('temp_cut.vtk')
def calc_surface_eigenmodes_nomask(surface_input_filename, output_eval_filename, output_emode_filename, num_modes):
"""Main function to calculate the eigenmodes of a cortical surface without application of a mask.
Parameters
----------
surface_input_filename : str
Filename of input surface
output_eval_filename : str
Filename of text file where the output eigenvalues will be stored
output_emode_filename : str
Filename of text file where the output eigenmodes will be stored
num_modes : int
Number of eigenmodes to be calculated
"""
# load surface (as a lapy object)
tria = TriaIO.import_vtk(surface_input_filename)
# calculate eigenvalues and eigenmodes
evals, emodes = calc_eig(tria, num_modes)
# save eigenmode results
np.savetxt(output_eval_filename, evals)
np.savetxt(output_emode_filename, emodes)
def main(raw_args=None):
parser = ArgumentParser(epilog="surface_eigenmodes.py -- A function to calculate the eigenmodes of a cortical surface. James Pang, Monash University, 2022 <[email protected]>")
parser.add_argument("surface_input_filename", help="An input surface in vtk format", metavar="surface_input.vtk")
parser.add_argument("output_eval_filename", help="An output text file where the eigenvalues will be stored", metavar="evals.txt")
parser.add_argument("output_emode_filename", help="An output text file where the eigenmodes will be stored", metavar="emodes.txt")
parser.add_argument("-save_cut", dest="save_cut", default=0, help="Logical value to decide whether to write the masked version of the input surface", metavar="0")
parser.add_argument("-N", dest="num_modes", default=20, help="Number of eigenmodes to be calculated, default=20", metavar="20")
parser.add_argument("-is_mask", dest="is_mask", default=1, help="Logical value to decide whether to apply the mask", metavar="1")
parser.add_argument("-mask", dest="mask_input_filename", help="An input mask text or gifti file", metavar="mask.txt")
#-------------------- Parsing the inputs from terminal: -------------------
args = parser.parse_args()
surface_input_filename = args.surface_input_filename
output_eval_filename = args.output_eval_filename
output_emode_filename = args.output_emode_filename
save_cut = int(args.save_cut)
num_modes = int(args.num_modes)
is_mask = int(args.is_mask)
mask_input_filename = args.mask_input_filename
#-------------------------------------------------------------------------------
if is_mask == 0:
calc_surface_eigenmodes_nomask(surface_input_filename, output_eval_filename, output_emode_filename, num_modes)
else:
if mask_input_filename == None:
print('ERROR: You need to provide a mask file')
else:
calc_surface_eigenmodes(surface_input_filename, mask_input_filename, output_eval_filename, output_emode_filename, save_cut, num_modes)
if __name__ == '__main__':
# running via commandline
main()
# # running within python
# surface_interest = 'fsLR_32k'
# structure = 'midthickness'
# hemispheres = ['lh', 'rh']
# num_modes = 200
# save_cut = 0
# for hemisphere in hemispheres:
# print('Processing ' + hemisphere)
# surface_input_filename = 'data/template_surfaces_volumes/' + surface_interest + '_' + structure + '-' + hemisphere + '.vtk'
# mask_input_filename = 'data/template_surfaces_volumes/' + surface_interest + '_cortex-' + hemisphere + '_mask.txt'
# # with cortex mask (remove medial wall)
# # this is the advisable way
# output_eval_filename = 'data/template_eigenmodes/' + surface_interest + '_' + structure + '-' + hemisphere + '_eval_' + str(num_modes) + '.txt'
# output_emode_filename = 'data/template_eigenmodes/' + surface_interest + '_' + structure + '-' + hemisphere + '_emode_' + str(num_modes) + '.txt'
# calc_surface_eigenmodes(surface_input_filename, mask_input_filename, output_eval_filename, output_emode_filename, save_cut, num_modes)
# # without cortex mask
# output_eval_filename = 'data/template_eigenmodes/' + 'nomask_' + surface_interest + '_' + structure + '-' + hemisphere + '_eval_' + str(num_modes) + '.txt'
# output_emode_filename = 'data/template_eigenmodes/' + 'nomask_' + surface_interest + '_' + structure + '-' + hemisphere + '_emode_' + str(num_modes) + '.txt'
# calc_surface_eigenmodes_nomask(surface_input_filename, output_eval_filename, output_emode_filename, num_modes)