gradient.h

Go to the documentation of this file.

/* Copyright (c) 2008-2022 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/.
 */
 
#ifndef __dwi_gradient_h__
#define __dwi_gradient_h__
 
#include "app.h"
#include "header.h"
#include "dwi/shells.h"
#include "file/config.h"
#include "file/path.h"
#include "math/condition_number.h"
#include "math/sphere.h"
#include "math/SH.h"
 
 
namespace MR
{
  namespace App { class OptionGroup; }
 
  namespace DWI
  {
 
    App::OptionGroup GradImportOptions();
    App::OptionGroup GradExportOptions();
 
    extern App::Option bvalue_scaling_option;
    extern const char* const bvalue_scaling_description;
 
 
 
 
    template <class MatrixType>
      inline void check_DW_scheme (const Header& header, const MatrixType& grad)
      {
        if (!grad.rows())
          throw Exception ("no valid diffusion gradient table found");
        if (grad.cols() < 4)
          throw Exception ("unexpected diffusion gradient table matrix dimensions");
 
        if (header.ndim() >= 4) {
          if (header.size (3) != (int) grad.rows())
            throw Exception ("number of studies in base image (" + str(header.size(3))
                + ") does not match number of rows in diffusion gradient table (" + str(grad.rows()) + ")");
        }
        else if (grad.rows() != 1)
          throw Exception ("For images with less than four dimensions, gradient table can have one row only");
      }
 
 
 
 
 
 
    template <class MatrixType, class IndexVectorType>
      inline Eigen::MatrixXd gen_direction_matrix (
          const MatrixType& grad,
          const IndexVectorType& dwi)
      {
        Eigen::MatrixXd dirs (dwi.size(),2);
        for (size_t i = 0; i < dwi.size(); i++) {
          dirs (i,0) = std::atan2 (grad (dwi[i],1), grad (dwi[i],0));
          auto z = grad (dwi[i],2) / grad.row (dwi[i]).template head<3>().norm();
          if (z >= 1.0)
            dirs(i,1) = 0.0;
          else if (z <= -1.0)
            dirs (i,1) = Math::pi;
          else
            dirs (i,1) = std::acos (z);
        }
        return dirs;
      }
 
 
 
 
 
    template <class MatrixType>
    default_type condition_number_for_lmax (const MatrixType& dirs, int lmax)
    {
      Eigen::MatrixXd g;
      if (dirs.cols() == 2) // spherical coordinates:
        g = dirs;
      else // Cartesian to spherical:
        g = Math::Sphere::cartesian2spherical (dirs).leftCols(2);
 
      return Math::condition_number (Math::SH::init_transform (g, lmax));
    }
 
 
 
 
    Eigen::MatrixXd load_bvecs_bvals (const Header& header, const std::string& bvecs_path, const std::string& bvals_path);
 
 
 
 
    void save_bvecs_bvals (const Header&, const std::string&, const std::string&);
 
 
 
    namespace
    {
      template <class MatrixType>
      std::string scheme2str (const MatrixType& G)
      {
        std::string dw_scheme;
        for (ssize_t row = 0; row < G.rows(); ++row) {
          std::string line = str(G(row,0), 10);
          for (ssize_t col = 1; col < G.cols(); ++col)
            line += "," + str(G(row,col), 10);
          add_line (dw_scheme, line);
        }
        return dw_scheme;
      }
    }
 
 
 
 
    template <class MatrixType>
      void set_DW_scheme (Header& header, const MatrixType& G)
      {
        if (!G.rows()) {
          auto it = header.keyval().find ("dw_scheme");
          if (it != header.keyval().end())
            header.keyval().erase (it);
          return;
        }
        try {
          check_DW_scheme (header, G);
          header.keyval()["dw_scheme"] = scheme2str (G);
        } catch (Exception&) {
          WARN ("attempt to add non-matching DW scheme to header - ignored");
        }
      }
 
 
 
 
    Eigen::MatrixXd parse_DW_scheme (const Header& header);
 
 
 
 
    Eigen::MatrixXd get_raw_DW_scheme (const Header& header);
 
 
 
    void clear_DW_scheme (Header&);
 
 
 
 
    template <class MatrixType>
    void stash_DW_scheme (Header& header, const MatrixType& grad)
    {
      clear_DW_scheme (header);
      if (grad.rows())
        header.keyval()["prior_dw_scheme"] = scheme2str (grad);
    }
 
 
 
 
 
    enum class BValueScalingBehaviour { Auto, UserOn, UserOff };
    BValueScalingBehaviour get_cmdline_bvalue_scaling_behaviour ();
 
 
    Eigen::MatrixXd get_DW_scheme (const Header& header, BValueScalingBehaviour bvalue_scaling = BValueScalingBehaviour::Auto);
 
 
 
 
 
    void export_grad_commandline (const Header& header);
 
 
 
    template <class MatrixType>
      Eigen::MatrixXd compute_SH2amp_mapping (
          const MatrixType& directions,
          bool lmax_from_command_line = true,
          int default_lmax = 8)
      {
        int lmax = -1;
        int lmax_from_ndir = Math::SH::LforN (directions.rows());
        bool lmax_set_from_commandline = false;
        if (lmax_from_command_line) {
          auto opt = App::get_options ("lmax");
          if (opt.size()) {
            lmax_set_from_commandline = true;
            lmax = to<int> (opt[0][0]);
            if (lmax % 2)
              throw Exception ("lmax must be an even number");
            if (lmax < 0)
              throw Exception ("lmax must be a non-negative number");
            if (lmax > lmax_from_ndir) {
              WARN ("not enough directions for lmax = " + str(lmax) + " - dropping down to " + str(lmax_from_ndir));
              lmax = lmax_from_ndir;
            }
          }
        }
 
        if (lmax < 0) {
          lmax = lmax_from_ndir;
          if (lmax > default_lmax)
            lmax = default_lmax;
        }
 
        INFO ("computing SH transform using lmax = " + str (lmax));
 
        int lmax_prev = lmax;
        Eigen::MatrixXd mapping;
        do {
          mapping = Math::SH::init_transform (directions, lmax);
          const default_type cond = Math::condition_number (mapping);
          if (cond < 10.0)
            break;
          WARN ("directions are poorly distributed for lmax = " + str(lmax) + " (condition number = " + str (cond) + ")");
          if (cond < 100.0 || lmax_set_from_commandline)
            break;
          lmax -= 2;
        } while (lmax >= 0);
 
        if (lmax_prev != lmax)
          WARN ("reducing lmax to " + str(lmax) + " to improve conditioning");
 
        return mapping;
      }
 
 
 
 
 
 
      inline size_t lmax_for_directions (const Eigen::MatrixXd& directions,
          const bool lmax_from_command_line = true,
          const int default_lmax = 8)
      {
        const auto mapping = compute_SH2amp_mapping (directions, lmax_from_command_line, default_lmax);
        return Math::SH::LforN (mapping.cols());
      }
 
 
  }
}
 
#endif