Skip to content

Commit

Permalink
Updated documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
VoodooCode14 committed Aug 25, 2023
1 parent 7744cfe commit 6e199d3
Show file tree
Hide file tree
Showing 5 changed files with 27 additions and 160 deletions.
2 changes: 1 addition & 1 deletion finnpy/source_reconstruction/bem_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@ def calc_skull_and_skin_models(subject_path, subject_id, preflood_height = 25, o
"""
Employs freesufers watershed algorithm to calculate skull and skin models.
Parameters
----------
subject_path : string
Expand All @@ -34,6 +33,7 @@ def calc_skull_and_skin_models(subject_path, subject_id, preflood_height = 25, o
overwrite : boolean
Flag to overwrite if files are already present. Defaults to False.
"""

if (overwrite == True or os.path.exists(subject_path + "bem/watershed/" + subject_id + "_inner_skull_surface") == False):

cmd = ["mri_watershed", "-h", str(preflood_height), "-useSRAS", "-surf", subject_path + "bem/watershed/" + subject_id, subject_path + "mri/T1.mgz", subject_path + "bem/watershed/ws.mgz"]
Expand Down
1 change: 0 additions & 1 deletion finnpy/source_reconstruction/coregistration_meg_mri.py
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,6 @@ def _registrate_3d_points_free(src_pts, tgt_pts, weights = [1., 10., 1.], initia
Initial transformation guess,
defaults to None for no translation/rotation/scaling.
Returns
-------
est_rotors : numpy.ndarray, shape(9,)
Expand Down
30 changes: 22 additions & 8 deletions finnpy/source_reconstruction/forward_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -93,15 +93,29 @@ def _apply_coreg_to_vertices(white_vert, coreg_trans_mat):

return trans_white_vert

def _update_invalid_vertices(approx_surface, geom_white_vert, white_valid_vert):
def _update_invalid_vertices(approx_surface, white_vert, valid_vert):
"""
Updates invalid vertices
Parameters
----------
approx_surface : scipy.spatial._qhull.Delaunay
Approximate in/out skull/skin surface.
white_vert : numpy.ndarray, shape(white_vtx_cnt, 3)
White matter surface vertices.
valid_vert : numpy.ndarray, shape(white_vtx_cnt, 3)
Binary list of valid vertices.
Returns
-------
trans_white_vert : numpy.ndarray, shape(white_vtx_cnt, 3)
Transformed vertices.
"""
data = geom_white_vert[np.asarray(white_valid_vert, dtype=bool)]
data = white_vert[np.asarray(valid_vert, dtype=bool)]
inside_check = (approx_surface.find_simplex(data) != -1)
white_valid_vert[np.where(white_valid_vert)[0][~inside_check]] = False
valid_vert[np.where(valid_vert)[0][~inside_check]] = False

return white_valid_vert
return valid_vert

def _calc_magnetic_fields(white_vertices_mri, reduced_in_skull_vert, meg_to_mri_trans, pre_fwd_solution):
"""
Expand Down Expand Up @@ -441,15 +455,15 @@ def _calc_acc_hem_normals(white_vert, white_faces):

return acc_normals

def _find_vertex_clusters(white_vert, white_valid_vert, white_faces):
def _find_vertex_clusters(white_vert, valid_vert, white_faces):
"""
Identifies which input vertices (white_vert) are presented by which valid vortex.
Parameters
----------
white_vert : numpy.ndarray, shape(lh_vert_cnt, 3)
MRI model vertices.
white_valid_vert : numpy.ndarray, shape(lh_vert_cnt,)
valid_vert : numpy.ndarray, shape(lh_vert_cnt,)
Valid/supporting vertices.
white_faces : numpy.ndarray, shape(lh_face_cnt, 3)
MRI model faces.
Expand All @@ -472,7 +486,7 @@ def _find_vertex_clusters(white_vert, white_valid_vert, white_faces):
edges = edges.tocoo()
edges_dists = np.linalg.norm(white_vert[edges.row, :] - white_vert[edges.col, :], axis = 1)
edges_adjacency = scipy.sparse.csr_matrix((edges_dists, (edges.row, edges.col)), shape = edges.shape)
_, _, min_idx = scipy.sparse.csgraph.dijkstra(edges_adjacency, indices = np.where(white_valid_vert)[0], min_only = True, return_predecessors = True)
_, _, min_idx = scipy.sparse.csgraph.dijkstra(edges_adjacency, indices = np.where(valid_vert)[0], min_only = True, return_predecessors = True)

#Accumulates the clusters
sort_near_idx = np.argsort(min_idx)
Expand All @@ -484,7 +498,7 @@ def _find_vertex_clusters(white_vert, white_valid_vert, white_faces):
for cluster_idx in range(len(starts)):
cluster_grp.append(np.sort(sort_near_idx[starts[cluster_idx]:ends[cluster_idx]]))
pre_cluster_indices = sort_min_idx[breaks - 1]
cluster_indices = np.searchsorted(pre_cluster_indices, np.where(white_valid_vert)[0])
cluster_indices = np.searchsorted(pre_cluster_indices, np.where(valid_vert)[0])

return (cluster_grp, cluster_indices)

Expand Down
3 changes: 1 addition & 2 deletions finnpy/source_reconstruction/sensor_covariance.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@

def _get_bio_channel_type_idx(raw_file, mask = None):
"""
Identifies bio channel types.
Parameters
Expand Down Expand Up @@ -83,7 +82,7 @@ def _empirically_estimate_cov(cov_data, meg_ch_indices, grad_ch_indices, valid_c
grad_ch_indices : list, int
Indices of gradiometer channels.
valid_ch_indices : numpy.ndarray, shape(ch_cnt,)
Binary list identifying channels as valid/invalid.
Binary list identifying channels as valid/invalid.
Returns
-------
Expand Down
151 changes: 3 additions & 148 deletions finnpy/source_reconstruction/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,7 @@ def run_subprocess_in_custom_working_directory(subject_id, cmd):
for c in iter(lambda: process.stderr.read(1), b""):
print(c.decode(), end = "")

shutil.rmtree(path_to_tmp_cwd) # And removed alter on
shutil.rmtree(path_to_tmp_cwd) # And removed later on

def init_fs_paths(subjects_path):
"""
Expand Down Expand Up @@ -274,8 +274,8 @@ def apply_inv_transformation(data, trans):
Results
-------
data : np.ndarray, shape(n, 3)
Transformed data.
trans_data : np.ndarray, shape(n, 3)
Transformed data.
"""
inv_trans = np.linalg.inv(trans)
tmp = np.dot(inv_trans[:3, :3], data.T).T
Expand Down Expand Up @@ -396,151 +396,6 @@ def get_eigenbasis(vortex_normals, valid_vert, cluster_grp, cluster_indices,
evec = evec.reshape(-1, 3)

return evec

def calc_acc_hem_normals(white_vert, geom_white_faces):
"""
Generatoes vertex normals from accumulated face normals.
See https://en.wikipedia.org/wiki/Vertex_normal.
Parameters
----------
white_vert : numpy.ndarray, shape(vert_cnt, 3)
MRI model vertices.
geom_white_faces : numpy.ndarray, shape(face_cnt, 3)
MRI model faces.
Returns
-------
acc_normals : numpy.ndarray, shape(vert_cnt, 3)
Array of surface normals for the input faces/vertices.
"""

#Calculate surface normals
vert0 = white_vert[geom_white_faces[:, 0], :]
vert1 = white_vert[geom_white_faces[:, 1], :]
vert2 = white_vert[geom_white_faces[:, 2], :]
vortex_normals = fast_cross_product(vert1 - vert0, vert2 - vert1)
len_vortex_normals = np.linalg.norm(vortex_normals, axis = 1)
vortex_normals[len_vortex_normals > 0] /= np.expand_dims(len_vortex_normals[len_vortex_normals > 0], axis = 1)

#Accumulate face normals
acc_normals = np.zeros((white_vert.shape[0], 3))
for outer_idx in range(3):
for inner_idx in range(3):
value = np.zeros((white_vert.shape[0],))
for vertex_idx in range(geom_white_faces.shape[0]):
value[geom_white_faces[vertex_idx, outer_idx]] += vortex_normals[vertex_idx, inner_idx]
acc_normals[:, inner_idx] += value

#Normalize face normals
len_acc_normals = np.linalg.norm(acc_normals, axis = 1)
acc_normals[len_acc_normals > 0] /= np.expand_dims(len_acc_normals[len_acc_normals > 0], axis = 1)

return acc_normals

def find_vertex_clusters(white_vert, white_valid_vert, geom_white_faces):
"""
Identifies which input vertices (white_vert) are presented by which valid vortex.
Parameters
----------
white_vert : numpy.ndarray, shape(lh_vert_cnt, 3)
MRI model vertices.
white_valid_vert : numpy.ndarray, shape(lh_vert_cnt,)
Valid/supporting vertices.
geom_white_faces : numpy.ndarray, shape(lh_face_cnt, 3)
MRI model faces.
Returns
-------
cluster_grp : list, len(n,)
Transformed surface model.
cluster_indices : list, len(n,)
Transformed surface model.
"""

#Get adjacency matrix to identify nearest valid vortex
edges = scipy.sparse.coo_matrix((np.ones(3 * geom_white_faces.shape[0]),
(np.concatenate((geom_white_faces[:, 0], geom_white_faces[:, 1], geom_white_faces[:, 2])),
np.concatenate((geom_white_faces[:, 1], geom_white_faces[:, 2], geom_white_faces[:, 0])))),
shape = (white_vert.shape[0], white_vert.shape[0]))
edges = edges.tocsr()
edges += edges.T
edges = edges.tocoo()
edges_dists = np.linalg.norm(white_vert[edges.row, :] - white_vert[edges.col, :], axis = 1)
edges_adjacency = scipy.sparse.csr_matrix((edges_dists, (edges.row, edges.col)), shape = edges.shape)
_, _, min_idx = scipy.sparse.csgraph.dijkstra(edges_adjacency, indices = np.where(white_valid_vert)[0], min_only = True, return_predecessors = True)

#Accumulates the clusters
sort_near_idx = np.argsort(min_idx)
sort_min_idx = min_idx[sort_near_idx]
breaks = np.where(sort_min_idx[1:] != sort_min_idx[:-1])[0] + 1
starts = [0] + breaks.tolist()
ends = breaks.tolist() + [len(min_idx)]
cluster_grp = list()
for cluster_idx in range(len(starts)):
cluster_grp.append(np.sort(sort_near_idx[starts[cluster_idx]:ends[cluster_idx]]))
pre_cluster_indices = sort_min_idx[breaks - 1]
cluster_indices = np.searchsorted(pre_cluster_indices, np.where(white_valid_vert)[0])

return (cluster_grp, cluster_indices)

def cortical_surface_reorient_fwd_model(lh_white_vert, lh_geom_white_faces, lh_valid_vert,
rh_white_vert, rh_geom_white_faces, rh_valid_vert,
fwd_sol, mri_to_head_trans):
"""
Transforms a fwd model into surface orientation (orthogonal to the respective surface cluster;
allowing for cortically constrained inverse modeling)
and shrinks it by making closeby channels project to the same destination.
This is different from a 'default' 3D transformation.
Parameters
----------
lh_white_vert : numpy.ndarray, shape(lh_vert_cnt, 3)
Lh vertices.
lh_geom_white_faces : numpy.ndarray, shape(lh_face_cnt, 3)
Lh faces.
lh_valid_vert : numpy.ndarray, shape(lh_vert_cnt,)
Valid/supporting lh vertices.
rh_white_vert : numpy.ndarray, shape(rh_vert_cnt, 3)
Rh vertices.
rh_geom_white_faces : numpy.ndarray, shape(rh_face_cnt, 3)
Rh faces.
rh_valid_vert : numpy.ndarray, shape(rh_vert_cnt,)
Valid/supporting rh vertices.
fwd_sol : numpy.ndarray, shape(meg_ch_cnt, valid_vtx_cnt * 3)
Forward solution with default orientation.
mri_to_head_trans : numpy.ndarray, shape(4, 4)
MRI to head transformation matrix.
Returns
-------
fwd_sol * rot : numpy.ndarray, shape(meg_ch_cnt, valid_vtx_cnt)
Transformed surface model.
"""

#Computes vertex normals
lh_acc_normals = calc_acc_hem_normals(lh_white_vert, lh_geom_white_faces)
rh_acc_normals = calc_acc_hem_normals(rh_white_vert, rh_geom_white_faces)

# Figures out which real vertices are represented by the same model vortex
(lh_cluster_grp, lh_cluster_indices) = find_vertex_clusters(lh_white_vert, lh_valid_vert, lh_geom_white_faces)
(rh_cluster_grp, rh_cluster_indices) = find_vertex_clusters(rh_white_vert, rh_valid_vert, rh_geom_white_faces)

#Determines the orientation
lh_orient = get_eigenbasis(lh_acc_normals, lh_valid_vert, lh_cluster_grp, lh_cluster_indices, mri_to_head_trans)
rh_orient = get_eigenbasis(rh_acc_normals, rh_valid_vert, rh_cluster_grp, rh_cluster_indices, mri_to_head_trans)

#Concatenates the eigenvector bases
orient = np.concatenate((lh_orient, rh_orient), axis = 0)

#Transforms eigenvector matrix into block matrix for ease of use
rot = orient_mat_to_block_format(orient[2::3, :])

#Applies rotation matrix to the fwd solution
return fwd_sol * rot

def calc_mri_maps(subj_path, fs_avg_path, hemisphere, overwrite):
"""
Expand Down

0 comments on commit 6e199d3

Please sign in to comment.