Simplify requirements for dmri input. Add support of external grad di…

…rs file as a main table

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<precisionType, MatrixType, VectorType> 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<seconds>(stop - start);
 		auto dm = duration_cast<minutes>(stop - start);

Ridgelets/Ridgelets.vcxproj

@@ -83,8 +83,8 @@
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
     <LinkIncremental>false</LinkIncremental>
-    <IncludePath>C:\Users\mukho\source\repos\SphericalRidgelet\Ridgelets\SphericalRidgelets\include;$(IncludePath)</IncludePath>
-    <SourcePath>C:\Users\mukho\source\repos\SphericalRidgelet\Ridgelets\SphericalRidgelets\src;$(SourcePath)</SourcePath>
+    <IncludePath>D:\VSProjects\SphericalRidgelets\Ridgelets\SphericalRidgelets\include;$(IncludePath)</IncludePath>
+    <SourcePath>D:\VSProjects\SphericalRidgelets\Ridgelets\SphericalRidgelets\src;$(SourcePath)</SourcePath>
   </PropertyGroup>
   <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
     <ClCompile>
@@ -197,12 +197,12 @@
   </ItemGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
   <ImportGroup Label="ExtensionTargets">
-    <Import Project="..\packages\Eigen.3.3.3\build\native\Eigen.targets" Condition="Exists('..\packages\Eigen.3.3.3\build\native\Eigen.targets')" />
+    <Import Project="..\packages\Eigen3.3.3.9\build\native\Eigen3.targets" Condition="Exists('..\packages\Eigen3.3.3.9\build\native\Eigen3.targets')" />
   </ImportGroup>
   <Target Name="EnsureNuGetPackageBuildImports" BeforeTargets="PrepareForBuild">
     <PropertyGroup>
       <ErrorText>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}.</ErrorText>
     </PropertyGroup>
-    <Error Condition="!Exists('..\packages\Eigen.3.3.3\build\native\Eigen.targets')" Text="$([System.String]::Format('$(ErrorText)', '..\packages\Eigen.3.3.3\build\native\Eigen.targets'))" />
+    <Error Condition="!Exists('..\packages\Eigen3.3.3.9\build\native\Eigen3.targets')" Text="$([System.String]::Format('$(ErrorText)', '..\packages\Eigen3.3.3.9\build\native\Eigen3.targets'))" />
   </Target>
 </Project>

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);

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;

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)

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<string>::const_iterator itKey = imgMetaKeys.begin();
 	string metaString;
 
-	for (; itKey != imgMetaKeys.end(); ++itKey)
-	{
-		itk::ExposeMetaData<string>(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<string>(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<string>(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);
+}

Ridgelets/packages.config

@@ -1,4 +1,4 @@
 <?xml version="1.0" encoding="utf-8"?>
 <packages>
-  <package id="Eigen" version="3.3.3" targetFramework="native" />
+  <package id="Eigen3" version="3.3.9" targetFramework="native" />
 </packages>