From 3368e0e478d6b1dc404f4b0a70b956544e78c5b5 Mon Sep 17 00:00:00 2001
From: Saad Jbabdi <saad@fmrib.ox.ac.uk>
Date: Tue, 13 Oct 2009 11:04:48 +0000
Subject: [PATCH] diffusion models
---
diffmodels.cc | 1154 +++++++++++++++++++++++++++++++++++++++++++++++++
diffmodels.h | 402 +++++++++++++++++
2 files changed, 1556 insertions(+)
create mode 100644 diffmodels.cc
create mode 100644 diffmodels.h
diff --git a/diffmodels.cc b/diffmodels.cc
new file mode 100644
index 0000000..669ba89
--- /dev/null
+++ b/diffmodels.cc
@@ -0,0 +1,1154 @@
+/* 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);
+}
+
diff --git a/diffmodels.h b/diffmodels.h
new file mode 100644
index 0000000..25aac74
--- /dev/null
+++ b/diffmodels.h
@@ -0,0 +1,402 @@
+/* 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
--
GitLab