diff --git a/bin/aseg2srf b/bin/aseg2srf new file mode 100755 index 0000000..4263dbb --- /dev/null +++ b/bin/aseg2srf @@ -0,0 +1,119 @@ +#!/bin/bash + +# Print usage if no argument is given +if [ -z "$1" ]; then +cat < [-l ] [-d] + +Options: +-s : Specify a list of subjects between quotes, + e.g. -s "john bill mary mark" or a text file + containing one subject per line. +-l : Specify a list of labels between quotes, + e.g. "10 11 14 52 253", or a text file + containing one label per line, or ignore + this option to convert all labels. +-d Debug mode. Leave all temporary files. + +Requirements: +FreeSurfer must have been configured and the variables +FREESURFER_HOME and SUBJECTS_DIR must have been correctly set. + +_____________________________________ +Anderson M. Winkler +Institute of Living / Yale University +Jul/2009 +http://brainder.org +EOU +exit +fi + +# List of labels to be converted if no list is specified +LABLIST="4 5 7 8 10 11 12 13 14 15 16 17 18 26 28 43 44 46 47 49 50 51 52 53 54 58 60 251 252 253 254 255" + +# Check and accept arguments +SBJLIST="" +DEBUG=N +while getopts 's:l:d' OPTION +do + case ${OPTION} in + s) SBJLIST=$( [[ -f ${OPTARG} ]] && cat ${OPTARG} || echo "${OPTARG}" ) ;; + l) LABLIST=$( [[ -f ${OPTARG} ]] && cat ${OPTARG} || echo "${OPTARG}" ) ;; + d) DEBUG=Y ;; + esac +done + +# Prepare a random string to save temporary files +#if hash md5 2> /dev/null ; then +# RND0=$(head -n 1 /dev/random | md5) +#elif hash md5sum 2> /dev/null ; then +# RND0=$(head -n 1 /dev/random | md5sum) +#fi +#RNDSTR=${RND0:0:12} +RNDSTR=${RANDOM:0:12} + +# Define a function for Ctrl+C as soon as the RNDSTR is defined +trap bashtrap INT +bashtrap() +{ + break ; break + [[ "${s}" != "" ]] && rm -rf ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR} # try to delete temp files + exit 1 +} + +# For each subject +for s in ${SBJLIST} ; do + + # Create directories for temp files and results + mkdir -p ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR} + mkdir -p ${SUBJECTS_DIR}/${s}/ascii + + # For each label + for lab in ${LABLIST} ; do + + # Label string + lab0=$(printf %03d ${lab}) + + # Pre-tessellate + echo "==> Pre-tessellating: ${s}, ${lab0}" + ${FREESURFER_HOME}/bin/mri_pretess \ + ${SUBJECTS_DIR}/${s}/mri/aseg.mgz ${lab} \ + ${SUBJECTS_DIR}/${s}/mri/norm.mgz \ + ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}_filled.mgz + + # Tessellate + echo "==> Tessellating: ${s}, ${lab0}" + ${FREESURFER_HOME}/bin/mri_tessellate \ + ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}_filled.mgz \ + ${lab} ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}_notsmooth + + # Smooth + echo "==> Smoothing: ${s}, ${lab0}" + ${FREESURFER_HOME}/bin/mris_smooth -nw \ + ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}_notsmooth \ + ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0} + + # Convert to ASCII + echo "==> Converting to ASCII: ${s}, ${lab0}" + ${FREESURFER_HOME}/bin/mris_convert \ + ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0} \ + ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}.asc + mv ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}/aseg_${lab0}.asc \ + ${SUBJECTS_DIR}/${s}/ascii/aseg_${lab0}.srf + done + + # Get rid of temp files + if [ "${DEBUG}" == "Y" ] ; then + echo "==> Temporary files for ${s} saved at:" + echo "${SUBJECTS_DIR}/${s}/tmp/${RNDSTR}" + else + echo "==> Removing temporary files for ${s}" + rm -rf ${SUBJECTS_DIR}/${s}/tmp/${RNDSTR} + fi + echo "==> Done: ${s}" +done +exit 0 diff --git a/bin/main.sh b/bin/main.sh index 4a8eef3..fe2bb19 100644 --- a/bin/main.sh +++ b/bin/main.sh @@ -13,6 +13,9 @@ sh $CODE/reconALL.sh source $CODE/postreconALL.sh #<-It returns all surfaces and volumes... +# Generate surfaces for subcortical structures +$CODE/aseg2srf -s ${SUBJECT} + #Now, for Steps 2 we split to three parallel streams: diff --git a/bin/run_surfaces_to_struct_dataset.py b/bin/run_surfaces_to_struct_dataset.py new file mode 100644 index 0000000..b2bc15e --- /dev/null +++ b/bin/run_surfaces_to_struct_dataset.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python3 + +import logging +import os +import sys + +from tvb.recon.flow.surfaces_to_structural_datasets import SurfacesToStructuralDataset +from tvb.recon.cli.runner import SimpleRunner + +def main(): + subjects_dir, subjid, source_lut, target_lut, weights_file, tract_lengths_file, out_file, out_surf = sys.argv[1:] + + logging.basicConfig(level=logging.INFO) + runner = SimpleRunner() + + flow = SurfacesToStructuralDataset( + os.path.join(subjects_dir, subjid, "surf"), + os.path.join(subjects_dir, subjid, "label"), + os.path.join(subjects_dir, subjid, "ascii"), + source_lut, + target_lut, + weights_file, + tract_lengths_file, + out_file, + out_surf + ) + + flow.run(runner) + + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/data/freesurfer_fsaverage/label/rh.aparc.annot b/data/freesurfer_fsaverage/label/rh.aparc.annot new file mode 100755 index 0000000..dc2e8aa Binary files /dev/null and b/data/freesurfer_fsaverage/label/rh.aparc.annot differ diff --git a/tvb/recon/algo/structural_dataset.py b/tvb/recon/algo/structural_dataset.py new file mode 100644 index 0000000..0dbdb21 --- /dev/null +++ b/tvb/recon/algo/structural_dataset.py @@ -0,0 +1,70 @@ + +import numpy as np +import os +import tempfile +from typing import List +from zipfile import ZipFile + + +class StructuralDataset: + def __init__(self, + orientations: np.ndarray, + areas: np.ndarray, + centers: np.ndarray, + cortical: np.ndarray, + weights_ut: np.ndarray, + tract_lengths_ut: np.ndarray, + names: List[str]): + + nregions = len(names) + + assert orientations.shape == (nregions, 3) + assert areas.shape == (nregions,) + assert centers.shape == (nregions, 3) + assert cortical.shape == (nregions,) + assert weights_ut.shape == (nregions, nregions) + assert tract_lengths_ut.shape == (nregions, nregions) + + # Upper triangular -> symmetric matrices + assert np.sum(np.tril(weights_ut, -1)) == 0 + assert np.sum(np.tril(tract_lengths_ut, -1)) == 0 + self.weights = weights_ut + weights_ut.transpose() - np.diag(np.diag(weights_ut)) + self.tract_lengths = tract_lengths_ut + tract_lengths_ut.transpose() - np.diag(np.diag(tract_lengths_ut)) + + self.orientations = orientations + self.areas = areas + self.centers = centers + self.cortical = cortical + self.names = names + + def save_to_txt_zip(self, filename: os.PathLike): + + tmpdir = tempfile.TemporaryDirectory() + + file_areas = os.path.join(tmpdir.name, 'areas.txt') + file_orientations = os.path.join(tmpdir.name, 'average_orientations.txt') + file_centres = os.path.join(tmpdir.name, 'centres.txt') + file_cortical = os.path.join(tmpdir.name, 'cortical.txt') + file_weights = os.path.join(tmpdir.name, 'weights.txt') + file_tract_lengths = os.path.join(tmpdir.name, 'tract_lengths.txt') + + np.savetxt(file_areas, self.areas, fmt='%.2f') + np.savetxt(file_orientations, self.orientations, fmt='%.2f %.2f %.2f') + np.savetxt(file_cortical, self.cortical, fmt='%d') + np.savetxt(file_weights, self.weights, fmt='%d') + np.savetxt(file_tract_lengths, self.tract_lengths, fmt='%.3f') + + with open(file_centres, 'w') as f: + for i, name in enumerate(self.names): + f.write('%s %.4f %.4f %.4f\n' % (name, self.centers[i, 0], self.centers[i, 1], self.centers[i, 2])) + + with ZipFile(filename, 'w') as zip_file: + zip_file.write(file_areas, os.path.basename(file_areas)) + zip_file.write(file_orientations, os.path.basename(file_orientations)) + zip_file.write(file_centres, os.path.basename(file_centres)) + zip_file.write(file_cortical, os.path.basename(file_cortical)) + zip_file.write(file_weights, os.path.basename(file_weights)) + zip_file.write(file_tract_lengths, os.path.basename(file_tract_lengths)) + + def save_to_h5(self, filename): + pass diff --git a/tvb/recon/cli/core.py b/tvb/recon/cli/core.py index 0267dbd..991ebef 100644 --- a/tvb/recon/cli/core.py +++ b/tvb/recon/cli/core.py @@ -31,8 +31,7 @@ def __new__(cls: type, name: str, bases: Tuple[type], attrs: Dict[str, object]) -> type: for key, value in list(attrs.items()): if not key.startswith('_'): - value: - str + value: str attrs[key] = Flag(key, value) return super().__new__(cls, name, bases, attrs) @@ -55,8 +54,7 @@ def __new__(cls: type, name: str, bases: Tuple[type], attrs: Dict[str, object]) -> type: for key, value in list(attrs.items()): if not key.startswith('_'): - value: - str + value: str attrs[key] = EnvVar(key, value) return super().__new__(cls, name, bases, attrs) diff --git a/tvb/recon/cli/fs.py b/tvb/recon/cli/fs.py index 58d2af5..ed6a35c 100644 --- a/tvb/recon/cli/fs.py +++ b/tvb/recon/cli/fs.py @@ -204,6 +204,14 @@ class mris_decimate(BaseFsCLI): exe = 'mris_decimate' +class mris_info(BaseFsCLI): + """ + The mris_info command from the fs package + + """ + exe = 'mris_info' + + class mris_smooth(BaseFsCLI): """ The mris_smooth command from the fs package. @@ -276,17 +284,17 @@ def convert(in_, out, # TODO in/out typed Union[Volume, VolFile] def binarize(in_file: str, out_file: str, - min_val: Union[int, float, NoneType]=None, + min_val: Union[int, float, None]=None, erode: int=0, dilate: int=0, - mask_file: Union[str, NoneType]=None, + mask_file: Union[str, None]=None, ): "Binarize image, possibly with dilation or erosion." cmd = mri_binarize args = [ cmd.exe, cmd.Flags.in_file, in_file, - cmd.Flags.out_file, out_file, + cmd.Flags.out_file, out_file] if min_val is not None: args += [cmd.Flags.min_val, str(min_val)] @@ -298,3 +306,10 @@ def binarize(in_file: str, out_file: str, args += [cmd.Flags.mask_file, mask_file] return args + +def mris_convert_run(in_: os.PathLike, out: os.PathLike): + return [mris_convert.exe, in_, out] + + +def mris_info_run(in_: os.PathLike): + return [mris_info.exe, in_] \ No newline at end of file diff --git a/tvb/recon/cli/runner.py b/tvb/recon/cli/runner.py index 3b0b0f1..eac119c 100644 --- a/tvb/recon/cli/runner.py +++ b/tvb/recon/cli/runner.py @@ -45,7 +45,7 @@ def tmp_fname(self, fname: str): pass def fname(self, fname: str): pass @abc.abstractmethod - def run(self, args): pass + def run(self, args, capture_output: bool): pass class SimpleRunner(Runner): @@ -110,10 +110,18 @@ def stringify_args(self, args) -> Sequence[str]: str_args.append(val) return str_args - def run(self, args, **kwargs): + def run(self, args, capture_output=False, **kwargs): str_args = self.stringify_args(args) for i, arg in enumerate(str_args): self.logger.debug('args[%d]: %r', i, arg) tic = time.time() - subprocess.check_call(str_args, **kwargs) + + if not capture_output: + subprocess.check_call(str_args, **kwargs) + output = None + else: + output = subprocess.check_output(str_args, stderr=subprocess.STDOUT, universal_newlines=True, **kwargs) + self.logger.debug('elapsed %.2fs' % (time.time() - tic, )) + return output + diff --git a/tvb/recon/flow/surfaces_to_structural_datasets.py b/tvb/recon/flow/surfaces_to_structural_datasets.py new file mode 100644 index 0000000..1497965 --- /dev/null +++ b/tvb/recon/flow/surfaces_to_structural_datasets.py @@ -0,0 +1,470 @@ +import enum +import numpy as np +import os +import os.path +import logging +import re +import tempfile +import time +from typing import List, Optional +from zipfile import ZipFile + +from tvb.recon.cli.runner import Runner, File +from tvb.recon.cli import fs +from tvb.recon.flow.core import Flow +from tvb.recon.algo.structural_dataset import StructuralDataset +from tvb.recon.io.factory import IOUtils + + +SUBCORTICAL_REG_INDS = [8, 10, 11, 12, 13, 16, 17, 18, 26, 47, 49, 50, 51, 52, 53, 54, 58] +FS_LUT_LH_SHIFT = 1000 +FS_LUT_RH_SHIFT = 2000 + + +class Hemisphere(enum.Enum): + rh = 'rh' + lh = 'lh' + + +class Surface: + def __init__(self, vertices: np.array, triangles: np.array, region_mapping: np.array): + assert vertices.ndim == 2 + assert triangles.ndim == 2 + assert region_mapping.ndim == 1 + + assert vertices.shape[1] == 3 + assert triangles.shape[1] == 3 + assert region_mapping.shape[0] == vertices.shape[0] + + self.vertices = vertices + self.triangles = triangles + self.region_mapping = region_mapping + + self.nverts = self.vertices.shape[0] + self.ntriangs = self.triangles.shape[0] + + self.vertex_triangles = compute_vertex_triangles(self.nverts, self.ntriangs, self.triangles) + self.triangle_normals = compute_triangle_normals(self.triangles, self.vertices) + self.triangle_angles = compute_triangle_angles(self.vertices, self.ntriangs, self.triangles) + self.vertex_normals = compute_vertex_normals(self.nverts, self.vertex_triangles, self.triangles, + self.triangle_angles, self.triangle_normals, self.vertices) + self.triangle_areas = compute_triangle_areas(self.vertices, self.triangles) + + def save_surf_zip(self, filename): + tmpdir = tempfile.TemporaryDirectory() + + file_vertices = os.path.join(tmpdir.name, 'vertices.txt') + file_triangles = os.path.join(tmpdir.name, 'triangles.txt') + file_normals = os.path.join(tmpdir.name, 'normals.txt') + + np.savetxt(file_vertices, self.vertices, fmt='%.6f %.6f %.6f') + np.savetxt(file_triangles, self.triangles, fmt='%d %d %d') + np.savetxt(file_normals, self.vertex_normals, fmt='%.6f %.6f %.6f') + + with ZipFile(filename, 'w') as zip_file: + zip_file.write(file_vertices, os.path.basename(file_vertices)) + zip_file.write(file_triangles, os.path.basename(file_triangles)) + zip_file.write(file_normals, os.path.basename(file_normals)) + + def save_region_mapping_txt(self, filename): + np.savetxt(filename, self.region_mapping, fmt="%d") + + +class ColorLut: + def __init__(self, filename: os.PathLike): + table = np.genfromtxt(os.fspath(filename), dtype=None) + + if len(table.dtype) == 6: + # id name R G B A + self.inds = table[table.dtype.names[0]] + self.names = table[table.dtype.names[1]].astype('U') + self.r = table[table.dtype.names[2]] + self.g = table[table.dtype.names[3]] + self.b = table[table.dtype.names[4]] + self.a = table[table.dtype.names[5]] + self.shortnames = np.zeros(len(self.names), dtype='U') + + elif len(table.dtype) == 7: + # id shortname name R G B A + self.inds = table[table.dtype.names[0]] + self.shortnames = table[table.dtype.names[1]].astype('U') + self.names = table[table.dtype.names[2]].astype('U') + self.r = table[table.dtype.names[3]] + self.g = table[table.dtype.names[4]] + self.b = table[table.dtype.names[5]] + self.a = table[table.dtype.names[6]] + + +class RegionIndexMapping: + + def __init__(self, color_lut_src_file: os.PathLike, color_lut_trg_file: os.PathLike): + self.src_table = ColorLut(color_lut_src_file) + self.trg_table = ColorLut(color_lut_trg_file) + + names_to_trg = dict(zip(self.trg_table.names, self.trg_table.inds)) + + self.src_to_trg = dict() + for src_ind, src_name in zip(self.src_table.inds, self.src_table.names): + trg_ind = names_to_trg.get(src_name, None) + if trg_ind is not None: + self.src_to_trg[src_ind] = trg_ind + + self.unknown_ind = names_to_trg.get('Unknown', 0) # zero as the default unknown area + + def source_to_target(self, index): + return self.src_to_trg.get(index, self.unknown_ind) + + +def merge_surfaces(surfaces: List[Surface]) -> Surface: + offsets = np.cumsum([0] + [vs.shape[0] for vs in [surf.vertices for surf in surfaces]][:-1]) + vertices = np.vstack([surf.vertices for surf in surfaces]) + triangles = np.vstack([ts + offset for ts, offset in zip([surf.triangles for surf in surfaces], offsets)]) + region_mappings = np.hstack([surf.region_mapping for surf in surfaces]) + return Surface(vertices, triangles, region_mappings) + + +def compute_vertex_triangles(number_of_vertices, number_of_triangles, triangles): + vertex_triangles = [[] for _ in range(number_of_vertices)] + for k in range(number_of_triangles): + vertex_triangles[triangles[k, 0]].append(k) + vertex_triangles[triangles[k, 1]].append(k) + vertex_triangles[triangles[k, 2]].append(k) + return vertex_triangles + + +def compute_triangle_normals(triangles, vertices): + """Calculates triangle normals.""" + tri_u = vertices[triangles[:, 1], :] - vertices[triangles[:, 0], :] + tri_v = vertices[triangles[:, 2], :] - vertices[triangles[:, 0], :] + tri_norm = np.cross(tri_u, tri_v) + + try: + triangle_normals = tri_norm / np.sqrt(np.sum(tri_norm ** 2, axis=1))[:, np.newaxis] + except FloatingPointError: + # TODO: NaN generation would stop execution, however for normals this case could maybe be + # handled in a better way. + triangle_normals = tri_norm + return triangle_normals + + +def compute_triangle_angles(vertices, number_of_triangles, triangles): + """ + Calculates the inner angles of all the triangles which make up a surface + """ + verts = vertices + # TODO: Should be possible with arrays, ie not nested loops... + # A short profile indicates this function takes 95% of the time to compute normals + # (this was a direct translation of some old matlab code) + angles = np.zeros((number_of_triangles, 3)) + for tt in range(number_of_triangles): + triangle = triangles[tt, :] + for ta in range(3): + ang = np.roll(triangle, -ta) + angles[tt, ta] = np.arccos(np.dot( + (verts[ang[1], :] - verts[ang[0], :]) / + np.sqrt(np.sum((verts[ang[1], :] - verts[ang[0], :]) ** 2, axis=0)), + (verts[ang[2], :] - verts[ang[0], :]) / + np.sqrt(np.sum((verts[ang[2], :] - verts[ang[0], :]) ** 2, axis=0)))) + + return angles + + +def compute_vertex_normals(number_of_vertices, vertex_triangles, triangles, + triangle_angles, triangle_normals, vertices): + """ + Estimates vertex normals, based on triangle normals weighted by the + angle they subtend at each vertex... + """ + vert_norms = np.zeros((number_of_vertices, 3)) + bad_normal_count = 0 + for k in range(number_of_vertices): + try: + tri_list = list(vertex_triangles[k]) + angle_mask = triangles[tri_list, :] == k + angles = triangle_angles[tri_list, :] + angles = angles[angle_mask][:, np.newaxis] + angle_scaling = angles / np.sum(angles, axis=0) + vert_norms[k, :] = np.mean(angle_scaling * triangle_normals[tri_list, :], axis=0) + # Scale by angle subtended. + vert_norms[k, :] = vert_norms[k, :] / np.sqrt(np.sum(vert_norms[k, :] ** 2, axis=0)) + # Normalise to unit vectors. + except (ValueError, FloatingPointError): + # If normals are bad, default to position vector + # A nicer solution would be to detect degenerate triangles and ignore their + # contribution to the vertex normal + vert_norms[k, :] = vertices[k] / np.sqrt(vertices[k].dot(vertices[k])) + bad_normal_count += 1 + if bad_normal_count: + print(" %d vertices have bad normals" % bad_normal_count) + return vert_norms + + +def compute_triangle_areas(vertices, triangles): + """Calculates the area of triangles making up a surface.""" + tri_u = vertices[triangles[:, 1], :] - vertices[triangles[:, 0], :] + tri_v = vertices[triangles[:, 2], :] - vertices[triangles[:, 0], :] + tri_norm = np.cross(tri_u, tri_v) + triangle_areas = np.sqrt(np.sum(tri_norm ** 2, axis=1)) / 2.0 + triangle_areas = triangle_areas[:, np.newaxis] + return triangle_areas + + +def compute_region_orientations(regions, vertex_normals, region_mapping): + """Compute the orientation of given regions from vertex_normals and region mapping""" + + average_orientation = np.zeros((regions.size, 3)) + # Average orientation of the region + for i, k in enumerate(regions): + orient = vertex_normals[region_mapping == k, :] + if orient.shape[0] > 0: + avg_orient = np.mean(orient, axis=0) + average_orientation[i, :] = avg_orient / np.sqrt(np.sum(avg_orient ** 2)) + + return average_orientation + + +def compute_region_areas(regions, triangle_areas, vertex_triangles, region_mapping): + """Compute the areas of given regions""" + + region_surface_area = np.zeros(regions.size) + avt = np.array(vertex_triangles) + # NOTE: Slightly overestimates as it counts overlapping border triangles, + # but, not really a problem provided triangle-size << region-size. + for i, k in enumerate(regions): + regs = map(set, avt[region_mapping == k]) + region_triangles = set.union(*regs) + if region_triangles: + region_surface_area[i] = triangle_areas[list(region_triangles)].sum() + + return region_surface_area + + +def compute_region_centers(regions, vertices, region_mapping): + + region_centers = np.zeros((regions.size, 3)) + for i, k in enumerate(regions): + vert = vertices[region_mapping == k, :] + if vert.shape[0] > 0: + region_centers[i, :] = np.mean(vert, axis=0) + + return region_centers + + +def compute_region_params(surface: Surface, subcortical: bool=False)\ + -> (np.ndarray, np.ndarray, np.ndarray, np.ndarray): + verts, triangs, region_mapping = surface.vertices, surface.triangles, surface.region_mapping + + regions = np.unique(region_mapping) + areas = compute_region_areas(regions, surface.triangle_areas, surface.vertex_triangles, region_mapping) + orientations = compute_region_orientations(regions, surface.vertex_normals, region_mapping) + centers = compute_region_centers(regions, verts, region_mapping) + + return regions, areas, orientations, centers + +def extract_vector(string: str, name: str) -> Optional[List[float]]: + r""" + Extract numerical vector from a block of text. The vector has to be on a single line with the format: + : (x0, x1, x2 [,...]) + If the vector in the correct format is missing, return None. + + >>> extract_vector("EXAMPLE\na: (1.0, 2.0, 3.0)\nb: (0.0, 0.0, 0.0)", "a") + [1.0, 2.0, 3.0] + + >>> extract_vector("EMPTY", "a") is None + True + """ + + for line in string.split("\n"): + match = re.match(r"""^\s* + (.+?) # name + \s*:\s* # separator + \(([0-9.,\s-]+)\) # vector: (x0, x1, ....) + \s*$""", + line, re.X) + + if match and match.group(1) == name: + try: + vector = [float(x) for x in match.group(2).split(",")] + return vector + except ValueError: + pass + + return None + + +def pial_to_verts_and_triangs(runner: Runner, pial_surf: File) -> (np.ndarray, np.ndarray): + pial_asc = runner.tmp_fname(os.path.basename(pial_surf.path + ".asc")) + runner.run(fs.mris_convert_run(pial_surf, pial_asc)) + + with open(pial_asc, 'r') as f: + f.readline() + nverts, ntriangs = [int(n) for n in f.readline().strip().split(' ')] + + vertices = np.genfromtxt(pial_asc.path, dtype=float, skip_header=2, skip_footer=ntriangs, usecols=(0, 1, 2)) + triangles = np.genfromtxt(pial_asc.path, dtype=int, skip_header=2+nverts, usecols=(0, 1, 2)) + assert vertices.shape == (nverts, 3) + assert triangles.shape == (ntriangs, 3) + + mris_info = runner.run(fs.mris_info_run(pial_surf), capture_output=True) + c_ras_list = extract_vector(mris_info, "c_(ras)") + assert c_ras_list is not None + vertices[:, 0:3] += np.array(c_ras_list) + + return vertices, triangles + + +def read_cortical_region_mapping(label_direc: os.PathLike, hemisphere: Hemisphere, fs_to_conn: RegionIndexMapping)\ + -> np.ndarray: + filename = os.path.join(label_direc, hemisphere.value + ".aparc.annot") + annot = IOUtils.read_annotation(filename) + region_mapping = annot.region_mapping + + region_mapping[region_mapping == -1] = 0 # Unknown regions in hemispheres + + # $FREESURFER_HOME/FreeSurferColorLUT.txt describes the shift + if hemisphere == Hemisphere.lh: + region_mapping += FS_LUT_LH_SHIFT + else: + region_mapping += FS_LUT_RH_SHIFT + + fs_to_conn_fun = np.vectorize(lambda n: fs_to_conn.source_to_target(n)) + region_mapping = fs_to_conn_fun(region_mapping) + + return region_mapping + + +def get_cortical_surfaces(runner: Runner, cort_surf_direc: os.PathLike, label_direc: os.PathLike, + region_index_mapping: RegionIndexMapping) -> Surface: + verts_l, triangs_l = pial_to_verts_and_triangs(runner, + File(os.path.join(cort_surf_direc, Hemisphere.lh.value + ".pial"))) + verts_r, triangs_r = pial_to_verts_and_triangs(runner, + File(os.path.join(cort_surf_direc, Hemisphere.rh.value + ".pial"))) + + region_mapping_l = read_cortical_region_mapping(label_direc, Hemisphere.lh, region_index_mapping) + region_mapping_r = read_cortical_region_mapping(label_direc, Hemisphere.rh, region_index_mapping) + + surface = merge_surfaces([Surface(verts_l, triangs_l, region_mapping_l), + Surface(verts_r, triangs_r, region_mapping_r)]) + + return surface + + +def get_subcortical_surfaces(subcort_surf_direc: os.PathLike, region_index_mapping: RegionIndexMapping) -> Surface: + surfaces = [] + + for fs_idx in SUBCORTICAL_REG_INDS: + conn_idx = region_index_mapping.source_to_target(fs_idx) + filename = os.path.join(subcort_surf_direc, 'aseg_%03d.srf' % fs_idx) + with open(filename, 'r') as f: + f.readline() + nverts, ntriangs = [int(n) for n in f.readline().strip().split(' ')] + + vertices = np.genfromtxt(filename, dtype=float, skip_header=2, skip_footer=ntriangs, usecols=(0, 1, 2)) + triangles = np.genfromtxt(filename, dtype=int, skip_header=2 + nverts, usecols=(0, 1, 2)) + region_mapping = conn_idx * np.ones(nverts, dtype=int) + surfaces.append(Surface(vertices, triangles, region_mapping)) + + surface = merge_surfaces(surfaces) + return surface + + +class SurfacesToStructuralDataset(Flow): + + def __init__(self, + cort_surf_direc: os.PathLike, + label_direc: os.PathLike, + subcort_surf_direc: os.PathLike, + source_lut: os.PathLike, + target_lut: os.PathLike, + weights_file: os.PathLike, + tract_lengths_file: os.PathLike, + struct_zip_file: os.PathLike, + out_surfaces_dir: os.PathLike=None): + """ + + Parameters + ---------- + cort_surf_direc: Directory that should contain: + rh.pial + lh.pial + label_direc: Directory that should contain: + rh.aparc.annot + lh.aparc.annot + subcort_surf_direc: Directory that should contain: + aseg_.srf + for each in SUBCORTICAL_REG_INDS + source_lut: File with the color look-up table used for the original parcellation + target_lut: File with the color look-up table used for the connectome generation + weights_file: text file with weights matrix (which should be upper triangular) + tract_lengths_file: text file with tract length matrix (which should be upper triangular) + struct_zip_file: zip file containing TVB structural dataset to be created + out_surfaces_dir: directory where to put the surfaces and region mappings in TVB format + """ + self.cort_surf_direc = cort_surf_direc + self.subcort_surf_direc = subcort_surf_direc + self.label_direc = label_direc + self.source_lut = source_lut + self.target_lut = target_lut + self.weights_file = weights_file + self.tract_lenghts_file = tract_lengths_file + self.struct_zip_file = struct_zip_file + self.out_surfaces_dir = out_surfaces_dir + + def run(self, runner: Runner, include_unknown=False): + + log = logging.getLogger('SurfacesToStructuralDataset') + tic = time.time() + + region_index_mapping = RegionIndexMapping(self.source_lut, self.target_lut) + + surf_subcort = get_subcortical_surfaces(self.subcort_surf_direc, region_index_mapping) + surf_cort = get_cortical_surfaces(runner, self.cort_surf_direc, self.label_direc, region_index_mapping) + + if self.out_surfaces_dir: + surf_subcort.save_region_mapping_txt(os.path.join(self.out_surfaces_dir, "region_mapping_subcort.txt")) + surf_subcort.save_surf_zip(os.path.join(self.out_surfaces_dir, "surface_subcort.zip")) + surf_cort.save_region_mapping_txt(os.path.join(self.out_surfaces_dir, "region_mapping_cort.txt")) + surf_cort.save_surf_zip(os.path.join(self.out_surfaces_dir, "surface_cort.zip")) + + region_params_subcort = compute_region_params(surf_subcort, True) + region_params_cort = compute_region_params(surf_cort, False) + + nregions = max(region_index_mapping.trg_table.inds) + 1 + orientations = np.zeros((nregions, 3)) + areas = np.zeros(nregions) + centers = np.zeros((nregions, 3)) + cortical = np.zeros(nregions, dtype=bool) + + for region_params, is_cortical in [(region_params_subcort, False), (region_params_cort, True)]: + regions, reg_areas, reg_orientations, reg_centers = region_params + orientations[regions, :] = reg_orientations + areas[regions] = reg_areas + centers[regions, :] = reg_centers + cortical[regions] = is_cortical + + weights = np.genfromtxt(os.fspath(self.weights_file)) + tract_lengths = np.genfromtxt(os.fspath(self.tract_lenghts_file)) + + if not include_unknown: + # Remove the region from orientations, areas and centers + indices = list(range(0, region_index_mapping.unknown_ind)) \ + + list(range(region_index_mapping.unknown_ind + 1, nregions)) + + names = region_index_mapping.trg_table.names[indices] + orientations = orientations[indices] + areas = areas[indices] + centers = centers[indices] + cortical = cortical[indices] + else: + # Add the region to weights and tract lengths + names = region_index_mapping.trg_table.names + + np.insert(weights, region_index_mapping.unknown_ind, 0.0, axis=0) + np.insert(weights, region_index_mapping.unknown_ind, 0.0, axis=1) + np.insert(tract_lengths, region_index_mapping.unknown_ind, 0.0, axis=0) + np.insert(tract_lengths, region_index_mapping.unknown_ind, 0.0, axis=1) + + dataset = StructuralDataset(orientations, areas, centers, cortical, weights, tract_lengths, names) + dataset.save_to_txt_zip(self.struct_zip_file) + + log.info('complete in %0.2fs', time.time() - tic) diff --git a/tvb/recon/tests/flow/test_surfaces_to_structural_datasets.py b/tvb/recon/tests/flow/test_surfaces_to_structural_datasets.py new file mode 100644 index 0000000..38ec526 --- /dev/null +++ b/tvb/recon/tests/flow/test_surfaces_to_structural_datasets.py @@ -0,0 +1,64 @@ +# -*- coding: utf-8 -*- + +import numpy as np +from tvb.recon.tests.base import BaseTest + +import tvb.recon.flow.surfaces_to_structural_datasets as stsd +from tvb.recon.flow.surfaces_to_structural_datasets import Surface + +FLOAT_TOL = 1e-16 + +class MinimalSurfaceTest(BaseTest): + + def setUp(self): + super().setUp() + + self.surf1 = Surface(np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0]]), + np.array([[0, 1, 2], [1, 3, 2]]), + np.array([1, 1, 1, 1])) + self.surf2 = Surface(np.array([[0, 0, 0], [0, 1, 0], [0, 0, 1], [0, 1, 1]]), + np.array([[0, 1, 2], [1, 3, 2]]), + np.array([1, 1, 1, 1])) + + self.nverts1 = self.surf1.vertices.shape[0] + self.nverts2 = self.surf2.vertices.shape[0] + self.ntri1 = self.surf1.triangles.shape[0] + self.ntri2 = self.surf2.triangles.shape[0] + + def test_merge_surfaces(self): + surf_merged = stsd.merge_surfaces([self.surf1, self.surf2]) + + self.assertEqual(surf_merged.vertices.shape, (self.nverts1 + self.nverts2, 3)) + self.assertEqual(surf_merged.triangles.shape, (self.ntri1 + self.ntri2, 3)) + self.assertEqual(surf_merged.region_mapping.shape, (self.nverts1 + self.nverts2,)) + + def test_compute_triangle_normals(self): + normals = stsd.compute_triangle_normals(self.surf1.triangles, self.surf1.vertices) + self.assertTrue(np.allclose(normals, + np.array([[0., 0., 1.], [0., 0., 1.]]), + atol=FLOAT_TOL)) + + def test_compute_vertex_normals(self): + vert_triangles = stsd.compute_vertex_triangles(self.nverts1, self.ntri1, self.surf1.triangles) + tri_angles = stsd.compute_triangle_angles(self.surf1.vertices, self.ntri1, self.surf1.triangles) + tri_normals = stsd.compute_triangle_normals(self.surf1.triangles, self.surf1.vertices) + vert_normals = stsd.compute_vertex_normals(self.nverts1, vert_triangles, self.surf1.triangles, + tri_angles, tri_normals, self.surf1.vertices) + self.assertTrue(np.allclose(vert_normals, + np.array([[0., 0., 1.], [0., 0., 1.], [0., 0., 1.], [0., 0., 1.]]), + atol=FLOAT_TOL)) + + def test_compute_triangle_areas(self): + areas = stsd.compute_triangle_areas(self.surf1.vertices, self.surf1.triangles) + self.assertTrue(np.allclose(areas, + np.array([0.5, 0.5]), + atol=FLOAT_TOL)) + + def test_compute_region_params(self): + surf = stsd.merge_surfaces([self.surf1, self.surf2]) + regions, areas, orientations, centers = stsd.compute_region_params(surf, False) + + self.assertTrue(np.equal(regions, [1])) + self.assertTrue(np.allclose(areas, [2.0], atol=FLOAT_TOL)) + self.assertTrue(np.allclose(orientations, np.array([1, 0, 1])*np.sqrt(2)/2, atol=FLOAT_TOL)) + self.assertTrue(np.allclose(centers, [0.25, 0.5, 0.25], atol=FLOAT_TOL))