Interp::Masked: Change implementation

In the original implementation, scanning through the dataset in order to produce an explicit binary voxel mask resulted in a noticeable delay in the commencement of tractography. With this change, the data in the voxel in which the point resides are instead checked for zero-filling at the point of the voxel() call.

core/interp/base.h

@@ -203,6 +203,11 @@ namespace MR
                                     pos[2] <= -0.5 || pos[2] >= bounds[2])));
         }
 
+        template <>
+        bool set_out_of_bounds (const bool& value) {
+          return ((out_of_bounds = value));
+        }
+
         template <class VectorType>
         Eigen::Vector3 intravoxel_offset (const VectorType& pos) {
           if (set_out_of_bounds (pos))

core/interp/masked.h

@@ -18,6 +18,7 @@
 #define __interp_masked_h__
 
 #include "datatype.h"
+#include "algo/loop.h"
 #include "interp/base.h"
 
 namespace MR
@@ -36,58 +37,38 @@ namespace MR
      *
      * (NaN values are permitted in order to be compatible with 3-vectors
      * images, i.e. sets of XYZ triplets; but there needs to be at least
-     * one non-NaN value in the voxel)
+     * one non-NaN & non-zero value in the voxel)
      */
-
     template <class InterpType>
     class Masked : public InterpType { MEMALIGN(Masked<InterpType>)
       public:
         using typename InterpType::value_type;
 
         Masked (const typename InterpType::image_type& parent,
                 const value_type value_when_out_of_bounds = Base<typename InterpType::image_type>::default_out_of_bounds_value()) :
-            InterpType (parent, value_when_out_of_bounds),
-            mask (Image<bool>::scratch (make_mask_header (parent), "scratch binary mask for masked interpolator"))
-        {
-          typename InterpType::image_type temp (parent);
-          compute_mask (temp);
-        }
-
-        //! Alternative constructor that uses a separate image to construct the mask
-        /* If the cost of accessing the underlying image is high (e.g. bootstrapping),
-         * this constructor enables use of the original plain data in the sampling of
-         * image values to construct the implicit mask */
-        template <class ImageType>
-        Masked (const typename InterpType::image_type& parent,
-                const ImageType& masking_data,
-                const value_type value_when_out_of_bounds = Base<typename InterpType::image_type>::default_out_of_bounds_value()) :
-            InterpType (parent, value_when_out_of_bounds),
-            mask (Image<bool>::scratch (make_mask_header (masking_data), "scratch binary mask for masked interpolator"))
-        {
-          assert (dimensions_match (parent, masking_data, 0, 3));
-          ImageType temp (masking_data);
-          if (masking_data.ndim() == 3)
-            copy (temp, mask);
-          else
-            compute_mask (temp);
-        }
+            InterpType (parent, value_when_out_of_bounds) { }
 
         //! Set the current position to <b>voxel space</b> position \a pos
-        /* Unlike other interpolators, this involves checking the value of
-         * an internally stored mask, which tracks those voxels that contain
-         * at least one finite & non-zero value.
+        /* Unlike other interpolators, this sets the .index() location of
+         * the parent image, and checks to see whether or not there is
+         * any finite & non-zero data present; if there is not, then the
+         * function returns false, as though the location is outside of the
+         * image FoV.
          *
          * See file interp/base.h for details. */
         template <class VectorType>
-        bool voxel (const VectorType& pos) {
+        bool voxel (const VectorType& pos)
+        {
           if (InterpType::set_out_of_bounds (pos))
             return false;
-          mask.index(0) = std::round (pos[0]);
-          mask.index(1) = std::round (pos[1]);
-          mask.index(2) = std::round (pos[2]);
-          if (!mask.value())
-            return false;
-          return InterpType::voxel (pos);
+          InterpType::image_type::index(0) = std::round (pos[0]);
+          InterpType::image_type::index(1) = std::round (pos[1]);
+          InterpType::image_type::index(2) = std::round (pos[2]);
+          for (auto l_inner = Loop(*this, 3) (*this); l_inner; ++l_inner) {
+            if (InterpType::image_type::value())
+              return InterpType::voxel (pos);
+          }
+          return InterpType::set_out_of_bounds (true);
         }
 
         //! Set the current position to <b>image space</b> position \a pos
@@ -104,35 +85,6 @@ namespace MR
           return voxel (Transform::scanner2voxel * pos.template cast<default_type>());
         }
 
-      protected:
-        Image<bool> mask;
-
-        template <class ImageType>
-        static Header make_mask_header (const ImageType& image)
-        {
-          Header H (image);
-          H.ndim() = 3;
-          H.datatype() = DataType::Bit;
-          return H;
-        }
-
-        template <class ImageType>
-        void compute_mask (ImageType& data)
-        {
-          // True if at least one finite & non-zero value
-          for (auto l_voxel = Loop("pre-computing implicit voxel mask for \"" + data.name() + "\"", mask) (mask); l_voxel; ++l_voxel) {
-            assign_pos_of (mask, 0, 3).to (data);
-            bool value = false;
-            for (auto l_inner = Loop(data, 3) (data); l_inner; ++l_inner) {
-              if (data.value()) {
-                value = true;
-                continue;
-              }
-            }
-            mask.value() = value;
-          }
-        }
-
     };
 
 

src/dwi/tractography/algorithms/fact.h

@@ -90,13 +90,10 @@ namespace MR
       bool init() override
       {
         if (!get_data (source)) return false;
-        if (!S.init_dir.allFinite()) {
-          if (!dir.allFinite())
-            dir = random_direction();
-        }
-        else
+        if (S.init_dir.allFinite())
           dir = S.init_dir;
-
+        else
+          dir = random_direction();
         return select_fixel (dir) >= S.threshold;
       }
 

src/dwi/tractography/algorithms/tensor_prob.h

@@ -120,7 +120,7 @@ namespace MR
               public:
                 Interp (const Bootstrap<Image<float>,WildBootstrap>& bootstrap_vox,
                         const Image<float>& vox) :
-                    Interpolator<Bootstrap<Image<float>,WildBootstrap> >::type (bootstrap_vox, vox)
+                    Interpolator<Bootstrap<Image<float>,WildBootstrap> >::type (bootstrap_vox)
                 {
                   for (size_t i = 0; i < 8; ++i)
                     raw_signals.push_back (Eigen::VectorXf (size(3)));
@@ -129,15 +129,13 @@ namespace MR
                 vector<Eigen::VectorXf> raw_signals;
 
                 bool get (const Eigen::Vector3f& pos, Eigen::VectorXf& data) {
-                  scanner (pos);
-                  if (out_of_bounds) {
+                  if (!scanner (pos)) {
                     data.fill (NaN);
                     return false;
                   }
 
                   data.setZero();
 
-                  // Modified to be consistent with the new Interp::Linear implementation
                   size_t i = 0;
                   for (ssize_t z = 0; z < 2; ++z) {
                     index(2) = clamp (P[2]+z, size(2));

testing/binaries/tests/tckgen

@@ -11,4 +11,4 @@ tckgen SIFT_phantom/fods.mif -algo sd_stream -seed_image SIFT_phantom/mask.mif -
 tckgen SIFT_phantom/dwi.mif -algo tensor_prob -seed_image SIFT_phantom/mask.mif -mask SIFT_phantom/mask.mif -minlength 4 -select 100 tmp.tck -force
 tckgen SIFT_phantom/fods.mif -algo ifod1 -seed_image SIFT_phantom/mask.mif -act SIFT_phantom/5tt.mif -backtrack -select 100 tmp.tck -force
 tckgen dwi.mif -algo tensor_det -seed_grid_per_voxel mrcrop/mask.mif 3 -nthread 0 tmp.tck -force && testing_diff_tck tmp.tck tckgen/tensor_det.tck
-tckgen dwi.mif -algo tensor_det -seed_grid_per_voxel mrcrop/mask.mif 3 tmp.tck -force && testing_diff_tck tmp.tck tckgen/tensor_det.tck -unordered
+tckgen dwi.mif -algo tensor_det -seed_grid_per_voxel mrcrop/mask.mif 3 tmp.tck -force && testing_diff_tck tmp.tck tckgen/tensor_det.tck -unordered && testing_diff_tck tckgen/tensor_det.tck tmp.tck -unordered