Merge pull request dipy#1443 from Garyfallidis/rb_rb

RecoBundles - recognition of bundles

dipy/align/bundlemin.pyx

@@ -132,8 +132,8 @@ def _bundle_minimum_distance_matrix(double [:, ::1] static,
     return np.asarray(D)
 
 
-def _bundle_minimum_distance(double [:, ::1] stat,
-                             double [:, ::1] mov,
+def _bundle_minimum_distance(double [:, ::1] static,
+                             double [:, ::1] moving,
                              cnp.npy_intp static_size,
                              cnp.npy_intp moving_size,
                              cnp.npy_intp rows,
@@ -207,8 +207,8 @@ def _bundle_minimum_distance(double [:, ::1] stat,
 
             for j in range(moving_size):
 
-                tmp = min_direct_flip_dist(&stat[i * rows, 0],
-                                       &mov[j * rows, 0], rows)
+                tmp = min_direct_flip_dist(&static[i * rows, 0],
+                                       &moving[j * rows, 0], rows)
 
                 if have_openmp:
                     openmp.omp_set_lock(&lock)
@@ -242,6 +242,98 @@ def _bundle_minimum_distance(double [:, ::1] stat,
     return dist
 
 
+
+def _bundle_minimum_distance_asymmetric(double [:, ::1] static,
+                                        double [:, ::1] moving,
+                                        cnp.npy_intp static_size,
+                                        cnp.npy_intp moving_size,
+                                        cnp.npy_intp rows):
+    """ MDF-based pairwise distance optimization function
+
+    We minimize the distance between moving streamlines of the same number of
+    points as they align with the static streamlines.
+
+    Parameters
+    -----------
+    static : array
+        Static streamlines
+    moving : array
+        Moving streamlines
+    static_size : int
+        Number of static streamlines
+    moving_size : int
+        Number of moving streamlines
+    rows : int
+        Number of points per streamline
+
+    Returns
+    -------
+    cost : double
+
+    Notes
+    -----
+    The difference with ``_bundle_minimum_distance`` is that we sum the
+    minimum values only for the static. Therefore, this is an asymetric
+    distance metric. This means that we are weighting only one direction of the
+    registration. Not both directions. This can be very useful when we want
+    to register a big set of bundles to a small set of bundles.
+    See [Wanyan17]_.
+
+    References
+    ----------
+    .. [Wanyan17] Wanyan and Garyfallidis, Important new insights for the
+    reduction of false positives in tractograms emerge from streamline-based
+    registration and pruning, International Society for Magnetic Resonance in
+    Medicine, Honolulu, Hawai, 2017.
+
+    """
+
+    cdef:
+        cnp.npy_intp i=0, j=0
+        double sum_i=0, sum_j=0, tmp=0
+        double inf = np.finfo('f8').max
+        double dist=0
+        double * min_j
+        openmp.omp_lock_t lock
+
+    with nogil:
+
+        if have_openmp:
+            openmp.omp_init_lock(&lock)
+
+        min_j = <double *> malloc(static_size * sizeof(double))
+
+        for i in range(static_size):
+            min_j[i] = inf
+
+        for i in prange(static_size):
+
+            for j in range(moving_size):
+
+                tmp = min_direct_flip_dist(&static[i * rows, 0],
+                                           &moving[j * rows, 0], rows)
+
+                if have_openmp:
+                    openmp.omp_set_lock(&lock)
+                if tmp < min_j[i]:
+                    min_j[i] = tmp
+
+                if have_openmp:
+                    openmp.omp_unset_lock(&lock)
+
+        if have_openmp:
+            openmp.omp_destroy_lock(&lock)
+
+        for i in range(static_size):
+            sum_i += min_j[i]
+
+        free(min_j)
+
+        dist = sum_i / <double>static_size
+
+    return dist
+
+
 def distance_matrix_mdf(streamlines_a, streamlines_b):
     r""" Minimum direct flipped distance matrix between two streamline sets