Python IOUtils.read_surface方法代码示例

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
def test_write_fs_surface():
    file_path = get_data_file(TEST_FS_SUBJECT, TEST_SURFACE_FOLDER, "lh.pial")
    original_surface = IOUtils.read_surface(file_path, False)
    triangles_number = len(original_surface.triangles)

    output_file_path = get_temporary_files_path("lh-test.pial")
    IOUtils.write_surface(output_file_path, original_surface)

    new_surface = IOUtils.read_surface(output_file_path, False)
    assert triangles_number == len(new_surface.triangles) == 327680

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
    def overlap_volume_surfaces(self, volume_background: os.PathLike, surfaces_path: os.PathLike,
                                use_center_surface: bool, use_cc_point: bool, snapshot_name: str=SNAPSHOT_NAME):
        volume = IOUtils.read_volume(volume_background)

        if use_cc_point:
            ras = self.generic_io.get_ras_coordinates(
                self.read_t1_affine_matrix())
        else:
            ras = volume.get_center_point()

        surfaces = [IOUtils.read_surface(os.path.expandvars(surface), use_center_surface) for surface in
                    surfaces_path]

        for projection in PROJECTIONS:
            try:
                x, y, background_matrix = volume.slice_volume(projection, ras)
            except IndexError:
                ras = volume.get_center_point()
                x, y, background_matrix = volume.slice_volume(projection, ras)
                self.logger.info("The volume center point has been used for %s snapshot of %s and %s.", projection,
                                 volume_background, surfaces_path)

            clear_flag = True
            for surface_index, surface in enumerate(surfaces):
                surf_x_array, surf_y_array = surface.cut_by_plane(
                    projection, ras)
                self.writer.write_matrix_and_surfaces(x, y, background_matrix, surf_x_array, surf_y_array,
                                                      surface_index, clear_flag)
                clear_flag = False
            self.writer.save_figure(
                self.generate_file_name(projection, snapshot_name))

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
    def test_merge_surfaces(self,):
        h5_surface_path = get_data_file("head2", "SurfaceCortical.h5")
        h5_surface = IOUtils.read_surface(h5_surface_path, False)

        vtx = h5_surface.vertices
        tri = h5_surface.triangles
        nv = vtx.size
        nt = tri.size
        nv2 = int(nv / 2)
        nt2 = int(nt / 2)

        lh_surface = Surface(vtx[:nv2], tri[:nt2])
        rh_surface = Surface(vtx[nv2:], tri[nt2:])

        # h5_region_mapping_path = get_data_file("head2", "RegionMapping.h5")
        # annotation = IOUtils.read_annotation(h5_region_mapping_path)
        #
        # lh_region_mapping = annotation.region_mapping[:len(annotation.region_mapping) / 2]
        # rh_region_mapping = annotation.region_mapping[len(annotation.region_mapping) / 2:]

        out_surface = self.service.merge_surfaces([lh_surface, rh_surface])
        self.assertEqual(len(out_surface.vertices),
                         len(lh_surface.vertices) + len(rh_surface.vertices))
        self.assertEqual(len(out_surface.triangles),
                         len(lh_surface.triangles) + len(rh_surface.triangles))

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
def test_write_write_brain_visa_surf():
    surface_path = get_data_file(
        TEST_FS_SUBJECT, TEST_SURFACE_FOLDER, "lh.pial")
    out_path = get_temporary_files_path("lh.pial.tri")

    surface = IOUtils.read_surface(surface_path, False)
    IOUtils.write_surface(out_path, surface)

    assert os.path.exists(out_path)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
 def test_overlap_surface_annotation(self):
     writer = ImageWriter(SNAPSHOTS_DIRECTORY)
     surface_path = get_data_file(self.head2, "SurfaceCortical.h5")
     surface = IOUtils.read_surface(surface_path, False)
     annot_path = get_data_file(self.head2, "RegionMapping.h5")
     annot = IOUtils.read_annotation(annot_path)
     annot.region_names = ['reg1', 'reg2']
     annot.regions_color_table = numpy.array(
         [[200, 200, 200, 255, 30567], [100, 150, 200, 255, 30568]])
     resulted_file_name = self.processor.generate_file_name(
         'surface_annotation', SNAPSHOT_NAME)
     writer.write_surface_with_annotation(surface, annot, resulted_file_name)
     fname = '%s0' % (resulted_file_name, )
     self._assert_writer_path_exists(fname)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
    def compute_gdist_mat(self, surf_name: str='pial', max_distance: float=40.0) -> numpy.ndarray:
        max_distance = float(max_distance)  # in case passed from sys.argv
        for h in 'rl':
            subjects_dir = os.environ['SUBJECTS_DIR']
            subject = os.environ['SUBJECT']
            surf_path = '%s/%s/surf/%sh.%s' % (subjects_dir,
                                               subject, h, surf_name)
            surface = IOUtils.read_surface(surf_path, False)
            mat_path = '%s/%s/surf/%sh.%s.gdist.mat' % (
                subjects_dir, subject, h, surf_name)
            mat = gdist.local_gdist_matrix(
                surface.vertices, surface.triangles.astype('<i4'), max_distance=max_distance)
            scipy.io.savemat(mat_path, {'gdist': mat})

            return mat

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
    def aseg_surf_conc_annot(self, surf_path: str, out_surf_path: str, annot_path: str,
                             label_indices: Union[numpy.ndarray, list], lut_path: Optional[str]=None) -> Surface:
        """
        Concatenate surfaces of one specific label of interest each, to create a single annotated surface.
        """

        lut_path = lut_path or default_lut_path()

        label_names, color_table = self.annotation_service.lut_to_annot_names_ctab(lut_path=lut_path,
                                                                                   labels=label_indices)
        label_indices = numpy.array(label_indices.split()).astype('i')

        #                  verts tri area_mask cras
        surfaces = []
        out_annotation = Annotation([], [], [])
        label_number = -1

        for label_index in label_indices:
            # TODO: This is hardcoded: /aseg-%06d here and also in pegasus dax generator
            this_surf_path = surf_path + "/aseg-%06d" % int(label_index)

            if os.path.exists(this_surf_path):
                ind_l, = numpy.where(label_indices == label_index)
                out_annotation.add_region_names_and_colors(
                    label_names[int(ind_l)],
                    color_table[ind_l, :])
                label_number += 1
                surfaces.append(IOUtils.read_surface(this_surf_path, False))
                out_annotation.add_region_mapping(
                    label_number * numpy.ones((surfaces[-1].n_vertices,), dtype='int64'))
        out_surface = self.merge_surfaces(surfaces)
        # out_annotation.regions_color_table = numpy.squeeze(numpy.array(out_annotation.regions_color_table).astype('i'))

        IOUtils.write_surface(out_surf_path, out_surface)
        IOUtils.write_annotation(annot_path, out_annotation)

        return out_surface

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
def merge_surfs(surf_lh, surf_rh, out_surf_path):
    s_lh = IOUtils.read_surface(surf_lh, False)
    s_rh = IOUtils.read_surface(surf_rh, False)
    surf = surfaceService.merge_surfaces([s_lh, s_rh])
    IOUtils.write_surface(out_surf_path, surf)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
def test_parse_h5_surface():
    h5_path = get_data_file('head2', 'SurfaceCortical.h5')
    surface = IOUtils.read_surface(h5_path, False)
    assert len(surface.vertices) == 16
    assert len(surface.triangles) == 24

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
def test_parse_not_existing_gifti_surface():
    file_path = "not_existing_surface.gii"
    with pytest.raises(FileNotFoundError):
        IOUtils.read_surface(file_path, False)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
def test_parse_gifti_centered_surface():
    file_path = get_data_file(
        TEST_MODIF_SUBJECT, TEST_SURFACE_FOLDER, "lh.pial.gii")
    surf = IOUtils.read_surface(file_path, True)
    assert len(surf.triangles) == 327680

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
def test_parse_centered_fs_surface():
    file_path = get_data_file(
        TEST_MODIF_SUBJECT, TEST_SURFACE_FOLDER, "lh-centered.pial")
    surf = IOUtils.read_surface(file_path, False)
    assert surf.center_ras == [0, 0, 0]

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
def test_parse_fs_surface():
    file_path = get_data_file(TEST_FS_SUBJECT, TEST_SURFACE_FOLDER, "lh.pial")
    surf = IOUtils.read_surface(file_path, False)
    assert len(surf.triangles) == 327680

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
 def convert_fs_to_brain_visa(self, in_surf_path: str, out_surf_path: Optional[str]=None):
     surface = IOUtils.read_surface(in_surf_path, False)
     if out_surf_path is None:
         out_surf_path = in_surf_path + '.tri'
     IOUtils.write_surface(out_surf_path, surface)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import read_surface [as 别名]
    def sample_vol_on_surf(self, surf_path: str, vol_path: str, annot_path: str, out_surf_path: str,
                           cras_path: str, add_string: str='', vertex_neighbourhood: int=1,
                           add_lbl: list=[], lut_path: Optional[str]=None) -> (Surface, Annotation):
        """
        Sample a volume of a specific label on a surface, by keeping only those surface vertices, the nearest voxel of
        which is of the given label (+ of possibly additional target labels, such as white matter).
        Allow optionally for vertices within a given voxel distance vn from the target voxels.
        """

        lut_path = lut_path or default_lut_path()

        # Read the inputs
        surface = IOUtils.read_surface(surf_path, False)

        annotation = IOUtils.read_annotation(annot_path)
        labels = self.annotation_service.annot_names_to_labels(annotation.region_names,
                                                               add_string=add_string, lut_path=lut_path)
        region_mapping_indexes = numpy.unique(annotation.region_mapping)

        volume_parser = VolumeIO()
        volume = volume_parser.read(vol_path)
        ras2vox_affine_matrix = numpy.linalg.inv(volume.affine_matrix)

        cras = numpy.loadtxt(cras_path)

        grid, n_grid = self.__prepare_grid(vertex_neighbourhood)

        # Initialize the output mask:
        verts_out_mask = numpy.repeat([False], surface.vertices.shape[0])
        for label_index in range(len(region_mapping_indexes)):

            self.logger.info("%s", add_string +
                             annotation.region_names[label_index])

            # Get the indexes of the vertices corresponding to this label:
            verts_indices_of_label, = numpy.where(
                annotation.region_mapping[:] == region_mapping_indexes[label_index])
            verts_indices_of_label_size = verts_indices_of_label.size
            if verts_indices_of_label_size == 0:
                continue

            # Add any additional labels
            all_labels = [labels[label_index]] + add_lbl

            # get the vertices for current label and add cras to take them to
            # scanner ras
            verts_of_label = surface.vertices[verts_indices_of_label, :]
            verts_of_label += numpy.repeat(numpy.expand_dims(
                cras, 1).T, verts_indices_of_label_size, axis=0)

            # Compute the nearest voxel coordinates using the affine transform
            ijk = numpy.round(
                ras2vox_affine_matrix.dot(numpy.c_[verts_of_label, numpy.ones(verts_indices_of_label_size)].T)[:3].T) \
                .astype('i')

            # Get the labels of these voxels:
            surf_vxls = volume.data[ijk[:, 0], ijk[:, 1], ijk[:, 2]]

            # Vertex mask to keep: those that correspond to voxels of one of
            # the target labels
            # surf_vxls==lbl if only one target label
            verts_keep, = numpy.where(numpy.in1d(surf_vxls, all_labels))
            verts_out_mask[verts_indices_of_label[verts_keep]] = True

            if vertex_neighbourhood > 0:
                # These are now the remaining indexes to be checked for
                # neighboring voxels
                verts_indices_of_label = numpy.delete(
                    verts_indices_of_label, verts_keep)
                ijk = numpy.delete(ijk, verts_keep, axis=0)

                for vertex_index in range(verts_indices_of_label.size):
                    # Generate the specific grid centered at the voxel ijk
                    ijk_grid = grid + \
                               numpy.tile(ijk[vertex_index, :], (n_grid, 1))

                    # Remove voxels outside the volume
                    indexes_within_limits = numpy.all([(ijk_grid[:, 0] >= 0), (ijk_grid[:, 0] < volume.dimensions[0]),
                                                       (ijk_grid[:, 1] >= 0), (ijk_grid[
                                                                               :, 1] < volume.dimensions[1]),
                                                       (ijk_grid[:, 2] >= 0), (ijk_grid[:, 2] < volume.dimensions[2])],
                                                      axis=0)
                    ijk_grid = ijk_grid[indexes_within_limits, :]
                    surf_vxls = volume.data[
                        ijk_grid[:, 0], ijk_grid[:, 1], ijk_grid[:, 2]]

                    # If any of the neighbors is of the target labels include
                    # the current vertex
                    # surf_vxls==lbl if only one target label
                    if numpy.any(numpy.in1d(surf_vxls, all_labels)):
                        verts_out_mask[
                            verts_indices_of_label[vertex_index]] = True

        # Vertex indexes and vertices to keep:
        verts_out_indices, = numpy.where(verts_out_mask)
        verts_out = surface.vertices[verts_out_indices]

        # TODO maybe: make sure that all voxels of this label correspond to at least one vertex.
        # Create a similar mask for faces by picking only triangles of which
        # all 3 vertices are included
#.........这里部分代码省略.........