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,
@@ -156,9 +156,20 @@ def _build_arg_parser():
                         '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 "
+                        choices=['None', 'continue', 'quack'],
+                        help="Region-Adaptive Propagation tractography method." \
+                        "To use option quack, you must install Quacktography " \
+                        "continue : goes straight in the RAP mask, " \
+                        "quack : uses a graph solution by quantum approach. " \
                         " [%(default)s]")
+    track_g.add_argument('--reps', type=int, default=2,
+                         help='Number of repetitions for the RAP method. '
+                             'Default: 2. This is only used for the quack '
+                             'method, not for continue.')
+    track_g.add_argument('--alpha', type=float, default=1.5,
+                         help='Alpha parameter for the RAP method. '
+                              '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 +303,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/rap.py

@@ -61,10 +61,41 @@ 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.")
+
+        prev_dir = np.array(prev_direction).round().astype(int)
+        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)
+        return line, prev_dir, is_line_valid