From 2f7a449f6062d9ed391a61e820d8e9edaf399bf3 Mon Sep 17 00:00:00 2001 From: Rinat Mukhometzianov <25082858+rmukh@users.noreply.github.com> Date: Thu, 21 Oct 2021 21:26:22 -0400 Subject: [PATCH] Simplify requirements for dmri input. Add support of external grad dirs file as a main table --- Ridgelets/Ridgelets.cpp | 19 ++-- Ridgelets/Ridgelets.vcxproj | 8 +- .../SphericalRidgelets/include/DATA_SOURCE.h | 2 +- .../include/SIGNAL_GENERATOR.h | 1 + .../SphericalRidgelets/src/DATA_SOURCE.cpp | 20 ++++- .../src/SIGNAL_GENERATOR.cpp | 87 ++++++++++++++----- Ridgelets/packages.config | 2 +- 7 files changed, 97 insertions(+), 42 deletions(-) diff --git a/Ridgelets/Ridgelets.cpp b/Ridgelets/Ridgelets.cpp index dcbb766..1d62e36 100644 --- a/Ridgelets/Ridgelets.cpp +++ b/Ridgelets/Ridgelets.cpp @@ -21,7 +21,7 @@ int main(int argc, char* argv[]) { #if defined(USE_FLOAT) || defined(UKF_USE_FLOAT) - std::cout << "The package was compiled for float precision computations" << std::endl; + cout << "The package was compiled for float precision computations" << endl; #endif Eigen::initParallel(); @@ -60,17 +60,12 @@ int main(int argc, char* argv[]) // 4D dMRI image to Eigen 2D Matrix MatrixType signal; MatrixType GradientDirections; // Matrix with dMRI image gradient directions - int res_dmri = data.DWI2Matrix(input_args.input_dmri, mask, signal, GradientDirections); + int res_dmri = data.DWI2Matrix(&input_args, mask, signal, GradientDirections); if (res_dmri) return EXIT_SUCCESS; - // Beginning of the main computational part SPH_RIDG ridg(input_args.sph_J, 1 / input_args.sph_rho); - MatrixType A; - cout << "Start computing R basis..." << endl; - ridg.RBasis(A, GradientDirections); - ridg.normBasis(A); - + // Read external gradients MatrixType ext_g; MatrixType ext_A; @@ -80,6 +75,12 @@ int main(int argc, char* argv[]) ridg.normBasis(ext_A); } + // Beginning of the main computational part + MatrixType A; + cout << "Start computing R basis..." << endl; + ridg.RBasis(A, GradientDirections); + ridg.normBasis(A); + if (!input_args.output_A.empty()) // -A { cout << "Saving A basis..." << endl; @@ -108,7 +109,7 @@ int main(int argc, char* argv[]) cout << "Start computing ridgelets coefficients..." << endl; auto start = high_resolution_clock::now(); //have a potentinal for optimization - slv.FISTA(C, input_args.n_splits, input_args.fista_iterations, input_args.fista_tolerance); + slv.FISTA(C, input_args.n_splits, input_args.fista_iterations, input_args.fista_tolerance); auto stop = high_resolution_clock::now(); auto ds = duration_cast(stop - start); auto dm = duration_cast(stop - start); diff --git a/Ridgelets/Ridgelets.vcxproj b/Ridgelets/Ridgelets.vcxproj index 64490a0..14cbdf3 100644 --- a/Ridgelets/Ridgelets.vcxproj +++ b/Ridgelets/Ridgelets.vcxproj @@ -83,8 +83,8 @@ false - C:\Users\mukho\source\repos\SphericalRidgelet\Ridgelets\SphericalRidgelets\include;$(IncludePath) - C:\Users\mukho\source\repos\SphericalRidgelet\Ridgelets\SphericalRidgelets\src;$(SourcePath) + D:\VSProjects\SphericalRidgelets\Ridgelets\SphericalRidgelets\include;$(IncludePath) + D:\VSProjects\SphericalRidgelets\Ridgelets\SphericalRidgelets\src;$(SourcePath) @@ -197,12 +197,12 @@ - + This project references NuGet package(s) that are missing on this computer. Use NuGet Package Restore to download them. For more information, see http://go.microsoft.com/fwlink/?LinkID=322105. The missing file is {0}. - + \ No newline at end of file diff --git a/Ridgelets/SphericalRidgelets/include/DATA_SOURCE.h b/Ridgelets/SphericalRidgelets/include/DATA_SOURCE.h index bbf995a..518de6b 100644 --- a/Ridgelets/SphericalRidgelets/include/DATA_SOURCE.h +++ b/Ridgelets/SphericalRidgelets/include/DATA_SOURCE.h @@ -32,7 +32,7 @@ class DATA_SOURCE void data_saving_info_out(unsigned long int coef_size, string name); void estimate_memory(MatrixType & s, MatrixType & A, const input_parse & params); int compute_splits(unsigned s_size); - int DWI2Matrix(string & dmri_file, MaskImagePointer & mask, MatrixType & signal, MatrixType & grad_dirs); + int DWI2Matrix(input_parse* input_args, MaskImagePointer & mask, MatrixType & signal, MatrixType & grad_dirs); void Matrix2DWI(DiffusionImagePointer & img, MaskImagePointer & mask, MatrixType & arr); void matrixToFile(const string & fname, MatrixType & matrix); void fileToMatrix(const string & fname, MatrixType & arr); diff --git a/Ridgelets/SphericalRidgelets/include/SIGNAL_GENERATOR.h b/Ridgelets/SphericalRidgelets/include/SIGNAL_GENERATOR.h index 1e8b2a1..0121a93 100644 --- a/Ridgelets/SphericalRidgelets/include/SIGNAL_GENERATOR.h +++ b/Ridgelets/SphericalRidgelets/include/SIGNAL_GENERATOR.h @@ -11,6 +11,7 @@ class SIGNAL_GENERATOR bool is_path_exists(const string & s); int readVolume(MatrixType & GradientDirections, DiffusionImagePointer & image); int ExtractMatrix(MaskImagePointer & mask, MatrixType & signal, MatrixType & grad_dirs); + void remove_row(MatrixType& a, Eigen::Index del); unsigned nGradImgs; // Number of gradient images unsigned nOfImgs; // Total number of images (including b0) dMRI_h_info h; diff --git a/Ridgelets/SphericalRidgelets/src/DATA_SOURCE.cpp b/Ridgelets/SphericalRidgelets/src/DATA_SOURCE.cpp index a0439d7..c1800da 100644 --- a/Ridgelets/SphericalRidgelets/src/DATA_SOURCE.cpp +++ b/Ridgelets/SphericalRidgelets/src/DATA_SOURCE.cpp @@ -217,7 +217,7 @@ int DATA_SOURCE::readMask(string inputMask, MaskImagePointer& image) { if (ext_vol.compare("nhdr")) { if (ext_vol.compare("nrrd")) { - cout << "NDHR or NRRD file formats only! Please, check file type." << endl; + cerr << "NDHR or NRRD file formats only! Please, check file type." << endl; return EXIT_FAILURE; } } @@ -351,9 +351,23 @@ int DATA_SOURCE::compute_splits(unsigned s_size) { return s; } -int DATA_SOURCE::DWI2Matrix(string &dmri_file, MaskImagePointer &mask, MatrixType &signal, MatrixType &grad_dirs) +int DATA_SOURCE::DWI2Matrix(input_parse* input_args, MaskImagePointer &mask, MatrixType &signal, MatrixType &grad_dirs) { - SIGNAL_GENERATOR sg(dmri_file); + if (!input_args->external_gradients.empty()) { + char answer; + cout << "External gradient file is found. Do you want to use gradients from it instead of the image's one? (y/n)" << endl; + cin >> answer; + if (answer == 'y') + DATA_SOURCE::fileGradientsToMatrix(input_args->external_gradients, grad_dirs); + else if (answer == 'n') + cout << "Volume's gradients will be used." << endl; + else { + cerr << "Incorrect answer. Rerun the program." << endl; + return EXIT_FAILURE; + } + } + + SIGNAL_GENERATOR sg(input_args->input_dmri); int res_dmri = sg.ExtractMatrix(mask, signal, grad_dirs); h = sg.h; if (res_dmri) diff --git a/Ridgelets/SphericalRidgelets/src/SIGNAL_GENERATOR.cpp b/Ridgelets/SphericalRidgelets/src/SIGNAL_GENERATOR.cpp index 9d4ee32..08ef743 100644 --- a/Ridgelets/SphericalRidgelets/src/SIGNAL_GENERATOR.cpp +++ b/Ridgelets/SphericalRidgelets/src/SIGNAL_GENERATOR.cpp @@ -9,7 +9,12 @@ bool SIGNAL_GENERATOR::is_path_exists(const string &s) } int SIGNAL_GENERATOR::readVolume(MatrixType & GradientDirections, DiffusionImagePointer & image) { - GradientDirections.resize(0, 3); + bool ext_grad_use = true; + if (GradientDirections.size() == 0) { + GradientDirections.resize(0, 3); + ext_grad_use = false; + } + bool is_b0 = false; precisionType b0 = 0; precisionType x, y, z; @@ -26,12 +31,12 @@ int SIGNAL_GENERATOR::readVolume(MatrixType & GradientDirections, DiffusionImage if (ext_vol.compare("nhdr")) { if (ext_vol.compare("nrrd")) { - cout << "NDHR or NRRD file formats only! Please, check file type." << endl; + cerr << "NDHR or NRRD file formats only! Please, check file type." << endl; return EXIT_FAILURE; } } if (!is_path_exists(inputVolume)) { - cout << "Input dMRI image is not exists! Please, check the path and file name." << endl; + cerr << "Input dMRI image is not exists! Please, check the path and file name." << endl; return EXIT_FAILURE; } @@ -60,43 +65,65 @@ int SIGNAL_GENERATOR::readVolume(MatrixType & GradientDirections, DiffusionImage vector::const_iterator itKey = imgMetaKeys.begin(); string metaString; - for (; itKey != imgMetaKeys.end(); ++itKey) - { - itk::ExposeMetaData(imgMetaDictionary, *itKey, metaString); - if (itKey->find("DWMRI_gradient") != string::npos) + // Get diffusion-encoding direction from dmri volume if an external file is not provided + if (!ext_grad_use) { + for (; itKey != imgMetaKeys.end(); ++itKey) { - //cout << *itKey << " -> " << metaString << endl; - istringstream ss(metaString.c_str()); - ss >> x; ss >> y; ss >> z; + itk::ExposeMetaData(imgMetaDictionary, *itKey, metaString); + if (itKey->find("DWMRI_gradient") != string::npos) + { + //cout << *itKey << " -> " << metaString << endl; + istringstream ss(metaString.c_str()); + ss >> x; ss >> y; ss >> z; + + ++nOfImgs; + // If the direction is 0.0, this is a reference image + if (x == sphZero && y == sphZero && z == sphZero) + continue; + GradientDirections.conservativeResize(GradientDirections.rows() + 1, GradientDirections.cols()); + GradientDirections.row(GradientDirections.rows() - 1) << x, y, z; + ++nGradImgs; + } + } + } + else { + for (int i = 0; i < GradientDirections.rows(); ++i) { ++nOfImgs; - // If the direction is 0.0, this is a reference image - if (x == sphZero && y == sphZero && z == sphZero) + if (GradientDirections(i, 0) == sphZero && GradientDirections(i, 1) == sphZero && GradientDirections(i, 2) == sphZero) continue; - - GradientDirections.conservativeResize(GradientDirections.rows() + 1, GradientDirections.cols()); - GradientDirections.row(GradientDirections.rows() - 1) << x, y, z; ++nGradImgs; } - else if (itKey->find("DWMRI_b-value") != string::npos) + for (int i = 0; i < GradientDirections.rows(); ++i) { + if (GradientDirections(i, 0) == sphZero && GradientDirections(i, 1) == sphZero && GradientDirections(i, 2) == sphZero) + remove_row(GradientDirections, i); + } + } + + // Get b-value from the dmri volume + for (; itKey != imgMetaKeys.end(); ++itKey) + { + itk::ExposeMetaData(imgMetaDictionary, *itKey, metaString); + if (itKey->find("DWMRI_b-value") != string::npos) { - //cout << *itKey << " -> " << metaString << endl; + // cout << *itKey << " -> " << metaString << endl; is_b0 = true; b0 = stod(metaString.c_str()); } } + // Normalize gradients GradientDirections.rowwise().normalize(); itk::Size<3> img_size = image->GetLargestPossibleRegion().GetSize(); + unsigned short int n_refs = nOfImgs - nGradImgs; cout << "Image size: (" << img_size[0] << " " << img_size[1] << " " << img_size[2] << "). "; - cout << "Number of gradient images: " << nGradImgs << ". Number of reference images: " << nOfImgs - nGradImgs << endl; - cout << "b-value " << b0 << endl; - - if (!is_b0) - { - cerr << "b-value not specified in file's header." << endl; - return EXIT_FAILURE; + cout << "Number of gradient images: " << nGradImgs << ". Number of reference images: " << n_refs << endl; + if (!is_b0) { + cout << "Be cautious! b-value is not specified in file's header. Sometimes it indicates that the input is corrupted." << endl; + } + else { + cout << "b-value " << b0 << endl; } return EXIT_SUCCESS; @@ -138,6 +165,7 @@ int SIGNAL_GENERATOR::ExtractMatrix(MaskImagePointer &mask, MatrixType &signal, unsigned first_grad_image_index = nOfImgs - nGradImgs; if (first_grad_image_index == 0) { cout << "Warning! There is no baseline image in the input file!" << endl; + cout << "Make sure that all diffusion-encoded volumes were normalized by an average b0!" << endl; avrg_b0 = 0.0; } signal = MatrixType::Zero(nGradImgs, N_of_voxels); @@ -246,3 +274,14 @@ int SIGNAL_GENERATOR::ExtractMatrix(MaskImagePointer &mask, MatrixType &signal, img = NULL; return EXIT_SUCCESS; } + +void SIGNAL_GENERATOR::remove_row(MatrixType& a, Eigen::Index del) +{ + unsigned cols = a.cols(); + unsigned rows = a.rows() - 1; + + if (del < rows) + a.block(del, 0, rows - del, cols) = a.block(del + 1, 0, rows - del, cols); + + a.conservativeResize(rows, cols); +} \ No newline at end of file diff --git a/Ridgelets/packages.config b/Ridgelets/packages.config index caa5e56..5e4558c 100644 --- a/Ridgelets/packages.config +++ b/Ridgelets/packages.config @@ -1,4 +1,4 @@  - + \ No newline at end of file