kurtosis v.0.0

kurtosis.cc

+/*  Copyright (C) 2008 University of Oxford  */
+
+/* S.Jbabdi */
+
+/*  CCOPYRIGHT  */
+
+#include <iostream>
+#include <cmath>
+#include "miscmaths/miscmaths.h"
+#include "miscmaths/minimize.h"
+#include "newmat.h"
+#include "newimage/newimageall.h"
+#include "dtifitOptions.h"
+
+using namespace std;
+using namespace NEWMAT;
+using namespace MISCMATHS;
+using namespace NEWIMAGE;
+using namespace DTIFIT;
+
+const float maxfloat=1e10;
+const float minfloat=1e-10;
+const float maxlogfloat=23;
+const float minlogfloat=-23;
+const int maxint=1000000000; 
+
+
+ReturnMatrix form_Kmat(const Matrix& r){
+  Matrix K(r.Ncols(),15);
+
+  for(int j=1;j<=r.Ncols();j++){
+    float x=r(1,j),y=r(2,j),z=r(3,j);
+    K(j,1) = MISCMATHS::pow(x,4);
+    K(j,2) = MISCMATHS::pow(y,4);
+    K(j,3) = MISCMATHS::pow(z,4);
+    K(j,4) = 4*MISCMATHS::pow(x,3)*y;
+    K(j,5) = 4*MISCMATHS::pow(x,3)*z;
+    K(j,6) = 4*MISCMATHS::pow(y,3)*x;
+    K(j,7) = 4*MISCMATHS::pow(y,3)*z;
+    K(j,8) = 4*MISCMATHS::pow(z,3)*x;
+    K(j,9) = 4*MISCMATHS::pow(z,3)*y;
+    K(j,10) = 6*MISCMATHS::pow(x,2)*MISCMATHS::pow(y,2);
+    K(j,11) = 6*MISCMATHS::pow(x,2)*MISCMATHS::pow(z,2);
+    K(j,12) = 6*MISCMATHS::pow(y,2)*MISCMATHS::pow(z,2);
+    K(j,13) = 12*MISCMATHS::pow(x,2)*y*z;
+    K(j,14) = 12*MISCMATHS::pow(y,2)*x*z;
+    K(j,15) = 12*MISCMATHS::pow(z,2)*x*y;
+    j+=1;
+  }
+  K.Release();
+  return K;
+}
+
+
+// note the order of the variable parameters
+// D11,D12,D13,D22,D23,D33,logS0
+// W1111,W2222,W333,W1112,W1113,W1222,W2223,W1333,
+// W2333,W1122,W1133,W2233,W1123,W1223,W1233
+class KurtosisNonlinCF : public gEvalFunction
+{
+protected:
+  ColumnVector m_A;
+  ColumnVector m_B;
+  Matrix       m_C;
+  Matrix       m_D;
+  int          m_n;
+
+public:
+  KurtosisNonlinCF(const ColumnVector& data,const Matrix& bvals,const Matrix& bvecs):gEvalFunction()
+  {
+
+    m_n = data.Nrows();
+
+    m_A.ReSize(m_n);
+    m_B.ReSize(m_n);
+    m_C.ReSize(m_n,6);
+    m_D.ReSize(m_n,15);
+
+    Matrix K = form_Kmat(bvecs);
+
+    for (int i=1;i<=m_n;i++){
+      if(data(i)>0){
+	m_A(i)=-log(data(i));
+      }
+      else{
+	m_A(i)=0;
+      }
+      m_B(i) = 1.0;
+      
+      m_C(i,1) = -bvals(1,i)*bvecs(1,i)*bvecs(1,i);
+      m_C(i,2) = -2*bvals(1,i)*bvecs(1,i)*bvecs(2,i);
+      m_C(i,3) = -2*bvals(1,i)*bvecs(1,i)*bvecs(3,i);
+      m_C(i,4) = -bvals(1,i)*bvecs(2,i)*bvecs(2,i);
+      m_C(i,5) = -2*bvals(1,i)*bvecs(2,i)*bvecs(3,i);
+      m_C(i,6) = -bvals(1,i)*bvecs(3,i)*bvecs(3,i);
+
+      for(int j=1;j<=15;j++)
+	m_D(i,j) = (bvals(1,i)*bvals(1,i)/6) * K(i,j);
+
+    }
+    
+  }
+
+  virtual ~KurtosisNonlinCF(){};
+
+  float evaluate(const ColumnVector& x) const{
+    float res=0;
+    res = ( m_A + m_B*x(7) + m_C*x.SubMatrix(1,6,1,1)
+	    + m_D*x.SubMatrix(8,22,1,1)*(x(1)+x(4)+x(6))*(x(1)+x(4)+x(6))/9 ).SumSquare();
+
+    return res;
+  }
+  ReturnMatrix g_evaluate(const ColumnVector& x) const{
+    ColumnVector sj_g(x.Nrows());
+    //    ColumnVector sj_gg;
+    
+    //    sj_gg = MISCMATHS::gradient(x,*this,1e-4);
+ 
+    ColumnVector sj_d(6);
+    ColumnVector sj_w(15);
+    sj_d = x.SubMatrix(1,6,1,1);
+    sj_w = x.SubMatrix(8,22,1,1);
+    double sj_t = x(1)+x(4)+x(6);
+    double sj_t2=sj_t*sj_t;
+
+    ColumnVector sj_func(m_n);
+    sj_func = m_A + m_B*x(7) + m_C*sj_d + m_D*sj_w*sj_t*sj_t/9;
+
+    sj_g(1) = 2*NEWMAT::SP(sj_func,m_C.SubMatrix(1,m_n,1,1)+2*sj_t/9*m_D*sj_w).Sum();
+    sj_g(2) = 2*NEWMAT::SP(sj_func,m_C.SubMatrix(1,m_n,2,2)).Sum();
+    sj_g(3) = 2*NEWMAT::SP(sj_func,m_C.SubMatrix(1,m_n,3,3)).Sum();
+    sj_g(4) = 2*NEWMAT::SP(sj_func,m_C.SubMatrix(1,m_n,4,4)+2*sj_t/9*m_D*sj_w).Sum();
+    sj_g(5) = 2*NEWMAT::SP(sj_func,m_C.SubMatrix(1,m_n,5,5)).Sum();
+    sj_g(6) = 2*NEWMAT::SP(sj_func,m_C.SubMatrix(1,m_n,6,6)+2*sj_t/9*m_D*sj_w).Sum();
+
+    sj_g(7) = 2*NEWMAT::SP(sj_func,m_B).Sum();
+ 
+    for(int sj_i=1,sj_j=8;sj_j<=x.Nrows();sj_i++,sj_j++)
+      sj_g(sj_j) = 2*NEWMAT::SP(sj_func,sj_t2/9*m_D.SubMatrix(1,m_n,sj_i,sj_i)).Sum();
+
+    sj_g.Release();
+    return sj_g;
+
+  }
+
+  const KurtosisNonlinCF& operator=(const KurtosisNonlinCF& par)
+  {
+    m_A   = par.m_A;
+    m_B   = par.m_B;
+    m_C   = par.m_C;
+    m_D   = par.m_D;
+    m_n = par.m_n;
+
+    return *this;
+    
+  }
+  KurtosisNonlinCF(const KurtosisNonlinCF& rhs):
+    m_A(rhs.m_A),m_B(rhs.m_B),m_C(rhs.m_C),m_D(rhs.m_D),m_n(rhs.m_n){
+    *this=rhs;
+  }
+
+};
+
+inline float PI() { return  3.14159265358979;}
+inline float min(float a,float b){
+  return a<b ? a:b;}
+inline float max(float a,float b){
+  return a>b ? a:b;}
+inline Matrix Anis()
+{ 
+  Matrix A(3,3);
+  A << 1 << 0 << 0
+    << 0 << 0 << 0
+    << 0 << 0 << 0;
+  return A;
+}
+
+inline Matrix Is()
+{ 
+  Matrix I(3,3);
+  I << 1 << 0 << 0
+    << 0 << 1 << 0
+    << 0 << 0 << 1;
+  return I;
+}
+
+inline ColumnVector Cross(const ColumnVector& A,const ColumnVector& B)
+{
+  ColumnVector res(3);
+  res << A(2)*B(3)-A(3)*B(2)
+      << A(3)*B(1)-A(1)*B(3)
+      << A(1)*B(2)-B(1)*A(2);
+  return res;
+}
+
+inline Matrix Cross(const Matrix& A,const Matrix& B)
+{
+  Matrix res(3,1);
+  res << A(2,1)*B(3,1)-A(3,1)*B(2,1)
+      << A(3,1)*B(1,1)-A(1,1)*B(3,1)
+      << A(1,1)*B(2,1)-B(1,1)*A(2,1);
+  return res;
+}
+
+float mod(float a, float b){
+  while(a>b){a=a-b;}
+  while(a<0){a=a+b;}
+  return a;
+}
+
+
+
+Matrix form_Amat(const Matrix& r,const Matrix& b)
+{
+  Matrix A(r.Ncols(),7);
+  Matrix tmpvec(3,1), tmpmat;
+  
+  for( int i = 1; i <= r.Ncols(); i++){
+    tmpvec << r(1,i) << r(2,i) << r(3,i);
+    tmpmat = tmpvec*tmpvec.t()*b(1,i);
+    A(i,1) = tmpmat(1,1);
+    A(i,2) = 2*tmpmat(1,2);
+    A(i,3) = 2*tmpmat(1,3);
+    A(i,4) = tmpmat(2,2);
+    A(i,5) = 2*tmpmat(2,3);
+    A(i,6) = tmpmat(3,3);
+    A(i,7) = 1;
+  }
+  return A;
+}
+inline SymmetricMatrix vec2tens(ColumnVector& Vec){
+  SymmetricMatrix tens(3);
+  tens(1,1)=Vec(1);
+  tens(2,1)=Vec(2);
+  tens(3,1)=Vec(3);
+  tens(2,2)=Vec(4);
+  tens(3,2)=Vec(5);
+  tens(3,3)=Vec(6);
+  return tens;
+}
+
+void tensorfit(DiagonalMatrix& Dd,ColumnVector& evec1,ColumnVector& evec2,ColumnVector& evec3,float& f,float& s0,ColumnVector& Dvec, const Matrix& Amat,const ColumnVector& S)
+{
+  //Initialise the parameters using traditional DTI analysis
+  ColumnVector logS(S.Nrows());
+  SymmetricMatrix tens;   //Basser's Diffusion Tensor;
+  //  DiagonalMatrix Dd;   //eigenvalues
+  Matrix Vd;   //eigenvectors
+  DiagonalMatrix Ddsorted(3);
+  float mDd, fsquared;
+
+  for ( int i = 1; i <= S.Nrows(); i++)
+    {
+      if(S(i)>0){
+	logS(i)=log(S(i));
+      }
+      else{
+	logS(i)=0;
+      }
+      //      logS(i)=(S(i)/S0)>0.01 ? log(S(i))-log(S0):log(0.01);
+    }
+
+  Dvec = -pinv(Amat)*logS;
+  if(  Dvec(7) >  -maxlogfloat ){
+    s0=exp(-Dvec(7));
+  }
+  else{
+    s0=S.MaximumAbsoluteValue();
+  }
+  for ( int i = 1; i <= S.Nrows(); i++)
+    {
+      if(s0<S.Sum()/S.Nrows()){ s0=S.MaximumAbsoluteValue();  }
+      logS(i)=(S(i)/s0)>0.01 ? log(S(i)):log(0.01*s0);
+    }
+  Dvec = -pinv(Amat)*logS;
+  s0=exp(-Dvec(7));
+  if(s0<S.Sum()/S.Nrows()){ s0=S.Sum()/S.Nrows();  }
+  tens = vec2tens(Dvec);
+  EigenValues(tens,Dd,Vd);
+  mDd = Dd.Sum()/Dd.Nrows();
+  int maxind = Dd(1) > Dd(2) ? 1:2;   //finding max,mid and min eigenvalues
+  maxind = Dd(maxind) > Dd(3) ? maxind:3;
+  int midind;
+  if( (Dd(1)>=Dd(2) && Dd(2)>=Dd(3)) || (Dd(1)<=Dd(2) && Dd(2)<=Dd(3)) ){midind=2;}
+  else if( (Dd(2)>=Dd(1) && Dd(1)>=Dd(3)) || (Dd(2)<=Dd(1) && Dd(1)<=Dd(3)) ){midind=1;}
+  else {midind=3;}
+  int minind = Dd(1) < Dd(2) ? 1:2;   //finding maximum eigenvalue
+  minind = Dd(minind) < Dd(3) ? minind:3;
+  Ddsorted << Dd(maxind) << Dd(midind) << Dd(minind);
+  Dd=Ddsorted;
+  evec1 << Vd(1,maxind) << Vd(2,maxind) << Vd(3,maxind);
+  evec2 << Vd(1,midind) << Vd(2,midind) << Vd(3,midind);
+  evec3 << Vd(1,minind) << Vd(2,minind) << Vd(3,minind);
+
+  float numer=1.5*((Dd(1)-mDd)*(Dd(1)-mDd)+(Dd(2)-mDd)*(Dd(2)-mDd)+(Dd(3)-mDd)*(Dd(3)-mDd));
+  float denom=(Dd(1)*Dd(1)+Dd(2)*Dd(2)+Dd(3)*Dd(3));
+ 
+  if(denom>0) fsquared=numer/denom;
+  else fsquared=0;
+  if(fsquared>0){f=sqrt(fsquared);}
+  else{f=0;}
+
+}
+
+void kurtosisfit(DiagonalMatrix& Dd,ColumnVector& evec1,ColumnVector& evec2, ColumnVector& evec3,
+		 float& f,float& s0,ColumnVector& Dvec, float& mk, ColumnVector& tens4, 
+		 const Matrix& Amat,const Matrix& Kmat,const ColumnVector& S,const Matrix& bvals,const Matrix& bvecs){
+
+  // initialise second-order tensor with simple tensor fit
+  // tensorfit(Dd,evec1,evec2,evec3,f,s0,Dvec,Amat,S);
+  //SymmetricMatrix tens;
+  //tens = vec2tens(Dvec);
+
+//   // initialise Kurtosis using Linear fit
+//   ColumnVector v(S.Nrows());
+//   for(int i=1;i<=S.Nrows();i++){
+//   float bDi = bvals(1,i)*(bvecs.Column(i).t()*tens*bvecs.Column(i)).AsScalar();
+//   if(bDi>0)
+//     v(i) = 6*(log(S(i)/s0)+bDi)/(bDi*bDi);
+//   else
+//     v(i) = 0;
+//   }
+//   tens4 = pinv(Kmat) * v;
+
+  // calculate DT and KT using non-linear fitting
+  KurtosisNonlinCF KNL(S,bvals,bvecs);
+  ColumnVector xmin(22);
+  KNL.minimize(xmin);
+
+
+  Dvec.SubMatrix(1,6,1,1) = xmin.SubMatrix(1,6,1,1);
+  tens4 = xmin.SubMatrix(8,22,1,1);
+  Dvec(7) = exp(xmin(7));  
+  s0 = Dvec(7);
+
+  // Tensor Stuff
+  float mDd, fsquared;
+  SymmetricMatrix tens;
+  DiagonalMatrix Ddsorted(3);
+  Matrix Vd;
+  tens = vec2tens(Dvec);
+  EigenValues(tens,Dd,Vd);
+  mDd = Dd.Sum()/Dd.Nrows();
+  int maxind = Dd(1) > Dd(2) ? 1:2;   //finding max,mid and min eigenvalues
+  maxind = Dd(maxind) > Dd(3) ? maxind:3;
+  int midind;
+  if( (Dd(1)>=Dd(2) && Dd(2)>=Dd(3)) || (Dd(1)<=Dd(2) && Dd(2)<=Dd(3)) ){midind=2;}
+  else if( (Dd(2)>=Dd(1) && Dd(1)>=Dd(3)) || (Dd(2)<=Dd(1) && Dd(1)<=Dd(3)) ){midind=1;}
+  else {midind=3;}
+  int minind = Dd(1) < Dd(2) ? 1:2;   //finding maximum eigenvalue
+  minind = Dd(minind) < Dd(3) ? minind:3;
+  Ddsorted << Dd(maxind) << Dd(midind) << Dd(minind);
+  Dd=Ddsorted;
+  evec1 << Vd(1,maxind) << Vd(2,maxind) << Vd(3,maxind);
+  evec2 << Vd(1,midind) << Vd(2,midind) << Vd(3,midind);
+  evec3 << Vd(1,minind) << Vd(2,minind) << Vd(3,minind);
+  float numer=1.5*((Dd(1)-mDd)*(Dd(1)-mDd)+(Dd(2)-mDd)*(Dd(2)-mDd)+(Dd(3)-mDd)*(Dd(3)-mDd));
+  float denom=(Dd(1)*Dd(1)+Dd(2)*Dd(2)+Dd(3)*Dd(3));
+  if(denom>0) fsquared=numer/denom;
+  else fsquared=0;
+  if(fsquared>0){f=sqrt(fsquared);}
+  else{f=0;}
+
+  // Kurtosis Stuff
+  mk = 0;
+  ColumnVector vec(S.Nrows());
+  vec = Kmat*tens4;
+  for(int i=1;i<=S.Nrows();i++){
+    if(bvals(1,i)>0)
+      mk+=vec(i)/(bvecs.Column(i).t()*tens*bvecs.Column(i)).AsScalar();
+  }
+  mk *= mDd*mDd;
+  //OUT(mk);
+}
+
+int main(int argc, char** argv)
+{
+  //parse command line
+  dtifitOptions& opts = dtifitOptions::getInstance();
+  int success=opts.parse_command_line(argc,argv);
+  if(!success) return 1;
+
+  if(opts.verbose.value()){
+    cout<<"data file "<<opts.dtidatafile.value()<<endl;
+    cout<<"mask file "<<opts.maskfile.value()<<endl;
+    cout<<"bvecs     "<<opts.bvecsfile.value()<<endl;
+    cout<<"bvals     "<<opts.bvalsfile.value()<<endl;
+    if(opts.littlebit.value()){
+      cout<<"min z     "<<opts.z_min.value()<<endl;
+      cout<<"max z     "<<opts.z_max.value()<<endl;
+      cout<<"min y     "<<opts.y_min.value()<<endl;
+      cout<<"max y     "<<opts.y_max.value()<<endl;
+      cout<<"min x     "<<opts.x_min.value()<<endl;
+      cout<<"max x     "<<opts.x_max.value()<<endl;
+    }
+  }
+  /////////////////////////////////////////
+  // read bvecs and bvals
+  // correct transpose and normalise bvecs
+  Matrix r = read_ascii_matrix(opts.bvecsfile.value());
+  if(r.Nrows()>3) r=r.t();
+  for(int i=1;i<=r.Ncols();i++){
+    float tmpsum=sqrt(r(1,i)*r(1,i)+r(2,i)*r(2,i)+r(3,i)*r(3,i));
+    if(tmpsum!=0){
+      r(1,i)=r(1,i)/tmpsum;
+      r(2,i)=r(2,i)/tmpsum;
+      r(3,i)=r(3,i)/tmpsum;
+    }  
+  }
+  Matrix b = read_ascii_matrix(opts.bvalsfile.value());
+  if(b.Nrows()>1) b=b.t();
+  //////////////////////////////////////////
+
+  volume4D<float> data;
+  volume<int> mask;
+  volumeinfo tempinfo;
+  if(opts.verbose.value()) cout<<"reading data"<<endl;
+  read_volume4D(data,opts.dtidatafile.value(),tempinfo);
+  if(opts.verbose.value()) cout<<"reading mask"<<endl;
+  read_volume(mask,opts.maskfile.value());
+  if(opts.verbose.value()) cout<<"ok"<<endl;
+  int minx=opts.littlebit.value() ? opts.x_min.value():0;
+  int maxx=opts.littlebit.value() ? opts.x_max.value():mask.xsize();
+  int miny=opts.littlebit.value() ? opts.y_min.value():0;
+  int maxy=opts.littlebit.value() ? opts.y_max.value():mask.ysize();
+  int minz=opts.littlebit.value() ? opts.z_min.value():0;
+  int maxz=opts.littlebit.value() ? opts.z_max.value():mask.zsize();
+  cout<<minx<<" "<<maxx<<" "<<miny<<" "<<maxy<<" "<<minz<<" "<<maxz<<endl;
+  if(opts.verbose.value()) cout<<"setting up vols"<<endl;
+
+  volume<float>   l1(maxx-minx,maxy-miny,maxz-minz);
+  volume<float>   l2(maxx-minx,maxy-miny,maxz-minz);
+  volume<float>   l3(maxx-minx,maxy-miny,maxz-minz);
+  volume<float>   MD(maxx-minx,maxy-miny,maxz-minz);
+  volume<float>   FA(maxx-minx,maxy-miny,maxz-minz);
+  volume<float>   S0(maxx-minx,maxy-miny,maxz-minz);
+  volume4D<float> V1(maxx-minx,maxy-miny,maxz-minz,3);
+  volume4D<float> V2(maxx-minx,maxy-miny,maxz-minz,3);
+  volume4D<float> V3(maxx-minx,maxy-miny,maxz-minz,3);
+  volume4D<float> Delements(maxx-minx,maxy-miny,maxz-minz,6);
+  volume<float>   MK(maxx-minx,maxy-miny,maxz-minz);
+  volume4D<float> KurtTens(maxx-minx,maxy-miny,maxz-minz,15);
+
+  if(opts.verbose.value()) cout<<"copying input properties to output volumes"<<endl;
+  copybasicproperties(data[0],l1);
+  copybasicproperties(data[0],l2);
+  copybasicproperties(data[0],l3);
+  copybasicproperties(data[0],MD);
+  copybasicproperties(data[0],FA);
+  copybasicproperties(data[0],S0);
+  copybasicproperties(data[0],V1[0]);
+  copybasicproperties(data[0],V2[0]);
+  copybasicproperties(data[0],V3[0]);
+  copybasicproperties(data[0],Delements[0]);
+  copybasicproperties(data[0],MK);
+  copybasicproperties(data[0],KurtTens[0]);
+
+  if(opts.verbose.value()) cout<<"zeroing output volumes"<<endl;
+  l1=0;l2=0;l3=0;MD=0;FA=0;S0=0;V1=0;V2=0;V3=0;Delements=0;MK=0;KurtTens=0;
+  if(opts.verbose.value()) cout<<"ok"<<endl;
+
+
+  DiagonalMatrix evals(3);
+  ColumnVector evec1(3),evec2(3),evec3(3);
+  ColumnVector tens4(15);
+  ColumnVector S(data.tsize());
+  float fa,s0,mk;
+  if(opts.verbose.value()) cout<<"Forming A matrix"<<endl;
+  Matrix Amat = form_Amat(r,b);
+  Matrix Kmat = form_Kmat(r);
+  if(opts.verbose.value()) cout<<"starting the fits"<<endl;
+  ColumnVector Dvec(7); Dvec=0;
+  for(int k = minz; k < maxz; k++){
+    cout<<k<<" slices processed"<<endl;
+      for(int j=miny; j < maxy; j++){
+	for(int i =minx; i< maxx; i++){
+	
+	  if(mask(i,j,k)>0){
+	    
+	    for(int t=0;t < data.tsize();t++){
+	      S(t+1)=data(i,j,k,t);
+	    }
+	    //tensorfit(evals,evec1,evec2,evec3,fa,s0,Dvec,Amat,S);
+	    kurtosisfit(evals,evec1,evec2,evec3,fa,s0,Dvec,mk,tens4,Amat,Kmat,S,b,r);
+
+
+	    l1(i-minx,j-miny,k-minz)=evals(1);
+	    l2(i-minx,j-miny,k-minz)=evals(2);
+	    l3(i-minx,j-miny,k-minz)=evals(3);
+	    MD(i-minx,j-miny,k-minz)=(evals(1)+evals(2)+evals(3))/3;
+	    FA(i-minx,j-miny,k-minz)=fa;
+	    S0(i-minx,j-miny,k-minz)=s0;
+	    V1(i-minx,j-miny,k-minz,0)=evec1(1);
+	    V1(i-minx,j-miny,k-minz,1)=evec1(2);
+	    V1(i-minx,j-miny,k-minz,2)=evec1(3);
+	    V2(i-minx,j-miny,k-minz,0)=evec2(1);
+	    V2(i-minx,j-miny,k-minz,1)=evec2(2);
+	    V2(i-minx,j-miny,k-minz,2)=evec2(3);
+	    V3(i-minx,j-miny,k-minz,0)=evec3(1);
+	    V3(i-minx,j-miny,k-minz,1)=evec3(2);
+	    V3(i-minx,j-miny,k-minz,2)=evec3(3);
+	    Delements(i-minx,j-miny,k-minz,0)=Dvec(1);
+	    Delements(i-minx,j-miny,k-minz,1)=Dvec(2);
+	    Delements(i-minx,j-miny,k-minz,2)=Dvec(3);
+	    Delements(i-minx,j-miny,k-minz,3)=Dvec(4);
+	    Delements(i-minx,j-miny,k-minz,4)=Dvec(5);
+	    Delements(i-minx,j-miny,k-minz,5)=Dvec(6);
+
+	    MK(i-minx,j-miny,k-minz)=mk;
+	    for(int iii=0;iii<15;iii++)
+	      KurtTens(i-minx,j-miny,k-minz,iii) = tens4(iii+1);
+
+	  }
+	}
+      }
+  }
+  
+    string fafile=opts.ofile.value()+"_FA";
+    string s0file=opts.ofile.value()+"_S0";
+    string l1file=opts.ofile.value()+"_L1";
+    string l2file=opts.ofile.value()+"_L2";
+    string l3file=opts.ofile.value()+"_L3";
+    string v1file=opts.ofile.value()+"_V1";
+    string v2file=opts.ofile.value()+"_V2";
+    string v3file=opts.ofile.value()+"_V3";
+    string MDfile=opts.ofile.value()+"_MD";
+    string tensfile=opts.ofile.value()+"_tensor";
+    string MKfile=opts.ofile.value()+"_MK";
+    string kurtosisfile=opts.ofile.value()+"_kurtosis";
+    if(opts.littlebit.value()){
+      fafile+="littlebit";
+      s0file+="littlebit";
+      l1file+="littlebit";
+      l2file+="littlebit";
+      l3file+="littlebit";
+      v1file+="littlebit";
+      v2file+="littlebit";
+      v3file+="littlebit";
+      MDfile+="littlebit";
+      tensfile+="littlebit";
+      MKfile+="littlebit";
+      kurtosisfile+="littlebit";
+    }
+  
+    save_volume(FA,fafile,tempinfo);
+    save_volume(S0,s0file,tempinfo);
+    save_volume(l1,l1file,tempinfo);
+    save_volume(l2,l2file,tempinfo);
+    save_volume(l3,l3file,tempinfo);
+    save_volume(MD,MDfile,tempinfo);
+    save_volume4D(V1,v1file,tempinfo);
+    save_volume4D(V2,v2file,tempinfo);
+    save_volume4D(V3,v3file,tempinfo);
+    save_volume(MK,MKfile,tempinfo);
+
+    if(opts.savetensor.value()){
+      save_volume4D(Delements,tensfile,tempinfo);
+      save_volume4D(KurtTens,kurtosisfile,tempinfo);
+    }
+      
+    
+  return 0;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+