Python IOUtils.write_volume方法代码示例

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import write_volume [as 别名]
    def label_with_dilation(self, to_label_nii_fname: os.PathLike, dilated_nii_fname: os.PathLike,
                            out_nii_fname: os.PathLike):
        """
        Labels a volume using its labeled dilation. The dilated volume is labeled using scipy.ndimage.label function.
        :param to_label_nii_fname: usually a CT-mask.nii.gz
        :param dilated_nii_fname: dilated version of the to_label_nii_fname volume
        """

        # TODO could make dilation with ndimage also.
        mask = IOUtils.read_volume(to_label_nii_fname)
        dil_mask = IOUtils.read_volume(dilated_nii_fname)

        lab, n = scipy.ndimage.label(dil_mask.data)

        # TODO: this change is from tvb-make. Keep it or not? It returns a different result than the old version.
        lab_xyz = list(self.compute_label_volume_centers(lab, dil_mask.affine_matrix))
        lab_sort = numpy.r_[:n + 1]
        # sort labels along AP axis
        for i, (val, _) in enumerate(sorted(lab_xyz, key=lambda t: t[1][1])):
            lab_sort[val] = i
        lab = lab_sort[lab]

        mask.data *= lab
        self.logger.info(
            '%d objects found when labeling the dilated volume.', n)

        IOUtils.write_volume(out_nii_fname, mask)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import write_volume [as 别名]
def test_write_volume():
    in_file_path = get_data_file(
        TEST_MODIF_SUBJECT, TEST_VOLUME_FOLDER, "T1.nii.gz")
    volume = IOUtils.read_volume(in_file_path)
    out_file_path = get_temporary_files_path('T1-out.nii.gz')
    IOUtils.write_volume(out_file_path, volume)
    assert os.path.exists(out_file_path)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import write_volume [as 别名]
def test_remove_zero_connectivity():
    service = VolumeService()

    data = numpy.array([[[0, 0, 1], [2, 3, 0]], [[4, 0, 0], [0, 0, 0]]])
    volume = Volume(data, [[1, 0, 0, 0], [0, 1, 0, 0],
                           [0, 0, 1, 0], [0, 0, 0, 1]], None)
    volume_path = get_temporary_files_path("tdi_lbl.nii.gz")

    IOUtils.write_volume(volume_path, volume)

    in_connectivity = numpy.array(
        [[10, 1, 0, 3], [0, 10, 0, 2], [0, 0, 0, 0], [0, 0, 0, 10]])
    connectivity_path = get_temporary_files_path("conn.csv")
    numpy.savetxt(connectivity_path, in_connectivity, fmt='%1d')

    tract_lengths_path = get_temporary_files_path("tract_lengths.csv")
    numpy.savetxt(tract_lengths_path, in_connectivity, fmt='%1d')

    service.remove_zero_connectivity_nodes(
        volume_path, connectivity_path, tract_lengths_path)

    assert os.path.exists(os.path.splitext(connectivity_path)[0] + ".npy")
    assert os.path.exists(os.path.splitext(tract_lengths_path)[0] + ".npy")

    vol = IOUtils.read_volume(volume_path)
    assert len(numpy.unique(vol.data)) == 4

    conn = numpy.array(numpy.genfromtxt(connectivity_path, dtype='int64'))
    assert numpy.array_equal(conn, [[20, 1, 3], [1, 20, 2], [3, 2, 20]])

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import write_volume [as 别名]
def test_label_with_dilation():
    service = VolumeService()

    ct_mask_data = numpy.array(
        [[[0, 0, 0], [0, 1, 0], [0, 1, 0]], [[1, 1, 1], [0, 0, 0], [0, 0, 0]], [[0, 0, 1], [0, 0, 0], [0, 0, 1]]])
    ct_mask_volume = Volume(ct_mask_data, [[1, 0, 0, 0], [0, 1, 0, 0], [
                            0, 0, 1, 0], [0, 0, 0, 1]], None)
    ct_mask_path = get_temporary_files_path("ct_mask.nii.gz")
    IOUtils.write_volume(ct_mask_path, ct_mask_volume)

    ct_dil_mask_data = numpy.array(
        [[[0, 0, 0], [1, 1, 1], [0, 1, 0]], [[1, 1, 1], [0, 0, 0], [0, 0, 0]], [[0, 1, 1], [0, 0, 0], [0, 1, 1]]])
    ct_dil_mask_volume = Volume(ct_dil_mask_data, [[1, 0, 0, 0], [
                                0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]], None)
    ct_dil_mask_path = get_temporary_files_path("ct_dil_mask.nii.gz")
    IOUtils.write_volume(ct_dil_mask_path, ct_dil_mask_volume)

    ct_result = get_temporary_files_path("ct_res.nii.gz")

    service.label_with_dilation(ct_mask_path, ct_dil_mask_path, ct_result)

    assert os.path.exists(ct_mask_path)
    assert os.path.exists(ct_dil_mask_path)
    assert os.path.exists(ct_result)

    vol = IOUtils.read_volume(ct_result)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import write_volume [as 别名]
    def simple_label_config(self, in_aparc_path: os.PathLike, out_volume_path: os.PathLike):
        """
        Relabel volume to have contiguous values like Mrtrix' labelconfig.
        :param in_aparc_path: volume voxel value is the index of the region it belongs to.
        :return: writes the labeled volume to out_volume_path.
        """

        aparc = IOUtils.read_volume(in_aparc_path)
        aparc = self._label_config(aparc)
        IOUtils.write_volume(out_volume_path, aparc)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import write_volume [as 别名]
    def label_vol_from_tdi(self, tdi_volume_path: os.PathLike, out_volume_path: os.PathLike, lo: float=0.5):
        """
        Creates a mask of the voxels with tract ends > lo and any other voxels become 0.
        Labels each voxel different from 0 with integer labels starting from 1.
        :param tdi_volume_path: volume voxel value is the sum of tract ends. Voxel without tract ends has value 0.
        :param lo: tract ends threshold used for masking.
        :return: writes labeled volume to :ut_volume_path.
        """

        nii_volume = IOUtils.read_volume(tdi_volume_path)
        tdi_volume = self._label_volume(nii_volume, lo)
        IOUtils.write_volume(out_volume_path, tdi_volume)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import write_volume [as 别名]
    def remove_zero_connectivity_nodes(self, node_volume_path: os.PathLike, connectivity_matrix_path: os.PathLike,
                                       tract_length_path: Optional[str]=None):
        """
        It removes network nodes with zero connectivity from the volume and connectivity matrices.
        The zero connectivity nodes will be labeled with 0 in the volume and the remaining labels will be updated.
        The connectivity matrices will be symmetric.
        :param node_volume_path: tdi_lbl.nii volume path
        :param connectivity_matrix_path: .csv file, output of Mrtrix3 tck2connectome
        :param tract_length_path: optional .csv tract lengths matrix
        :return: overwrites the input volume and matrices with the processed ones. Also saves matrices as .npy.
        """

        node_volume = IOUtils.read_volume(node_volume_path)

        connectivity = numpy.array(numpy.genfromtxt(
            connectivity_matrix_path, dtype='int64'))
        connectivity = connectivity + connectivity.T
        connectivity_row_sum = numpy.sum(connectivity, axis=0)

        nodes_to_keep_indices = connectivity_row_sum > 0
        connectivity = connectivity[nodes_to_keep_indices, :][
                       :, nodes_to_keep_indices]

        numpy.save(os.path.splitext(connectivity_matrix_path)
                   [0] + NPY_EXTENSION, connectivity)
        numpy.savetxt(connectivity_matrix_path, connectivity, fmt='%1d')

        if os.path.exists(str(tract_length_path)):
            connectivity = numpy.array(numpy.genfromtxt(
                tract_length_path, dtype='int64'))
            connectivity = connectivity[nodes_to_keep_indices, :][
                           :, nodes_to_keep_indices]

            numpy.save(os.path.splitext(tract_length_path)
                       [0] + NPY_EXTENSION, connectivity)
            numpy.savetxt(tract_length_path, connectivity, fmt='%1d')

        else:
            self.logger.warning("Path %s is not valid.", tract_length_path)

        nodes_to_remove_indices, = numpy.where(~nodes_to_keep_indices)
        nodes_to_remove_indices += 1

        for node_index in nodes_to_remove_indices:
            node_volume.data[node_volume.data == node_index] = 0

        node_volume.data[node_volume.data > 0] = numpy.r_[
                                                 1:(connectivity.shape[0] + 1)]

        IOUtils.write_volume(node_volume_path, node_volume)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import write_volume [as 别名]
def create_tvb_dataset(atlas_suffix: AtlasSuffix, mri_direc: os.PathLike,
                       region_details_direc: os.PathLike,
                       weights_file: os.PathLike,
                       tracts_file: os.PathLike,
                       out_dir: os.PathLike,
                       bring_t1=False):
    weights_matrix = numpy.loadtxt(str(weights_file), dtype='i', delimiter=' ')
    weights_matrix += weights_matrix.T

    tracts_matrix = numpy.loadtxt(str(tracts_file), dtype='f', delimiter=' ')
    tracts_matrix += tracts_matrix.T

    is_cortical_rm = numpy.genfromtxt(
        os.path.join(region_details_direc, AsegFiles.CORTICAL_TXT.value.replace("%s", atlas_suffix)), usecols=[0],
        dtype='i')
    region_names = numpy.genfromtxt(
        os.path.join(region_details_direc, AsegFiles.CENTERS_TXT.value.replace("%s", atlas_suffix)), usecols=[0],
        dtype="str")
    region_centers = numpy.genfromtxt(
        os.path.join(region_details_direc, AsegFiles.CENTERS_TXT.value.replace("%s", atlas_suffix)), usecols=[1, 2, 3])
    region_areas = numpy.genfromtxt(
        os.path.join(region_details_direc, AsegFiles.AREAS_TXT.value.replace("%s", atlas_suffix)), usecols=[0])
    region_orientations = numpy.genfromtxt(
        os.path.join(region_details_direc, AsegFiles.ORIENTATIONS_TXT.value.replace("%s", atlas_suffix)),
        usecols=[0, 1, 2])
    rm_idx = numpy.genfromtxt(
        os.path.join(region_details_direc, AsegFiles.RM_TO_APARC_ASEG_TXT.value.replace("%s", atlas_suffix)),
        usecols=[0, 1], dtype='i')
    rm_index_dict = dict(zip(rm_idx[:, 0], rm_idx[:, 1]))
    print(rm_index_dict)

    genericIO = GenericIO()
    genericIO.write_connectivity_zip(out_dir, weights_matrix, tracts_matrix, is_cortical_rm, region_names,
                                     region_centers, region_areas, region_orientations, atlas_suffix)

    aparc_aseg_file = os.path.join(mri_direc, T1Files.APARC_ASEG_NII_GZ.value.replace("%s", atlas_suffix))
    aparc_aseg_volume = IOUtils.read_volume(aparc_aseg_file)

    volume_service = VolumeService()
    aparc_aseg_cor_volume = volume_service.change_labels_of_aparc_aseg(atlas_suffix, aparc_aseg_volume, rm_index_dict,
                                                                       weights_matrix.shape[0])
    IOUtils.write_volume(os.path.join(out_dir, OutputConvFiles.APARC_ASEG_COR_NII_GZ.value.replace("%s", atlas_suffix)),
                         aparc_aseg_cor_volume)

    if bring_t1:
        shutil.copy2(os.path.join(mri_direc, "T1.nii.gz"), out_dir)

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import write_volume [as 别名]

#.........这里部分代码省略.........
            ijk2ijk = volume.affine_matrix.dot(
                numpy.dot(xyz2xyz, numpy.linalg.inv(mask_vol.affine_matrix)))

        # Construct a grid template of voxels +/- vn voxels around each ijk
        # voxel, sharing at least a corner
        grid = numpy.meshgrid(list(range(-vn, vn + 1, 1)), list(
            range(-vn, vn + 1, 1)), list(range(-vn, vn + 1, 1)), indexing='ij')
        grid = numpy.c_[numpy.array(grid[0]).flatten(), numpy.array(
            grid[1]).flatten(), numpy.array(grid[2]).flatten()]
        n_grid = grid.shape[0]

        out_volume = Volume(numpy.array(volume.data),
                            volume.affine_matrix, volume.header)

        # Initialize output indexes
        out_ijk = []

        # For each target label:
        for label_index in range(number_of_labels):
            current_label = labels[label_index]
            # Get the indexes of all voxels of this label:
            label_voxels_i, label_voxels_j, label_voxels_k = numpy.where(
                volume.data == current_label)

            for voxel_index in range(label_voxels_i.size):
                current_voxel_i, current_voxel_j, current_voxel_k = \
                    label_voxels_i[voxel_index], label_voxels_j[
                        voxel_index], label_voxels_k[voxel_index]
                # TODO if necessary: deal with voxels at the edge of the image, such as brain stem ones...
                #     if any([(i==0), (i==mask_shape[0]-1),(j==0), (j==mask_shape[0]-1),(k==0), (k==mask_shape[0]-1)]):
                #               mask_shape[i,j,k]=0
                #               continue

                # ...get the corresponding voxel in the mask volume:
                ijk = numpy.round(ijk2ijk.dot(numpy.array(
                    [current_voxel_i, current_voxel_j, current_voxel_k, 1]))[:3]).astype('i')

                # Make sure this point is within image limits
                for cc in range(3):
                    if ijk[cc] < 0:
                        ijk[cc] = 0

                    elif ijk[cc] >= mask_vol.dimensions[cc]:
                        ijk[cc] = mask_vol.dimensions[cc] - 1

                # If this is a voxel to keep, set it so...
                if mask_vol.data[ijk[0], ijk[1], ijk[2]] >= th:
                    out_volume.data[current_voxel_i, current_voxel_j,
                                    current_voxel_k] = labels_mask[label_index]
                    out_ijk.append(
                        [current_voxel_i, current_voxel_j, current_voxel_k])

                elif vn > 0:
                    # If not, and as long as vn>0 check whether any of its vn neighbors is a mask voxel.
                    # Generate the specific grid centered at the vertex ijk
                    ijk_grid = grid + numpy.tile(ijk, (n_grid, 1))

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

                    try:
                        # If none of these points is a mask point:
                        if (mask_vol.data[ijk_grid[:, 0], ijk_grid[
                                                          :, 1], ijk_grid[:, 2]] < th).all():
                            out_volume.data[
                                current_voxel_i, current_voxel_j, current_voxel_k] = labels_nomask[label_index]

                        else:  # if any of them is a mask point:
                            out_volume.data[
                                current_voxel_i, current_voxel_j, current_voxel_k] = labels_mask[label_index]
                            out_ijk.append(
                                [current_voxel_i, current_voxel_j, current_voxel_k])

                    except ValueError:  # empty grid
                        self.logger.error("Error at voxel ( %s, %s, %s ): It appears to have no common-face neighbors "
                                          "inside the image!", str(
                            current_voxel_i), str(current_voxel_j),
                                          str(current_voxel_k))
                        return

                else:
                    out_volume.data[current_voxel_i, current_voxel_j,
                                    current_voxel_k] = labels_nomask[label_index]

        if out_vol_path is None:
            out_vol_path = in_vol_path

        IOUtils.write_volume(out_vol_path, out_volume)

        # Save the output indexes that survived masking
        out_ijk = numpy.vstack(out_ijk)
        filepath = os.path.splitext(out_vol_path)[0]
        numpy.save(filepath + "-idx.npy", out_ijk)
        numpy.savetxt(filepath + "-idx.txt", out_ijk, fmt='%d')

# 需要导入模块: from tvb.recon.io.factory import IOUtils [as 别名]
# 或者: from tvb.recon.io.factory.IOUtils import write_volume [as 别名]

#.........这里部分代码省略.........
            labels=labels, ctx=ctx)
        number_of_labels = len(labels)
        # Set the labels for the surfaces
        if labels_surf is None:
            labels_surf = labels
        else:
            # Read the surface labels and make sure there is one for each label
            labels_surf = numpy.array(labels_surf.split()).astype('i')
            if len(labels_surf) == 1:
                labels_surf = numpy.repeat(
                    labels_inner, number_of_labels).tolist()
            elif len(labels_surf) != number_of_labels:
                self.logger.warning(
                    "Output labels for surface voxels are neither of length "
                    "1 nor of length equal to the one of target labels.")
                return
            else:
                labels_surf = labels_surf.tolist()
        # Read the inner, non-surface labels
        labels_inner = numpy.array(labels_inner.split()).astype('i')
        # ...and make sure there is one for each label
        if len(labels_inner) == 1:
            labels_inner = numpy.repeat(
                labels_inner, number_of_labels).tolist()
        elif len(labels_inner) != number_of_labels:
            self.logger.warning(
                "Output labels for inner voxels are neither of length 1 nor "
                "of length equal to the one of the target labels.")
            return
        else:
            labels_inner = labels_inner.tolist()

        # Read the input volume...
        volume = IOUtils.read_volume(in_vol_path)

        # Neigbors' grid sharing a face
        eye3 = numpy.identity(3)
        border_grid = numpy.c_[eye3, -eye3].T.astype('i')
        n_border = 6

        out_volume = Volume(numpy.array(volume.data),
                            volume.affine_matrix, volume.header)

        # Initialize output indexes
        out_ijk = []

        for label_index in range(number_of_labels):
            current_label = labels[label_index]
            # Get the indexes of all voxels of this label:
            label_volxels_i, label_voxels_j, label_voxels_k = numpy.where(
                volume.data == current_label)
            # and for each voxel
            for voxel_index in range(label_volxels_i.size):
                # indexes of this voxel:
                current_voxel_i, current_voxel_j, current_voxel_k = \
                    label_volxels_i[voxel_index], label_voxels_j[
                        voxel_index], label_voxels_k[voxel_index]
                # Create the neighbors' grid sharing a face
                ijk_grid = border_grid + \
                           numpy.tile(numpy.array(
                               [current_voxel_i, current_voxel_j, current_voxel_k]), (n_border, 1))
                # Remove voxels outside the image
                indices_inside_image = 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[indices_inside_image, :]
                try:
                    # If all face neighbors are of the same label...
                    if numpy.all(volume.data[ijk_grid[:, 0], ijk_grid[:, 1], ijk_grid[:, 2]] == numpy.tile(
                            volume.data[current_voxel_i,
                                        current_voxel_j, current_voxel_k],
                            (n_border, 1))):
                        # ...set this voxel to the corresponding inner target label
                        out_volume.data[current_voxel_i, current_voxel_j,
                                        current_voxel_k] = labels_inner[label_index]
                    else:
                        # ...set this voxel to the corresponding surface target label
                        out_volume.data[current_voxel_i, current_voxel_j,
                                        current_voxel_k] = labels_surf[label_index]
                        out_ijk.append(
                            [current_voxel_i, current_voxel_j, current_voxel_k])
                except ValueError:  # empty grid
                    self.logger.error("Error at voxel ( %s, %s, %s ) of label %s: It appears to have no common-face "
                                      "neighbors inside the image!", str(
                        current_voxel_i), str(current_voxel_j),
                                      str(current_voxel_k), str(current_label))
                    return

        if out_vol_path is None:
            out_vol_path = in_vol_path

        IOUtils.write_volume(out_vol_path, out_volume)

        # save the output indexes that survived masking
        out_ijk = numpy.vstack(out_ijk)
        filepath = os.path.splitext(out_vol_path)[0]
        numpy.save(filepath + "-idx.npy", out_ijk)
        numpy.savetxt(filepath + "-idx.txt", out_ijk, fmt='%d')