1import logging
2import abc
3import numpy as np
4from dipy.core.optimize import Optimizer
5from dipy.align.bundlemin import (_bundle_minimum_distance,
6                                  _bundle_minimum_distance_asymmetric,
7                                  distance_matrix_mdf)
8from dipy.tracking.streamline import (transform_streamlines,
9                                      unlist_streamlines,
10                                      center_streamlines,
11                                      set_number_of_points,
12                                      select_random_set_of_streamlines,
13                                      length,
14                                      Streamlines)
15from dipy.segment.clustering import qbx_and_merge
16from dipy.core.geometry import (compose_transformations,
17                                compose_matrix,
18                                decompose_matrix)
19from time import time
20
21DEFAULT_BOUNDS = [(-35, 35), (-35, 35), (-35, 35),
22                  (-45, 45), (-45, 45), (-45, 45),
23                  (0.6, 1.4), (0.6, 1.4), (0.6, 1.4),
24                  (-10, 10), (-10, 10), (-10, 10)]
25
26logger = logging.getLogger(__name__)
27
28
29class StreamlineDistanceMetric(object, metaclass=abc.ABCMeta):
30
31    def __init__(self, num_threads=None):
32        """ An abstract class for the metric used for streamline registration
33
34        If the two sets of streamlines match exactly then method ``distance``
35        of this object should be minimum.
36
37        Parameters
38        ----------
39        num_threads : int, optional
40            Number of threads to be used for OpenMP parallelization. If None
41            (default) the value of OMP_NUM_THREADS environment variable is used
42            if it is set, otherwise all available threads are used. If < 0 the
43            maximal number of threads minus |num_threads + 1| is used (enter -1
44            to use as many threads as possible). 0 raises an error. Only
45            metrics using OpenMP will use this variable.
46        """
47        self.static = None
48        self.moving = None
49        self.num_threads = num_threads
50
51    @abc.abstractmethod
52    def setup(self, static, moving):
53        pass
54
55    @abc.abstractmethod
56    def distance(self, xopt):
57        """ calculate distance for current set of parameters
58        """
59        pass
60
61
62class BundleMinDistanceMetric(StreamlineDistanceMetric):
63    """ Bundle-based Minimum Distance aka BMD
64
65    This is the cost function used by the StreamlineLinearRegistration
66
67    Methods
68    -------
69    setup(static, moving)
70    distance(xopt)
71
72    References
73    ----------
74    .. [Garyfallidis14] Garyfallidis et al., "Direct native-space fiber
75                        bundle alignment for group comparisons", ISMRM,
76                        2014.
77    """
78
79    def setup(self, static, moving):
80        """ Setup static and moving sets of streamlines
81
82        Parameters
83        ----------
84        static : streamlines
85            Fixed or reference set of streamlines.
86        moving : streamlines
87            Moving streamlines.
88
89        Notes
90        -----
91        Call this after the object is initiated and before distance.
92        """
93
94        self._set_static(static)
95        self._set_moving(moving)
96
97    def _set_static(self, static):
98        static_centered_pts, st_idx = unlist_streamlines(static)
99        self.static_centered_pts = np.ascontiguousarray(static_centered_pts,
100                                                        dtype=np.float64)
101        self.block_size = st_idx[0]
102
103    def _set_moving(self, moving):
104        self.moving_centered_pts, _ = unlist_streamlines(moving)
105
106    def distance(self, xopt):
107        """ Distance calculated from this Metric
108
109        Parameters
110        ----------
111        xopt : sequence
112            List of affine parameters as an 1D vector,
113
114        """
115        return bundle_min_distance_fast(xopt,
116                                        self.static_centered_pts,
117                                        self.moving_centered_pts,
118                                        self.block_size,
119                                        self.num_threads)
120
121
122class BundleMinDistanceMatrixMetric(StreamlineDistanceMetric):
123    """ Bundle-based Minimum Distance aka BMD
124
125    This is the cost function used by the StreamlineLinearRegistration
126
127    Methods
128    -------
129    setup(static, moving)
130    distance(xopt)
131
132    Notes
133    -----
134    The difference with BundleMinDistanceMetric is that this creates
135    the entire distance matrix and therefore requires more memory.
136
137    """
138
139    def setup(self, static, moving):
140        """ Setup static and moving sets of streamlines
141
142        Parameters
143        ----------
144        static : streamlines
145            Fixed or reference set of streamlines.
146        moving : streamlines
147            Moving streamlines.
148
149        Notes
150        -----
151        Call this after the object is initiated and before distance.
152
153        Num_threads is not used in this class. Use ``BundleMinDistanceMetric``
154        for a faster, threaded and less memory hungry metric
155        """
156        self.static = static
157        self.moving = moving
158
159    def distance(self, xopt):
160        """ Distance calculated from this Metric
161
162        Parameters
163        ----------
164        xopt : sequence
165            List of affine parameters as an 1D vector
166        """
167        return bundle_min_distance(xopt, self.static, self.moving)
168
169
170class BundleMinDistanceAsymmetricMetric(BundleMinDistanceMetric):
171    """ Asymmetric Bundle-based Minimum distance
172
173    This is a cost function that can be used by the
174    StreamlineLinearRegistration class.
175
176    """
177
178    def distance(self, xopt):
179        """ Distance calculated from this Metric
180
181        Parameters
182        ----------
183        xopt : sequence
184            List of affine parameters as an 1D vector
185        """
186        return bundle_min_distance_asymmetric_fast(xopt,
187                                                   self.static_centered_pts,
188                                                   self.moving_centered_pts,
189                                                   self.block_size)
190
191
192class BundleSumDistanceMatrixMetric(BundleMinDistanceMatrixMetric):
193    """ Bundle-based Sum Distance aka BMD
194
195    This is a cost function that can be used by the
196    StreamlineLinearRegistration class.
197
198    Methods
199    -------
200    setup(static, moving)
201    distance(xopt)
202
203    Notes
204    -----
205    The difference with BundleMinDistanceMatrixMetric is that it uses
206    uses the sum of the distance matrix and not the sum of mins.
207    """
208
209    def distance(self, xopt):
210        """ Distance calculated from this Metric
211
212        Parameters
213        ----------
214        xopt : sequence
215            List of affine parameters as an 1D vector
216        """
217        return bundle_sum_distance(xopt, self.static, self.moving)
218
219
220class StreamlineLinearRegistration(object):
221
222    def __init__(self, metric=None, x0="rigid", method='L-BFGS-B',
223                 bounds=None, verbose=False, options=None, evolution=False,
224                 num_threads=None):
225        r""" Linear registration of 2 sets of streamlines [Garyfallidis15]_.
226
227        Parameters
228        ----------
229        metric : StreamlineDistanceMetric,
230            If None and fast is False then the BMD distance is used. If fast
231            is True then a faster implementation of BMD is used. Otherwise,
232            use the given distance metric.
233
234        x0 : array or int or str
235            Initial parametrization for the optimization.
236
237            If 1D array with:
238                a) 6 elements then only rigid registration is performed with
239                the 3 first elements for translation and 3 for rotation.
240                b) 7 elements also isotropic scaling is performed (similarity).
241                c) 12 elements then translation, rotation (in degrees),
242                scaling and shearing is performed (affine).
243
244                Here is an example of x0 with 12 elements:
245                ``x0=np.array([0, 10, 0, 40, 0, 0, 2., 1.5, 1, 0.1, -0.5, 0])``
246
247                This has translation (0, 10, 0), rotation (40, 0, 0) in
248                degrees, scaling (2., 1.5, 1) and shearing (0.1, -0.5, 0).
249
250            If int:
251                a) 6
252                    ``x0 = np.array([0, 0, 0, 0, 0, 0])``
253                b) 7
254                    ``x0 = np.array([0, 0, 0, 0, 0, 0, 1.])``
255                c) 12
256                    ``x0 = np.array([0, 0, 0, 0, 0, 0, 1., 1., 1, 0, 0, 0])``
257
258            If str:
259                a) "rigid"
260                    ``x0 = np.array([0, 0, 0, 0, 0, 0])``
261                b) "similarity"
262                    ``x0 = np.array([0, 0, 0, 0, 0, 0, 1.])``
263                c) "affine"
264                    ``x0 = np.array([0, 0, 0, 0, 0, 0, 1., 1., 1, 0, 0, 0])``
265
266        method : str,
267            'L_BFGS_B' or 'Powell' optimizers can be used. Default is
268            'L_BFGS_B'.
269
270        bounds : list of tuples or None,
271            If method == 'L_BFGS_B' then we can use bounded optimization.
272            For example for the six parameters of rigid rotation we can set
273            the bounds = [(-30, 30), (-30, 30), (-30, 30),
274                          (-45, 45), (-45, 45), (-45, 45)]
275            That means that we have set the bounds for the three translations
276            and three rotation axes (in degrees).
277
278        verbose : bool, optional.
279            If True, if True then information about the optimization is shown.
280            Default: False.
281
282        options : None or dict,
283            Extra options to be used with the selected method.
284
285        evolution : boolean
286            If True save the transformation for each iteration of the
287            optimizer. Default is False. Supported only with Scipy >= 0.11.
288
289        num_threads : int, optional
290            Number of threads to be used for OpenMP parallelization. If None
291            (default) the value of OMP_NUM_THREADS environment variable is used
292            if it is set, otherwise all available threads are used. If < 0 the
293            maximal number of threads minus |num_threads + 1| is used (enter -1
294            to use as many threads as possible). 0 raises an error. Only
295            metrics using OpenMP will use this variable.
296
297        References
298        ----------
299        .. [Garyfallidis15] Garyfallidis et al. "Robust and efficient linear
300           registration of white-matter fascicles in the space of streamlines",
301           NeuroImage, 117, 124--140, 2015
302
303        .. [Garyfallidis14] Garyfallidis et al., "Direct native-space fiber
304           bundle alignment for group comparisons", ISMRM, 2014.
305
306        .. [Garyfallidis17] Garyfallidis et al. Recognition of white matter
307           bundles using local and global streamline-based
308           registration and clustering, Neuroimage, 2017.
309        """
310
311        self.x0 = self._set_x0(x0)
312        self.metric = metric
313
314        if self.metric is None:
315            self.metric = BundleMinDistanceMetric(num_threads=num_threads)
316
317        self.verbose = verbose
318        self.method = method
319        if self.method not in ['Powell', 'L-BFGS-B']:
320            raise ValueError('Only Powell and L-BFGS-B can be used')
321        self.bounds = bounds
322        self.options = options
323        self.evolution = evolution
324
325    def optimize(self, static, moving, mat=None):
326        """ Find the minimum of the provided metric.
327
328        Parameters
329        ----------
330        static : streamlines
331            Reference or fixed set of streamlines.
332        moving : streamlines
333            Moving set of streamlines.
334        mat : array
335            Transformation (4, 4) matrix to start the registration. ``mat``
336            is applied to moving. Default value None which means that initial
337            transformation will be generated by shifting the centers of moving
338            and static sets of streamlines to the origin.
339
340        Returns
341        -------
342        map : StreamlineRegistrationMap
343
344        """
345
346        msg = 'need to have the same number of points. Use '
347        msg += 'set_number_of_points from dipy.tracking.streamline'
348
349        if not np.all(np.array(list(map(len, static))) == static[0].shape[0]):
350            raise ValueError('Static streamlines ' + msg)
351
352        if not np.all(np.array(list(map(len, moving))) == moving[0].shape[0]):
353            raise ValueError('Moving streamlines ' + msg)
354
355        if not np.all(np.array(list(map(len, moving))) == static[0].shape[0]):
356            raise ValueError('Static and moving streamlines ' + msg)
357
358        if mat is None:
359            static_centered, static_shift = center_streamlines(static)
360            moving_centered, moving_shift = center_streamlines(moving)
361            static_mat = compose_matrix44([static_shift[0], static_shift[1],
362                                           static_shift[2], 0, 0, 0])
363
364            moving_mat = compose_matrix44([-moving_shift[0], -moving_shift[1],
365                                           -moving_shift[2], 0, 0, 0])
366        else:
367            static_centered = static
368            moving_centered = transform_streamlines(moving, mat)
369            static_mat = np.eye(4)
370            moving_mat = mat
371
372        self.metric.setup(static_centered, moving_centered)
373
374        distance = self.metric.distance
375
376        if self.method == 'Powell':
377
378            if self.options is None:
379                self.options = {'xtol': 1e-6, 'ftol': 1e-6, 'maxiter': 1e6}
380
381            opt = Optimizer(distance, self.x0.tolist(),
382                            method=self.method, options=self.options,
383                            evolution=self.evolution)
384
385        if self.method == 'L-BFGS-B':
386
387            if self.options is None:
388                self.options = {'maxcor': 10, 'ftol': 1e-7,
389                                'gtol': 1e-5, 'eps': 1e-8,
390                                'maxiter': 100}
391
392            opt = Optimizer(distance, self.x0.tolist(),
393                            method=self.method,
394                            bounds=self.bounds, options=self.options,
395                            evolution=self.evolution)
396        if self.verbose:
397            opt.print_summary()
398
399        opt_mat = compose_matrix44(opt.xopt)
400
401        mat = compose_transformations(moving_mat, opt_mat, static_mat)
402
403        mat_history = []
404
405        if opt.evolution is not None:
406            for vecs in opt.evolution:
407                mat_history.append(
408                    compose_transformations(moving_mat,
409                                            compose_matrix44(vecs),
410                                            static_mat))
411
412        srm = StreamlineRegistrationMap(mat, opt.xopt, opt.fopt,
413                                        mat_history, opt.nfev, opt.nit)
414        del opt
415        return srm
416
417    def _set_x0(self, x0):
418        """ check if input is of correct type"""
419
420        if hasattr(x0, 'ndim'):
421
422            if len(x0) not in [3, 6, 7, 9, 12]:
423                m_ = 'Only 1D arrays of 3, 6, 7, 9 and 12 elements are allowed'
424                raise ValueError(m_)
425            if x0.ndim != 1:
426                raise ValueError("Array should have only one dimension")
427            return x0
428
429        if isinstance(x0, str):
430
431            if x0.lower() == 'translation':
432                return np.zeros(3)
433
434            if x0.lower() == 'rigid':
435                return np.zeros(6)
436
437            if x0.lower() == 'similarity':
438                return np.array([0, 0, 0, 0, 0, 0, 1.])
439
440            if x0.lower() == 'scaling':
441                return np.array([0, 0, 0, 0, 0, 0, 1., 1., 1.])
442
443            if x0.lower() == 'affine':
444                return np.array([0, 0, 0, 0, 0, 0, 1., 1., 1., 0, 0, 0])
445
446        if isinstance(x0, int):
447            if x0 not in [3, 6, 7, 9, 12]:
448                msg = 'Only 3, 6, 7, 9 and 12 are accepted as integers'
449                raise ValueError(msg)
450            else:
451                if x0 == 3:
452                    return np.zeros(3)
453                if x0 == 6:
454                    return np.zeros(6)
455                if x0 == 7:
456                    return np.array([0, 0, 0, 0, 0, 0, 1.])
457                if x0 == 9:
458                    return np.array([0, 0, 0, 0, 0, 0, 1., 1., 1.])
459                if x0 == 12:
460                    return np.array([0, 0, 0, 0, 0, 0, 1., 1., 1., 0, 0, 0])
461
462        raise ValueError('Wrong input')
463
464
465class StreamlineRegistrationMap(object):
466
467    def __init__(self, matopt, xopt, fopt, matopt_history, funcs, iterations):
468        r""" A map holding the optimum affine matrix and some other parameters
469        of the optimization
470
471        Parameters
472        ----------
473
474        matrix : array,
475            4x4 affine matrix which transforms the moving to the static
476            streamlines
477
478        xopt : array,
479            1d array with the parameters of the transformation after centering
480
481        fopt : float,
482            final value of the metric
483
484        matrix_history : array
485            All transformation matrices created during the optimization
486
487        funcs : int,
488            Number of function evaluations of the optimizer
489
490        iterations : int
491            Number of iterations of the optimizer
492        """
493
494        self.matrix = matopt
495        self.xopt = xopt
496        self.fopt = fopt
497        self.matrix_history = matopt_history
498        self.funcs = funcs
499        self.iterations = iterations
500
501    def transform(self, moving):
502        """ Transform moving streamlines to the static.
503
504        Parameters
505        ----------
506        moving : streamlines
507
508        Returns
509        -------
510        moved : streamlines
511
512        Notes
513        -----
514
515        All this does is apply ``self.matrix`` to the input streamlines.
516        """
517
518        return transform_streamlines(moving, self.matrix)
519
520
521def bundle_sum_distance(t, static, moving, num_threads=None):
522    """ MDF distance optimization function (SUM)
523
524    We minimize the distance between moving streamlines as they align
525    with the static streamlines.
526
527    Parameters
528    -----------
529    t : ndarray
530        t is a vector of affine transformation parameters with
531        size at least 6.
532        If the size is 6, t is interpreted as translation + rotation.
533        If the size is 7, t is interpreted as translation + rotation +
534        isotropic scaling.
535        If size is 12, t is interpreted as translation + rotation +
536        scaling + shearing.
537
538    static : list
539        Static streamlines
540
541    moving : list
542        Moving streamlines. These will be transformed to align with
543        the static streamlines
544
545    num_threads : int, optional
546        Number of threads. If -1 then all available threads will be used.
547
548    Returns
549    -------
550    cost: float
551
552    """
553
554    aff = compose_matrix44(t)
555    moving = transform_streamlines(moving, aff)
556    d01 = distance_matrix_mdf(static, moving)
557    return np.sum(d01)
558
559
560def bundle_min_distance(t, static, moving):
561    """ MDF-based pairwise distance optimization function (MIN)
562
563    We minimize the distance between moving streamlines as they align
564    with the static streamlines.
565
566    Parameters
567    -----------
568    t : ndarray
569        t is a vector of affine transformation parameters with
570        size at least 6.
571        If size is 6, t is interpreted as translation + rotation.
572        If size is 7, t is interpreted as translation + rotation +
573        isotropic scaling.
574        If size is 12, t is interpreted as translation + rotation +
575        scaling + shearing.
576
577    static : list
578        Static streamlines
579
580    moving : list
581        Moving streamlines.
582
583    Returns
584    -------
585    cost: float
586
587    """
588    aff = compose_matrix44(t)
589    moving = transform_streamlines(moving, aff)
590    d01 = distance_matrix_mdf(static, moving)
591
592    rows, cols = d01.shape
593    return 0.25 * (np.sum(np.min(d01, axis=0)) / float(cols) +
594                   np.sum(np.min(d01, axis=1)) / float(rows)) ** 2
595
596
597def bundle_min_distance_fast(t, static, moving, block_size, num_threads=None):
598    """ MDF-based pairwise distance optimization function (MIN)
599
600    We minimize the distance between moving streamlines as they align
601    with the static streamlines.
602
603    Parameters
604    -----------
605    t : array
606        1D array. t is a vector of affine transformation parameters with
607        size at least 6.
608        If the size is 6, t is interpreted as translation + rotation.
609        If the size is 7, t is interpreted as translation + rotation +
610        isotropic scaling.
611        If size is 12, t is interpreted as translation + rotation +
612        scaling + shearing.
613
614    static : array
615        N*M x 3 array. All the points of the static streamlines. With order of
616        streamlines intact. Where N is the number of streamlines and M
617        is the number of points per streamline.
618
619    moving : array
620        K*M x 3 array. All the points of the moving streamlines. With order of
621        streamlines intact. Where K is the number of streamlines and M
622        is the number of points per streamline.
623
624    block_size : int
625        Number of points per streamline. All streamlines in static and moving
626        should have the same number of points M.
627
628    num_threads : int, optional
629        Number of threads to be used for OpenMP parallelization. If None
630        (default) the value of OMP_NUM_THREADS environment variable is used
631        if it is set, otherwise all available threads are used. If < 0 the
632        maximal number of threads minus |num_threads + 1| is used (enter -1 to
633        use as many threads as possible). 0 raises an error.
634
635    Returns
636    -------
637    cost: float
638
639    Notes
640    -----
641    This is a faster implementation of ``bundle_min_distance``, which requires
642    that all the points of each streamline are allocated into an ndarray
643    (of shape N*M by 3, with N the number of points per streamline and M the
644    number of streamlines). This can be done by calling
645    `dipy.tracking.streamlines.unlist_streamlines`.
646
647    """
648
649    aff = compose_matrix44(t)
650    moving = np.dot(aff[:3, :3], moving.T).T + aff[:3, 3]
651    moving = np.ascontiguousarray(moving, dtype=np.float64)
652
653    rows = static.shape[0] // block_size
654    cols = moving.shape[0] // block_size
655
656    return _bundle_minimum_distance(static, moving,
657                                    rows,
658                                    cols,
659                                    block_size,
660                                    num_threads)
661
662
663def bundle_min_distance_asymmetric_fast(t, static, moving, block_size):
664    """ MDF-based pairwise distance optimization function (MIN)
665
666    We minimize the distance between moving streamlines as they align
667    with the static streamlines.
668
669    Parameters
670    -----------
671    t : array
672        1D array. t is a vector of affine transformation parameters with
673        size at least 6.
674        If the size is 6, t is interpreted as translation + rotation.
675        If the size is 7, t is interpreted as translation + rotation +
676        isotropic scaling.
677        If size is 12, t is interpreted as translation + rotation +
678        scaling + shearing.
679
680    static : array
681        N*M x 3 array. All the points of the static streamlines. With order of
682        streamlines intact. Where N is the number of streamlines and M
683        is the number of points per streamline.
684
685    moving : array
686        K*M x 3 array. All the points of the moving streamlines. With order of
687        streamlines intact. Where K is the number of streamlines and M
688        is the number of points per streamline.
689
690    block_size : int
691        Number of points per streamline. All streamlines in static and moving
692        should have the same number of points M.
693
694    Returns
695    -------
696    cost: float
697    """
698
699    aff = compose_matrix44(t)
700    moving = np.dot(aff[:3, :3], moving.T).T + aff[:3, 3]
701    moving = np.ascontiguousarray(moving, dtype=np.float64)
702
703    rows = static.shape[0] // block_size
704    cols = moving.shape[0] // block_size
705
706    return _bundle_minimum_distance_asymmetric(static, moving,
707                                               rows,
708                                               cols,
709                                               block_size)
710
711
712def remove_clusters_by_size(clusters, min_size=0):
713    ob = filter(lambda c: len(c) >= min_size, clusters)
714
715    centroids = Streamlines()
716    for cluster in ob:
717        centroids.append(cluster.centroid)
718
719    return centroids
720
721
722def progressive_slr(static, moving, metric, x0, bounds, method='L-BFGS-B',
723                    verbose=False, num_threads=None):
724    """ Progressive SLR
725
726    This is an utility function that allows for example to do affine
727    registration using Streamline-based Linear Registration (SLR)
728    [Garyfallidis15]_ by starting with translation first, then rigid,
729    then similarity, scaling and finally affine.
730
731    Similarly, if for example, you want to perform rigid then you start with
732    translation first. This progressive strategy can helps with finding the
733    optimal parameters of the final transformation.
734
735    Parameters
736    ----------
737    static : Streamlines
738    moving : Streamlines
739    metric : StreamlineDistanceMetric
740    x0 : string
741        Could be any of 'translation', 'rigid', 'similarity', 'scaling',
742        'affine'
743    bounds : array
744        Boundaries of registration parameters. See variable `DEFAULT_BOUNDS`
745        for example.
746    method : string
747        L_BFGS_B' or 'Powell' optimizers can be used. Default is 'L_BFGS_B'.
748    verbose :  bool, optional.
749        If True, log messages. Default:
750    num_threads : int, optional
751        Number of threads to be used for OpenMP parallelization. If None
752        (default) the value of OMP_NUM_THREADS environment variable is used
753        if it is set, otherwise all available threads are used. If < 0 the
754        maximal number of threads minus |num_threads + 1| is used (enter -1 to
755        use as many threads as possible). 0 raises an error. Only metrics
756        using OpenMP will use this variable.
757
758    References
759    ----------
760    .. [Garyfallidis15] Garyfallidis et al. "Robust and efficient linear
761        registration of white-matter fascicles in the space of streamlines",
762        NeuroImage, 117, 124--140, 2015
763    """
764    if verbose:
765        logger.info('Progressive Registration is Enabled')
766
767    if x0 == 'translation' or x0 == 'rigid' or \
768       x0 == 'similarity' or x0 == 'scaling' or x0 == 'affine':
769        if verbose:
770            logger.info(' Translation  (3 parameters)...')
771        slr_t = StreamlineLinearRegistration(metric=metric,
772                                             x0='translation',
773                                             bounds=bounds[:3],
774                                             method=method)
775
776        slm_t = slr_t.optimize(static, moving)
777
778    if x0 == 'rigid' or x0 == 'similarity' or \
779       x0 == 'scaling' or x0 == 'affine':
780
781        x_translation = slm_t.xopt
782        x = np.zeros(6)
783        x[:3] = x_translation
784        if verbose:
785            logger.info(' Rigid  (6 parameters) ...')
786        slr_r = StreamlineLinearRegistration(metric=metric,
787                                             x0=x,
788                                             bounds=bounds[:6],
789                                             method=method)
790        slm_r = slr_r.optimize(static, moving)
791
792    if x0 == 'similarity' or x0 == 'scaling' or x0 == 'affine':
793
794        x_rigid = slm_r.xopt
795        x = np.zeros(7)
796        x[:6] = x_rigid
797        x[6] = 1.
798        if verbose:
799            logger.info(' Similarity (7 parameters) ...')
800        slr_s = StreamlineLinearRegistration(metric=metric,
801                                             x0=x,
802                                             bounds=bounds[:7],
803                                             method=method)
804        slm_s = slr_s.optimize(static, moving)
805
806    if x0 == 'scaling' or x0 == 'affine':
807
808        x_similarity = slm_s.xopt
809        x = np.zeros(9)
810        x[:6] = x_similarity[:6]
811        x[6:] = np.array((x_similarity[6],) * 3)
812        if verbose:
813            logger.info(' Scaling (9 parameters) ...')
814
815        slr_c = StreamlineLinearRegistration(metric=metric,
816                                             x0=x,
817                                             bounds=bounds[:9],
818                                             method=method)
819        slm_c = slr_c.optimize(static, moving)
820
821    if x0 == 'affine':
822
823        x_scaling = slm_c.xopt
824        x = np.zeros(12)
825        x[:9] = x_scaling[:9]
826        x[9:] = np.zeros(3)
827        if verbose:
828            logger.info(' Affine (12 parameters) ...')
829
830        slr_a = StreamlineLinearRegistration(metric=metric,
831                                             x0=x,
832                                             bounds=bounds[:12],
833                                             method=method)
834        slm_a = slr_a.optimize(static, moving)
835
836    if x0 == 'translation':
837        slm = slm_t
838    elif x0 == 'rigid':
839        slm = slm_r
840    elif x0 == 'similarity':
841        slm = slm_s
842    elif x0 == 'scaling':
843        slm = slm_c
844    elif x0 == 'affine':
845        slm = slm_a
846    else:
847        raise ValueError('Incorrect SLR transform')
848
849    return slm
850
851
852def slr_with_qbx(static, moving,
853                 x0='affine',
854                 rm_small_clusters=50,
855                 maxiter=100,
856                 select_random=None,
857                 verbose=False,
858                 greater_than=50,
859                 less_than=250,
860                 qbx_thr=[40, 30, 20, 15],
861                 nb_pts=20,
862                 progressive=True, rng=None, num_threads=None):
863    """ Utility function for registering large tractograms.
864
865    For efficiency, we apply the registration on cluster centroids and remove
866    small clusters.
867
868    Parameters
869    ----------
870    static : Streamlines
871    moving : Streamlines
872
873    x0 : str, optional.
874        rigid, similarity or affine transformation model (default affine)
875
876    rm_small_clusters : int, optional
877        Remove clusters that have less than `rm_small_clusters` (default 50)
878
879    select_random : int, optional.
880        If not, None selects a random number of streamlines to apply clustering
881        Default None.
882
883    verbose : bool, optional
884        If True, logs information about optimization. Default: False
885
886    greater_than : int, optional
887            Keep streamlines that have length greater than
888            this value (default 50)
889
890    less_than : int, optional
891            Keep streamlines have length less than this value (default 250)
892
893    qbx_thr : variable int
894            Thresholds for QuickBundlesX (default [40, 30, 20, 15])
895
896    np_pts : int, optional
897            Number of points for discretizing each streamline (default 20)
898
899    progressive : boolean, optional
900            (default True)
901
902    rng : RandomState
903        If None creates RandomState in function.
904
905    num_threads : int, optional
906        Number of threads to be used for OpenMP parallelization. If None
907        (default) the value of OMP_NUM_THREADS environment variable is used
908        if it is set, otherwise all available threads are used. If < 0 the
909        maximal number of threads minus |num_threads + 1| is used (enter -1 to
910        use as many threads as possible). 0 raises an error. Only metrics
911        using OpenMP will use this variable.
912
913    Notes
914    -----
915    The order of operations is the following. First short or long streamlines
916    are removed. Second, the tractogram or a random selection of the tractogram
917    is clustered with QuickBundles. Then SLR [Garyfallidis15]_ is applied.
918
919    References
920    ----------
921    .. [Garyfallidis15] Garyfallidis et al. "Robust and efficient linear
922    registration of white-matter fascicles in the space of streamlines",
923    NeuroImage, 117, 124--140, 2015
924    .. [Garyfallidis14] Garyfallidis et al., "Direct native-space fiber
925            bundle alignment for group comparisons", ISMRM, 2014.
926    .. [Garyfallidis17] Garyfallidis et al. Recognition of white matter
927    bundles using local and global streamline-based registration and
928    clustering, Neuroimage, 2017.
929    """
930    if rng is None:
931        rng = np.random.RandomState()
932
933    if verbose:
934        logger.info('Static streamlines size {}'.format(len(static)))
935        logger.info('Moving streamlines size {}'.format(len(moving)))
936
937    def check_range(streamline, gt=greater_than, lt=less_than):
938
939        if (length(streamline) > gt) & (length(streamline) < lt):
940            return True
941        else:
942            return False
943
944    streamlines1 = Streamlines(static[np.array([check_range(s)
945                                                for s in static])])
946    streamlines2 = Streamlines(moving[np.array([check_range(s)
947                                                for s in moving])])
948    if verbose:
949        logger.info('Static streamlines after length reduction {}'
950                    .format(len(streamlines1)))
951        logger.info('Moving streamlines after length reduction {}'
952                    .format(len(streamlines2)))
953
954    if select_random is not None:
955        rstreamlines1 = select_random_set_of_streamlines(streamlines1,
956                                                         select_random,
957                                                         rng=rng)
958    else:
959        rstreamlines1 = streamlines1
960
961    rstreamlines1 = set_number_of_points(rstreamlines1, nb_pts)
962
963    rstreamlines1._data.astype('f4')
964
965    cluster_map1 = qbx_and_merge(rstreamlines1, thresholds=qbx_thr, rng=rng)
966    qb_centroids1 = remove_clusters_by_size(cluster_map1, rm_small_clusters)
967
968    if select_random is not None:
969        rstreamlines2 = select_random_set_of_streamlines(streamlines2,
970                                                         select_random,
971                                                         rng=rng)
972    else:
973        rstreamlines2 = streamlines2
974
975    rstreamlines2 = set_number_of_points(rstreamlines2, nb_pts)
976    rstreamlines2._data.astype('f4')
977
978    cluster_map2 = qbx_and_merge(rstreamlines2, thresholds=qbx_thr, rng=rng)
979
980    qb_centroids2 = remove_clusters_by_size(cluster_map2, rm_small_clusters)
981
982    if verbose:
983        t = time()
984
985    if not progressive:
986        slr = StreamlineLinearRegistration(x0=x0,
987                                           options={'maxiter': maxiter},
988                                           num_threads=num_threads)
989        slm = slr.optimize(qb_centroids1, qb_centroids2)
990    else:
991        bounds = DEFAULT_BOUNDS
992
993        slm = progressive_slr(qb_centroids1, qb_centroids2,
994                              x0=x0, metric=None,
995                              bounds=bounds, num_threads=num_threads)
996
997    if verbose:
998        logger.info('QB static centroids size %d' % len(qb_centroids1,))
999        logger.info('QB moving centroids size %d' % len(qb_centroids2,))
1000        duration = time() - t
1001        logger.info('SLR finished in  %0.3f seconds.' % (duration,))
1002        if slm.iterations is not None:
1003            logger.info('SLR iterations: %d ' % (slm.iterations,))
1004
1005    moved = slm.transform(moving)
1006
1007    return moved, slm.matrix, qb_centroids1, qb_centroids2
1008
1009
1010# In essence whole_brain_slr can be thought as a combination of
1011# SLR on QuickBundles centroids and some thresholding see
1012# Garyfallidis et al. Recognition of white matter
1013# bundles using local and global streamline-based registration and
1014# clustering, NeuroImage, 2017.
1015whole_brain_slr = slr_with_qbx
1016
1017
1018def _threshold(x, th):
1019    return np.maximum(np.minimum(x, th), -th)
1020
1021
1022def compose_matrix44(t, dtype=np.double):
1023    """ Compose a 4x4 transformation matrix
1024
1025    Parameters
1026    -----------
1027    t : ndarray
1028        This is a 1D vector of affine transformation parameters with
1029        size at least 3.
1030        If the size is 3, t is interpreted as translation.
1031        If the size is 6, t is interpreted as translation + rotation.
1032        If the size is 7, t is interpreted as translation + rotation +
1033        isotropic scaling.
1034        If the size is 9, t is interpreted as translation + rotation +
1035        anisotropic scaling.
1036        If size is 12, t is interpreted as translation + rotation +
1037        scaling + shearing.
1038
1039    Returns
1040    -------
1041    T : ndarray
1042        Homogeneous transformation matrix of size 4x4.
1043
1044    """
1045    if isinstance(t, list):
1046        t = np.array(t)
1047    size = t.size
1048
1049    if size not in [3, 6, 7, 9, 12]:
1050        raise ValueError('Accepted number of parameters is 3, 6, 7, 9 and 12')
1051
1052    MAX_DIST = 1e10
1053    scale, shear, angles, translate = (None, ) * 4
1054    translate = _threshold(t[0:3], MAX_DIST)
1055    if size in [6, 7, 9, 12]:
1056        angles = np.deg2rad(t[3:6])
1057    if size == 7:
1058        scale = np.array((t[6],) * 3)
1059    if size in [9, 12]:
1060        scale = t[6: 9]
1061    if size == 12:
1062        shear = t[9: 12]
1063    return compose_matrix(scale=scale, shear=shear,
1064                          angles=angles,
1065                          translate=translate)
1066
1067
1068def decompose_matrix44(mat, size=12):
1069    """ Given a 4x4 homogeneous matrix return the parameter vector
1070
1071    Parameters
1072    -----------
1073    mat : array
1074        Homogeneous 4x4 transformation matrix
1075    size : int
1076        Size of the output vector. 3, for translation, 6 for rigid,
1077        7 for similarity, 9 for scaling and 12 for affine. Default is 12.
1078
1079    Returns
1080    -------
1081    t : ndarray
1082        One dimensional ndarray of 3, 6, 7, 9 or 12 affine parameters.
1083
1084    """
1085    scale, shear, angles, translate, _ = decompose_matrix(mat)
1086
1087    t = np.zeros(12)
1088    t[:3] = translate
1089    if size == 3:
1090        return t[:3]
1091    t[3: 6] = np.rad2deg(angles)
1092    if size == 6:
1093        return t[:6]
1094    if size == 7:
1095        t[6] = np.mean(scale)
1096        return t[:7]
1097    if size == 9:
1098        t[6:9] = scale
1099        return t[:9]
1100    if size == 12:
1101        t[6: 9] = scale
1102        t[9: 12] = shear
1103        return t
1104
1105    raise ValueError('Size can be 3, 6, 7, 9 or 12')
1106