diffusion models

diffmodels.cc

+/*  Diffusion model fitting
+
+    Timothy Behrens, Saad Jbabdi  - FMRIB Image Analysis Group
+ 
+    Copyright (C) 2005 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include "diffmodels.h"
+
+
+
+////////////////////////////////////////////////
+//       DIFFUSION TENSOR MODEL
+////////////////////////////////////////////////
+void DTI::linfit(){
+  ColumnVector logS(npts);
+  ColumnVector Dvec(7);
+  for (int i=1;i<=npts; i++){
+    if(Y(i)>0)
+      logS(i)=log(Y(i));
+    else
+      logS(i)=0;
+  }
+  Dvec = -iAmat*logS;
+  if(Dvec(7)>-23)
+    m_s0=exp(-Dvec(7));
+  else
+    m_s0=Y.MaximumAbsoluteValue();
+  for (int i=1;i<=Y.Nrows();i++){
+    if(m_s0<Y.Sum()/Y.Nrows()){ m_s0=Y.MaximumAbsoluteValue();  }
+    logS(i)=(Y(i)/m_s0)>0.01 ? log(Y(i)):log(0.01*m_s0);
+  }
+  Dvec = -iAmat*logS;
+  m_sse = (Amat*Dvec+logS).SumSquare();
+  m_s0=exp(-Dvec(7));
+  if(m_s0<Y.Sum()/Y.Nrows()){ m_s0=Y.Sum()/Y.Nrows();  }
+  vec2tens(Dvec);
+  calc_tensor_parameters();
+
+  m_covar.ReSize(7);
+  float dof=logS.Nrows()-7;
+  float sig2=m_sse/dof;
+  m_covar << sig2*(Amat.t()*Amat).i();
+}
+ColumnVector DTI::calc_md_grad(const ColumnVector& _tens)const{
+  ColumnVector g(6);
+  g = 0;
+  g(1) = 1/3.0;
+  g(4) = 1/3.0;
+  g(6) = 1/3.0;
+  return g;
+}
+// this will only work if the determinant is strictly positive
+ReturnMatrix DTI::calc_fa_grad(const ColumnVector& _intens)const{
+  ColumnVector gradv(6),ik(6),k(6);
+  float m = (_intens(1)+_intens(4)+_intens(6))/3.0;
+  SymmetricMatrix K(3),iK(3),M(3);
+
+  // rescale input matrix
+  vec2tens(_intens,M);
+  //M /=m;
+
+  m = M.Trace()/3.0;
+
+  K = M - m*IdentityMatrix(3);
+  tens2vec(K,k);
+  iK << K.i();
+  tens2vec(iK,ik);
+
+  float p   = K.SumSquare()/6.0;
+  float q   = K.Determinant()/2.0;
+  float h   = std::sqrt(p*p*p-q*q)/q;
+  float phi = std::atan(h)/3.0;
+  if(q<0)phi+=M_PI;
+
+
+  float _l1 = m + 2.0*std::sqrt(p)*std::cos(phi);
+  float _l2 = m - std::sqrt(p)*(std::cos(phi)+std::sqrt(3.0)*std::sin(phi));
+  float _l3 = m - std::sqrt(p)*(std::cos(phi)-std::sqrt(3.0)*std::sin(phi));
+
+  float t  = 6.0/9.0*(_l1*_l1+_l2*_l2+_l3*_l3 - _l1*_l2-_l1*_l3-_l2*_l3);
+  float b  = _l1*_l1+_l2*_l2+_l3*_l3;
+
+  float _fa  = std::sqrt(3.0/2.0)*std::sqrt(t/b);
+  
+
+  float dfadl1 = 3.0/4.0/_fa * ( 6.0/9.0*(2.0*_l1-_l2-_l3)/b - t/b/b*2.0*_l1 );
+  float dfadl2 = 3.0/4.0/_fa * ( 6.0/9.0*(2.0*_l2-_l1-_l3)/b - t/b/b*2.0*_l2 );
+  float dfadl3 = 3.0/4.0/_fa * ( 6.0/9.0*(2.0*_l3-_l1-_l2)/b - t/b/b*2.0*_l3 );
+
+  
+
+  // determine dkdx
+  ColumnVector dkdx(6);
+  dkdx << 2.0/3.0 << 1.0 << 1.0 << 2.0/3.0 << 1.0 << 2.0/3.0;
+
+  for(int i=1;i<=6;i++){
+    float dL1dx=0,dL2dx=0,dL3dx=0;
+    if(i==1||i==4||i==6){
+      dL1dx=1.0/3.0;dL2dx=1.0/3.0;dL3dx=1.0/3.0;
+    }
+    //
+    float p_p = k(i)/3.0 * dkdx(i);
+    float q_p = q*ik(i) * dkdx(i);
+    float h_p = (3.0*p*p*p_p - 2.0*q_p*q*(1+h*h))/2.0/h/q/q;
+
+    float phi_p = h_p/(1+h*h)/3.0;
+
+    dL1dx += p_p/std::sqrt(p)*std::cos(phi) - 2.0*std::sqrt(p)*phi_p*std::sin(phi);
+    dL2dx -= std::sqrt(p)*(.5*p_p*(m-_l2)+phi_p*(std::sin(phi)-std::sqrt(3.0)*std::cos(phi)));
+    dL3dx -= std::sqrt(p)*(.5*p_p*(m-_l3)+phi_p*(std::sin(phi)+std::sqrt(3.0)*std::cos(phi)));
+
+    //
+    gradv(i) = dfadl1*dL1dx + dfadl2*dL2dx + dfadl3*dL3dx;
+  }
+  gradv.Release();
+  return gradv;
+}
+float DTI::calc_fa_var()const{
+  ColumnVector grd;
+  ColumnVector vtens;
+  tens2vec(m_tens,vtens);
+  grd = calc_fa_grad(vtens);
+  ColumnVector g(7);
+  g.SubMatrix(1,6,1,1) = grd;
+  g(7) = 0;
+  
+  return((g.t()*m_covar*g).AsScalar());
+}
+
+void DTI::rot2angles(const Matrix& rot,float& th1,float& th2,float& th3)const{
+  if(rot(3,1)!=1 && rot(3,1)!=-1){
+    th2 = -asin(rot(3,1));
+    float c=std::cos(th2);
+    th1 = atan2(rot(3,2)/c,rot(3,3)/c);
+    th3 = atan2(rot(2,1)/c,rot(1,1)/c);
+  }
+  else{
+    th1 = atan2(rot(1,2),rot(1,3));
+    th2 = -sign(rot(3,1))*M_PI/2;
+    th3 = 0;
+  }
+}
+void DTI::angles2rot(const float& th1,const float& th2,const float& th3,Matrix& rot)const{
+  float c1=std::cos(th1),s1=std::sin(th1);
+  float c2=std::cos(th2),s2=std::sin(th2);
+  float c3=std::cos(th3),s3=std::sin(th3);
+
+  rot(1,1) = c2*c3;    rot(1,2) = s1*s2*c3-c1*s3;    rot(3,1) = c1*s2*c3+s1*s3;
+  rot(2,1) = c2*s3;    rot(2,2) = s1*s2*s3+c1*c3;    rot(3,2) = c1*s2*s3-s1*c3;
+  rot(3,1) = -s2;      rot(3,2) = s1*c2;             rot(3,3) = c1*c2;
+}
+
+
+// nonlinear tensor fitting 
+void DTI::nonlinfit(){
+  // initialise with linear fitting
+  linfit();
+
+  print();
+
+  // set starting parameters
+  // params = s0, log(l1),log(l2), log(l3), th1, th2, th3
+  ColumnVector start(nparams);
+
+  start(1) = m_s0;
+  // eigenvalues
+  start(2) = m_l1>0?std::log(m_l1):std::log(1e-5);
+  start(3) = m_l2>0?std::log(m_l2):std::log(1e-5);
+  start(4) = m_l3>0?std::log(m_l3):std::log(1e-5);
+  // angles
+  float th1,th2,th3;
+  Matrix rot(3,3);
+  rot.Row(1) = m_v1.t();
+  rot.Row(2) = m_v2.t();
+  rot.Row(3) = m_v3.t();
+  rot2angles(rot,th1,th2,th3);
+  start(5) = th1;
+  start(6) = th2;
+  start(7) = th3;
+
+
+  // do the fit
+  NonlinParam  lmpar(start.Nrows(),NL_LM); 
+  lmpar.SetGaussNewtonType(LM_L);
+  lmpar.SetStartingEstimate(start);
+
+
+  NonlinOut status;
+  status = nonlin(lmpar,(*this));
+  ColumnVector final_par(nparams);
+  final_par = lmpar.Par();
+
+
+  // finalise parameters
+  m_s0 = final_par(1);
+  m_l1 = exp(final_par(2));
+  m_l2 = exp(final_par(3));
+  m_l3 = exp(final_par(4));
+
+  angles2rot(final_par(5),final_par(6),final_par(7),rot);
+  m_v1 = rot.Row(1).t();
+  m_v2 = rot.Row(2).t();
+  m_v3 = rot.Row(3).t();
+
+  sort();
+
+  m_tens << m_l1*m_v1*m_v1.t() + m_l2*m_v2*m_v2.t() + m_l3*m_v3*m_v3.t();
+  calc_tensor_parameters();
+
+  print();
+  //exit(1);
+
+}
+void DTI::sort(){
+  vector< pair<float,int> > ls(3);
+  vector<ColumnVector> vs(3);
+  ls[0].first=m_l1;
+  ls[0].second=0;
+  ls[1].first=m_l2;
+  ls[1].second=1;
+  ls[2].first=m_l3;
+  ls[2].second=2;
+  vs[0]=m_v1;vs[1]=m_v2;vs[2]=m_v3;
+  
+  std::sort(ls.begin(),ls.end());
+
+  m_l1 = ls[2].first;
+  m_v1 = vs[ ls[2].second ];
+  m_l2 = ls[1].first;
+  m_v2 = vs[ ls[1].second ];
+  m_l3 = ls[0].first;
+  m_v3 = vs[ ls[0].second ];
+  
+}
+void DTI::calc_tensor_parameters(){
+  Matrix Vd;
+  DiagonalMatrix Dd(3);
+  // mean, eigenvalues and eigenvectors
+  EigenValues(m_tens,Dd,Vd);
+  m_md = Dd.Sum()/Dd.Nrows();
+  m_l1 = Dd(3,3);
+  m_l2 = Dd(2,2);
+  m_l3 = Dd(1,1);
+  m_v1 = Vd.Column(3);
+  m_v2 = Vd.Column(2);
+  m_v3 = Vd.Column(1);
+  // mode
+  float e1=m_l1-m_md, e2=m_l2-m_md, e3=m_l3-m_md;
+  float n = (e1 + e2 - 2*e3)*(2*e1 - e2 - e3)*(e1 - 2*e2 + e3);
+  float d = (e1*e1 + e2*e2 + e3*e3 - e1*e2 - e2*e3 - e1*e3);
+  d = sqrt(bigger(0, d));
+  d = 2*d*d*d;
+  m_mo = smaller(bigger(d ? n/d : 0.0, -1),1);
+  // fa
+  float numer=1.5*((m_l1-m_md)*(m_l1-m_md)+(m_l2-m_md)*(m_l2-m_md)+(m_l3-m_md)*(m_l3-m_md));
+  float denom=(m_l1*m_l1+m_l2*m_l2+m_l3*m_l3);
+  if(denom>0) m_fa=numer/denom;
+  else m_fa=0;
+  if(m_fa>0){m_fa=sqrt(m_fa);}
+  else{m_fa=0;}
+}
+// now the nonlinear fitting routines
+double DTI::cf(const NEWMAT::ColumnVector& p)const{
+  //cout << "CF" << endl;
+  //OUT(p.t());
+  double cfv = 0.0;
+  double err = 0.0;
+  ////////////////////////////////////
+  ColumnVector ls(3);
+  Matrix rot(3,3);
+  angles2rot(p(5),p(6),p(7),rot);
+  for(int k=2;k<=4;k++){
+    ls(k-1) = exp(p(k));
+  }
+  ////////////////////////////////////
+  for(int i=1;i<=Y.Nrows();i++){
+    err = p(1)*anisoterm(i,ls,rot) - Y(i); 
+    cfv += err*err; 
+  }  
+  //OUT(cfv);
+  //cout<<"--------"<<endl;
+  return(cfv);
+}
+
+NEWMAT::ReturnMatrix DTI::grad(const NEWMAT::ColumnVector& p)const{
+  NEWMAT::ColumnVector gradv(p.Nrows());
+
+  cout<<"grad"<<endl;
+  OUT(p.t());
+
+  gradv = 0.0;
+  ////////////////////////////////////
+  ColumnVector ls(3);
+  Matrix rot(3,3);
+  Matrix rot1(3,3),rot2(3,3),rot3(3,3);
+  angles2rot(p(5),p(6),p(7),rot);
+
+  angles2rot(p(5)+M_PI/2.0,p(6),p(7),rot1);
+  angles2rot(p(5),p(6)+M_PI/2.0,p(7),rot2);
+  angles2rot(p(5),p(6),p(7)+M_PI/2.0,rot3);
+  for(int k=2;k<=4;k++){
+    ls(k-1) = exp(p(k));
+  }
+  ////////////////////////////////////
+  Matrix J(npts,nparams);
+  ColumnVector x(3);
+  ColumnVector diff(npts);
+  float sig;
+  for(int i=1;i<=Y.Nrows();i++){
+    sig = p(1)*anisoterm(i,ls,rot);
+    
+    J(i,1) = sig/p(1);
+
+    x = rotproduct(bvecs.Column(i),rot);
+    J(i,2) = -bvals(1,i)*x(1)*sig*ls(1);
+    J(i,3) = -bvals(1,i)*x(2)*sig*ls(2);
+    J(i,4) = -bvals(1,i)*x(3)*sig*ls(3);
+
+    x = rotproduct(bvecs.Column(i),rot1,rot);
+    J(i,5) = -2.0*bvals(1,i)*(ls(1)*x(1)+ls(2)*x(2)+ls(3)*x(3))*sig;
+    x = rotproduct(bvecs.Column(i),rot2,rot);
+    J(i,6) = -2.0*bvals(1,i)*(ls(1)*x(1)+ls(2)*x(2)+ls(3)*x(3))*sig;
+    x = rotproduct(bvecs.Column(i),rot3,rot);
+    J(i,7) = -2.0*bvals(1,i)*(ls(1)*x(1)+ls(2)*x(2)+ls(3)*x(3))*sig;
+
+    diff(i) = sig - Y(i);
+  }
+
+  OUT(diff.t());
+  OUT(J.t());
+  
+  gradv = 2.0*J.t()*diff;
+
+  OUT(gradv.t());
+  cout<<"------"<<endl;
+
+  gradv.Release();
+  return gradv;
+}
+
+//this uses Gauss-Newton approximation
+boost::shared_ptr<BFMatrix> DTI::hess(const NEWMAT::ColumnVector& p,boost::shared_ptr<BFMatrix> iptr)const{
+  boost::shared_ptr<BFMatrix>   hessm;
+  if (iptr && iptr->Nrows()==(unsigned int)p.Nrows() && iptr->Ncols()==(unsigned int)p.Nrows()) hessm = iptr;
+  else hessm = boost::shared_ptr<BFMatrix>(new FullBFMatrix(p.Nrows(),p.Nrows()));
+
+  cout<<"hess"<<endl;
+  OUT(p.t());
+  
+  ////////////////////////////////////
+  ColumnVector ls(3);
+  Matrix rot(3,3);
+  Matrix rot1(3,3),rot2(3,3),rot3(3,3);
+  angles2rot(p(5),p(6),p(7),rot);
+
+  angles2rot(p(5)+M_PI/2,p(6),p(7),rot1);
+  angles2rot(p(5),p(6)+M_PI/2,p(7),rot2);
+  angles2rot(p(5),p(6),p(7)+M_PI/2,rot3);
+  for(int k=2;k<=4;k++){
+    ls(k-1) = exp(p(k));
+  }
+  ////////////////////////////////////
+  Matrix J(npts,nparams);
+  ColumnVector x(3);
+  float sig;
+  for(int i=1;i<=Y.Nrows();i++){
+    sig = p(1)*anisoterm(i,ls,rot);
+    
+    J(i,1) = sig/p(1);
+
+    x = rotproduct(bvecs.Column(i),rot);
+    J(i,2) = -bvals(1,i)*x(1)*sig*ls(1);
+    J(i,3) = -bvals(1,i)*x(2)*sig*ls(2);
+    J(i,4) = -bvals(1,i)*x(3)*sig*ls(3);
+
+    x = rotproduct(bvecs.Column(i),rot1,rot);
+    J(i,5) = -2.0*bvals(1,i)*(ls(1)*x(1)+ls(2)*x(2)+ls(3)*x(3))*sig;
+    x = rotproduct(bvecs.Column(i),rot2,rot);
+    J(i,6) = -2.0*bvals(1,i)*(ls(1)*x(1)+ls(2)*x(2)+ls(3)*x(3))*sig;
+    x = rotproduct(bvecs.Column(i),rot3,rot);
+    J(i,7) = -2.0*bvals(1,i)*(ls(1)*x(1)+ls(2)*x(2)+ls(3)*x(3))*sig;
+  }
+  
+
+  for (int i=1; i<=p.Nrows(); i++){
+    for (int j=i; j<=p.Nrows(); j++){
+      sig = 0.0;
+      for(int k=1;k<=J.Nrows();k++)
+	sig += 2.0*(J(k,i)*J(k,j));
+      hessm->Set(i,j,sig);
+    }
+  }
+  for (int j=1; j<=p.Nrows(); j++) {
+    for (int i=j+1; i<=p.Nrows(); i++) {
+      hessm->Set(i,j,hessm->Peek(j,i));
+    }
+  }
+
+  hessm->Print();
+  cout<<"------"<<endl;
+
+  return(hessm);
+}
+
+ColumnVector DTI::rotproduct(const ColumnVector& x,const Matrix& R)const{
+  ColumnVector ret(3);
+  
+  for(int i=1;i<=3;i++)
+    ret(i) = x(1)*x(1)*R(1,i)*R(1,i)+x(2)*x(2)*R(2,i)*R(2,i)+x(3)*x(3)*R(3,i)*R(3,i)
+      +2.0*( x(1)*R(1,i)*(x(2)*R(2,i)+x(3)*R(3,i)) +x(2)*x(3)*R(2,i)*R(3,i) );   
+  
+  return ret;
+}
+ColumnVector DTI::rotproduct(const ColumnVector& x,const Matrix& R1,const Matrix& R2)const{
+  ColumnVector ret(3);
+  
+  for(int i=1;i<=3;i++)
+    ret(i) = x(1)*x(1)*R1(1,i)*R2(1,i)+x(2)*x(2)*R1(2,i)*R2(2,i)+x(3)*x(3)*R1(3,i)*R2(3,i)
+      +( x(1)*R1(1,i)*(x(2)*R2(2,i)+x(3)*R2(3,i)) +x(2)*x(3)*R1(2,i)*R2(3,i) )
+      +( x(1)*R2(1,i)*(x(2)*R1(2,i)+x(3)*R1(3,i)) +x(2)*x(3)*R2(2,i)*R1(3,i) );   
+  
+  return ret;
+}
+
+float DTI::anisoterm(const int& pt,const ColumnVector& ls,const Matrix& rot)const{
+  ColumnVector x(3);
+  x = rotproduct(bvecs.Column(pt),rot);
+
+  return exp(-bvals(1,pt)*(ls(1)*x(1)+ls(2)*x(2)+ls(3)*x(3)));
+}
+
+
+////////////////////////////////////////////////
+//       PARTIAL VOLUME MODEL - SINGLE SHELL
+////////////////////////////////////////////////
+
+void PVM_single::fit(){
+
+  // initialise with a tensor
+  DTI dti(Y,bvecs,bvals);
+  dti.linfit();
+
+  // set starting parameters for nonlinear fitting
+  float _th,_ph;
+  cart2sph(dti.get_v1(),_th,_ph);
+
+  ColumnVector start(nparams);
+  start(1) = dti.get_s0();
+  start(2) = dti.get_md()>0?dti.get_md()*2:0.001; // empirically found that d~2*MD
+  start(3) = f2x(dti.get_fa()); // first pvf = FA 
+  start(4) = _th;
+  start(5) = _ph;
+  float sumf=x2f(start(2));
+  float tmpsumf=sumf;
+  for(int ii=2,i=6;ii<=nfib;ii++,i+=3){
+    float denom=2;
+    do{
+      start(i) = f2x(x2f(start(i-3))/denom);
+      denom *= 2;
+      tmpsumf = sumf + x2f(start(i));
+    }while(tmpsumf>=1);
+    sumf += x2f(start(i));
+    cart2sph(dti.get_v(ii),_th,_ph);
+    start(i+1) = _th;
+    start(i+2) = _ph;
+  }
+
+
+  // do the fit
+  NonlinParam  lmpar(start.Nrows(),NL_LM); 
+  lmpar.SetGaussNewtonType(LM_L);
+  lmpar.SetStartingEstimate(start);
+
+  NonlinOut status;
+  status = nonlin(lmpar,(*this));
+  ColumnVector final_par(nparams);
+  final_par = lmpar.Par();
+
+
+  // finalise parameters
+  m_s0 = final_par(1);
+  m_d  = final_par(2);
+  for(int k=1;k<=nfib;k++){
+    int kk = 3 + 3*(k-1);
+    m_f(k)  = x2f(final_par(kk));
+    m_th(k) = final_par(kk+1);
+    m_ph(k) = final_par(kk+2);
+  }
+  sort();
+  fix_fsum();
+}
+void PVM_single::sort(){
+  vector< pair<float,int> > fvals(nfib);
+  ColumnVector ftmp(nfib),thtmp(nfib),phtmp(nfib);
+  ftmp=m_f;thtmp=m_th;phtmp=m_ph;
+  for(int i=1;i<=nfib;i++){
+    pair<float,int> p(m_f(i),i);
+    fvals[i-1] = p;
+  }
+  std::sort(fvals.begin(),fvals.end());
+  for(int i=1,ii=nfib-1;ii>=0;i++,ii--){
+    m_f(i)  = ftmp(fvals[ii].second);
+    m_th(i) = thtmp(fvals[ii].second);
+    m_ph(i) = phtmp(fvals[ii].second);
+  }
+}
+void PVM_single::fix_fsum(){
+  float sumf=0;
+  for(int i=1;i<=nfib;i++){
+    sumf+=m_f(i);
+    if(sumf>=1){for(int j=i;j<=nfib;j++)m_f(i)=0;break;}
+  }
+}
+ReturnMatrix PVM_single::get_prediction()const{
+  ColumnVector pred(npts);
+  ColumnVector p(nparams);
+  p(1) = m_s0;
+  p(2) = m_d;
+  for(int i=3,ii=1;ii<=nfib;i+=3,ii++){
+    p(i)   = f2x(m_f(ii));
+    p(i+1) = m_th(ii);
+    p(i+2) = m_ph(ii);
+  }
+  pred = forwardModel(p);
+
+  pred.Release();
+  return pred;
+}
+NEWMAT::ReturnMatrix PVM_single::forwardModel(const NEWMAT::ColumnVector& p)const{
+  //cout<<"FORWARD"<<endl;
+  //OUT(p.t());
+  ColumnVector pred(npts);
+  pred = 0;
+  float val;
+  float _d = p(2);
+  ////////////////////////////////////
+  ColumnVector fs(nfib);
+  Matrix x(nfib,3);
+  float sumf=0;
+  for(int k=1;k<=nfib;k++){
+    int kk = 3+3*(k-1);
+    fs(k) = x2f(p(kk));
+    sumf += fs(k);
+    x(k,1) = sin(p(kk+1))*cos(p(kk+2));
+    x(k,2) = sin(p(kk+1))*sin(p(kk+2));
+    x(k,3) = cos(p(kk+1));
+  }
+  ////////////////////////////////////
+  for(int i=1;i<=Y.Nrows();i++){
+    val = 0.0;
+    for(int k=1;k<=nfib;k++){
+      val += fs(k)*anisoterm(i,_d,x.Row(k).t());
+    }
+    pred(i) = p(1)*((1-sumf)*isoterm(i,_d)+val); 
+  }  
+  pred.Release();
+  //cout<<"----"<<endl;
+  return pred;
+}
+double PVM_single::cf(const NEWMAT::ColumnVector& p)const{
+  //cout<<"CF"<<endl;
+  //OUT(p.t());
+  double cfv = 0.0;
+  double err;
+  float _d = p(2);
+  ////////////////////////////////////
+  ColumnVector fs(nfib);
+  Matrix x(nfib,3);
+  float sumf=0;
+  for(int k=1;k<=nfib;k++){
+    int kk = 3+3*(k-1);
+    fs(k) = x2f(p(kk));
+    sumf += fs(k);
+    x(k,1) = sin(p(kk+1))*cos(p(kk+2));
+    x(k,2) = sin(p(kk+1))*sin(p(kk+2));
+    x(k,3) = cos(p(kk+1));
+  }
+  ////////////////////////////////////
+  for(int i=1;i<=Y.Nrows();i++){
+    err = 0.0;
+    for(int k=1;k<=nfib;k++){
+      err += fs(k)*anisoterm(i,_d,x.Row(k).t());
+    }
+    err = (p(1)*((1-sumf)*isoterm(i,_d)+err) - Y(i)); 
+    cfv += err*err; 
+  }  
+  //OUT(cfv);
+  //cout<<"----"<<endl;
+  return(cfv);
+}
+
+NEWMAT::ReturnMatrix PVM_single::grad(const NEWMAT::ColumnVector& p)const{
+  //cout<<"GRAD"<<endl;
+  //OUT(p.t());
+  NEWMAT::ColumnVector gradv(p.Nrows());
+  gradv = 0.0;
+  float _d = p(2);
+  ////////////////////////////////////
+  ColumnVector fs(nfib);
+  Matrix x(nfib,3);ColumnVector xx(3);
+  float sumf=0;
+  for(int k=1;k<=nfib;k++){
+    int kk = 3+3*(k-1);
+    fs(k) = x2f(p(kk));
+    sumf += fs(k);
+    x(k,1) = sin(p(kk+1))*cos(p(kk+2));
+    x(k,2) = sin(p(kk+1))*sin(p(kk+2));
+    x(k,3) = cos(p(kk+1));
+  }
+  ////////////////////////////////////
+  Matrix J(npts,nparams);
+  ColumnVector diff(npts);
+  float sig;
+  for(int i=1;i<=Y.Nrows();i++){
+    sig = 0;
+    J.Row(i)=0;
+    for(int k=1;k<=nfib;k++){
+      int kk = 3+3*(k-1);
+      xx = x.Row(k).t();
+      sig += fs(k)*anisoterm(i,_d,xx);
+      // other stuff for derivatives
+      // d
+      J(i,2) += p(1)*fs(k)*anisoterm_d(i,_d,xx);
+      // f
+      J(i,kk) = p(1)*(anisoterm(i,_d,xx)-isoterm(i,_d)) * two_pi*sign(p(kk))*1/(1+p(kk)*p(kk));
+      // th
+      J(i,kk+1) = p(1)*fs(k)*anisoterm_th(i,_d,xx,p(kk+1),p(kk+2));
+      // ph
+      J(i,kk+2) = p(1)*fs(k)*anisoterm_ph(i,_d,xx,p(kk+1),p(kk+2));
+    }
+    sig = p(1)*((1-sumf)*isoterm(i,_d)+sig);
+    diff(i) = sig - Y(i);
+    // other stuff for derivatives
+    J(i,1) = sig/p(1);
+    J(i,2) += p(1)*(1-sumf)*isoterm_d(i,_d);
+  }
+  
+  gradv = 2*J.t()*diff;
+  //OUT(gradv.t());
+  //cout<<"----"<<endl;
+  gradv.Release();
+  return gradv;
+
+
+}
+//this uses Gauss-Newton approximation
+boost::shared_ptr<BFMatrix> PVM_single::hess(const NEWMAT::ColumnVector& p,boost::shared_ptr<BFMatrix> iptr)const{
+  //cout<<"HESS"<<endl;
+  //OUT(p.t());
+  boost::shared_ptr<BFMatrix>   hessm;
+  if (iptr && iptr->Nrows()==(unsigned int)p.Nrows() && iptr->Ncols()==(unsigned int)p.Nrows()) hessm = iptr;
+  else hessm = boost::shared_ptr<BFMatrix>(new FullBFMatrix(p.Nrows(),p.Nrows()));
+
+  float _d = p(2);
+  ////////////////////////////////////
+  ColumnVector fs(nfib);
+  Matrix x(nfib,3);ColumnVector xx(3);
+  float sumf=0;
+  for(int k=1;k<=nfib;k++){
+    int kk = 3+3*(k-1);
+    fs(k) = x2f(p(kk));
+    sumf += fs(k);
+    x(k,1) = sin(p(kk+1))*cos(p(kk+2));
+    x(k,2) = sin(p(kk+1))*sin(p(kk+2));
+    x(k,3) = cos(p(kk+1));
+  }
+  ////////////////////////////////////
+  Matrix J(npts,nparams);
+  vector<Matrix> J2(npts);
+  ColumnVector diff(npts);
+  float sig;
+  for(int i=1;i<=Y.Nrows();i++){
+    sig = 0;
+    J.Row(i)=0;
+    J2[i-1].ReSize(nparams,nparams);
+    J2[i-1] = 0;
+    for(int k=1;k<=nfib;k++){
+      int kk = 3+3*(k-1);
+      xx = x.Row(k).t();
+      sig += fs(k)*anisoterm(i,_d,xx);
+      // other stuff for derivatives
+      // change of variable
+      float cov = two_pi*sign(p(kk))*1/(1+p(kk)*p(kk));
+      // d
+      J(i,2) += p(1)*fs(k)*anisoterm_d(i,_d,xx);
+      // f
+      J(i,kk) = p(1)*(anisoterm(i,_d,xx)-isoterm(i,_d)) * cov;
+      // th
+      J(i,kk+1) = p(1)*fs(k)*anisoterm_th(i,_d,xx,p(kk+1),p(kk+2));
+      // ph
+      J(i,kk+2) = p(1)*fs(k)*anisoterm_ph(i,_d,xx,p(kk+1),p(kk+2));
+      // 2nd order derivatives
+      J2[i-1](1,kk)   = J(i,kk)/p(1) *cov;                  // s0,f
+      J2[i-1](1,kk+1) = J(i,kk+1)/p(1);                // s0,th
+      J2[i-1](1,kk+2) = J(i,kk+2)/p(1);                // s0,ph
+      J2[i-1](2,2)    += p(1)*fs(k)*anisoterm_dd(i,_d,xx); // d,d
+      J2[i-1](2,kk)   = p(1)*(anisoterm_d(i,_d,xx) - isoterm_d(i,_d)) *cov; // d,f
+      J2[i-1](2,kk+1) = p(1)*fs(k)*(anisoterm_dth(i,_d,xx,p(kk+1),p(k+2))); // d,th
+      J2[i-1](2,kk+2) = p(1)*fs(k)*(anisoterm_dph(i,_d,xx,p(kk+1),p(k+2))); // d,ph
+
+      J2[i-1](kk,kk)   = J(i,kk)*(-2*p(kk)*cov); // f,f 
+      J2[i-1](kk,kk+1) = J(i,kk+1)/fs(k) *cov;  // f,th
+      J2[i-1](kk,kk+2) = J(i,kk+2)/fs(k) *cov;  // f,ph
+      J2[i-1](kk+1,kk+1) = p(1)*fs(k)*anisoterm_thth(i,_d,xx,p(kk+1),p(k+2)); // th,th
+      J2[i-1](kk+1,kk+2) = p(1)*fs(k)*anisoterm_thph(i,_d,xx,p(kk+1),p(k+2)); // th,ph
+      J2[i-1](kk+2,kk+2) = p(1)*fs(k)*anisoterm_phph(i,_d,xx,p(kk+1),p(k+2)); // ph,ph
+
+    }
+    sig = p(1)*((1-sumf)*isoterm(i,_d)+sig);
+    diff(i) = sig - Y(i);
+    // other stuff for derivatives
+    J(i,1) = sig/p(1);
+    J(i,2) += p(1)*(1-sumf)*isoterm_d(i,_d);
+
+    // 2nd order derivatives
+    J2[i-1](1,1) = 0;                                      // s0,s0
+    J2[i-1](1,2) = J(i,2)/p(1);                            // s0,d
+    J2[i-1](2,2) += p(1)*(1-sumf)*isoterm_dd(i,_d);        // d,d
+    
+  }
+  
+
+  for (int i=1; i<=p.Nrows(); i++){
+    for (int j=i; j<=p.Nrows(); j++){
+      sig = 0.0;
+      for(int k=1;k<=J.Nrows();k++)
+	sig += 2*(J(k,i)*J(k,j));// + J2[k-1](i,j)*diff(k));
+      hessm->Set(i,j,sig);
+    }
+  }
+  for (int j=1; j<=p.Nrows(); j++) {
+    for (int i=j+1; i<=p.Nrows(); i++) {
+      hessm->Set(i,j,hessm->Peek(j,i));
+    }
+  }
+  //hessm->Print();
+  //cout<<"----"<<endl;
+  return(hessm);
+}
+
+
+
+
+////////////////////////////////////////////////
+//       PARTIAL VOLUME MODEL - MULTIPLE SHELLS
+////////////////////////////////////////////////
+
+void PVM_multi::fit(){
+
+  // initialise with a tensor
+  DTI dti(Y,bvecs,bvals);
+  dti.linfit();
+
+  // set starting parameters for nonlinear fitting
+  float _th,_ph;
+  cart2sph(dti.get_v1(),_th,_ph);
+
+  float _a,_b;
+  _a = 1; // start with d=d_std
+  _b = dti.get_md()>0?dti.get_md()*2:0.001;
+
+  ColumnVector start(nparams);
+  start(1) = dti.get_s0();
+  start(2) = _a;
+  start(3) = _b;
+  start(4) = f2x(dti.get_fa());                   // first pvf = FA 
+  start(5) = _th;
+  start(6) = _ph;
+  float sumf=x2f(start(2));
+  float tmpsumf=sumf;
+  for(int ii=2,i=7;ii<=nfib;ii++,i+=3){
+    float denom=2;
+    do{
+      start(i) = f2x(x2f(start(i-3))/denom);
+      denom *= 2;
+      tmpsumf = sumf + x2f(start(i));
+    }while(tmpsumf>=1);
+    sumf += x2f(start(i));
+    cart2sph(dti.get_v(ii),_th,_ph);
+    start(i+1) = _th;
+    start(i+2) = _ph;
+  }
+  
+  // do the fit
+  NonlinParam  lmpar(start.Nrows(),NL_LM); 
+  lmpar.SetGaussNewtonType(LM_L);
+  lmpar.SetStartingEstimate(start);
+
+  NonlinOut status;
+  status = nonlin(lmpar,(*this));
+  ColumnVector final_par(nparams);
+  final_par = lmpar.Par();
+
+
+  // finalise parameters
+  m_s0     = final_par(1);
+  m_d      = final_par(2)*final_par(3);
+  m_d_std  = std::sqrt(final_par(2))*final_par(3);
+  for(int i=4,k=1;k<=nfib;i+=3,k++){
+    m_f(k)  = x2f(final_par(i));
+    m_th(k) = final_par(i+1);
+    m_ph(k) = final_par(i+2);
+  }
+  sort();
+  fix_fsum();
+}
+void PVM_multi::sort(){
+  vector< pair<float,int> > fvals(nfib);
+  ColumnVector ftmp(nfib),thtmp(nfib),phtmp(nfib);
+  ftmp=m_f;thtmp=m_th;phtmp=m_ph;
+  for(int i=1;i<=nfib;i++){
+    pair<float,int> p(m_f(i),i);
+    fvals[i-1] = p;
+  }
+  std::sort(fvals.begin(),fvals.end());
+  for(int i=1,ii=nfib-1;ii>=0;i++,ii--){
+    m_f(i)  = ftmp(fvals[ii].second);
+    m_th(i) = thtmp(fvals[ii].second);
+    m_ph(i) = phtmp(fvals[ii].second);
+  }
+}
+void PVM_multi::fix_fsum(){
+  float sumf=0;
+  for(int i=1;i<=nfib;i++){
+    sumf+=m_f(i);
+    if(sumf>=1){for(int j=i;j<=nfib;j++)m_f(i)=0;break;}
+  }
+}
+ReturnMatrix PVM_multi::get_prediction()const{
+  ColumnVector pred(npts);
+  ColumnVector p(nparams);
+  p(1) = m_s0;
+  p(2) = m_d*m_d/m_d_std/m_d_std;
+  p(3) = m_d_std*m_d_std/m_d;
+  for(int k=1;k<=nfib;k++){
+    int kk = 4+3*(k-1);
+    p(kk)   = f2x(m_f(k));
+    p(kk+1) = m_th(k);
+    p(kk+2) = m_ph(k);
+  }
+  pred = forwardModel(p);
+  pred.Release();
+  return pred;
+}
+NEWMAT::ReturnMatrix PVM_multi::forwardModel(const NEWMAT::ColumnVector& p)const{
+  ColumnVector pred(npts);
+  pred = 0;
+  float val;
+  float _a = p(2);
+  float _b = p(3);
+  ////////////////////////////////////
+  ColumnVector fs(nfib);
+  Matrix x(nfib,3);
+  float sumf=0;
+  for(int k=1;k<=nfib;k++){
+    int kk = 4+3*(k-1);
+    fs(k) = x2f(p(kk));
+    sumf += fs(k);
+    x(k,1) = sin(p(kk+1))*cos(p(kk+2));
+    x(k,2) = sin(p(kk+1))*sin(p(kk+2));
+    x(k,3) = cos(p(kk+1));
+  }
+  ////////////////////////////////////
+  for(int i=1;i<=Y.Nrows();i++){
+    val = 0.0;
+    for(int k=1;k<=nfib;k++){
+      val += fs(k)*anisoterm(i,_a,_b,x.Row(k).t());
+    }
+    pred(i) = p(1)*((1-sumf)*isoterm(i,_a,_b)+val); 
+  }  
+  pred.Release();
+  return pred;
+}
+double PVM_multi::cf(const NEWMAT::ColumnVector& p)const{
+  //cout<<"CF"<<endl;
+  //OUT(p.t());
+  double cfv = 0.0;
+  double err;
+  float _a = p(2);
+  float _b = p(3);
+  ////////////////////////////////////
+  ColumnVector fs(nfib);
+  Matrix x(nfib,3);
+  float sumf=0;
+  for(int k=1;k<=nfib;k++){
+    int kk = 4+3*(k-1);
+    fs(k) = x2f(p(kk));
+    sumf += fs(k);
+    x(k,1) = sin(p(kk+1))*cos(p(kk+2));
+    x(k,2) = sin(p(kk+1))*sin(p(kk+2));
+    x(k,3) = cos(p(kk+1));
+  }
+  ////////////////////////////////////
+  for(int i=1;i<=Y.Nrows();i++){
+    err = 0.0;
+    for(int k=1;k<=nfib;k++){
+      err += fs(k)*anisoterm(i,_a,_b,x.Row(k).t());
+    }
+    err = (p(1)*((1-sumf)*isoterm(i,_a,_b)+err) - Y(i)); 
+    cfv += err*err; 
+  }  
+  //OUT(cfv);
+  //cout<<"----"<<endl;
+  return(cfv);
+}
+
+NEWMAT::ReturnMatrix PVM_multi::grad(const NEWMAT::ColumnVector& p)const{
+  //cout<<"GRAD"<<endl;
+  //OUT(p.t());
+  NEWMAT::ColumnVector gradv(p.Nrows());
+  gradv = 0.0;
+  float _a = p(2);
+  float _b = p(3);
+  ////////////////////////////////////
+  ColumnVector fs(nfib);
+  Matrix x(nfib,3);ColumnVector xx(3);
+  float sumf=0;
+  for(int k=1;k<=nfib;k++){
+    int kk = 4+3*(k-1);
+    fs(k) = x2f(p(kk));
+    sumf += fs(k);
+    x(k,1) = sin(p(kk+1))*cos(p(kk+2));
+    x(k,2) = sin(p(kk+1))*sin(p(kk+2));
+    x(k,3) = cos(p(kk+1));
+  }
+  ////////////////////////////////////
+  Matrix J(npts,nparams);
+  ColumnVector diff(npts);
+  float sig;
+  for(int i=1;i<=Y.Nrows();i++){
+    sig = 0;
+    J.Row(i)=0;
+    for(int k=1;k<=nfib;k++){
+      int kk = 4+3*(k-1);
+      xx = x.Row(k).t();
+      sig += fs(k)*anisoterm(i,_a,_b,xx);
+      // other stuff for derivatives
+      // alpha
+      J(i,2) += p(1)*fs(k)*anisoterm_a(i,_a,_b,xx);
+      // beta
+      J(i,2) += p(1)*fs(k)*anisoterm_b(i,_a,_b,xx) * (-p(3)*p(3)); // change of variable beta=1/beta
+      // f
+      J(i,kk) = p(1)*(anisoterm(i,_a,_b,xx)-isoterm(i,_a,_b)) * two_pi*sign(p(kk))*1/(1+p(kk)*p(kk));
+      // th
+      J(i,kk+1) = p(1)*fs(k)*anisoterm_th(i,_a,_b,xx,p(kk+1),p(kk+2));
+      // ph
+      J(i,kk+2) = p(1)*fs(k)*anisoterm_ph(i,_a,_b,xx,p(kk+1),p(kk+2));
+    }
+    sig = p(1)*((1-sumf)*isoterm(i,_a,_b)+sig);
+    diff(i) = sig - Y(i);
+    // other stuff for derivatives
+    J(i,1) = sig/p(1);
+    J(i,2) += p(1)*(1-sumf)*isoterm_a(i,_a,_b);
+    J(i,3) += p(1)*(1-sumf)*isoterm_b(i,_a,_b);
+  }
+  
+  gradv = 2*J.t()*diff;
+  //OUT(gradv.t());
+  //cout<<"----"<<endl;
+  gradv.Release();
+  return gradv;
+
+
+}
+//this uses Gauss-Newton approximation
+boost::shared_ptr<BFMatrix> PVM_multi::hess(const NEWMAT::ColumnVector& p,boost::shared_ptr<BFMatrix> iptr)const{
+  //cout<<"HESS"<<endl;
+  //OUT(p.t());
+  boost::shared_ptr<BFMatrix>   hessm;
+  if (iptr && iptr->Nrows()==(unsigned int)p.Nrows() && iptr->Ncols()==(unsigned int)p.Nrows()) hessm = iptr;
+  else hessm = boost::shared_ptr<BFMatrix>(new FullBFMatrix(p.Nrows(),p.Nrows()));
+
+  float _a = p(2);
+  float _b = p(3);
+  ////////////////////////////////////
+  ColumnVector fs(nfib);
+  Matrix x(nfib,3);ColumnVector xx(3);
+  float sumf=0;
+  for(int k=1;k<=nfib;k++){
+    int kk = 4+3*(k-1);
+    fs(k) = x2f(p(kk));
+    sumf += fs(k);
+    x(k,1) = sin(p(kk+1))*cos(p(kk+2));
+    x(k,2) = sin(p(kk+1))*sin(p(kk+2));
+    x(k,3) = cos(p(kk+1));
+  }
+  ////////////////////////////////////
+  Matrix J(npts,nparams);
+  ColumnVector diff(npts);
+  float sig;
+  for(int i=1;i<=Y.Nrows();i++){
+    sig = 0;
+    J.Row(i)=0;
+    for(int k=1;k<=nfib;k++){
+      int kk = 4+3*(k-1);
+      xx = x.Row(k).t();
+      sig += fs(k)*anisoterm(i,_a,_b,xx);
+      // other stuff for derivatives
+      // change of variable
+      float cov = two_pi*sign(p(kk))*1/(1+p(kk)*p(kk));
+      // alpha
+      J(i,2) += p(1)*fs(k)*anisoterm_a(i,_a,_b,xx);
+      // beta
+      J(i,3) += p(1)*fs(k)*anisoterm_b(i,_a,_b,xx);
+      // f
+      J(i,kk) = p(1)*(anisoterm(i,_a,_b,xx)-isoterm(i,_a,_b)) * cov;
+      // th
+      J(i,kk+1) = p(1)*fs(k)*anisoterm_th(i,_a,_b,xx,p(kk+1),p(kk+2));
+      // ph
+      J(i,kk+2) = p(1)*fs(k)*anisoterm_ph(i,_a,_b,xx,p(kk+1),p(kk+2));
+    }
+    sig = p(1)*((1-sumf)*isoterm(i,_a,_b)+sig);
+    diff(i) = sig - Y(i);
+    // other stuff for derivatives
+    J(i,1) = sig/p(1);
+    J(i,2) += p(1)*(1-sumf)*isoterm_a(i,_a,_b);
+    J(i,3) += p(1)*(1-sumf)*isoterm_b(i,_a,_b);
+
+  }
+  
+
+  for (int i=1; i<=p.Nrows(); i++){
+    for (int j=i; j<=p.Nrows(); j++){
+      sig = 0.0;
+      for(int k=1;k<=J.Nrows();k++)
+	sig += 2*(J(k,i)*J(k,j));
+      hessm->Set(i,j,sig);
+    }
+  }
+  for (int j=1; j<=p.Nrows(); j++) {
+    for (int i=j+1; i<=p.Nrows(); i++) {
+      hessm->Set(i,j,hessm->Peek(j,i));
+    }
+  }
+  //hessm->Print();
+  //cout<<"----"<<endl;
+  return(hessm);
+}
+
+
+
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////
+//               USEFUL FUNCTIONS TO CALCULATE DERIVATIVES
+///////////////////////////////////////////////////////////////////////////////////////////////
+// functions
+float PVM_single::isoterm(const int& pt,const float& _d)const{
+  return(std::exp(-bvals(1,pt)*_d));
+}
+float PVM_single::anisoterm(const int& pt,const float& _d,const ColumnVector& x)const{
+  float dp = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  return(std::exp(-bvals(1,pt)*_d*dp*dp));
+}
+float PVM_single::bvecs_fibre_dp(const int& pt,const float& _th,const float& _ph)const{
+  float angtmp = cos(_ph-beta(pt))*sinalpha(pt)*sin(_th) + cosalpha(pt)*cos(_th);
+  return(angtmp*angtmp);
+}
+float PVM_single::bvecs_fibre_dp(const int& pt,const ColumnVector& x)const{
+  float dp = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  return(dp*dp);
+}
+// 1st order derivatives
+float PVM_single::isoterm_d(const int& pt,const float& _d)const{
+  return(-bvals(1,pt)*std::exp(-bvals(1,pt)*_d));
+}
+float PVM_single::anisoterm_d(const int& pt,const float& _d,const ColumnVector& x)const{
+  float dp = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  return(-bvals(1,pt)*dp*dp*std::exp(-bvals(1,pt)*_d*dp*dp));
+}
+float PVM_single::anisoterm_th(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const{
+  float dp  = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  float dp1 = cos(_th)*(bvecs(1,pt)*cos(_ph) + bvecs(2,pt)*sin(_ph)) - bvecs(3,pt)*sin(_th);
+  return(-2*bvals(1,pt)*_d*dp*dp1*std::exp(-bvals(1,pt)*_d*dp*dp));
+}
+float PVM_single::anisoterm_ph(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const{
+  float dp  = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  float dp1 = sin(_th)*(-bvecs(1,pt)*sin(_ph) + bvecs(2,pt)*cos(_ph));
+  return(-2*bvals(1,pt)*_d*dp*dp1*std::exp(-bvals(1,pt)*_d*dp*dp));
+}
+// 2nd order derivatives
+float PVM_single::isoterm_dd(const int& pt,const float& _d)const{
+  return(bvals(1,pt)*bvals(1,pt)*std::exp(-bvals(1,pt)*_d));
+}
+float PVM_single::anisoterm_dd(const int& pt,const float& _d,const ColumnVector& x)const{
+  float dp = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  dp *= dp;
+  return(bvals(1,pt)*dp*bvals(1,pt)*dp*std::exp(-bvals(1,pt)*_d*dp));
+}
+float PVM_single::anisoterm_dth(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const{
+  float dp  = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  float dp1 = cos(_th)*(bvecs(1,pt)*cos(_ph) + bvecs(2,pt)*sin(_ph)) - bvecs(3,pt)*sin(_th);
+  return( -2*bvals(1,pt)*dp*dp1*(1-bvals(1,pt)*_d*dp*dp)*std::exp(-bvals(1,pt)*_d*dp*dp) );
+}
+float PVM_single::anisoterm_dph(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const{
+  float dp  = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  float dp1 = sin(_th)*(-bvecs(1,pt)*sin(_ph) + bvecs(2,pt)*cos(_ph));
+  return( -2*bvals(1,pt)*dp*dp1*(1-bvals(1,pt)*_d*dp*dp)*std::exp(-bvals(1,pt)*_d*dp*dp) );
+}
+float PVM_single::anisoterm_thth(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const{
+  float dp  = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  return( -2*bvals(1,pt)*_d*std::exp(-bvals(1,pt)*_d*dp*dp)* ( (1-2*bvals(1,pt)*dp*dp) -dp*dp ) );
+}
+float PVM_single::anisoterm_phph(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const{
+  float dp  = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  float dp1 = (1-cos(2*_th))/2.0;
+  float dp2 = -bvecs(1,pt)*x(1) - bvecs(2,pt)*x(2);
+  return( -2*bvals(1,pt)*_d*std::exp(-bvals(1,pt)*_d*dp*dp)* ( (1-2*bvals(1,pt)*dp*dp)*dp1 +dp*dp2 ) );
+}
+float PVM_single::anisoterm_thph(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const{
+  float dp  = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  float dp2 = cos(_th)*(-bvecs(1,pt)*sin(_ph) + bvecs(2,pt)*cos(_ph));
+  return( -2*bvals(1,pt)*_d*std::exp(-bvals(1,pt)*_d*dp*dp)* ( dp*dp2 ) );
+}
+
+
+
+////// NOW FOR MULTISHELL
+// functions
+float PVM_multi::isoterm(const int& pt,const float& _a,const float& _b)const{
+  return(std::exp(-_a*std::log(1+bvals(1,pt)*_b)));
+}
+float PVM_multi::anisoterm(const int& pt,const float& _a,const float& _b,const ColumnVector& x)const{
+  float dp = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  return(std::exp(-_a*std::log(1+bvals(1,pt)*_b*dp*dp)));
+}
+// 1st order derivatives
+float PVM_multi::isoterm_a(const int& pt,const float& _a,const float& _b)const{
+    return(-std::log(1+bvals(1,pt)*_b)*std::exp(-_a*std::log(1+bvals(1,pt)*_b)));
+}
+float PVM_multi::isoterm_b(const int& pt,const float& _a,const float& _b)const{
+      return(-_a*bvals(1,pt)/(1+bvals(1,pt)*_b)*std::exp(-_a*std::log(1+bvals(1,pt)*_b)));
+}
+float PVM_multi::anisoterm_a(const int& pt,const float& _a,const float& _b,const ColumnVector& x)const{
+  float dp = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  return(-std::log(1+bvals(1,pt)*dp*dp*_b)*std::exp(-_a*std::log(1+bvals(1,pt)*dp*dp*_b)));
+}
+float PVM_multi::anisoterm_b(const int& pt,const float& _a,const float& _b,const ColumnVector& x)const{
+  float dp = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  return(-_a*bvals(1,pt)*dp*dp/(1+bvals(1,pt)*dp*dp*_b)*std::exp(-_a*std::log(1+bvals(1,pt)*dp*dp*_b)));
+}
+float PVM_multi::anisoterm_th(const int& pt,const float& _a,const float& _b,const ColumnVector& x,const float& _th,const float& _ph)const{
+  float dp = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  float dp1 = cos(_th)*(bvecs(1,pt)*cos(_ph) + bvecs(2,pt)*sin(_ph)) - bvecs(3,pt)*sin(_th);
+  return(-_a*_b*bvals(1,pt)/(1+bvals(1,pt)*dp*dp*_b)*std::exp(-_a*std::log(1+bvals(1,pt)*dp*dp*_b))*2*dp*dp1);
+}
+float PVM_multi::anisoterm_ph(const int& pt,const float& _a,const float& _b,const ColumnVector& x,const float& _th,const float& _ph)const{
+  float dp = bvecs(1,pt)*x(1)+bvecs(2,pt)*x(2)+bvecs(3,pt)*x(3);
+  float dp1 = sin(_th)*(-bvecs(1,pt)*sin(_ph) + bvecs(2,pt)*cos(_ph));
+  return(-_a*_b*bvals(1,pt)/(1+bvals(1,pt)*dp*dp*_b)*std::exp(-_a*std::log(1+bvals(1,pt)*dp*dp*_b))*2*dp*dp1);
+}
+

diffmodels.h

+/*  Diffusion model fitting
+
+    Timothy Behrens, Saad Jbabdi  - FMRIB Image Analysis Group
+ 
+    Copyright (C) 2005 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#if !defined (diffmodels_h)
+#define diffmodels_h
+
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#define WANT_STREAM
+#define WANT_MATH
+#include <string>
+#include "utils/log.h"
+#include "utils/tracer_plus.h"
+#include "miscmaths/miscmaths.h"
+#include "miscmaths/nonlin.h"
+#include "stdlib.h"
+
+
+
+using namespace NEWMAT;
+using namespace MISCMATHS;
+
+
+#define two_pi 0.636619772
+#define f2x(x) (std::tan((x)/two_pi))
+#define x2f(x) (std::abs(two_pi*std::atan((x))))
+#define bigger(a,b) ((a)>(b)?(a):(b))
+#define smaller(a,b) ((a)>(b)?(b):(a))
+
+
+////////////////////////////////////////////////
+//       DIFFUSION TENSOR MODEL
+////////////////////////////////////////////////
+
+class DTI : public NonlinCF{
+public: 
+  DTI(const ColumnVector& iY,
+      const Matrix& ibvecs,const Matrix& ibvals):bvecs(ibvecs),bvals(ibvals){
+    Y = iY;
+    npts = Y.Nrows();
+    m_v1.ReSize(3);
+    m_v2.ReSize(3);
+    m_v3.ReSize(3);
+    form_Amat();
+    nparams=7;
+  }
+  ~DTI(){}
+  void linfit();
+  void nonlinfit();
+  void calc_tensor_parameters();
+  void sort();
+  void set_data(const ColumnVector& data){Y=data;}
+  float get_fa()const{return m_fa;}
+  float get_md()const{return m_md;}
+  float get_s0()const{return m_s0;}
+  float get_mo()const{return m_mo;}
+  ColumnVector get_v1()const{return m_v1;}
+  ColumnVector get_v2()const{return m_v2;}
+  ColumnVector get_v3()const{return m_v3;}
+  float get_l1()const{return m_l1;}
+  float get_l2()const{return m_l2;}
+  float get_l3()const{return m_l3;}
+  ColumnVector get_eigen()const{ColumnVector x(3);x<<m_l1<<m_l2<<m_l3;return x;}
+  ColumnVector get_tensor()const{
+    ColumnVector x(6);
+    x << m_tens(1,1)
+      << m_tens(2,1)
+      << m_tens(3,1)
+      << m_tens(2,2)
+      << m_tens(3,2)
+      << m_tens(3,3);
+    return x;
+  }
+  ColumnVector get_v(const int& i)const{if(i==1)return m_v1;else if(i==2)return m_v2;else return m_v3;}
+  ColumnVector get_prediction()const;
+  SymmetricMatrix get_covar()const{return m_covar;}
+  ColumnVector get_data()const{return Y;}
+  Matrix get_Amat()const{return Amat;}
+
+  // derivatives of tensor functions w.r.t. tensor parameters
+  ReturnMatrix calc_fa_grad(const ColumnVector& _tens)const;
+  float calc_fa_var()const;
+  ColumnVector calc_md_grad(const ColumnVector& _tens)const;
+  ColumnVector calc_mo_grad(const ColumnVector& _tens)const;
+
+  // conversion between rotation matrix and angles
+  void rot2angles(const Matrix& rot,float& th1,float& th2,float& th3)const;
+  void angles2rot(const float& th1,const float& th2,const float& th3,Matrix& rot)const;
+
+  void print()const{
+    cout << "DTI FIT RESULTS " << endl;
+    cout << "S0   :" << m_s0 << endl;
+    cout << "MD   :" << m_md << endl;
+    cout << "FA   :" << m_fa << endl;
+    cout << "MO   :" << m_mo << endl;
+    ColumnVector x(3);
+    x=m_v1;
+    if(x(3)<0)x=-x;
+    float _th,_ph;cart2sph(x,_th,_ph);
+    cout << "TH   :" << _th*180.0/M_PI << " deg" << endl; 
+    cout << "PH   :" << _ph*180.0/M_PI << " deg" << endl; 
+    cout << "V1   : " << x(1) << " " << x(2) << " " << x(3) << endl;
+  }
+  void form_Amat(){
+    Amat.ReSize(bvecs.Ncols(),7);
+    Matrix tmpvec(3,1), tmpmat;
+    for( int i = 1; i <= bvecs.Ncols(); i++){
+      tmpvec << bvecs(1,i) << bvecs(2,i) << bvecs(3,i);
+      tmpmat = tmpvec*tmpvec.t()*bvals(1,i);
+      Amat(i,1) = tmpmat(1,1);
+      Amat(i,2) = 2*tmpmat(1,2);
+      Amat(i,3) = 2*tmpmat(1,3);
+      Amat(i,4) = tmpmat(2,2);
+      Amat(i,5) = 2*tmpmat(2,3);
+      Amat(i,6) = tmpmat(3,3);
+      Amat(i,7) = 1;
+    }
+    iAmat = pinv(Amat);
+  }
+  void vec2tens(const ColumnVector& Vec){
+    m_tens.ReSize(3);
+    m_tens(1,1)=Vec(1);
+    m_tens(2,1)=Vec(2);
+    m_tens(3,1)=Vec(3);
+    m_tens(2,2)=Vec(4);
+    m_tens(3,2)=Vec(5);
+    m_tens(3,3)=Vec(6);
+  }
+  void vec2tens(const ColumnVector& Vec,SymmetricMatrix& Tens)const{
+    Tens.ReSize(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);
+  }
+  void tens2vec(const SymmetricMatrix& Tens,ColumnVector& Vec)const{
+    Vec.ReSize(6);
+    Vec<<Tens(1,1)<<Tens(2,1)<<Tens(3,1)<<Tens(2,2)<<Tens(3,2)<<Tens(3,3);
+  }
+
+  // nonlinear fitting routines
+  NEWMAT::ReturnMatrix grad(const NEWMAT::ColumnVector& p)const;
+  boost::shared_ptr<BFMatrix> hess(const NEWMAT::ColumnVector&p,boost::shared_ptr<BFMatrix> iptr)const;
+  double cf(const NEWMAT::ColumnVector& p)const;
+  NEWMAT::ReturnMatrix forwardModel(const NEWMAT::ColumnVector& p)const;
+  
+  ColumnVector rotproduct(const ColumnVector& x,const Matrix& R)const;
+  ColumnVector rotproduct(const ColumnVector& x,const Matrix& R1,const Matrix& R2)const;
+  float anisoterm(const int& pt,const ColumnVector& ls,const Matrix& xx)const;
+  
+private:
+  const Matrix& bvecs;
+  const Matrix& bvals;
+  ColumnVector Y;
+  Matrix Amat,iAmat;
+  int npts,nparams;
+  ColumnVector m_v1,m_v2,m_v3;
+  float m_l1,m_l2,m_l3;
+  float m_fa,m_s0,m_md,m_mo;
+  float m_sse;
+  SymmetricMatrix m_tens;
+  SymmetricMatrix m_covar;
+};
+
+////////////////////////////////////////////////
+//       Partial Volume Models
+////////////////////////////////////////////////
+
+// Generic class
+class PVM {
+public:
+  PVM(const ColumnVector& iY,
+      const Matrix& ibvecs, const Matrix& ibvals,
+      const int& nfibres):Y(iY),bvecs(ibvecs),bvals(ibvals){
+    
+    npts    = Y.Nrows();
+    nfib    = nfibres;
+    
+    cart2sph(ibvecs,alpha,beta);
+    
+    cosalpha.ReSize(npts);
+    sinalpha.ReSize(npts);
+    for(int i=1;i<=npts;i++){
+      sinalpha(i) = sin(alpha(i));
+      cosalpha(i) = cos(alpha(i));
+    }
+    
+  }
+  virtual ~PVM(){}
+  
+  // PVM virtual routines
+  virtual void fit()  = 0;
+  virtual void sort() = 0;
+  virtual void print()const = 0;
+  virtual void print(const ColumnVector& p)const = 0;
+  
+  virtual ReturnMatrix get_prediction()const = 0;
+
+protected:
+  const ColumnVector& Y;
+  const Matrix& bvecs;
+  const Matrix& bvals;
+  ColumnVector alpha;
+  ColumnVector sinalpha;
+  ColumnVector cosalpha;
+  ColumnVector beta;
+  
+  float npts;
+  int   nfib;
+
+};
+
+// Model 1 : mono-exponential (for single shell)
+class PVM_single : public PVM, public NonlinCF {
+public:
+  PVM_single(const ColumnVector& iY,
+	     const Matrix& ibvecs, const Matrix& ibvals,
+	     const int& nfibres):PVM(iY,ibvecs,ibvals,nfibres){
+
+    nparams = nfib*3 + 2;
+
+    m_f.ReSize(nfib);
+    m_th.ReSize(nfib);
+    m_ph.ReSize(nfib);
+  }
+  ~PVM_single(){}
+
+  // routines from NonlinCF
+  NEWMAT::ReturnMatrix grad(const NEWMAT::ColumnVector& p)const;
+  boost::shared_ptr<BFMatrix> hess(const NEWMAT::ColumnVector&p,boost::shared_ptr<BFMatrix> iptr)const;
+  double cf(const NEWMAT::ColumnVector& p)const;
+  NEWMAT::ReturnMatrix forwardModel(const NEWMAT::ColumnVector& p)const;
+
+  // other routines
+  void fit();
+  void sort();
+  void fix_fsum();
+  void print()const{
+    cout << "PVM (Single) FIT RESULTS " << endl;
+    cout << "S0   :" << m_s0 << endl;
+    cout << "D    :" << m_d << endl;
+    for(int i=1;i<=nfib;i++){
+      cout << "F" << i << "   :" << m_f(i) << endl;
+      ColumnVector x(3);
+      x << sin(m_th(i))*cos(m_ph(i)) << sin(m_th(i))*sin(m_ph(i)) << cos(m_th(i));
+      if(x(3)<0)x=-x;
+      float _th,_ph;cart2sph(x,_th,_ph);
+      cout << "TH" << i << "  :" << _th*180.0/M_PI << " deg" << endl; 
+      cout << "PH" << i << "  :" << _ph*180.0/M_PI << " deg" << endl; 
+      cout << "DIR" << i << "   : " << x(1) << " " << x(2) << " " << x(3) << endl;
+    }
+  }
+  void print(const ColumnVector& p)const{
+    cout << "PARAMETER VALUES " << endl;
+    cout << "S0   :" << p(1) << endl;
+    cout << "D    :" << p(2) << endl;
+    for(int i=3,ii=1;ii<=nfib;i+=3,ii++){
+      cout << "F" << ii << "   :" << x2f(p(i)) << endl;
+      cout << "TH" << ii << "  :" << p(i+1)*180.0/M_PI << " deg" << endl; 
+      cout << "PH" << ii << "  :" << p(i+2)*180.0/M_PI << " deg" << endl; 
+    }
+  }
+
+  float get_s0()const{return m_s0;}
+  float get_d()const{return m_d;}
+  ColumnVector get_f()const{return m_f;}
+  ColumnVector get_th()const{return m_th;}
+  ColumnVector get_ph()const{return m_ph;}
+  float get_f(const int& i)const{return m_f(i);}
+  float get_th(const int& i)const{return m_th(i);}
+  float get_ph(const int& i)const{return m_ph(i);}
+  ReturnMatrix get_prediction()const;
+
+  // useful functions for calculating signal and its derivatives
+  // functions
+  float isoterm(const int& pt,const float& _d)const;
+  float anisoterm(const int& pt,const float& _d,const ColumnVector& x)const;
+  float bvecs_fibre_dp(const int& pt,const float& _th,const float& _ph)const;
+  float bvecs_fibre_dp(const int& pt,const ColumnVector& x)const;
+  // 1st order derivatives
+  float isoterm_d(const int& pt,const float& _d)const;
+  float anisoterm_d(const int& pt,const float& _d,const ColumnVector& x)const;
+  float anisoterm_th(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const;
+  float anisoterm_ph(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const;
+  // 2nd order derivatives
+  float isoterm_dd(const int& pt,const float& _d)const;
+  float anisoterm_dd(const int& pt,const float& _d,const ColumnVector& x)const;
+  float anisoterm_dth(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const;
+  float anisoterm_dph(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const;
+  float anisoterm_thth(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const;
+  float anisoterm_phph(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const;
+  float anisoterm_thph(const int& pt,const float& _d,const ColumnVector& x,const float& _th,const float& _ph)const;
+
+
+private:  
+  int   nparams;
+  float m_s0;
+  float m_d;
+  ColumnVector m_f;
+  ColumnVector m_th;
+  ColumnVector m_ph;
+};
+
+
+
+////////////////////////////////////////////////
+//       Partial Volume Models
+////////////////////////////////////////////////
+
+// Model 2 : non-mono-exponential (for multiple shells)
+class PVM_multi : public PVM, public NonlinCF {
+public:
+  PVM_multi(const ColumnVector& iY,
+	    const Matrix& ibvecs, const Matrix& ibvals,
+	    const int& nfibres):PVM(iY,ibvecs,ibvals,nfibres){
+
+    nparams = nfib*3 + 3;
+
+    m_f.ReSize(nfib);
+    m_th.ReSize(nfib);
+    m_ph.ReSize(nfib);
+  }
+  ~PVM_multi(){}
+
+  // routines from NonlinCF
+  NEWMAT::ReturnMatrix grad(const NEWMAT::ColumnVector& p)const;
+  boost::shared_ptr<BFMatrix> hess(const NEWMAT::ColumnVector&p,boost::shared_ptr<BFMatrix> iptr)const;
+  double cf(const NEWMAT::ColumnVector& p)const;
+  NEWMAT::ReturnMatrix forwardModel(const NEWMAT::ColumnVector& p)const;
+
+  // other routines
+  void fit();
+  void sort();
+  void fix_fsum();
+  void print()const{
+    cout << "PVM (MULTI) FIT RESULTS " << endl;
+    cout << "S0    :" << m_s0 << endl;
+    cout << "D     :" << m_d << endl;
+    cout << "D_STD :" << m_d_std << endl;
+    for(int i=1;i<=nfib;i++){
+      cout << "F" << i << "    :" << m_f(i) << endl;
+      cout << "TH" << i << "   :" << m_th(i) << endl; 
+      cout << "PH" << i << "   :" << m_ph(i) << endl; 
+    }
+  }
+  void print(const ColumnVector& p)const{
+    cout << "PARAMETER VALUES " << endl;
+    cout << "S0    :" << p(1) << endl;
+    cout << "D     :" << p(2) << endl;
+    cout << "D_STD :" << p(3) << endl;
+    for(int i=3,ii=1;ii<=nfib;i+=3,ii++){
+      cout << "F" << ii << "    :" << x2f(p(i)) << endl;
+      cout << "TH" << ii << "   :" << p(i+1) << endl; 
+      cout << "PH" << ii << "   :" << p(i+2) << endl; 
+    }
+  }
+
+  float get_s0()const{return m_s0;}
+  float get_d()const{return m_d;}
+  float get_d_std()const{return m_d_std;}
+  ColumnVector get_f()const{return m_f;}
+  ColumnVector get_th()const{return m_th;}
+  ColumnVector get_ph()const{return m_ph;}
+  float get_f(const int& i)const{return m_f(i);}
+  float get_th(const int& i)const{return m_th(i);}
+  float get_ph(const int& i)const{return m_ph(i);}
+
+  ReturnMatrix get_prediction()const;
+
+  // useful functions for calculating signal and its derivatives
+  // functions
+  float isoterm(const int& pt,const float& _a,const float& _b)const;
+  float anisoterm(const int& pt,const float& _a,const float& _b,const ColumnVector& x)const;
+  // 1st order derivatives
+  float isoterm_a(const int& pt,const float& _a,const float& _b)const;
+  float anisoterm_a(const int& pt,const float& _a,const float& _b,const ColumnVector& x)const;
+  float isoterm_b(const int& pt,const float& _a,const float& _b)const;
+  float anisoterm_b(const int& pt,const float& _a,const float& _b,const ColumnVector& x)const;
+  float anisoterm_th(const int& pt,const float& _a,const float& _b,const ColumnVector& x,const float& _th,const float& _ph)const;
+  float anisoterm_ph(const int& pt,const float& _a,const float& _b,const ColumnVector& x,const float& _th,const float& _ph)const;
+  
+private:
+  int   nparams;
+  float m_s0;
+  float m_d;
+  float m_d_std;
+  ColumnVector m_f;
+  ColumnVector m_th;
+  ColumnVector m_ph;
+};
+
+
+
+#endif