Python QuickBundles.cluster方法代码示例

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def test_quickbundles_with_python_metric():

    class MDFpy(dipymetric.Metric):
        def are_compatible(self, shape1, shape2):
            return shape1 == shape2

        def dist(self, features1, features2):
            dist = np.sqrt(np.sum((features1 - features2)**2, axis=1))
            dist = np.sum(dist / len(features1))
            return dist

    rdata = streamline_utils.set_number_of_points(data, 10)
    qb = QuickBundles(threshold=2 * threshold, metric=MDFpy())

    clusters = qb.cluster(rdata)

    # By default `refdata` refers to data being clustered.
    assert_equal(clusters.refdata, rdata)
    # Set `refdata` to return indices instead of actual data points.
    clusters.refdata = None
    assert_array_equal(list(itertools.chain(*clusters)),
                       list(itertools.chain(*clusters_truth)))

    # Cluster read-only data
    for datum in rdata:
        datum.setflags(write=False)

    # Cluster data with different dtype (should be converted into float32)
    for datatype in [np.float64, np.int32, np.int64]:
        newdata = [datum.astype(datatype) for datum in rdata]
        clusters = qb.cluster(newdata)
        assert_equal(clusters.centroids[0].dtype, np.float32)

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def test_quickbundles_memory_leaks():
    qb = QuickBundles(threshold=2*threshold)

    type_name_pattern = "memoryview"
    initial_types_refcount = get_type_refcount(type_name_pattern)

    qb.cluster(data)
    # At this point, all memoryviews created during clustering should be freed.
    assert_equal(get_type_refcount(type_name_pattern), initial_types_refcount)

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def bench_quickbundles():
    dtype = "float32"
    repeat = 10
    nb_points = 12

    streams, hdr = nib.trackvis.read(get_fnames('fornix'))
    fornix = [s[0].astype(dtype) for s in streams]
    fornix = streamline_utils.set_number_of_points(fornix, nb_points)

    # Create eight copies of the fornix to be clustered (one in each octant).
    streamlines = []
    streamlines += [s + np.array([100, 100, 100], dtype) for s in fornix]
    streamlines += [s + np.array([100, -100, 100], dtype) for s in fornix]
    streamlines += [s + np.array([100, 100, -100], dtype) for s in fornix]
    streamlines += [s + np.array([100, -100, -100], dtype) for s in fornix]
    streamlines += [s + np.array([-100, 100, 100], dtype) for s in fornix]
    streamlines += [s + np.array([-100, -100, 100], dtype) for s in fornix]
    streamlines += [s + np.array([-100, 100, -100], dtype) for s in fornix]
    streamlines += [s + np.array([-100, -100, -100], dtype) for s in fornix]

    # The expected number of clusters of the fornix using threshold=10 is 4.
    threshold = 10.
    expected_nb_clusters = 4 * 8

    print("Timing QuickBundles 1.0 vs. 2.0")

    qb2 = QB_New(threshold)
    qb2_time = measure("clusters = qb2.cluster(streamlines)", repeat)
    print("QuickBundles2 time: {0:.4}sec".format(qb2_time))
    print("Speed up of {0}x".format(qb1_time / qb2_time))
    clusters = qb2.cluster(streamlines)
    sizes2 = map(len, clusters)
    indices2 = map(lambda c: c.indices, clusters)
    assert_equal(len(clusters), expected_nb_clusters)
    assert_array_equal(list(sizes2), sizes1)
    assert_arrays_equal(indices2, indices1)

    qb = QB_New(threshold, metric=MDFpy())
    qb3_time = measure("clusters = qb.cluster(streamlines)", repeat)
    print("QuickBundles2_python time: {0:.4}sec".format(qb3_time))
    print("Speed up of {0}x".format(qb1_time / qb3_time))
    clusters = qb.cluster(streamlines)
    sizes3 = map(len, clusters)
    indices3 = map(lambda c: c.indices, clusters)
    assert_equal(len(clusters), expected_nb_clusters)
    assert_array_equal(list(sizes3), sizes1)
    assert_arrays_equal(indices3, indices1)

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def main():
    parser = buildArgsParser()
    args = parser.parse_args()

    full_tfile = nib.streamlines.load(args.full_tfile)
    model_tfile = nib.streamlines.load(args.model_tfile)
    model_mask = nib.load(args.model_mask)

    # Bring streamlines to voxel space and where coordinate (0,0,0) represents the corner of a voxel.
    model_tfile.tractogram.apply_affine(np.linalg.inv(model_mask.affine))
    model_tfile.streamlines._data += 0.5  # Shift of half a voxel
    full_tfile.tractogram.apply_affine(np.linalg.inv(model_mask.affine))
    full_tfile.streamlines._data += 0.5  # Shift of half a voxel

    assert(model_mask.get_data().sum() == create_binary_map(model_tfile.streamlines, model_mask).sum())

    # Resample streamlines
    full_streamlines = set_number_of_points(full_tfile.streamlines, args.nb_points_resampling)
    model_streamlines = set_number_of_points(model_tfile.streamlines, args.nb_points_resampling)

    # Segment model
    rng = np.random.RandomState(42)
    indices = np.arange(len(model_streamlines))
    rng.shuffle(indices)
    qb = QuickBundles(args.qb_threshold)
    clusters = qb.cluster(model_streamlines, ordering=indices)

    # Try to find optimal assignment threshold
    best_threshold = None
    best_f1_score = -np.inf
    thresholds = np.arange(-2, 10, 0.2) + args.qb_threshold
    for threshold in thresholds:
        indices = qb.find_closest(clusters, full_streamlines, threshold=threshold)
        nb_assignments = np.sum(indices != -1)

        mask = create_binary_map(full_tfile.streamlines[indices != -1], model_mask)

        overlap_per_bundle = _compute_overlap(model_mask.get_data(), mask)
        overreach_per_bundle = _compute_overreach(model_mask.get_data(), mask)
        # overreach_norm_gt_per_bundle = _compute_overreach_normalize_gt(model_mask.get_data(), mask)
        f1_score = _compute_f1_score(overlap_per_bundle, overreach_per_bundle)
        if best_f1_score < f1_score:
            best_threshold = threshold
            best_f1_score = f1_score

        print("{}:\t {}/{} ({:.1%}) {:.1%}/{:.1%} ({:.1%}) {}/{}".format(
            threshold,
            nb_assignments, len(model_streamlines), nb_assignments/len(model_streamlines),
            overlap_per_bundle, overreach_per_bundle, f1_score,
            mask.sum(), model_mask.get_data().sum()))

        if overlap_per_bundle >= 1:
            break


    print("Best threshold: {} with F1-Score of {}".format(best_threshold, best_f1_score))

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def test_quickbundles_with_not_order_invariant_metric():
    metric = dipymetric.AveragePointwiseEuclideanMetric()
    qb = QuickBundles(threshold=np.inf, metric=metric)

    streamline = np.arange(10*3, dtype=dtype).reshape((-1, 3))
    streamlines = [streamline, streamline[::-1]]

    clusters = qb.cluster(streamlines)
    assert_equal(len(clusters), 1)
    assert_array_equal(clusters[0].centroid, streamline)

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def test_quickbundles_shape_uncompatibility():
    # QuickBundles' old default metric (AveragePointwiseEuclideanMetric, aka MDF)
    # requires that all streamlines have the same number of points.
    metric = dipymetric.AveragePointwiseEuclideanMetric()
    qb = QuickBundles(threshold=20., metric=metric)
    assert_raises(ValueError, qb.cluster, data)

    # QuickBundles' new default metric (AveragePointwiseEuclideanMetric, aka MDF
    #  combined with ResampleFeature) will automatically resample streamlines so
    #  they all have 18 points.
    qb = QuickBundles(threshold=20.)
    clusters1 = qb.cluster(data)

    feature = dipymetric.ResampleFeature(nb_points=18)
    metric = dipymetric.AveragePointwiseEuclideanMetric(feature)
    qb = QuickBundles(threshold=20., metric=metric)
    clusters2 = qb.cluster(data)

    assert_array_equal(list(itertools.chain(*clusters1)), list(itertools.chain(*clusters2)))

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def test_quickbundles_streamlines():
    rdata = streamline_utils.set_number_of_points(data, 10)
    qb = QuickBundles(threshold=2*threshold)

    clusters = qb.cluster(rdata)
    # By default `refdata` refers to data being clustered.
    assert_equal(clusters.refdata, rdata)
    # Set `refdata` to return indices instead of actual data points.
    clusters.refdata = None
    assert_array_equal(list(itertools.chain(*clusters)),
                       list(itertools.chain(*clusters_truth)))

    # Cluster read-only data
    for datum in rdata:
        datum.setflags(write=False)

    # Cluster data with different dtype (should be converted into float32)
    for datatype in [np.float64, np.int32, np.int64]:
        newdata = [datum.astype(datatype) for datum in rdata]
        clusters = qb.cluster(newdata)
        assert_equal(clusters.centroids[0].dtype, np.float32)

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def _prepare_gt_bundles_info(bundles_dir, bundles_masks_dir,
                             gt_bundles_attribs, ref_anat_fname):
    # Ref bundles will contain {'name': 'name_of_the_bundle',
    #                           'threshold': thres_value,
    #                           'streamlines': list_of_streamlines}

    dummy_attribs = {'orientation': 'LPS'}
    qb = QuickBundles(20, metric=AveragePointwiseEuclideanMetric())

    ref_bundles = []

    for bundle_idx, bundle_f in enumerate(sorted(os.listdir(bundles_dir))):
        bundle_name = os.path.splitext(os.path.basename(bundle_f))[0]

        bundle_attribs = gt_bundles_attribs.get(os.path.basename(bundle_f))
        if bundle_attribs is None:
            raise ValueError(
                "Missing basic bundle attribs for {0}".format(bundle_f))

        # Already resample to avoid doing it for each iteration of chunking
        orig_strl = [s for s in get_tracts_voxel_space_for_dipy(
                        os.path.join(bundles_dir, bundle_f),
                        ref_anat_fname, dummy_attribs)]

        resamp_bundle = set_number_of_points(orig_strl, NB_POINTS_RESAMPLE)
        resamp_bundle = [s.astype('f4') for s in resamp_bundle]

        bundle_cluster_map = qb.cluster(resamp_bundle)
        bundle_cluster_map.refdata = resamp_bundle

        bundle_mask = nib.load(os.path.join(bundles_masks_dir,
                                            bundle_name + '.nii.gz'))

        ref_bundles.append({'name': bundle_name,
                            'threshold': bundle_attribs['cluster_threshold'],
                            'cluster_map': bundle_cluster_map,
                            'mask': bundle_mask})

    return ref_bundles

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def bundle_analysis(model_bundle_folder, bundle_folder, orig_bundle_folder,
                    metric_folder, group, subject, no_disks=100,
                    out_dir=''):
    """
    Applies statistical analysis on bundles and saves the results
    in a directory specified by ``out_dir``.

    Parameters
    ----------
    model_bundle_folder : string
        Path to the input model bundle files. This path may contain
        wildcards to process multiple inputs at once.
    bundle_folder : string
        Path to the input bundle files in common space. This path may
        contain wildcards to process multiple inputs at once.
    orig_folder : string
        Path to the input bundle files in native space. This path may
        contain wildcards to process multiple inputs at once.
    metric_folder : string
        Path to the input dti metric or/and peak files. It will be used as
        metric for statistical analysis of bundles.
    group : string
        what group subject belongs to e.g. control or patient
    subject : string
        subject id e.g. 10001
    no_disks : integer, optional
        Number of disks used for dividing bundle into disks. (Default 100)
    out_dir : string, optional
        Output directory (default input file directory)

    References
    ----------
    .. [Chandio19] Chandio, B.Q., S. Koudoro, D. Reagan, J. Harezlak,
    E. Garyfallidis, Bundle Analytics: a computational and statistical
    analyses framework for tractometric studies, Proceedings of:
    International Society of Magnetic Resonance in Medicine (ISMRM),
    Montreal, Canada, 2019.

    """

    dt = dict()

    mb = os.listdir(model_bundle_folder)
    mb.sort()
    bd = os.listdir(bundle_folder)
    bd.sort()
    org_bd = os.listdir(orig_bundle_folder)
    org_bd.sort()
    n = len(org_bd)

    for io in range(n):
        mbundles, _ = load_trk(os.path.join(model_bundle_folder, mb[io]))
        bundles, _ = load_trk(os.path.join(bundle_folder, bd[io]))
        orig_bundles, _ = load_trk(os.path.join(orig_bundle_folder,
                                   org_bd[io]))

        mbundle_streamlines = set_number_of_points(mbundles,
                                                   nb_points=no_disks)

        metric = AveragePointwiseEuclideanMetric()
        qb = QuickBundles(threshold=25., metric=metric)
        clusters = qb.cluster(mbundle_streamlines)
        centroids = Streamlines(clusters.centroids)

        print('Number of centroids ', len(centroids.data))
        print('Model bundle ', mb[io])
        print('Number of streamlines in bundle in common space ',
              len(bundles))
        print('Number of streamlines in bundle in original space ',
              len(orig_bundles))

        _, indx = cKDTree(centroids.data, 1,
                          copy_data=True).query(bundles.data, k=1)

        metric_files_names = os.listdir(metric_folder)
        _, affine = load_nifti(os.path.join(metric_folder, "fa.nii.gz"))

        affine_r = np.linalg.inv(affine)
        transformed_orig_bundles = transform_streamlines(orig_bundles,
                                                         affine_r)

        for mn in range(0, len(metric_files_names)):

            ind = np.array(indx)
            fm = metric_files_names[mn][:2]
            bm = mb[io][:-4]
            dt = dict()
            metric_name = os.path.join(metric_folder,
                                       metric_files_names[mn])

            if metric_files_names[mn][2:] == '.nii.gz':
                metric, _ = load_nifti(metric_name)

                dti_measures(transformed_orig_bundles, metric, dt, fm,
                             bm, subject, group, ind, out_dir)

            else:
                fm = metric_files_names[mn][:3]
                metric = load_peaks(metric_name)
                peak_values(bundles, metric, dt, fm, bm, subject, group,
#.........这里部分代码省略.........

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
streams, hdr = tv.read(fname)

streamlines = [i[0] for i in streams]

"""
Perform QuickBundles clustering using the MDF metric and a 10mm distance
threshold. Keep in mind that since the MDF metric requires streamlines to have
the same number of points, the clustering algorithm will internally use a
representation of streamlines that have been automatically downsampled/upsampled
so they have only 12 points (To set manually the number of points,
see :ref:`clustering-examples-ResampleFeature`).
"""

qb = QuickBundles(threshold=10.)
clusters = qb.cluster(streamlines)

"""
`clusters` is a `ClusterMap` object which contains attributes that
provide information about the clustering result.
"""

print("Nb. clusters:", len(clusters))
print("Cluster sizes:", map(len, clusters))
print("Small clusters:", clusters < 10)
print("Streamlines indices of the first cluster:\n", clusters[0].indices)
print("Centroid of the last cluster:\n", clusters[-1].centroid)

"""

::

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def score_from_files(filename, masks_dir, bundles_dir,
                     tracts_attribs, basic_bundles_attribs,
                     save_segmented=False, save_IBs=False,
                     save_VBs=False, save_VCWPs=False,
                     segmented_out_dir='', segmented_base_name='',
                     verbose=False):
    """
    Computes all metrics in order to score a tractogram.

    Given a ``tck`` file of streamlines and a folder containing masks,
    compute the percent of: Valid Connections (VC), Invalid Connections (IC),
    Valid Connections but Wrong Path (VCWP), No Connections (NC),
    Average Bundle Coverage (ABC), Average ROIs Coverage (ARC),
    coverage per bundles and coverage per ROIs. It also provides the number of:
    Valid Bundles (VB), Invalid Bundles (IB) and streamlines per bundles.


    Parameters
    ------------
    filename : str
       name of a tracts file
    masks_dir : str
       name of the directory containing the masks
    save_segmented : bool
        if true, saves the segmented VC, IC, VCWP and NC

    Returns
    ---------
    scores : dict
        dictionnary containing a score for each metric
    indices : dict
        dictionnary containing the indices of streamlines composing VC, IC,
        VCWP and NC

    """
    if verbose:
        logging.basicConfig(level=logging.DEBUG)

    rois_dir = masks_dir + "rois/"
    bundles_masks_dir = masks_dir + "bundles/"
    wm_file = masks_dir + "wm.nii.gz"
    wm = nib.load(wm_file)

    streamlines = load_streamlines(filename, wm_file, tracts_attribs)

    ROIs = [nib.load(rois_dir + f) for f in sorted(os.listdir(rois_dir))]
    bundles_masks = [nib.load(bundles_masks_dir + f) for f in sorted(os.listdir(bundles_masks_dir))]
    ref_bundles = []

    # Ref bundles will contain {'name': 'name_of_the_bundle', 'threshold': thres_value,
    #                           'streamlines': list_of_streamlines}
    dummy_attribs = {'orientation': 'LPS'}
    qb = QuickBundles(threshold=REF_BUNDLES_THRESHOLD, metric=AveragePointwiseEuclideanMetric())

    out_centroids_dir = os.path.join(segmented_out_dir, os.path.pardir, "centroids")
    if not os.path.isdir(out_centroids_dir):
        os.mkdir(out_centroids_dir)

    rng = np.random.RandomState(42)

    for bundle_idx, bundle_f in enumerate(sorted(os.listdir(bundles_dir))):
        bundle_attribs = basic_bundles_attribs.get(os.path.basename(bundle_f))
        if bundle_attribs is None:
            raise ValueError("Missing basic bundle attribs for {0}".format(bundle_f))

        # # Already resample to avoid doing it for each iteration of chunking
        # orig_strl = [s for s in get_tracts_voxel_space_for_dipy(
        #                         os.path.join(bundles_dir, bundle_f),
        #                         wm_file, dummy_attribs)]
        orig_strl = load_streamlines(os.path.join(bundles_dir, bundle_f), wm_file, dummy_attribs)
        resamp_bundle = set_number_of_points(orig_strl, NB_POINTS_RESAMPLE)
        # resamp_bundle = [s.astype('f4') for s in resamp_bundle]

        indices = np.arange(len(resamp_bundle))
        rng.shuffle(indices)
        bundle_cluster_map = qb.cluster(resamp_bundle, ordering=indices)

        # bundle_cluster_map.refdata = resamp_bundle

        bundle_mask_inv = nib.Nifti1Image((1 - bundles_masks[bundle_idx].get_data()) * wm.get_data(),
                                          bundles_masks[bundle_idx].get_affine())

        ref_bundles.append({'name': os.path.basename(bundle_f).replace('.fib', '').replace('.tck', ''),
                            'threshold': bundle_attribs['cluster_threshold'],
                            'cluster_map': bundle_cluster_map,
                            'mask': bundles_masks[bundle_idx],
                            'mask_inv': bundle_mask_inv})

        logging.debug("{}: {} centroids".format(ref_bundles[-1]['name'], len(bundle_cluster_map)))
        nib.streamlines.save(nib.streamlines.Tractogram(bundle_cluster_map.centroids, affine_to_rasmm=np.eye(4)),
                             os.path.join(out_centroids_dir, ref_bundles[-1]['name'] + ".tck"))

    score_func = score_auto_extract_auto_IBs

    return score_func(streamlines, bundles_masks, ref_bundles, ROIs, wm,
                      save_segmented=save_segmented, save_IBs=save_IBs,
                      save_VBs=save_VBs, save_VCWPs=save_VCWPs,
                      out_segmented_strl_dir=segmented_out_dir,
                      base_out_segmented_strl=segmented_base_name,
                      ref_anat_fname=wm_file)

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
from dipy.segment.metric import (AveragePointwiseEuclideanMetric,
                                 ResampleFeature)
from dipy.segment.clustering import QuickBundles

feature = ResampleFeature(nb_points=100)
metric = AveragePointwiseEuclideanMetric(feature)

"""
Since we are going to include all of the streamlines in the single cluster
from the streamlines, we set the threshold to `np.inf`. We pull out the
centroid as the standard.
"""

qb = QuickBundles(np.inf, metric=metric)

cluster_cst_l = qb.cluster(model_cst_l)
standard_cst_l = cluster_cst_l.centroids[0]

cluster_af_l = qb.cluster(model_af_l)
standard_af_l = cluster_af_l.centroids[0]

"""
We use the centroid streamline for each atlas bundle as the standard to orient
all of the streamlines in each bundle from the individual subject. Here, the
affine used is the one from the transform between the atlas and individual
tractogram. This is so that the orienting is done relative to the space of the
individual, and not relative to the atlas space.
"""

import dipy.tracking.streamline as dts

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def auto_extract_VCs(streamlines, ref_bundles):
    # Streamlines = list of all streamlines

    VC = 0
    VC_idx = set()

    found_vbs_info = {}
    for bundle in ref_bundles:
        found_vbs_info[bundle['name']] = {'nb_streamlines': 0,
                                          'streamlines_indices': set()}

    # Need to bookkeep because we chunk for big datasets
    processed_strl_count = 0
    chunk_size = 5000
    chunk_it = 0

    nb_bundles = len(ref_bundles)
    bundles_found = [False] * nb_bundles

    logging.debug("Starting scoring VCs")

    qb = QuickBundles(threshold=20, metric=AveragePointwiseEuclideanMetric())

    # Start loop here for big datasets
    while processed_strl_count < len(streamlines):
        logging.debug("Starting chunk: {0}".format(chunk_it))

        strl_chunk = streamlines[chunk_it * chunk_size:
                                 (chunk_it + 1) * chunk_size]

        processed_strl_count += len(strl_chunk)
        cur_chunk_VC_idx, cur_chunk_IC_idx, cur_chunk_VCWP_idx = set(), set(), set()

        # Already resample and run quickbundles on the submission chunk,
        # to avoid doing it at every call of auto_extract
        rstreamlines = set_number_of_points(strl_chunk, NB_POINTS_RESAMPLE)

        # qb.cluster had problem with f8
        rstreamlines = [s.astype('f4') for s in rstreamlines]

        chunk_cluster_map = qb.cluster(rstreamlines)
        chunk_cluster_map.refdata = strl_chunk

        logging.debug("Starting VC identification through auto_extract")

        for bundle_idx, ref_bundle in enumerate(ref_bundles):
            # The selected indices are from [0, len(strl_chunk)]
            selected_streamlines_indices = auto_extract(ref_bundle['cluster_map'],
                                                        chunk_cluster_map,
                                                        clean_thr=ref_bundle['threshold'])

            # Remove duplicates, when streamlines are assigned to multiple VBs.
            selected_streamlines_indices = set(selected_streamlines_indices) - \
                                           cur_chunk_VC_idx
            cur_chunk_VC_idx |= selected_streamlines_indices

            nb_selected_streamlines = len(selected_streamlines_indices)

            if nb_selected_streamlines:
                bundles_found[bundle_idx] = True
                VC += nb_selected_streamlines

                # Shift indices to match the real number of streamlines
                global_select_strl_indices = set([v + chunk_it * chunk_size
                                                 for v in selected_streamlines_indices])
                vb_info = found_vbs_info.get(ref_bundle['name'])
                vb_info['nb_streamlines'] += nb_selected_streamlines
                vb_info['streamlines_indices'] |= global_select_strl_indices

                VC_idx |= global_select_strl_indices
            else:
                global_select_strl_indices = set()

        chunk_it += 1

    # Compute bundle overlap, overreach and f1_scores and update found_vbs_info
    for bundle_idx, ref_bundle in enumerate(ref_bundles):
        bundle_name = ref_bundle["name"]
        bundle_mask = ref_bundle["mask"]

        vb_info = found_vbs_info[bundle_name]

        # Streamlines are in voxel space since that's how they were
        # loaded in the scoring function.
        tractogram = Tractogram(streamlines=(streamlines[i] for i in vb_info['streamlines_indices']),
                                affine_to_rasmm=bundle_mask.affine)

        scores = {}
        if len(tractogram) > 0:
            scores = compute_bundle_coverage_scores(tractogram, bundle_mask)

        vb_info['overlap'] = scores.get("OL", 0)
        vb_info['overreach'] = scores.get("OR", 0)
        vb_info['overreach_norm'] = scores.get("ORn", 0)
        vb_info['f1_score'] = scores.get("F1", 0)

    return VC_idx, found_vbs_info

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
"""
Fiber clustering
----------------

Based on an agglomerative clustering, and a geometric distance.
"""

clustering_outdir = os.path.join(outdir, "clustering")
cluster_file = os.path.join(clustering_outdir, "clusters.json")
if not os.path.isdir(clustering_outdir):
    os.mkdir(clustering_outdir)
if not os.path.isfile(cluster_file):
    fibers_18 = [resample(track, nb_pol=18) for track in fibers]
    qb = QuickBundles(threshold=10.)
    clusters_ = qb.cluster(fibers_18)
    clusters = {}
    for cnt, cluster in enumerate(clusters_):
        clusters[str(cnt)] = {"indices": cluster.indices}
    with open(cluster_file, "w") as open_file:
        json.dump(clusters, open_file, indent=4)
else:
    with open(cluster_file) as open_file:
        clusters = json.load(open_file)

if 1: #use_vtk:
    ren = pvtk.ren()
    colors = numpy.ones((len(fibers),))
    nb_clusters = len(clusters)
    for clusterid, item in clusters.items():
        indices = item["indices"]

# 需要导入模块: from dipy.segment.clustering import QuickBundles [as 别名]
# 或者: from dipy.segment.clustering.QuickBundles import cluster [as 别名]
def slr_with_qb(static, moving,
                x0='affine',
                rm_small_clusters=50,
                maxiter=100,
                select_random=None,
                verbose=False,
                greater_than=50,
                less_than=250,
                qb_thr=15,
                nb_pts=20,
                progressive=True, num_threads=None):
    """ Utility function for registering large tractograms.

    For efficiency we apply the registration on cluster centroids and remove
    small clusters.

    Parameters
    ----------
    static : Streamlines
    moving : Streamlines
    x0 : str
        rigid, similarity or affine transformation model (default affine)

    rm_small_clusters : int
        Remove clusters that have less than `rm_small_clusters` (default 50)

    verbose : bool,
        If True then information about the optimization is shown.

    select_random : int
        If not None select a random number of streamlines to apply clustering
        Default None.

    options : None or dict,
        Extra options to be used with the selected method.

    num_threads : int
        Number of threads. If None (default) then all available threads
        will be used. Only metrics using OpenMP will use this variable.

    Notes
    -----
    The order of operations is the following. First short or long streamlines
    are removed. Second the tractogram or a random selection of the tractogram
    is clustered with QuickBundles. Then SLR [Garyfallidis15]_ is applied.

    References
    ----------
    .. [Garyfallidis15] Garyfallidis et al. "Robust and efficient linear
            registration of white-matter fascicles in the space of streamlines"
            , NeuroImage, 117, 124--140, 2015
    .. [Garyfallidis14] Garyfallidis et al., "Direct native-space fiber
            bundle alignment for group comparisons", ISMRM, 2014.
    .. [Garyfallidis17] Garyfallidis et al. Recognition of white matter
            bundles using local and global streamline-based registration and
            clustering, Neuroimage, 2017.
    """
    if verbose:
        print('Static streamlines size {}'.format(len(static)))
        print('Moving streamlines size {}'.format(len(moving)))

    def check_range(streamline, gt=greater_than, lt=less_than):

        if (length(streamline) > gt) & (length(streamline) < lt):
            return True
        else:
            return False

    # TODO change this to the new Streamlines API
    streamlines1 = [s for s in static if check_range(s)]
    streamlines2 = [s for s in moving if check_range(s)]

    if verbose:

        print('Static streamlines after length reduction {}'
              .format(len(streamlines1)))
        print('Moving streamlines after length reduction {}'
              .format(len(streamlines2)))

    if select_random is not None:
        rstreamlines1 = select_random_set_of_streamlines(streamlines1,
                                                         select_random)
    else:
        rstreamlines1 = streamlines1

    rstreamlines1 = set_number_of_points(rstreamlines1, nb_pts)
    qb1 = QuickBundles(threshold=qb_thr)
    rstreamlines1 = [s.astype('f4') for s in rstreamlines1]
    cluster_map1 = qb1.cluster(rstreamlines1)
    clusters1 = remove_clusters_by_size(cluster_map1, rm_small_clusters)
    qb_centroids1 = [cluster.centroid for cluster in clusters1]

    if select_random is not None:
        rstreamlines2 = select_random_set_of_streamlines(streamlines2,
                                                         select_random)
    else:
        rstreamlines2 = streamlines2

    rstreamlines2 = set_number_of_points(rstreamlines2, nb_pts)
    qb2 = QuickBundles(threshold=qb_thr)
#.........这里部分代码省略.........