estimate_response: simple script to get single-shell response function

cmd/sh2peaks.cpp

@@ -69,7 +69,11 @@ void usage ()
 
   + Option ("mask",
             "only perform computation within the specified binary brain mask image.")
-  + Argument ("image").type_image_in();
+  + Argument ("image").type_image_in()
+
+  + Option ("null", 
+            "specify value for non-peak directions (default: NaN).")
+  +   Argument ("value").type_float();
 }
 
 
@@ -80,7 +84,7 @@ typedef float value_type;
 class Direction
 {
   public:
-    Direction () : a (NAN) { }
+    Direction () : a (NaN) { }
     Direction (const Direction& d) : a (d.a), v (d.v) { }
     Direction (value_type phi, value_type theta) : a (1.0), v (std::cos (phi) *std::sin (theta), std::sin (phi) *std::sin (theta), std::cos (theta)) { }
     value_type a;
@@ -123,7 +127,7 @@ class DataLoader
           assign_pos_of(sh).to(*mask);
           if (!mask->value()) {
             for (auto l = Loop(3) (sh); l; ++l)
-              item.data[sh.index(3)] = NAN;
+              item.data[sh.index(3)] = NaN;
           }
         } else {
           // iterates over SH coefficients
@@ -156,15 +160,17 @@ class Processor
                int npeaks,
                std::vector<Direction> true_peaks,
                value_type threshold,
-               Image<value_type>* ipeaks_data) :
+               Image<value_type>* ipeaks_data, 
+               value_type value_if_not_peak = NaN) :
       dirs_vox (dirs_data),
       dirs (directions),
       lmax (lmax),
       npeaks (npeaks),
       true_peaks (true_peaks),
       threshold (threshold),
       peaks_out (npeaks),
-      ipeaks_vox (ipeaks_data) { }
+      ipeaks_vox (ipeaks_data),
+      value_if_not_peak (value_if_not_peak) { }
 
     bool operator() (const Item& item) {
 
@@ -174,7 +180,7 @@ class Processor
 
       if (check_input (item)) {
         for (auto l = Loop(3) (dirs_vox); l; ++l)
-          dirs_vox.value() = NAN;
+          dirs_vox.value() = NaN;
         return true;
       }
 
@@ -186,7 +192,7 @@ class Processor
         if (std::isfinite (p.a)) {
           for (size_t j = 0; j < all_peaks.size(); j++) {
             if (std::abs (p.v.dot (all_peaks[j].v)) > DOT_THRESHOLD) {
-              p.a = NAN;
+              p.a = NaN;
               break;
             }
           }
@@ -243,7 +249,8 @@ class Processor
         dirs_vox.value() = peaks_out[n].a*peaks_out[n].v[2];
         dirs_vox.index(3)++;
       }
-      for (; dirs_vox.index(3) < 3*npeaks; dirs_vox.index(3)++) dirs_vox.value() = NAN;
+      for (; dirs_vox.index(3) < 3*npeaks; dirs_vox.index(3)++) 
+        dirs_vox.value() = value_if_not_peak;
 
       return true;
     }
@@ -256,6 +263,7 @@ class Processor
     value_type threshold;
     std::vector<Direction> peaks_out;
     copy_ptr<Image<value_type> > ipeaks_vox;
+    const value_type value_if_not_peak;
 
     bool check_input (const Item& item) {
       if (ipeaks_vox) {
@@ -317,7 +325,9 @@ void run ()
   if (true_peaks.size()) 
     npeaks = true_peaks.size();
 
-  value_type threshold = get_option_value("threshold", -INFINITY);
+  value_type threshold = get_option_value ("threshold", -INFINITY);
+
+  value_type value_if_not_peak = get_option_value ("null", NaN);
 
   auto header = Header(SH_data);
   header.datatype() = DataType::Float32;
@@ -338,7 +348,7 @@ void run ()
 
   DataLoader loader (SH_data, mask_data.get());
   Processor processor (peaks, dirs, Math::SH::LforN (SH_data.size (3)),
-      npeaks, true_peaks, threshold, ipeaks_data.get());
+      npeaks, true_peaks, threshold, ipeaks_data.get(), value_if_not_peak);
 
   Thread::run_queue (loader, Thread::batch (Item()), Thread::multi (processor));
 }

scripts/__mrtrix_bash_helper_functions

@@ -0,0 +1,27 @@
+
+function get_config_parameter ()
+{
+  local value=$([ -f ~/.mrtrix.conf ] && sed -n -e 's/^[[:space:]]*'"$1"':[[:space:]]*//p' < ~/.mrtrix.conf)
+  [ "x$value" == "x" ] && [ -f /etc/mrtrix.conf ] && value=$(sed -n -e 's/^[[:space:]]*'"$1"':[[:space:]]*//p' < /etc/mrtrix.conf)
+  if [ "x$value" == "x" ]; then return 1; else echo "$value"; fi
+}
+
+function get_tmpdir ()
+{
+  local tmpdir=$(get_config_parameter TmpFileDir || 
+  ( [ "x$TMPDIR" != "x" ] && echo "$TMPDIR" ) ||
+  ( [ "x$TMP" != "x" ] && echo "$TMP" ) ||
+  ( [ "x$TEMP" != "x" ] && echo "$TEMP" ) ||
+  ( [ -d /tmp ] && echo /tmp ) ||
+  return 1 )
+  echo "$tmpdir"
+}
+
+
+function create_tmpdir ()
+{
+  tmpdir=$(get_tmpdir) || ( echo "could not find temporary folder location!"; return 1)
+  tmpdir=$(mktemp -d -p "$tmpdir" mrtrix.tmp-XXXXXX)
+  echo "$tmpdir"
+}
+

scripts/estimate_response

@@ -0,0 +1,136 @@
+#!/bin/bash
+
+if [ $# -lt 3 ]; then 
+  cat <<EOD
+
+ usage: estimate_response [OPTIONS] dwi mask response
+
+ where:
+
+    dwi         input DWI images
+    mask        input brain mask
+    response    output text file containing response function coefficients
+
+ OPTIONS:
+
+    -niter N        set the number of iterations to use (default: 3)
+    -shell b        set the shell to select from the DWI (default: max)
+    -lmax L         set the maximum harmonic order to use (default: 8).
+    -sf mask        output the final mask of single-fibre voxels 
+    -peaks image    output image of peak orientations
+
+EOD
+
+  exit 1
+fi
+
+niter=5
+
+while [[ $# > 3 ]]; do
+  case $1 in 
+    -niter)
+      niter="$2";
+      shift
+      ;;
+    -shell)
+      shell="-shell $2";
+      shift
+      ;;
+    -lmax)
+      lmax="-lmax $2";
+      shift
+      ;;
+    -sf)
+      sf_mask="$2";
+      shift
+      ;;
+    -peaks)
+      output_peak="$2";
+      shift
+      ;;
+    *)
+      ;;
+  esac
+  shift
+done
+
+
+
+# set fail on error:
+set -e
+
+# import helper functions:
+. "$(dirname "$(which "$0")")"/__mrtrix_bash_helper_functions
+
+tmpdir=$(create_tmpdir)
+echo 'NOTE: temporary files will be placed in "'$tmpdir'"'
+
+response="$3"
+dwi="$tmpdir"/dwi.mif
+sh="$tmpdir"/sh.mif
+fod="$tmpdir"/fod.mif
+peaks="$tmpdir"/peaks
+peak_ratio="$tmpdir"/peak_ratio.mif
+
+mask="$2"
+[ "x$sf_mask" == "x" ] && sf_mask="$tmpdir/sf_mask.mif"
+
+# start with sharpest possible response function (i.e. SH decomposition of a flat disc):
+echo "1,-1,1" > "$response"
+
+# preload DWI with optimal stride (should also allow -shell option in due course)
+dwiextract "$1" $shell "$dwi" 
+
+# compute SH fit to DWI:
+amp2sh "$dwi" $lmax "$sh"
+
+# set initial lmax=4:
+current_lmax="-lmax 4"
+
+for ((n=1; n<=niter; n++)); do
+
+  echo ''
+  echo "### iteration $n ###"
+
+  # estimate FOD within mask:
+  rm -f "$fod"
+  dwi2fod "$dwi" $current_lmax -mask "$mask" "$response" "$fod" 
+  # make sure lmax is set as per command-line in later iterations:
+  current_lmax="$lmax"
+
+  # extract 2 largest peaks into a 4D image split over volumes:
+  rm -f "${peaks}"-?.mif
+  sh2peaks "$fod" -mask "$mask" -num 2 -null 0 "${peaks}-[].mif" 
+
+  # compute metric of single-fibre-ness (admittedly not pretty, but clear enough):
+  # this computes: sqrt(peak1)*(1-peak2/peak1)^2 
+  mrmath "${peaks}-0.mif" "${peaks}-1.mif" "${peaks}-2.mif" rms "$tmpdir"/peak1.mif -force -quiet
+  mrmath "${peaks}-3.mif" "${peaks}-4.mif" "${peaks}-5.mif" rms "$tmpdir"/peak2.mif -force -quiet
+  mrcalc "$tmpdir"/peak1.mif -sqrt 1 "$tmpdir"/peak2.mif "$tmpdir"/peak1.mif -div -sub 2 -pow -mult "$peak_ratio" -force -quiet
+
+  # get mask of 300 voxels with largest metric:
+  rm -f "$sf_mask"
+  mrthreshold "$peak_ratio" -top 300 "$sf_mask" -quiet
+
+  # update estimate of response function within that mask, using FOD primary peak orientation:
+  rm -f "$response"
+  sh2response "$sh" "$sf_mask" "${peaks}-[0:2].mif" "$response" 
+
+
+  # refine mask to speed up later iterations:
+  if [ "$mask" == "$2" ]; then
+    mask="$tmpdir/mask.mif" 
+    mrthreshold "$peak_ratio" -top 3000 - | maskfilter - dilate "$mask"
+  fi
+
+
+  echo "current response coefficients: [ "$(cat "$response")" ]"
+
+done
+
+[ "x$output_peak" != "x" ] && mrconvert "${peaks}-[0:2].mif" "$output_peak" 
+
+# clean up:
+rm -rf "$tmpdir"
+
+