tckgen PTT: Differentiate probe length from step size

Co-authored-by: Dogu Baran Aydogan <baran.aydogan@gmail.com>

cmd/tckgen.cpp

@@ -226,6 +226,7 @@ void usage ()
 
   + DWI::Tractography::Algorithms::FODOptions
   + DWI::Tractography::Algorithms::SecondOrderOptions
+  + DWI::Tractography::Algorithms::PTTOptions
 
   + DWI::GradImportOptions();
 
@@ -258,7 +259,8 @@ void run ()
     Algorithms::load_FOD_options (properties);
   if (algorithm == 2 || algorithm == 3)
     Algorithms::load_2ndorder_options (properties);
-
+  if (algorithm == 3)
+    Algorithms::load_PTT_options (properties);
 
   //load ROIs and tractography specific options
   //NB must occur before seed check below due to -select option override

src/dwi/tractography/algorithms/ptt.cpp

@@ -0,0 +1,55 @@
+/* Copyright (c) 2008-2023 the MRtrix3 contributors.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ * Covered Software is provided under this License on an "as is"
+ * basis, without warranty of any kind, either expressed, implied, or
+ * statutory, including, without limitation, warranties that the
+ * Covered Software is free of defects, merchantable, fit for a
+ * particular purpose or non-infringing.
+ * See the Mozilla Public License v. 2.0 for more details.
+ *
+ * For more details, see http://www.mrtrix.org/.
+ */
+
+#include "dwi/tractography/algorithms/ptt.h"
+
+
+namespace MR
+{
+  namespace DWI
+  {
+    namespace Tractography
+    {
+      namespace Algorithms
+      {
+
+
+
+        using namespace App;
+
+
+
+        const OptionGroup PTTOptions = OptionGroup ("Options specific to the Parallel Transport Tractography (PTT) algorithm")
+        + Option ("probe_length", "length of probe used to sample FOD amplitudes; not necessarily equal to step size (default: " + str(Defaults::probelength_voxels_ptt) + " x voxel size)")
+          + Argument ("value").type_float (0.0);
+
+
+
+        void load_PTT_options (Tractography::Properties& properties)
+        {
+          auto opt = get_options ("probe_length");
+          if (opt.size())
+            properties["probe_length"] = str<float> (opt[0][0]);
+        }
+
+
+
+      }
+    }
+  }
+}
+
+

src/dwi/tractography/algorithms/ptt.h

@@ -36,6 +36,8 @@ namespace MR
 
 
         using namespace MR::DWI::Tractography::Tracking;
+        extern const App::OptionGroup PTTOptions;
+        void load_PTT_options (Tractography::Properties&);
 
 
 
@@ -53,7 +55,8 @@ namespace MR
                 max_trials_calibration_seeding (1000),
                 max_trials_sampling (Defaults::max_trials_per_step),
                 kmax (2.0f / vox()),
-                nsamples (Defaults::secondorder_nsamples)
+                nsamples (Defaults::secondorder_nsamples),
+                probe_length (Defaults::probelength_voxels_ptt * vox())
             {
               try {
                 Math::SH::check (source);
@@ -78,7 +81,7 @@ namespace MR
               // TODO In addition this needs to properly deal with potential differences between step size and probe length
 
               set_step_and_angle (Defaults::stepsize_voxels_ptt, Defaults::angle_ptt, true);
-
+              properties.set (probe_length, "probe_length");
               set_num_points();
 
               set_cutoff (Defaults::cutoff_fod * (is_act() ? Defaults::cutoff_act_multiplier : 1.0));
@@ -95,9 +98,9 @@ namespace MR
                 precomputer = new Math::SH::PrecomputedAL<float>(lmax);
 
               properties.set (nsamples, "samples_per_step");
-              ts.reserve (nsamples);
+              probe_t.reserve (nsamples);
               for (size_t i = 0; i != nsamples; ++i)
-                ts.push_back (step_size * float(i) / float(nsamples));
+                probe_t.push_back (probe_length * float(i) / float(nsamples-1));
 
             }
 
@@ -114,7 +117,8 @@ namespace MR
             size_t max_trials_calibration_tracking, max_trials_calibration_seeding, max_trials_sampling;
             float kmax;
             size_t nsamples;
-            vector<float> ts;
+            float probe_length;
+            vector<float> probe_t;
             Math::SH::PrecomputedAL<float> *precomputer;
           };
 
@@ -176,7 +180,8 @@ namespace MR
             if (!std::isfinite (sample.f))
               return MODEL;
 
-            F = probe.back();
+            const Eigen::Matrix<float, 4, 4> propagator = P (sample.k1, sample.k2, S.step_size);
+            F = propagator * F;
             pos = F.x();
             dir = F.T();
             return CONTINUE;
@@ -308,7 +313,7 @@ namespace MR
             probe[0] = F;
             for (size_t v = 1; v != S.nsamples; ++v)
             {
-              propagator = P (k1, k2, S.ts[v]);
+              propagator = P (k1, k2, S.probe_t[v]);
               probe[v].noalias() = propagator * F;
             }
           }

src/dwi/tractography/tracking/tractography.h

@@ -48,7 +48,7 @@ namespace MR
 
           constexpr size_t max_trials_per_step = 1000;
 
-
+          constexpr float probelength_voxels_ptt = 0.25f;
           constexpr float stepsize_voxels_firstorder = 0.1f;
           constexpr float stepsize_voxels_ifod2 = 0.5f;
           constexpr float stepsize_voxels_ptt = 0.025f;