Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def check_dist_to_pial_vertices(elc_name, subject_fol, threshold):
(electrodes, names, hemis, threshold) = utils.load(op.join(output_fol, '{}_electrodes.pkl'.format(int(threshold))))
hemi_indices, close_verts_indices, all_dists, dural_mask = fect.get_t1_voxels_inside_dural(electrodes, subject_fol)
elc_ind = names.index(elc_name)
t1_tkras_coords = np.array([electrodes[elc_ind]])
verts, faces, normals = {}, {}, {}
for hemi in ['lh', 'rh']:
verts[hemi], faces[hemi] = nib.freesurfer.read_geometry(op.join(subject_fol, 'surf', '{}.dural'.format(hemi)))
normals[hemi] = fect.calc_normals(verts[hemi], faces[hemi])
hemi = hemis[elc_ind]
dists = cdist(t1_tkras_coords, verts[hemi])
close_verts = np.argmin(dists, axis=1)
is_inside = fect.point_in_mesh(t1_tkras_coords[0], verts[hemi][close_verts[0]], normals[hemi][close_verts[0]])
# vert_norm = np.linalg.norm(vertices[close_verts][0])
# elc_norm = np.linalg.norm(t1_tkras_coords[0])
print(is_inside)
def load_surf(base_path, name):
return nibabel.freesurfer.read_geometry(
os.path.join(base_path, 'surf/' + name))
def compute_adjacency(hemi, min_dist, max_dist, projfrac,step_dist):
v, f = nib.freesurfer.read_geometry("%s/fsaverage/surf/%s.white" % ((os.environ["SUBJECTS_DIR"]),hemi))
v = v.astype(np.float32, order = "C")
f = f.astype(np.int32, order = "C")
v, f = mergeIdenticalVertices(v, f)
v, f = removeNonManifoldTriangles(v, f)
vn = computeNormals(v, f)
t = nib.freesurfer.read_morph_data("%s/fsaverage/surf/%s.thickness" % ((os.environ["SUBJECTS_DIR"]),hemi))
v_ = projNormFracThick(v, vn, t, projfrac) # project to midthickness
nib.freesurfer.io.write_geometry("%s.midthickness" % hemi, v_, f)
thresholds = np.arange(min_dist, max_dist, step=step_dist, dtype = np.float32)
adjacency = compute(v_, f, thresholds)
mean_full_lh = img_mean_lh.get_data()
mean_lh = np.squeeze(mean_full_lh)
img_mean_rh = nib.freesurfer.mghformat.load("rh.mean.%s.%s.mgh" % (surface,FWHM))
mean_full_rh = img_mean_rh.get_data()
mean_rh = np.squeeze(mean_full_rh)
#TFCE
if opts.triangularmesh:
print("Creating adjacency set")
if opts.inputsurfs:
# 3 Neighbour vertex connectity
v_lh, faces_lh = nib.freesurfer.read_geometry(opts.inputsurfs[0])
v_rh, faces_rh = nib.freesurfer.read_geometry(opts.inputsurfs[1])
else:
v_lh, faces_lh = nib.freesurfer.read_geometry("%s/fsaverage/surf/lh.sphere" % os.environ["SUBJECTS_DIR"])
v_rh, faces_rh = nib.freesurfer.read_geometry("%s/fsaverage/surf/rh.sphere" % os.environ["SUBJECTS_DIR"])
adjac_lh = create_adjac_vertex(v_lh,faces_lh)
adjac_rh = create_adjac_vertex(v_rh,faces_rh)
elif opts.adjfiles:
print("Loading prior adjacency set")
arg_adjac_lh = opts.adjfiles[0]
arg_adjac_rh = opts.adjfiles[1]
adjac_lh = np.load(arg_adjac_lh)
adjac_rh = np.load(arg_adjac_rh)
elif opts.dist:
print("Loading prior adjacency set for %s mm" % opts.dist[0])
adjac_lh = np.load("%s/adjacency_sets/lh_adjacency_dist_%s.0_mm.npy" % (scriptwd,str(opts.dist[0])))
adjac_rh = np.load("%s/adjacency_sets/rh_adjacency_dist_%s.0_mm.npy" % (scriptwd,str(opts.dist[0])))
else:
print("Error")
if opts.noweight or opts.triangularmesh:
vdensity_lh = 1
def convert_fs_to_brain_visa(fs_surf):
v, f = nibabel.freesurfer.read_geometry(fs_surf)
write_brain_visa_surf(fs_surf + '.tri', v, f)
else:
bin_mask_lh = mean_lh != 0
bin_mask_rh = mean_rh != 0
data_lh = data_lh[bin_mask_lh]
num_vertex_lh = data_lh.shape[0]
data_rh = data_rh[bin_mask_rh]
num_vertex_rh = data_rh.shape[0]
num_vertex = num_vertex_lh + num_vertex_rh
all_vertex = data_full_lh.shape[0]
#TFCE
if opts.triangularmesh:
print("Creating adjacency set")
if opts.inputsurfs:
# 3 Neighbour vertex connectity
v_lh, faces_lh = nib.freesurfer.read_geometry(opts.inputsurfs[0])
v_rh, faces_rh = nib.freesurfer.read_geometry(opts.inputsurfs[1])
else:
v_lh, faces_lh = nib.freesurfer.read_geometry("%s/fsaverage/surf/lh.sphere" % os.environ["SUBJECTS_DIR"])
v_rh, faces_rh = nib.freesurfer.read_geometry("%s/fsaverage/surf/rh.sphere" % os.environ["SUBJECTS_DIR"])
adjac_lh = create_adjac_vertex(v_lh,faces_lh)
adjac_rh = create_adjac_vertex(v_rh,faces_rh)
elif opts.adjfiles:
print("Loading prior adjacency set")
arg_adjac_lh = opts.adjfiles[0]
arg_adjac_rh = opts.adjfiles[1]
adjac_lh = np.load(arg_adjac_lh)
adjac_rh = np.load(arg_adjac_rh)
elif opts.dist:
print("Loading prior adjacency set for %s mm" % opts.dist[0])
adjac_lh = np.load("%s/adjacency_sets/lh_adjacency_dist_%s.0_mm.npy" % (scriptwd,str(opts.dist[0])))
adjac_rh = np.load("%s/adjacency_sets/rh_adjacency_dist_%s.0_mm.npy" % (scriptwd,str(opts.dist[0])))
if relative:
denominator = np.sum(abs(couple_signal))
euclidean = euclidean / denominator if denominator else 0
sum_tmp += euclidean
count += 1
return sum_tmp / float(count) if count else 0
if __name__ == '__main__':
from nibabel.freesurfer import read_geometry
from froi.io.surf_io import read_scalar_data
from networkx import Graph
from graph_tool import graph2parcel, node_attr2array
import nibabel as nib
coords, faces = read_geometry('/nfs/t1/nsppara/corticalsurface/fsaverage5/surf/rh.inflated')
scalar = read_scalar_data('/nfs/t3/workingshop/chenxiayu/data/region-growing-froi/S1/surf/'
'rh_zstat1_1w_fracavg.mgz')
# faces = np.array([[1, 2, 3], [0, 1, 3]])
# scalar = np.array([[1], [2], [3], [4]])
graph = mesh2graph(faces, vtx_signal=scalar, weight_normalization=True)
graph, parcel_neighbors = graph2parcel(graph, n=5000)
labels = [attrs['label'] for attrs in graph.node.values()]
print 'finish ncut!'
labels = np.unique(labels)
print len(labels)
print np.max(labels)
arr = node_attr2array(graph, ('label',))
# zero_idx = np.where(map(lambda x: x not in parcel_neighbors[800], arr))
def read_brain_mesh(surface_path):
u"""Read triangular format Freesurfer brain surface.
Parameters
----------
surface_path: str
Path to the brain surface.
Returns
-------
vertices : numpy.array
Array of vertex (x, y, z) coordinates, of size number_of_vertices x 3.
faces : numpy.array
Array defining mesh triangles, of size number_of_faces x 3.
"""
vertices, faces = freesurfer.read_geometry(surface_path)
faces = faces.astype(np.uint32)
return vertices, faces
def load_geometry(self):
"""Load geometry of the surface.
Parameters
----------
None
Returns
-------
None
"""
surf_path = path.join(self.data_path, 'surf',
'%s.%s' % (self.hemi, self.surf))
coords, faces = freesurfer.read_geometry(surf_path)
if self.units == 'm':
coords /= 1000.
if self.offset is not None:
if self.hemi == 'lh':
coords[:, 0] -= (np.max(coords[:, 0]) + self.offset)
else:
coords[:, 0] -= (np.min(coords[:, 0]) + self.offset)
nn = _compute_normals(coords, faces)
if self.coords is None:
self.coords = coords
self.faces = faces
self.nn = nn
else:
self.coords[:] = coords
self.faces[:] = faces
all_vertex = data_full_lh.shape[0]
else:
data_lh = np.squeeze(nib.freesurfer.mghformat.load("%s/lh.all.%s.%s.mgh" % (sfolder, surface,FWHM)).get_data())
data_rh = np.squeeze(nib.freesurfer.mghformat.load("%s/rh.all.%s.%s.mgh" % (sfolder, surface,FWHM)).get_data())
data.append(np.hstack((data_lh[mask_lh].T,data_rh[mask_rh].T)))
data_lh = data_rh = []
data = np.array(data)
#TFCE
if opts.triangularmesh:
print("Creating adjacency set")
if opts.inputsurfs:
# 3 Neighbour vertex connectity
v_lh, faces_lh = nib.freesurfer.read_geometry(opts.inputsurfs[0])
v_rh, faces_rh = nib.freesurfer.read_geometry(opts.inputsurfs[1])
else:
v_lh, faces_lh = nib.freesurfer.read_geometry("%s/fsaverage/surf/lh.sphere" % os.environ["SUBJECTS_DIR"])
v_rh, faces_rh = nib.freesurfer.read_geometry("%s/fsaverage/surf/rh.sphere" % os.environ["SUBJECTS_DIR"])
adjac_lh = create_adjac_vertex(v_lh,faces_lh)
adjac_rh = create_adjac_vertex(v_rh,faces_rh)
elif opts.adjfiles:
print("Loading prior adjacency set")
arg_adjac_lh = opts.adjfiles[0]
arg_adjac_rh = opts.adjfiles[1]
adjac_lh = np.load(arg_adjac_lh)
adjac_rh = np.load(arg_adjac_rh)
elif opts.dist:
print("Loading prior adjacency set for %s mm" % opts.dist[0])
adjac_lh = np.load("%s/adjacency_sets/lh_adjacency_dist_%s.0_mm.npy" % (scriptwd,str(opts.dist[0])))
adjac_rh = np.load("%s/adjacency_sets/rh_adjacency_dist_%s.0_mm.npy" % (scriptwd,str(opts.dist[0])))
else: