Quac rap

src/scilpy/cli/scil_tracking_local_dev.py

@@ -67,7 +67,7 @@
                              load_matrix_in_any_format)
 from scilpy.image.volume_space_management import DataVolume
 from scilpy.tracking.propagator import ODFPropagator
-from scilpy.tracking.rap import RAPContinue
+from scilpy.tracking.rap import RAPContinue, RAPGraph
 from scilpy.tracking.seed import SeedGenerator, CustomSeedsDispenser
 from scilpy.tracking.tracker import Tracker
 from scilpy.tracking.utils import (add_mandatory_options_tracking,
@@ -151,14 +151,27 @@ def _build_arg_parser():
                           "with -nt 1,000,000, \nyou can create tractogram_2 "
                           "with \n--skip 1,000,000.")
 
-    track_g.add_argument('--rap_mask', default=None,
+    rap_g = p.add_argument_group('RAP options')
+    rap_g.add_argument('--rap_mask', default=None,
                          help='Region-Adaptive Propagation mask (.nii.gz).\n'
-                              'Region-Adaptive Propagation tractography will start within '
-                              'this mask.')
-    track_g.add_argument('--rap_method', default='None',
-                         choices=['None', 'continue'],
-                         help="Region-Adaptive Propagation tractography method "
-                              " [%(default)s]")
+                        'Region-Adaptive Propagation tractography will start within '
+                        'this mask.')
+    rap_g.add_argument('--rap_method', default='None',
+                        choices=['None', 'continue', 'quack'],
+                        help="Region-Adaptive Propagation tractography method.\n" 
+                        "To use option quack, you must install Quacktography \n" 
+                        "continue : goes straight in the RAP mask, \n" 
+                        "quack : uses a graph solution by quantum approach. " 
+                        " [%(default)s]")
+    rap_g.add_argument('--reps', type=int, default=2,
+                         help='Number of repetitions for the RAP method. \n'
+                             'Default: 2. This is only used for the quack '
+                             'method, not for continue.')
+    rap_g.add_argument('--alpha', type=float, default=1.5,
+                         help='Alpha parameter for the RAP method. \n'
+                              'Default: 1.5. This is only used for the quack '
+                              'method, not for continue.')
+
 
     m_g = p.add_argument_group('Memory options')
     add_processes_arg(m_g)
@@ -292,6 +305,10 @@ def main():
     if args.rap_method == "continue":
         rap = RAPContinue(rap_mask, propagator, max_nbr_pts,
                           step_size=vox_step_size)
+    elif args.rap_method == "quack":
+        rap = RAPGraph(rap_mask, rap_img, propagator, max_nbr_pts, fodf=odf_sh_img, reps=args.reps,
+                       alpha=args.alpha)
+
     else:
         rap = None
 

src/scilpy/tracking/propagator.py

@@ -402,7 +402,7 @@ def __init__(self, datavolume, step_size,
             get_sh_order_and_fullness(self.datavolume.nb_coeffs)
         self.basis = basis
         self.is_legacy = is_legacy
-        self.B = sh_to_sf_matrix(self.sphere, sh_order, self.basis,
+        self.B = sh_to_sf_matrix(self.sphere, sh_order=sh_order, basis_type=self.basis,
                                  smooth=0.006, return_inv=False,
                                  full_basis=full_basis, legacy=self.is_legacy)
 

src/scilpy/tracking/rap.py

@@ -1,6 +1,8 @@
 # -*- coding: utf-8 -*-
 
 import numpy as np
+from scilpy.tracking.utils import TrackingDirection
+
 
 
 class RAP:
@@ -61,10 +63,44 @@ def rap_multistep_propagate(self, line, prev_direction):
 
 
 class RAPGraph(RAP):
-    def __init__(self, mask_rap, propagator, max_nbr_pts, neighboorhood_size):
+    def __init__(self, mask_rap, rap_img, propagator, max_nbr_pts, fodf, reps, alpha):
+        """
+        RAPGraph class for the quantum Graph solution for a region.
+
+        Parameters
+        ----------
+        fodf: DataVolume
+            The FODF volume used to compute the RAP.
+        reps: int
+            Number of repetitions used in the quantum circuit.
+        alpha: float
+            Initial paramater to search the cost landscape.
+        """
         super().__init__(mask_rap, propagator, max_nbr_pts)
-        self.neighboorhood_size = neighboorhood_size
+
+
+        self.fodf = fodf
+        self.rap_img = rap_img
+        self.reps = reps
+        self.alpha = alpha
 
 
     def rap_multistep_propagate(self, line, prev_direction):
-        raise NotImplementedError
+        try:
+            from quactography.solver.rap_tracking import quack_rap
+        except ImportError:
+            raise ImportError("quactography is not installed. "
+                              "Please install it to use RAPGraph.\n"
+                              "Add: Follow instructions here: https://github.com/scilus/quactography")
+
+        prev_dir = np.array(prev_direction)
+        seg, prev_dir, is_line_valid = quack_rap(self.rap_img, self.fodf, line[-1].round().astype(int),
+                                                  reps = self.reps,
+                                                  alpha = self.alpha, 
+                                                  prev_direction = prev_dir,
+                                                  theta = self.propagator.theta,
+                                                  threshold = self.propagator.sf_threshold)
+        line.extend(seg)
+        last_dir = TrackingDirection(prev_direction)
+
+        return line, last_dir, is_line_valid