From c38b202c8bcad821a1009838e5a1f9e633263a3f Mon Sep 17 00:00:00 2001
From: Moises Fernandez <moisesf@fmrib.ox.ac.uk>
Date: Thu, 28 Feb 2013 16:12:38 +0000
Subject: [PATCH] Changes in grad-nonlin functions
---
CUDA/xfibres_gpu.cu | 25 +++++--------------------
CUDA/xfibres_gpu.cuh | 2 --
CUDA/xfibres_gpu.h | 2 --
xfibres_gpu.cc | 26 ++------------------------
4 files changed, 7 insertions(+), 48 deletions(-)
diff --git a/CUDA/xfibres_gpu.cu b/CUDA/xfibres_gpu.cu
index 1829acb..f0ccf62 100644
--- a/CUDA/xfibres_gpu.cu
+++ b/CUDA/xfibres_gpu.cu
@@ -31,8 +31,6 @@ void xfibres_gpu( //INPUT
const Matrix bvecs,
const Matrix bvals,
const Matrix gradm,
- const Matrix Qform,
- const Matrix Qform_inv,
const NEWIMAGE::volume<int> vol2matrixkey,
const NEWMAT::Matrix matrix2volkey,
const NEWIMAGE::volume<float> mask,
@@ -64,7 +62,7 @@ void xfibres_gpu( //INPUT
vector<Matrix> bvecs_vec, bvals_vec;
///// FIT /////
- prepare_data_gpu_FIT(datam,bvecs,bvals,gradm,Qform,Qform_inv,datam_vec, bvecs_vec, bvals_vec, datam_host, bvecs_host, bvals_host, alpha_host, beta_host, params_host, tau_host);
+ prepare_data_gpu_FIT(datam,bvecs,bvals,gradm,datam_vec, bvecs_vec, bvals_vec, datam_host, bvecs_host, bvals_host, alpha_host, beta_host, params_host, tau_host);
thrust::device_vector<double> datam_gpu=datam_host;
thrust::device_vector<double> bvecs_gpu=bvecs_host;
@@ -125,7 +123,7 @@ void xfibres_gpu( //INPUT
// Correct bvals/bvecs accounting for Gradient Nonlinearities
// ColumnVector grad_nonlin has 9 entries, corresponding to the 3 components of each of the x,y and z gradient deviation
-void correct_bvals_bvecs(const Matrix& bvals,const Matrix& bvecs, const ColumnVector& grad_nonlin, const Matrix& Qform, const Matrix& Qform_inv, Matrix& bvals_c, Matrix& bvecs_c){
+void correct_bvals_bvecs(const Matrix& bvals,const Matrix& bvecs, const ColumnVector& grad_nonlin, Matrix& bvals_c, Matrix& bvecs_c){
bvals_c=bvals; bvecs_c=bvecs;
Matrix L(3,3); //gradient coil tensor
float mag;
@@ -137,22 +135,11 @@ void correct_bvals_bvecs(const Matrix& bvals,const Matrix& bvecs, const ColumnVe
for (int l=1; l<=bvals.Ncols(); l++){
if (bvals(1,l)>0){ //do not correct b0s
- //Rotate bvecs to scanner's coordinate system
- ColumnVector bvec_tmp(3);
- bvec_tmp=Qform*bvecs.Column(l);
- bvec_tmp(1)=-bvec_tmp(1); //Sign-flip X coordinate
-
- //Correct for grad-nonlin in scanner's coordinate system
- bvecs_c.Column(l)=(Id+L)*bvec_tmp;//bvecs.Column(l);
+ bvecs_c.Column(l)=(Id+L)*bvecs.Column(l);
mag=sqrt(bvecs_c(1,l)*bvecs_c(1,l)+bvecs_c(2,l)*bvecs_c(2,l)+bvecs_c(3,l)*bvecs_c(3,l));
if (mag!=0)
bvecs_c.Column(l)=bvecs_c.Column(l)/mag;
- bvals_c(1,l)=mag*mag*bvals(1,l);
- bvec_tmp=bvecs_c.Column(l);
-
- //Rotate corrected bvecs back to voxel coordinate system
- bvec_tmp(1)=-bvec_tmp(1); //Sign-flip X coordinate
- bvecs_c.Column(l)=Qform_inv*bvec_tmp;
+ bvals_c(1,l)=mag*mag*bvals(1,l);//mag^2 as b propto |G|^2
}
}
}
@@ -335,8 +322,6 @@ void prepare_data_gpu_FIT( //INPUT
const Matrix bvecs,
const Matrix bvals,
const Matrix gradm,
- const Matrix Qform,
- const Matrix Qform_inv,
//OUTPUT
vector<ColumnVector>& datam_vec,
vector<Matrix>& bvecs_vec,
@@ -373,7 +358,7 @@ void prepare_data_gpu_FIT( //INPUT
if (opts.grad_file.set()){
for(int vox=0;vox<nvox;vox++){
- correct_bvals_bvecs(bvals,bvecs, gradm.Column(vox+1),Qform,Qform_inv,bvals_vec[vox],bvecs_vec[vox]); //correct for gradient nonlinearities
+ correct_bvals_bvecs(bvals,bvecs, gradm.Column(vox+1),bvals_vec[vox],bvecs_vec[vox]); //correct for gradient nonlinearities
MISCMATHS::cart2sph(bvecs_vec[vox],alpha,beta);
for(int dir=0;dir<NDIRECTIONS;dir++){
diff --git a/CUDA/xfibres_gpu.cuh b/CUDA/xfibres_gpu.cuh
index 05dae84..4dce61a 100644
--- a/CUDA/xfibres_gpu.cuh
+++ b/CUDA/xfibres_gpu.cuh
@@ -35,8 +35,6 @@ void prepare_data_gpu_FIT( //INPUT
const Matrix bvecs,
const Matrix bvals,
const Matrix gradm,
- const Matrix Qform,
- const Matrix Qform_inv,
//OUTPUT
vector<ColumnVector>& datam_vec,
vector<Matrix>& bvecs_vec,
diff --git a/CUDA/xfibres_gpu.h b/CUDA/xfibres_gpu.h
index 309fcdd..6089189 100644
--- a/CUDA/xfibres_gpu.h
+++ b/CUDA/xfibres_gpu.h
@@ -17,8 +17,6 @@ void xfibres_gpu( //INPUT
const NEWMAT::Matrix bvals,
const NEWMAT::Matrix bvecs,
const NEWMAT::Matrix gradm,
- const NEWMAT::Matrix Qform,
- const NEWMAT::Matrix Qform_inv,
const NEWIMAGE::volume<int> vol2matrixkey,
const NEWMAT::Matrix matrix2volkey,
const NEWIMAGE::volume<float> mask,
diff --git a/xfibres_gpu.cc b/xfibres_gpu.cc
index c965b04..67b769d 100644
--- a/xfibres_gpu.cc
+++ b/xfibres_gpu.cc
@@ -6,10 +6,6 @@
/* CCOPYRIGHT */
-#ifndef EXPOSE_TREACHEROUS
-#define EXPOSE_TREACHEROUS
-#endif
-
#include "CUDA/xfibres_gpu.h"
#include "xfibresoptions.h"
#include "newmat.h"
@@ -22,9 +18,6 @@
using namespace Xfibres;
-
-void Return_Qform(const Matrix& qform_mat, Matrix& QMat, const float xdim, const float ydim, const float zdim);
-
//////////////////////////////////////////////////////////
// XFIBRES CPU PART. IT CALLS TO GPU PART
//////////////////////////////////////////////////////////
@@ -114,13 +107,10 @@ int xfibres_gpu(char *subjdir,int slice, int nextras, char** extras)
//Samples samples(vol2matrixkey,matrix2volkey,datam.Ncols(),datam.Nrows()); //in xfibres_gpu.cu
//Read Gradient Non_linearity Maps if provided
- NEWIMAGE::volume4D<float> grad; Matrix gradm, Qform, Qform_inv;
+ NEWIMAGE::volume4D<float> grad; Matrix gradm;
if (opts.grad_file.set()){
read_volume4D(grad,opts.grad_file.value());
gradm=grad.matrix(mask);
- //Get the scale-free Qform matrix to rotate bvecs back to scanner's coordinate system
- Return_Qform(data.qform_mat(), Qform, data.xdim(), data.ydim(), data.zdim());
- Qform_inv=Qform.i();
}
myfile << "Number of Voxels: " << datam.Ncols() << endl;
@@ -129,7 +119,7 @@ int xfibres_gpu(char *subjdir,int slice, int nextras, char** extras)
if(opts.rician.value() && !opts.nonlin.value())
cout<<"Rician noise model requested. Non-linear parameter initialization will be performed, overriding other initialization options!"<<endl;
- xfibres_gpu(datam,bvecs,bvals,gradm,Qform,Qform_inv,vol2matrixkey,matrix2volkey,mask,slice,subjdir);
+ xfibres_gpu(datam,bvecs,bvals,gradm,vol2matrixkey,matrix2volkey,mask,slice,subjdir);
gettimeofday(&t2,NULL);
time=timeval_diff(&t2,&t1);
@@ -139,15 +129,3 @@ int xfibres_gpu(char *subjdir,int slice, int nextras, char** extras)
return 0;
}
-
-
-//Get the scale-free Qform matrix to rotate bvecs back to scanner's coordinate system
-//After applying this matrix, make sure to sign flip the x coordinate of the resulted bvecs!
-//If you apply the inverse of the matrix, sign flip the x coordinate of the input bvecs
-void Return_Qform(const Matrix& qform_mat, Matrix& QMat, const float xdim, const float ydim, const float zdim){
- Matrix QMat_tmp; DiagonalMatrix Scale(3);
- QMat_tmp=qform_mat.SubMatrix(1,3,1,3);
- Scale(1)=xdim; Scale(2)=ydim; Scale(3)=zdim;
- QMat_tmp=Scale.i()*QMat_tmp;
- QMat=QMat_tmp;
-}
--
GitLab