From 160516d7ceccac0726b6756cba510b8d26252d1f Mon Sep 17 00:00:00 2001 From: Beatrice Schuster Date: Mon, 24 Jan 2022 22:59:21 +0000 Subject: [PATCH 1/7] Add skeleton for new response framework --- cmd/minimal.cpp | 23 +++++++++++++++++++++++ src/dwi/sdeconv/response.h | 2 ++ 2 files changed, 25 insertions(+) create mode 100644 cmd/minimal.cpp create mode 100644 src/dwi/sdeconv/response.h diff --git a/cmd/minimal.cpp b/cmd/minimal.cpp new file mode 100644 index 0000000000..ee905d5d8d --- /dev/null +++ b/cmd/minimal.cpp @@ -0,0 +1,23 @@ +#include "command.h" +// #include "dwi/sdeconv/response.h" + +using namespace MR; +using namespace App; + +void usage () { + + AUTHOR = "Joe Bloggs"; + SYNOPSIS = "example minimal command"; + DESCRIPTION + + "What the command does" + + "multiple paragraphs can be added like this"; + + ARGUMENTS + + Argument ("arg", "argument description").type_text(); +} + + +void run() +{ + +} \ No newline at end of file diff --git a/src/dwi/sdeconv/response.h b/src/dwi/sdeconv/response.h new file mode 100644 index 0000000000..73f1825e23 --- /dev/null +++ b/src/dwi/sdeconv/response.h @@ -0,0 +1,2 @@ +// #ifndef __response_h__ +// #define __response_h__ \ No newline at end of file From f147a4c64884f698864635ec77cca33e45fd4eb8 Mon Sep 17 00:00:00 2001 From: Beatrice Schuster Date: Tue, 25 Jan 2022 14:09:27 +0000 Subject: [PATCH 2/7] Add skeleton of new response class --- cmd/minimal.cpp | 36 ++++++++++++++++++++++++++++++++++-- src/dwi/sdeconv/response.h | 8 ++++++-- 2 files changed, 40 insertions(+), 4 deletions(-) diff --git a/cmd/minimal.cpp b/cmd/minimal.cpp index ee905d5d8d..a845ab3dea 100644 --- a/cmd/minimal.cpp +++ b/cmd/minimal.cpp @@ -1,5 +1,7 @@ #include "command.h" -// #include "dwi/sdeconv/response.h" +#include "dwi/sdeconv/response.h" +#include "dwi/shells.h" + using namespace MR; using namespace App; @@ -7,7 +9,7 @@ using namespace App; void usage () { AUTHOR = "Joe Bloggs"; - SYNOPSIS = "example minimal command"; + SYNOPSIS = "Estimate response function coefficients for a given b-value"; DESCRIPTION + "What the command does" + "multiple paragraphs can be added like this"; @@ -17,6 +19,36 @@ void usage () { } + +// Create a class to represent the response +// public attributes of the class: +// 1. Method to retrieve the coeffs at a given b-value +// 2. Method to query the max harmonic order (that it stores) + +class Response{ + using BValueList = decltype(std::declval().col(0)); + + public: + Eigen::VectorXd retrieve_coeffs(const BValueList bval, const Eigen::Vector3d& dir) + { + // this method retrieves the SH coefficients for a given b-value + Eigen::MatrixXd M, Y; + Eigen::VectorXd amplitudes, b; + + // coef_vec = S - M + // M[i] = Y[i] + + + // std::distance(x.begin(), std::min_element(x.begin(), x.end())); + } + + + int max_harmonic()b + { + + } +}; + void run() { diff --git a/src/dwi/sdeconv/response.h b/src/dwi/sdeconv/response.h index 73f1825e23..8b0b167cba 100644 --- a/src/dwi/sdeconv/response.h +++ b/src/dwi/sdeconv/response.h @@ -1,2 +1,6 @@ -// #ifndef __response_h__ -// #define __response_h__ \ No newline at end of file +#ifndef __dwi_sdeconv_response_h__ +#define __dwi_sdeconv_response_h__ + + + +#endif From 34a7c1bb99ffc1b895b7b6486b041f542afd7d2c Mon Sep 17 00:00:00 2001 From: Beatrice Schuster Date: Tue, 1 Feb 2022 14:09:53 +0000 Subject: [PATCH 3/7] Load matrix function --- src/dwi/sdeconv/response.h | 104 +++++++++++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) diff --git a/src/dwi/sdeconv/response.h b/src/dwi/sdeconv/response.h index 8b0b167cba..0aa7b32d58 100644 --- a/src/dwi/sdeconv/response.h +++ b/src/dwi/sdeconv/response.h @@ -1,6 +1,110 @@ #ifndef __dwi_sdeconv_response_h__ #define __dwi_sdeconv_response_h__ +#include "math/ZSH.h" +#include "math/math.h" +#include "mrtrix.h" +// Create a class to represent the response +// public attributes of the class: +// 1. Method to retrieve the coeffs at a given b-value +// 2. Method to query the max harmonic order (that it stores) + +// const MR::App::OptionGroup GradImportOptions = DWI::GradImportOptions(); + +namespace MR +{ + namespace DWI + { + namespace SDeconv + { + class Response + { + // add method to load matrixes and bvalues that correspond to it; + error handling () + + // to make sure that bvals.size() = original.coefs.rows() to make sure the size is the same. + // - if it's not the case, we throw an exception: assert it and check for it + + // constructor or method taking a const str as arg; load the data from that; + // - create constructor with the same arguments as the method + // get the comments, find the one with shells, strip the shells to get the nrs, parse_float for bvalues vector + + public: + // template + Eigen::VectorXd retrieve_coeffs(const double bval) + { + Eigen::VectorXd coefs(original_coefs.cols()); + + int i = 0; + vector dist_up, dist_down; + + while (bval < original_bvals[i]) + i++; + + // Linear interpolation + float ratio = (bval - original_bvals[i]) / (original_bvals[i + 1] - original_bvals[i]); + coefs = (1.0 - ratio) * original_coefs.row(i) + ratio * original_coefs.row(i + 1); + + if (i == original_coefs.cols()) + throw Exception("Iterator is out of boundaries"); + + // if (!(original_coefs.coeffs().minCoeff() < bval < original_coefs.coeffs().maxCoeff())) + // throw Exception("b-value is out of boundaries, cannot interpolate"); + + // // calculate difference between bval and elements of original_bvals + // if (bval == i) + // return bval; + // else if (i > bval) + // dist_up.push_back(abs(i - bval)); + // else + // dist_down.push_back(abs(i - bval)); + + // // find bval neighbours + // float up_bound = find(dist_up.begin(), dist_up.end(), min_element(dist_up.begin(), dist_up.end())), + // low_bound = find(dist_down.begin(), dist_down.end(), min_element(dist_down.begin(), dist_down.end())); + + // // interpolate + // float ratio = (bval - low_bound) / (up_bound - low_bound); + // coefs = (1.0 - ratio) * original_coefs.row(low_bound) + original_coefs.row(up_bound) * ratio; + + return coefs; + } + + int lmax() const + { + + return Math::ZSH::LforN(original_coefs.cols()); + } + + // LOAD MATRIX + Eigen::MatrixXd load_matrix(const std::string &filename, vector* comments) + { + auto vec = load_matrix_2D_vector(filename, comments); + Eigen::MatrixXd M(vec.size(), vec[0].size()); + + for (ssize_t r = 0; r < vec.size(); r++) + for (ssize_t c = 0; c < vec[0].size(); r++) + M(r, c) = vec[r][c]; + + return M; + } + + + // const bool shell_bvalues = App::get_options("shell_bvalues").size(); + // const bool shell_sizes = App::get_options("shell_sizes").size(); + // const bool shell_indices = App::get_options("shell_indices").size(); + + // Eigen::MatrixXd grad = DWI::get_DW_scheme(header, DWI::get_cmdline_bvalue_scaling_behaviour()); + + // if (shell_bvalues || shell_sizes || shell_indices) + // print_shells(grad, shell_bvalues, shell_sizes, shell_indices); + + private: + Eigen::MatrixXd original_coefs; + vector original_bvals; + }; + } + } +} #endif From b273db423cb2ce57ebae72bbcaf42eccc65568d7 Mon Sep 17 00:00:00 2001 From: Beatrice Schuster Date: Tue, 1 Feb 2022 15:38:54 +0000 Subject: [PATCH 4/7] Response function now functional --- cmd/minimal.cpp | 34 ++----------- src/dwi/sdeconv/response.h | 99 +++++++++++++------------------------- 2 files changed, 37 insertions(+), 96 deletions(-) diff --git a/cmd/minimal.cpp b/cmd/minimal.cpp index a845ab3dea..a8b935cfd5 100644 --- a/cmd/minimal.cpp +++ b/cmd/minimal.cpp @@ -1,6 +1,6 @@ #include "command.h" #include "dwi/sdeconv/response.h" -#include "dwi/shells.h" +// #include "dwi/shells.h" using namespace MR; @@ -20,36 +20,10 @@ void usage () { -// Create a class to represent the response -// public attributes of the class: -// 1. Method to retrieve the coeffs at a given b-value -// 2. Method to query the max harmonic order (that it stores) - -class Response{ - using BValueList = decltype(std::declval().col(0)); - - public: - Eigen::VectorXd retrieve_coeffs(const BValueList bval, const Eigen::Vector3d& dir) - { - // this method retrieves the SH coefficients for a given b-value - Eigen::MatrixXd M, Y; - Eigen::VectorXd amplitudes, b; - - // coef_vec = S - M - // M[i] = Y[i] - - - // std::distance(x.begin(), std::min_element(x.begin(), x.end())); - } - - - int max_harmonic()b - { - - } -}; - void run() { + DWI::SDeconv::Response response(argument[0]); + VAR(response.lmax()); + VAR(response.coeffs(3500)); } \ No newline at end of file diff --git a/src/dwi/sdeconv/response.h b/src/dwi/sdeconv/response.h index 0aa7b32d58..e93b81e437 100644 --- a/src/dwi/sdeconv/response.h +++ b/src/dwi/sdeconv/response.h @@ -4,13 +4,7 @@ #include "math/ZSH.h" #include "math/math.h" #include "mrtrix.h" - -// Create a class to represent the response -// public attributes of the class: -// 1. Method to retrieve the coeffs at a given b-value -// 2. Method to query the max harmonic order (that it stores) - -// const MR::App::OptionGroup GradImportOptions = DWI::GradImportOptions(); +#include "debug.h" namespace MR { @@ -20,84 +14,57 @@ namespace MR { class Response { - // add method to load matrixes and bvalues that correspond to it; + error handling () - - // to make sure that bvals.size() = original.coefs.rows() to make sure the size is the same. - // - if it's not the case, we throw an exception: assert it and check for it - - // constructor or method taking a const str as arg; load the data from that; - // - create constructor with the same arguments as the method - // get the comments, find the one with shells, strip the shells to get the nrs, parse_float for bvalues vector public: - // template - Eigen::VectorXd retrieve_coeffs(const double bval) + Response() {} + Response(const std::string &filename) { load(filename); } + + Eigen::VectorXd coeffs(const double bval) { - Eigen::VectorXd coefs(original_coefs.cols()); + if (bval < original_bvals[0]) + throw Exception("bvalue out of bounds"); - int i = 0; - vector dist_up, dist_down; + if (bval >= original_bvals.back()) + return original_coefs.row(original_coefs.rows() - 1); - while (bval < original_bvals[i]) + int i = 0; + while (bval > original_bvals[i + 1]) i++; - // Linear interpolation float ratio = (bval - original_bvals[i]) / (original_bvals[i + 1] - original_bvals[i]); - coefs = (1.0 - ratio) * original_coefs.row(i) + ratio * original_coefs.row(i + 1); - - if (i == original_coefs.cols()) - throw Exception("Iterator is out of boundaries"); - - // if (!(original_coefs.coeffs().minCoeff() < bval < original_coefs.coeffs().maxCoeff())) - // throw Exception("b-value is out of boundaries, cannot interpolate"); - - // // calculate difference between bval and elements of original_bvals - // if (bval == i) - // return bval; - // else if (i > bval) - // dist_up.push_back(abs(i - bval)); - // else - // dist_down.push_back(abs(i - bval)); - - // // find bval neighbours - // float up_bound = find(dist_up.begin(), dist_up.end(), min_element(dist_up.begin(), dist_up.end())), - // low_bound = find(dist_down.begin(), dist_down.end(), min_element(dist_down.begin(), dist_down.end())); - - // // interpolate - // float ratio = (bval - low_bound) / (up_bound - low_bound); - // coefs = (1.0 - ratio) * original_coefs.row(low_bound) + original_coefs.row(up_bound) * ratio; - - return coefs; + return (1.0 - ratio) * original_coefs.row(i) + ratio * original_coefs.row(i + 1); } - int lmax() const - { + int lmax() const { return Math::ZSH::LforN(original_coefs.cols()); } - return Math::ZSH::LforN(original_coefs.cols()); - } + const Eigen::VectorXd &bvalues() const { return original_bvals; } // LOAD MATRIX - Eigen::MatrixXd load_matrix(const std::string &filename, vector* comments) + void load(const std::string &filename) { - auto vec = load_matrix_2D_vector(filename, comments); - Eigen::MatrixXd M(vec.size(), vec[0].size()); + vector comments; + auto vec = load_matrix_2D_vector(filename, &comments); - for (ssize_t r = 0; r < vec.size(); r++) - for (ssize_t c = 0; c < vec[0].size(); r++) - M(r, c) = vec[r][c]; + //load bvalues + for (const auto &comment : comments) + { + auto n = comment.rfind("Shells:"); + if (n != comment.npos) + original_bvals = parse_floats(comment.substr(n + 7)); + } - return M; - } + if (original_bvals.size() != vec.size()) + throw Exception("Number of b-values does not match the number of rows in response"); - - // const bool shell_bvalues = App::get_options("shell_bvalues").size(); - // const bool shell_sizes = App::get_options("shell_sizes").size(); - // const bool shell_indices = App::get_options("shell_indices").size(); + if (!std::is_sorted(original_bvals.begin(), original_bvals.end())) + throw Exception("B-values are not sorted in ascending order"); - // Eigen::MatrixXd grad = DWI::get_DW_scheme(header, DWI::get_cmdline_bvalue_scaling_behaviour()); + original_coefs.resize(vec.size(), vec[0].size()); - // if (shell_bvalues || shell_sizes || shell_indices) - // print_shells(grad, shell_bvalues, shell_sizes, shell_indices); + for (ssize_t r = 0; r < vec.size(); r++) + for (ssize_t c = 0; c < vec[0].size(); c++) + original_coefs(r, c) = vec[r][c]; + } private: Eigen::MatrixXd original_coefs; From 98f8f895975c20c8f9423b88b2f3b6a2b18ea78b Mon Sep 17 00:00:00 2001 From: Beatrice Schuster Date: Tue, 8 Feb 2022 17:24:11 +0100 Subject: [PATCH 5/7] dwi2fod: first attempt at arbitrary dwi encoding fit --- cmd/dwi2fod.cpp | 1 - src/dwi/sdeconv/msmt_csd.h | 401 +++++++++++++++++-------------------- src/dwi/sdeconv/response.h | 17 +- 3 files changed, 195 insertions(+), 224 deletions(-) diff --git a/cmd/dwi2fod.cpp b/cmd/dwi2fod.cpp index 0473efd34f..6b00024209 100644 --- a/cmd/dwi2fod.cpp +++ b/cmd/dwi2fod.cpp @@ -320,7 +320,6 @@ void run () MSMT_Processor processor (shared, mask, odfs, dwi_modelled); auto dwi = header_in.get_image().with_direct_io (3); ThreadedLoop ("performing MSMT CSD (" - + str(shared.num_shells()) + " shell" + (shared.num_shells() > 1 ? "s" : "") + ", " + str(num_tissues) + " tissue" + (num_tissues > 1 ? "s" : "") + ")", dwi, 0, 3) .run (processor, dwi); diff --git a/src/dwi/sdeconv/msmt_csd.h b/src/dwi/sdeconv/msmt_csd.h index 409e3d2582..ad82082da0 100644 --- a/src/dwi/sdeconv/msmt_csd.h +++ b/src/dwi/sdeconv/msmt_csd.h @@ -20,11 +20,11 @@ #include "header.h" #include "types.h" #include "dwi/gradient.h" -#include "dwi/shells.h" #include "math/constrained_least_squares.h" #include "math/math.h" #include "math/SH.h" #include "math/ZSH.h" +#include "dwi/sdeconv/response.h" #include "dwi/directions/predefined.h" @@ -41,248 +41,217 @@ namespace MR extern const App::OptionGroup MSMT_CSD_options; - class MSMT_CSD { MEMALIGN(MSMT_CSD) + class MSMT_CSD + { + MEMALIGN(MSMT_CSD) + public: + class Shared + { + MEMALIGN(Shared) public: + Shared(const Header &dwi_header) : grad(DWI::get_DW_scheme(dwi_header)), + HR_dirs(DWI::Directions::electrostatic_repulsion_300()), + solution_min_norm_regularisation(DEFAULT_MSMTCSD_NORM_LAMBDA), + constraint_min_norm_regularisation(DEFAULT_MSMTCSD_NEG_LAMBDA) {} + + void parse_cmdline_options() + { + using namespace App; + auto opt = get_options("lmax"); + if (opt.size()) + lmax = parse_ints(opt[0][0]); + opt = get_options("directions"); + if (opt.size()) + HR_dirs = load_matrix(opt[0][0]); + opt = get_options("norm_lambda"); + if (opt.size()) + solution_min_norm_regularisation = opt[0][0]; + opt = get_options("neg_lambda"); + if (opt.size()) + constraint_min_norm_regularisation = opt[0][0]; + } - class Shared { MEMALIGN(Shared) - public: - Shared (const Header& dwi_header) : - grad (DWI::get_DW_scheme (dwi_header)), - shells (grad), - HR_dirs (DWI::Directions::electrostatic_repulsion_300()), - solution_min_norm_regularisation (DEFAULT_MSMTCSD_NORM_LAMBDA), - constraint_min_norm_regularisation (DEFAULT_MSMTCSD_NEG_LAMBDA) { shells.select_shells(false,false,false); } - - - void parse_cmdline_options() + void set_responses(const vector &files) + { + lmax_response.clear(); + for (const auto &s : files) + { + Response r; + try { - using namespace App; - auto opt = get_options ("lmax"); - if (opt.size()) - lmax = parse_ints (opt[0][0]); - opt = get_options ("directions"); - if (opt.size()) - HR_dirs = load_matrix (opt[0][0]); - opt = get_options ("norm_lambda"); - if (opt.size()) - solution_min_norm_regularisation = opt[0][0]; - opt = get_options ("neg_lambda"); - if (opt.size()) - constraint_min_norm_regularisation = opt[0][0]; + r.load(s); } - - - - void set_responses (const vector& files) + catch (Exception &e) { - lmax_response.clear(); - for (const auto& s : files) { - Eigen::MatrixXd r; - try { - r = load_matrix (s); - } catch (Exception& e) { - throw Exception (e, "File \"" + s + "\" is not a valid response function file"); - } - responses.push_back (std::move (r)); - } - prepare_responses(); - response_files = files; + throw Exception(e, "File \"" + s + "\" is not a valid response function file"); } + responses.push_back(std::move(r)); + } + response_files = files; + } - void set_responses (const vector& matrices) + void init() + { + if (lmax.empty()) + { + lmax = lmax_response; + for (size_t t = 0; t != num_tissues(); ++t) { - responses = matrices; - prepare_responses(); + lmax[t] = std::min(uint32_t(DEFAULT_MSMTCSD_LMAX), lmax[t]); } - - - - void init() + } + else + { + if (lmax.size() != num_tissues()) + throw Exception("Number of lmaxes specified (" + str(lmax.size()) + ") does not match number of tissues (" + str(num_tissues()) + ")"); + for (const auto i : lmax) { - if (lmax.empty()) { - lmax = lmax_response; - for (size_t t = 0; t != num_tissues(); ++t) { - lmax[t] = std::min (uint32_t(DEFAULT_MSMTCSD_LMAX), lmax[t]); - } - } else { - if (lmax.size() != num_tissues()) - throw Exception ("Number of lmaxes specified (" + str(lmax.size()) + ") does not match number of tissues (" + str(num_tissues()) + ")"); - for (const auto i : lmax) { - if (i % 2) - throw Exception ("Each value of lmax must be a non-negative even integer"); - } - } - - for (size_t t = 0; t != num_tissues(); ++t) { - if (size_t(responses[t].rows()) != num_shells()) - throw Exception ("number of rows in response functions must match number of b-value shells; " - "number of shells is " + str(num_shells()) + ", but file \"" + response_files[t] + "\" contains " + str(responses[t].rows()) + " rows"); - // Pad response functions out to the requested lmax for this tissue - responses[t].conservativeResizeLike (Eigen::MatrixXd::Zero (num_shells(), Math::ZSH::NforL (lmax[t]))); - } - - ////////////////////////////////////////////////// - // Set up the constrained least squares problem // - ////////////////////////////////////////////////// - - size_t nparams = 0; - uint32_t maxlmax = 0; - for (size_t i = 0; i < num_tissues(); i++) { - nparams += Math::SH::NforL (lmax[i]); - maxlmax = std::max (maxlmax, lmax[i]); - } - - INFO ("initialising multi-tissue CSD for " + str(num_tissues()) + " tissue types, with " + str (nparams) + " parameters"); - - Eigen::MatrixXd C = Eigen::MatrixXd::Zero (grad.rows(), nparams); - - vector dwilist; - for (size_t i = 0; i != size_t(grad.rows()); i++) - dwilist.push_back(i); - - Eigen::MatrixXd directions = DWI::gen_direction_matrix (grad, dwilist); - Eigen::MatrixXd SHT = Math::SH::init_transform (directions, maxlmax); - for (ssize_t i = 0; i < SHT.rows(); i++) - for (ssize_t j = 0; j < SHT.cols(); j++) - if (std::isnan (SHT(i,j))) - SHT(i,j) = 0.0; - - // TODO: is this just computing the Associated Legrendre polynomials...? - Eigen::MatrixXd delta(1,2); delta << 0, 0; - Eigen::MatrixXd DSH__ = Math::SH::init_transform (delta, maxlmax); - Eigen::VectorXd DSH_ = DSH__.row(0); - Eigen::VectorXd DSH (maxlmax/2+1); - size_t j = 0; - for (ssize_t i = 0; i < DSH_.size(); i++) - if (DSH_[i] != 0.0) { - DSH[j] = DSH_[i]; - j++; - } - - size_t pbegin = 0; - for (size_t tissue_idx = 0; tissue_idx < num_tissues(); ++tissue_idx) { - const size_t tissue_lmax = lmax[tissue_idx]; - const size_t tissue_n = Math::SH::NforL (tissue_lmax); - const size_t tissue_nmzero = tissue_lmax/2+1; - - for (size_t shell_idx = 0; shell_idx < num_shells(); ++shell_idx) { - Eigen::VectorXd response_ = responses[tissue_idx].row (shell_idx); - response_ = (response_.array() / DSH.head (tissue_nmzero).array()).matrix(); - Eigen::VectorXd fconv (tissue_n); - int li = 0; int mi = 0; - for (int l = 0; l <= int(tissue_lmax); l+=2) { - for (int m = -l; m <= l; m++) { - fconv[mi] = response_[li]; - mi++; - } - li++; - } - vector vols = shells[shell_idx].get_volumes(); - for (size_t idx = 0; idx < vols.size(); idx++) { - Eigen::VectorXd SHT_(SHT.row (vols[idx]).head (tissue_n)); - SHT_ = (SHT_.array()*fconv.array()).matrix(); - C.row (vols[idx]).segment (pbegin, tissue_n) = SHT_; - } - } - pbegin += tissue_n; - } - - vector m (num_tissues()); - vector n (num_tissues()); - size_t M = 0; - size_t N = 0; - - Eigen::MatrixXd HR_SHT = Math::SH::init_transform (HR_dirs, maxlmax); - - for (size_t i = 0; i != num_tissues(); i++) { - if (lmax[i] > 0) - m[i] = HR_dirs.rows(); - else - m[i] = 1; - M += m[i]; - n[i] = Math::SH::NforL (lmax[i]); - N += n[i]; - } - - Eigen::MatrixXd A (Eigen::MatrixXd::Zero (M, N)); - size_t b_m = 0; size_t b_n = 0; - for (size_t i = 0; i != num_tissues(); i++) { - A.block (b_m,b_n,m[i],n[i]) = HR_SHT.block (0,0,m[i],n[i]); - b_m += m[i]; - b_n += n[i]; - } - problem = Math::ICLS::Problem (C, A, Eigen::VectorXd(), 0, - solution_min_norm_regularisation, constraint_min_norm_regularisation); - - INFO ("Multi-shell, multi-tissue CSD initialised successfully"); + if (i % 2) + throw Exception("Each value of lmax must be a non-negative even integer"); + } + } + + for (size_t t = 0; t != num_tissues(); ++t) + { + // more appropriate test would be to check that the max bval is in the gradient table + } + + ////////////////////////////////////////////////// + // Set up the constrained least squares problem // + ////////////////////////////////////////////////// + + size_t nparams = 0; + uint32_t maxlmax = 0; + for (size_t i = 0; i < num_tissues(); i++) + { + nparams += Math::SH::NforL(lmax[i]); + maxlmax = std::max(maxlmax, lmax[i]); + } + + INFO("initialising multi-tissue CSD for " + str(num_tissues()) + " tissue types, with " + str(nparams) + " parameters"); + + Eigen::MatrixXd C = Eigen::MatrixXd::Zero(grad.rows(), nparams); + + vector dwilist; + for (size_t i = 0; i != size_t(grad.rows()); i++) + dwilist.push_back(i); + + Eigen::MatrixXd directions = DWI::gen_direction_matrix(grad, dwilist); + Eigen::MatrixXd SHT = Math::SH::init_transform(directions, maxlmax); + for (ssize_t i = 0; i < SHT.rows(); i++) + for (ssize_t j = 0; j < SHT.cols(); j++) + if (std::isnan(SHT(i, j))) + SHT(i, j) = 0.0; + + // TODO: is this just computing the Associated Legrendre polynomials...? + Eigen::MatrixXd delta(1, 2); + delta << 0, 0; + Eigen::MatrixXd DSH__ = Math::SH::init_transform(delta, maxlmax); + Eigen::VectorXd DSH_ = DSH__.row(0); + Eigen::VectorXd DSH(maxlmax / 2 + 1); + size_t j = 0; + for (ssize_t i = 0; i < DSH_.size(); i++) + if (DSH_[i] != 0.0) + { + DSH[j] = DSH_[i]; + j++; } + size_t pbegin = 0; + for (size_t tissue_idx = 0; tissue_idx < num_tissues(); ++tissue_idx) + { + const size_t tissue_lmax = lmax[tissue_idx]; + const size_t tissue_n = Math::SH::NforL(tissue_lmax); + const size_t tissue_nmzero = tissue_lmax / 2 + 1; - size_t num_shells() const { return shells.count(); } - size_t num_tissues() const { return responses.size(); } - - - const Eigen::MatrixXd grad; - DWI::Shells shells; - Eigen::MatrixXd HR_dirs; - vector lmax, lmax_response; - vector responses; - vector response_files; - Math::ICLS::Problem problem; - double solution_min_norm_regularisation, constraint_min_norm_regularisation; - - - private: - - void prepare_responses() { - for (size_t t = 0; t != num_tissues(); ++t) { - Eigen::MatrixXd& r (responses[t]); - size_t n = 0; - for (size_t row = 0; row < size_t(r.rows()); row++) { - for (size_t col = 0; col < size_t(r.cols()); col++) { - if (r(row,col)) - n = std::max (n, col+1); - } + // for (size_t shell_idx = 0; shell_idx < num_shells(); ++shell_idx) { + // loop over all volumes + for (size_t vol = 0; vol < C.rows(); ++vol) + { + Eigen::VectorXd response_ = responses[tissue_idx].coeffs(grad(vol, 3)); + response_ = (response_.array() / DSH.head(tissue_nmzero).array()).matrix(); + Eigen::VectorXd fconv(tissue_n); + int li = 0; + int mi = 0; + for (int l = 0; l <= int(tissue_lmax); l += 2) + { + for (int m = -l; m <= l; m++) + { + fconv[mi] = response_[li]; + mi++; } - // Clip off any empty columns, i.e. degrees containing zero coefficients for all shells - r.conservativeResize (r.rows(), n); - // Store the lmax for each tissue based on their response functions; - // if the user doesn't manually specify lmax, these will determine the - // lmax of each tissue ODF output, with a further default lmax=8 - // restriction at that stage - lmax_response.push_back (Math::ZSH::LforN (r.cols())); + li++; } + Eigen::VectorXd SHT_(SHT.row(vol).head(tissue_n)); + SHT_ = (SHT_.array() * fconv.array()).matrix(); + C.row(vol).segment(pbegin, tissue_n) = SHT_; } + pbegin += tissue_n; + } + + vector m(num_tissues()); + vector n(num_tissues()); + size_t M = 0; + size_t N = 0; + + Eigen::MatrixXd HR_SHT = Math::SH::init_transform(HR_dirs, maxlmax); + + for (size_t i = 0; i != num_tissues(); i++) + { + if (lmax[i] > 0) + m[i] = HR_dirs.rows(); + else + m[i] = 1; + M += m[i]; + n[i] = Math::SH::NforL(lmax[i]); + N += n[i]; + } + + Eigen::MatrixXd A(Eigen::MatrixXd::Zero(M, N)); + size_t b_m = 0; + size_t b_n = 0; + for (size_t i = 0; i != num_tissues(); i++) + { + A.block(b_m, b_n, m[i], n[i]) = HR_SHT.block(0, 0, m[i], n[i]); + b_m += m[i]; + b_n += n[i]; + } + problem = Math::ICLS::Problem(C, A, Eigen::VectorXd(), 0, + solution_min_norm_regularisation, constraint_min_norm_regularisation); + + INFO("Multi-shell, multi-tissue CSD initialised successfully"); + } - }; - - + size_t num_tissues() const { return responses.size(); } + const Eigen::MatrixXd grad; + Eigen::MatrixXd HR_dirs; + vector lmax, lmax_response; + vector responses; + vector response_files; + Math::ICLS::Problem problem; + double solution_min_norm_regularisation, constraint_min_norm_regularisation; - MSMT_CSD (const Shared& shared_data) : - niter (0), - shared (shared_data), - solver (shared.problem) { } + }; - void operator() (const Eigen::VectorXd& data, Eigen::VectorXd& output) { - niter = solver (output, data); - } + MSMT_CSD(const Shared &shared_data) : niter(0), + shared(shared_data), + solver(shared.problem) {} - size_t niter; - const Shared& shared; - - private: - Math::ICLS::Solver solver; + void operator()(const Eigen::VectorXd &data, Eigen::VectorXd &output) + { + niter = solver(output, data); + } + size_t niter; + const Shared &shared; + private: + Math::ICLS::Solver solver; }; - - } } } #endif - - diff --git a/src/dwi/sdeconv/response.h b/src/dwi/sdeconv/response.h index e93b81e437..082b35126c 100644 --- a/src/dwi/sdeconv/response.h +++ b/src/dwi/sdeconv/response.h @@ -18,6 +18,7 @@ namespace MR public: Response() {} Response(const std::string &filename) { load(filename); } + Response(Response &&other) = default; Eigen::VectorXd coeffs(const double bval) { @@ -25,19 +26,21 @@ namespace MR throw Exception("bvalue out of bounds"); if (bval >= original_bvals.back()) - return original_coefs.row(original_coefs.rows() - 1); + return original_coeffs.row(original_coeffs.rows() - 1); int i = 0; while (bval > original_bvals[i + 1]) i++; float ratio = (bval - original_bvals[i]) / (original_bvals[i + 1] - original_bvals[i]); - return (1.0 - ratio) * original_coefs.row(i) + ratio * original_coefs.row(i + 1); + return (1.0 - ratio) * original_coeffs.row(i) + ratio * original_coeffs.row(i + 1); } - int lmax() const { return Math::ZSH::LforN(original_coefs.cols()); } + int lmax() const { return Math::ZSH::LforN(original_coeffs.cols()); } - const Eigen::VectorXd &bvalues() const { return original_bvals; } + const vector &bvalues() const { return original_bvals; } + + // int size() const { return original_bvals.size(); } // LOAD MATRIX void load(const std::string &filename) @@ -59,15 +62,15 @@ namespace MR if (!std::is_sorted(original_bvals.begin(), original_bvals.end())) throw Exception("B-values are not sorted in ascending order"); - original_coefs.resize(vec.size(), vec[0].size()); + original_coeffs.resize(vec.size(), vec[0].size()); for (ssize_t r = 0; r < vec.size(); r++) for (ssize_t c = 0; c < vec[0].size(); c++) - original_coefs(r, c) = vec[r][c]; + original_coeffs(r, c) = vec[r][c]; } private: - Eigen::MatrixXd original_coefs; + Eigen::MatrixXd original_coeffs; vector original_bvals; }; } From 58aa1758f95bbb0c7c7af0b3839f3dae727c4d20 Mon Sep 17 00:00:00 2001 From: Beatrice Schuster Date: Wed, 23 Feb 2022 15:44:45 +0100 Subject: [PATCH 6/7] updated working code for msmt csd --- cmd/dwi2fod.cpp | 2 -- src/dwi/sdeconv/msmt_csd.h | 4 +--- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/cmd/dwi2fod.cpp b/cmd/dwi2fod.cpp index 6b00024209..382458f74e 100644 --- a/cmd/dwi2fod.cpp +++ b/cmd/dwi2fod.cpp @@ -284,7 +284,6 @@ void run () if (argument.size() % 2) throw Exception ("MSMT_CSD algorithm expects pairs of (input response function & output FOD image) to be provided"); - DWI::SDeconv::MSMT_CSD::Shared shared (header_in); shared.parse_cmdline_options(); @@ -301,7 +300,6 @@ void run () } catch (Exception& e) { throw Exception (e, "MSMT_CSD algorithm expects the first file in each argument pair to be an input response function file"); } - shared.init(); DWI::stash_DW_scheme (header_out, shared.grad); diff --git a/src/dwi/sdeconv/msmt_csd.h b/src/dwi/sdeconv/msmt_csd.h index ad82082da0..a6fb5e7d72 100644 --- a/src/dwi/sdeconv/msmt_csd.h +++ b/src/dwi/sdeconv/msmt_csd.h @@ -94,10 +94,9 @@ namespace MR { if (lmax.empty()) { - lmax = lmax_response; for (size_t t = 0; t != num_tissues(); ++t) { - lmax[t] = std::min(uint32_t(DEFAULT_MSMTCSD_LMAX), lmax[t]); + lmax.push_back(std::min(DEFAULT_MSMTCSD_LMAX, responses[t].lmax())); } } else @@ -110,7 +109,6 @@ namespace MR throw Exception("Each value of lmax must be a non-negative even integer"); } } - for (size_t t = 0; t != num_tissues(); ++t) { // more appropriate test would be to check that the max bval is in the gradient table From 9ed37c454cc84e49e9157befd55594e361c25f23 Mon Sep 17 00:00:00 2001 From: J-Donald Tournier Date: Wed, 23 Feb 2022 15:33:25 +0000 Subject: [PATCH 7/7] dwi2fod: fix minor glitch in SH->RH conversion --- src/dwi/sdeconv/msmt_csd.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dwi/sdeconv/msmt_csd.h b/src/dwi/sdeconv/msmt_csd.h index a6fb5e7d72..f382e9d4e7 100644 --- a/src/dwi/sdeconv/msmt_csd.h +++ b/src/dwi/sdeconv/msmt_csd.h @@ -162,12 +162,12 @@ namespace MR const size_t tissue_n = Math::SH::NforL(tissue_lmax); const size_t tissue_nmzero = tissue_lmax / 2 + 1; - // for (size_t shell_idx = 0; shell_idx < num_shells(); ++shell_idx) { // loop over all volumes for (size_t vol = 0; vol < C.rows(); ++vol) { Eigen::VectorXd response_ = responses[tissue_idx].coeffs(grad(vol, 3)); - response_ = (response_.array() / DSH.head(tissue_nmzero).array()).matrix(); + response_.conservativeResize (tissue_nmzero); + response_.array() = (response_.array() / DSH.head(tissue_nmzero).array()); Eigen::VectorXd fconv(tissue_n); int li = 0; int mi = 0;