1- def fit_spheres_to_mri (subjects_dir , subject , bem_surf , trans , n_spheres ,show_spheres = False ):
1+ def fit_spheres_to_mri (
2+ subjects_dir , subject , bem_surf , trans , n_spheres , show_spheres = False
3+ ):
24 """Fits two spheres to MRI using BEM, such that spheres fit while brain but
35 do not encroach on sensors. For use with Milti-SSS Maxwell Filtering
46
@@ -33,8 +35,9 @@ def fit_spheres_to_mri(subjects_dir, subject, bem_surf, trans, n_spheres,show_sp
3335 """
3436 ## --- required imports
3537
36- import nibabel as nib
3738 import os
39+
40+ import nibabel as nib
3841 import numpy as np
3942 import vedo
4043 from scipy .spatial import KDTree
@@ -49,58 +52,62 @@ def fit_spheres_to_mri(subjects_dir, subject, bem_surf, trans, n_spheres,show_sp
4952
5053 ## --- begin
5154 mindist = 2e-3
52- assert bem_surf [0 ]['id' ] == FIFF .FIFFV_BEM_SURF_ID_HEAD
53- assert bem_surf [2 ]['id' ] == FIFF .FIFFV_BEM_SURF_ID_BRAIN
55+ assert bem_surf [0 ]["id" ] == FIFF .FIFFV_BEM_SURF_ID_HEAD
56+ assert bem_surf [2 ]["id" ] == FIFF .FIFFV_BEM_SURF_ID_BRAIN
5457 scalp , _ , inner_skull = bem_surf
55- inside_scalp = _CheckInside (scalp , mode = ' pyvista' )
56- inside_skull = _CheckInside (inner_skull , mode = ' pyvista' )
57- m3_to_cc = 100 ** 3
58- assert inside_scalp (inner_skull ['rr' ]).all ()
59- assert not inside_skull (scalp ['rr' ]).any ()
60- b = vedo .Mesh ([inner_skull ['rr' ], inner_skull [' tris' ]])
61- s = vedo .Mesh ([scalp ['rr' ], scalp [' tris' ]])
62- s_tree = KDTree (scalp ['rr' ])
58+ inside_scalp = _CheckInside (scalp , mode = " pyvista" )
59+ inside_skull = _CheckInside (inner_skull , mode = " pyvista" )
60+ m3_to_cc = 100 ** 3
61+ assert inside_scalp (inner_skull ["rr" ]).all ()
62+ assert not inside_skull (scalp ["rr" ]).any ()
63+ b = vedo .Mesh ([inner_skull ["rr" ], inner_skull [" tris" ]])
64+ s = vedo .Mesh ([scalp ["rr" ], scalp [" tris" ]])
65+ s_tree = KDTree (scalp ["rr" ])
6366 brain_volume = b .volume ()
64- print (f' Brain vedo: { brain_volume * m3_to_cc :8.2f} cc' )
65- brain_vol = nib .load (os .path .join (subjects_dir ,subject ,' mri' , ' brainmask.mgz' ))
67+ print (f" Brain vedo: { brain_volume * m3_to_cc :8.2f} cc" )
68+ brain_vol = nib .load (os .path .join (subjects_dir , subject , " mri" , " brainmask.mgz" ))
6669 brain_rr = np .array (np .where (brain_vol .get_fdata ())).T
67- brain_rr = apply_trans (brain_vol .header .get_vox2ras_tkr (), brain_rr ) / 1000. #apply a transformation matrix
68- del brain_vol #delete brain volume
70+ brain_rr = (
71+ apply_trans (brain_vol .header .get_vox2ras_tkr (), brain_rr ) / 1000.0
72+ ) # apply a transformation matrix
73+ del brain_vol # delete brain volume
6974 brain_rr = brain_rr [inside_skull (brain_rr )]
7075 vox_to_m3 = 1e-9
7176 brain_volume_vox = len (brain_rr ) * vox_to_m3
7277
7378 def _print_q (title , got , want ):
74- title = f' { title } :' .ljust (15 )
75- print (f' { title } { got * m3_to_cc :8.2f} cc ({ (want - got ) / want * 100 :6.2f} %)' )
79+ title = f" { title } :" .ljust (15 )
80+ print (f" { title } { got * m3_to_cc :8.2f} cc ({ (want - got ) / want * 100 :6.2f} %)" )
7681
77- _print_q (' Brain vox' , brain_volume_vox , brain_volume_vox )
82+ _print_q (" Brain vox" , brain_volume_vox , brain_volume_vox )
7883
7984 # 1. Compute a naive sphere using the center of mass of brain surf verts
80- naive_c = np .mean (inner_skull ['rr' ], axis = 0 )
81- naive_r = np .min (np .linalg .norm (inner_skull ['rr' ] - naive_c , axis = 1 ))
82- naive_v = 4 / 3 * np .pi * naive_r ** 3
83- _print_q (' Naive sphere' , naive_v , brain_volume )
85+ naive_c = np .mean (inner_skull ["rr" ], axis = 0 )
86+ naive_r = np .min (np .linalg .norm (inner_skull ["rr" ] - naive_c , axis = 1 ))
87+ naive_v = 4 / 3 * np .pi * naive_r ** 3
88+ _print_q (" Naive sphere" , naive_v , brain_volume )
8489 s1 = vedo .Sphere (naive_c , naive_r , res = 100 )
85- _print_q (' Naive vedo' , s1 .volume (), brain_volume )
90+ _print_q (" Naive vedo" , s1 .volume (), brain_volume )
8691
8792 # 2. Now use the larger radius (to head) plus mesh arithmetic
8893 better_r = s_tree .query (naive_c )[0 ] - mindist
8994 s1 = vedo .Sphere (naive_c , better_r , res = 24 )
90- _print_q (' Better vedo' , s1 .boolean ("intersect" , b ).volume (), brain_volume )
95+ _print_q (" Better vedo" , s1 .boolean ("intersect" , b ).volume (), brain_volume )
9196 v = np .sum (np .linalg .norm (brain_rr - naive_c , axis = 1 ) <= better_r ) * vox_to_m3
92- _print_q (' Better vox' , v , brain_volume_vox )
97+ _print_q (" Better vox" , v , brain_volume_vox )
9398
9499 # 3. Now optimize one sphere
95- from scipy .optimize import fmin_cobyla #constrained optimization by linear approximation
100+ from scipy .optimize import (
101+ fmin_cobyla , # constrained optimization by linear approximation
102+ )
96103
97104 def _cost (c ):
98105 cs = c .reshape (- 1 , 3 )
99- rs = np .maximum (s_tree .query (cs )[0 ] - mindist , 0. )
106+ rs = np .maximum (s_tree .query (cs )[0 ] - mindist , 0.0 )
100107 resid = brain_volume
101108 mask = None
102109 for c , r in zip (cs , rs ):
103- if not (r and s .contains (c )): # was is_inside
110+ if not (r and s .contains (c )): # was is_inside
104111 continue
105112 m = np .linalg .norm (brain_rr - c , axis = 1 ) <= r
106113 if mask is None :
@@ -114,67 +121,72 @@ def _cost(c):
114121
115122 def _cons (c ):
116123 cs = c .reshape (- 1 , 3 )
117- sign = np .array ([2 * s .contains (c ) - 1 for c in cs ], float ) # was "is_inside"
124+ sign = np .array ([2 * s .contains (c ) - 1 for c in cs ], float ) # was "is_inside"
118125 cons = sign * s_tree .query (cs )[0 ] - mindist
119126 return cons
120127
121128 x = naive_c
122129 c_opt_1 = fmin_cobyla (_cost , x , _cons , rhobeg = 1e-2 , rhoend = 1e-4 )
123130 v_opt_1 = brain_volume_vox - _cost (c_opt_1 )
124- _print_q (' COBYLA 1' , v_opt_1 , brain_volume_vox )
131+ _print_q (" COBYLA 1" , v_opt_1 , brain_volume_vox )
125132
126133 # 4. Now optimize two spheres
127134 x = np .concatenate ([c_opt_1 , naive_c ])
128135 c_opt_2 = fmin_cobyla (_cost , x , _cons , rhobeg = 1e-2 , rhoend = 1e-4 )
129136 v_opt_2 = brain_volume_vox - _cost (c_opt_2 )
130- _print_q (' COBYLA 2' , v_opt_2 , brain_volume_vox )
137+ _print_q (" COBYLA 2" , v_opt_2 , brain_volume_vox )
131138
132139 # 4. Finally, three spheres (not perfect, not global opt)
133140 x = np .concatenate ([c_opt_2 , naive_c ])
134141 c_opt_3 = fmin_cobyla (_cost , x , _cons , rhobeg = 1e-2 , rhoend = 1e-4 )
135142 v_opt_3 = brain_volume_vox - _cost (c_opt_3 )
136- _print_q (' COBYLA 3' , v_opt_3 , brain_volume_vox )
143+ _print_q (" COBYLA 3" , v_opt_3 , brain_volume_vox )
137144
138145 if show_spheres :
146+ import matplotlib
139147 import pyvista as pv
140148 import pyvistaqt
141- import matplotlib
149+
142150 plotter = pyvistaqt .BackgroundPlotter (
143- shape = (1 , 2 ), window_size = (1200 , 300 ),
144- editor = False , menu_bar = False , toolbar = False )
145- plotter .background_color = 'w'
146- brain_mesh = pv .make_tri_mesh (inner_skull ['rr' ], inner_skull ['tris' ])
147- scalp_mesh = pv .make_tri_mesh (scalp ['rr' ], scalp ['tris' ])
148- colors = matplotlib .rcParams ['axes.prop_cycle' ].by_key ()['color' ]
151+ shape = (1 , 2 ),
152+ window_size = (1200 , 300 ),
153+ editor = False ,
154+ menu_bar = False ,
155+ toolbar = False ,
156+ )
157+ plotter .background_color = "w"
158+ brain_mesh = pv .make_tri_mesh (inner_skull ["rr" ], inner_skull ["tris" ])
159+ scalp_mesh = pv .make_tri_mesh (scalp ["rr" ], scalp ["tris" ])
160+ colors = matplotlib .rcParams ["axes.prop_cycle" ].by_key ()["color" ]
149161 mesh_kwargs = dict (render = False , reset_camera = False , smooth_shading = True )
150162 for ci , cs in enumerate ((c_opt_1 , c_opt_2 , c_opt_3 )):
151163 plotter .subplot (0 , ci )
152- plotter .camera .position = (0. , - 0.5 , 0 )
153- plotter .camera .focal_point = (0. , 0. , 0. )
164+ plotter .camera .position = (0.0 , - 0.5 , 0 )
165+ plotter .camera .focal_point = (0.0 , 0.0 , 0.0 )
154166 plotter .camera .azimuth = 90
155167 plotter .camera .elevation = 0
156- plotter .camera .up = (0. , 0. , 1. )
157- plotter .add_mesh (brain_mesh , opacity = 0.2 , color = 'k' , ** mesh_kwargs )
158- plotter .add_mesh (scalp_mesh , opacity = 0.1 , color = ' tan' , ** mesh_kwargs )
168+ plotter .camera .up = (0.0 , 0.0 , 1.0 )
169+ plotter .add_mesh (brain_mesh , opacity = 0.2 , color = "k" , ** mesh_kwargs )
170+ plotter .add_mesh (scalp_mesh , opacity = 0.1 , color = " tan" , ** mesh_kwargs )
159171 for c , color in zip (cs .reshape (- 1 , 3 ), colors ):
160172 sphere = pv .Sphere (s_tree .query (c )[0 ] - mindist , c )
161173 plotter .add_mesh (sphere , opacity = 0.5 , color = color , ** mesh_kwargs )
162174 plotter .show ()
163175
164176 # Ready centers to output, transform into device space
165- mri_head_t = invert_transform (read_trans (trans ))
177+ mri_head_t = invert_transform (read_trans (trans ))
166178 if mri_head_t ["from" ] == FIFF .FIFFV_COORD_HEAD :
167179 mri_head_t = invert_transform (mri_head_t )
168- assert mri_head_t [' from' ] == FIFF .FIFFV_COORD_MRI , mri_head_t [' from' ]
169- centers = []
170- for use in (c_opt_1 ,c_opt_2 ,c_opt_3 ):
180+ assert mri_head_t [" from" ] == FIFF .FIFFV_COORD_MRI , mri_head_t [" from" ]
181+ centers = []
182+ for use in (c_opt_1 , c_opt_2 , c_opt_3 ):
171183 centers .append (apply_trans (mri_head_t , use .reshape (- 1 , 3 )))
172- if n_spheres == 1 :
184+ if n_spheres == 1 :
173185 return centers [0 ]
174- if n_spheres == 2 :
186+ if n_spheres == 2 :
175187 return centers [1 ]
176- if n_spheres == 3 :
177- print ("Warning: use of mSSS with three origins and expansions is not tested or recommended" )
188+ if n_spheres == 3 :
189+ print (
190+ "Warning: use of mSSS with three origins and expansions is not tested or recommended"
191+ )
178192 return centers [2 ]
179-
180-
0 commit comments