fixelcfestats: New option -mask

docs/reference/commands/fixelcfestats.rst

@@ -62,6 +62,8 @@ Additional options for fixelcfestats
 
 -  **-angle value** the max angle threshold for assigning streamline tangents to fixels (Default: 45 degrees)
 
+-  **-mask file** provide a fixel data file containing a mask of those fixels to be used during processing
+
 Standard options
 ^^^^^^^^^^^^^^^^
 
@@ -92,7 +94,7 @@ Raffelt, D.; Smith, RE.; Ridgway, GR.; Tournier, JD.; Vaughan, DN.; Rose, S.; He
 
 
 
-**Author:** David Raffelt (david.raffelt@florey.edu.au)
+**Author:** David Raffelt (david.raffelt@florey.edu.au) and Robert E. Smith (robert.smith@florey.edu.au)
 
 **Copyright:** Copyright (c) 2008-2017 the MRtrix3 contributors.
 

lib/mrtrix3/app.py

@@ -7,13 +7,11 @@
 # - cmdline.addCitation(), cmdline.addDescription(), cmdline.setCopyright() as needed
 # - Add arguments and options to 'cmdline' as needed
 # - parse()
-# - checkOutputFile() as needed
+# - checkOutputPath() as needed
 # - makeTempDir() if the script requires a temporary directory
 # - gotoTempDir() if the script is using a temporary directory
 # - complete()
 #
-# checkOutputFile() can be called at any time after parse() to check if an intended output
-#   file already exists
 
 
 

src/gui/mrview/tool/fixel/base_fixel.cpp

@@ -118,35 +118,37 @@ namespace MR
                "void main() {\n";
 
           if (fixel.use_discard_lower())
-            source += "  if (v_threshold[0] < lower) return;\n";
+            source += "  if (v_threshold[0] < lower || isnan(v_threshold[0])) return;\n";
           if (fixel.use_discard_upper())
-            source += "  if (v_threshold[0] > upper) return;\n";
+            source += "  if (v_threshold[0] > upper || isnan(v_threshold[0])) return;\n";
 
           switch (scale_type) {
             case Unity:
-              source += "   vec4 line_offset = length_mult * vec4 (v_dir[0], 0);\n";
+              source += "  vec4 line_offset = length_mult * vec4 (v_dir[0], 0);\n";
               break;
             case Value:
-              source += "   vec4 line_offset = length_mult * v_scale[0] * vec4 (v_dir[0], 0);\n";
+              source += "  if (isnan(v_scale[0])) return;\n"
+                        "  vec4 line_offset = length_mult * v_scale[0] * vec4 (v_dir[0], 0);\n";
               break;
           }
 
           switch (color_type) {
             case CValue:
               if (!ColourMap::maps[colourmap].special) {
-                source += "    float amplitude = clamp (";
+                source += "  if (isnan(v_colour[0])) return;\n"
+                          "  float amplitude = clamp (";
                 if (fixel.scale_inverted())
                   source += "1.0 -";
                 source += " scale * (v_colour[0] - offset), 0.0, 1.0);\n";
               }
               source +=
-                std::string ("    vec3 color;\n") +
+                std::string ("  vec3 color;\n") +
                 ColourMap::maps[colourmap].glsl_mapping +
-                "   fColour = color;\n";
+                "  fColour = color;\n";
               break;
             case Direction:
               source +=
-                "   fColour = normalize (abs (v_dir[0]));\n";
+                "  fColour = normalize (abs (v_dir[0]));\n";
               break;
             default:
               break;

src/stats/cfe.cpp

@@ -23,13 +23,15 @@ namespace MR
 
 
 
-      TrackProcessor::TrackProcessor (Image<uint32_t>& fixel_indexer,
+      TrackProcessor::TrackProcessor (Image<index_type>& fixel_indexer,
                                       const vector<direction_type>& fixel_directions,
+                                      Image<bool>& fixel_mask,
                                       vector<uint16_t>& fixel_TDI,
-                                      vector<std::map<uint32_t, connectivity> >& connectivity_matrix,
+                                      init_connectivity_matrix_type& connectivity_matrix,
                                       const value_type angular_threshold) :
                                         fixel_indexer        (fixel_indexer) ,
                                         fixel_directions     (fixel_directions),
+                                        fixel_mask           (fixel_mask),
                                         fixel_TDI            (fixel_TDI),
                                         connectivity_matrix  (connectivity_matrix),
                                         angular_threshold_dp (std::cos (angular_threshold * (Math::pi/180.0))) { }
@@ -39,26 +41,30 @@ namespace MR
       bool TrackProcessor::operator() (const SetVoxelDir& in)
       {
         // For each voxel tract tangent, assign to a fixel
-        vector<int32_t> tract_fixel_indices;
+        vector<index_type> tract_fixel_indices;
         for (SetVoxelDir::const_iterator i = in.begin(); i != in.end(); ++i) {
           assign_pos_of (*i).to (fixel_indexer);
           fixel_indexer.index(3) = 0;
-          uint32_t num_fibres = fixel_indexer.value();
-          if (num_fibres > 0) {
+          const index_type num_fixels = fixel_indexer.value();
+          if (num_fixels > 0) {
             fixel_indexer.index(3) = 1;
-            uint32_t first_index = fixel_indexer.value();
-            uint32_t last_index = first_index + num_fibres;
-            uint32_t closest_fixel_index = 0;
+            const index_type first_index = fixel_indexer.value();
+            const index_type last_index = first_index + num_fixels;
+            // Note: Streamlines can still be assigned to a fixel that is outside the mask;
+            //   however this will not be permitted to contribute to the matrix
+            index_type closest_fixel_index = num_fixels;
             value_type largest_dp = 0.0;
             const direction_type dir (i->get_dir().normalized());
-            for (uint32_t j = first_index; j < last_index; ++j) {
+            for (index_type j = first_index; j < last_index; ++j) {
               const value_type dp = std::abs (dir.dot (fixel_directions[j]));
               if (dp > largest_dp) {
                 largest_dp = dp;
-                closest_fixel_index = j;
+                fixel_mask.index(0) = j;
+                if (fixel_mask.value())
+                  closest_fixel_index = j;
               }
             }
-            if (largest_dp > angular_threshold_dp) {
+            if (closest_fixel_index != num_fixels && largest_dp > angular_threshold_dp) {
               tract_fixel_indices.push_back (closest_fixel_index);
               fixel_TDI[closest_fixel_index]++;
             }
@@ -86,11 +92,11 @@ namespace MR
 
 
 
-      Enhancer::Enhancer (const vector<std::map<uint32_t, connectivity> >& connectivity_map,
+      Enhancer::Enhancer (const norm_connectivity_matrix_type& connectivity_matrix,
                           const value_type dh,
                           const value_type E,
                           const value_type H) :
-          connectivity_map (connectivity_map),
+          connectivity_matrix (connectivity_matrix),
           dh (dh),
           E (E),
           H (H) { }
@@ -101,13 +107,13 @@ namespace MR
       {
         enhanced_stats = vector_type::Zero (stats.size());
         value_type max_enhanced_stat = 0.0;
-        for (size_t fixel = 0; fixel < connectivity_map.size(); ++fixel) {
-          std::map<uint32_t, connectivity>::const_iterator connected_fixel;
+        vector<NormMatrixElement>::const_iterator connected_fixel;
+        for (size_t fixel = 0; fixel < connectivity_matrix.size(); ++fixel) {
           for (value_type h = this->dh; h < stats[fixel]; h +=  this->dh) {
             value_type extent = 0.0;
-            for (connected_fixel = connectivity_map[fixel].begin(); connected_fixel != connectivity_map[fixel].end(); ++connected_fixel)
-              if (stats[connected_fixel->first] > h)
-                extent += connected_fixel->second.value;
+            for (connected_fixel = connectivity_matrix[fixel].begin(); connected_fixel != connectivity_matrix[fixel].end(); ++connected_fixel)
+              if (stats[connected_fixel->index()] > h)
+                extent += connected_fixel->value();
             enhanced_stats[fixel] += std::pow (extent, E) * std::pow (h, H);
           }
           if (enhanced_stats[fixel] > max_enhanced_stat)

src/stats/cfe.h

@@ -17,6 +17,7 @@
 
 #include "image.h"
 #include "image_helpers.h"
+#include "types.h"
 #include "math/math.h"
 #include "math/stats/typedefs.h"
 
@@ -30,11 +31,11 @@ namespace MR
     namespace CFE
     {
 
+      using index_type = uint32_t;
       using value_type = Math::Stats::value_type;
       using vector_type = Math::Stats::vector_type;
       using connectivity_value_type = float;
       using direction_type = Eigen::Matrix<value_type, 3, 1>;
-      using connectivity_vector_type = Eigen::Array<connectivity_value_type, Eigen::Dynamic, 1>;
       using SetVoxelDir = DWI::Tractography::Mapping::SetVoxelDir;
 
 
@@ -43,14 +44,38 @@ namespace MR
       @{ */
 
 
-      class connectivity { MEMALIGN(connectivity)
+      class connectivity { NOMEMALIGN
         public:
           connectivity () : value (0.0) { }
           connectivity (const connectivity_value_type v) : value (v) { }
           connectivity_value_type value;
       };
 
 
+      // A class to store fixel index / connectivity value pairs
+      //   only after the connectivity matrix has been thresholded / normalised
+      class NormMatrixElement
+      { NOMEMALIGN
+        public:
+          NormMatrixElement (const index_type fixel_index,
+                             const connectivity_value_type connectivity_value) :
+              fixel_index (fixel_index),
+              connectivity_value (connectivity_value) { }
+          FORCE_INLINE index_type index() const { return fixel_index; }
+          FORCE_INLINE connectivity_value_type value() const { return connectivity_value; }
+          FORCE_INLINE void normalise (const connectivity_value_type norm_factor) { connectivity_value *= norm_factor; }
+        private:
+          const index_type fixel_index;
+          connectivity_value_type connectivity_value;
+      };
+
+
+
+      // Different types are used depending on whether the connectivity matrix
+      //   is in the process of being built, or whether it has been normalised
+      using init_connectivity_matrix_type = vector<std::map<index_type, connectivity>>;
+      using norm_connectivity_matrix_type = vector<vector<NormMatrixElement>>;
+
 
 
       /**
@@ -59,19 +84,21 @@ namespace MR
       class TrackProcessor { MEMALIGN(TrackProcessor)
 
         public:
-          TrackProcessor (Image<uint32_t>& fixel_indexer,
+          TrackProcessor (Image<index_type>& fixel_indexer,
                           const vector<direction_type>& fixel_directions,
+                          Image<bool>& fixel_mask,
                           vector<uint16_t>& fixel_TDI,
-                          vector<std::map<uint32_t, connectivity> >& connectivity_matrix,
+                          init_connectivity_matrix_type& connectivity_matrix,
                           const value_type angular_threshold);
 
           bool operator () (const SetVoxelDir& in);
 
         private:
-          Image<uint32_t> fixel_indexer;
+          Image<index_type> fixel_indexer;
           const vector<direction_type>& fixel_directions;
+          Image<bool> fixel_mask;
           vector<uint16_t>& fixel_TDI;
-          vector<std::map<uint32_t, connectivity> >& connectivity_matrix;
+          init_connectivity_matrix_type& connectivity_matrix;
           const value_type angular_threshold_dp;
       };
 
@@ -80,15 +107,15 @@ namespace MR
 
       class Enhancer : public Stats::EnhancerBase { MEMALIGN (Enhancer)
         public:
-          Enhancer (const vector<std::map<uint32_t, connectivity> >& connectivity_map,
+          Enhancer (const norm_connectivity_matrix_type& connectivity_matrix,
                     const value_type dh, const value_type E, const value_type H);
 
 
           value_type operator() (const vector_type& stats, vector_type& enhanced_stats) const override;
 
 
         protected:
-          const vector<std::map<uint32_t, connectivity> >& connectivity_map;
+          const norm_connectivity_matrix_type& connectivity_matrix;
           const value_type dh, E, H;
       };