From f4bfb9a46fb50ceec0d9c1578591f8f13df89fbc Mon Sep 17 00:00:00 2001
From: Matthew Webster <mwebster@fmrib.ox.ac.uk>
Date: Fri, 21 Jan 2011 13:23:22 +0000
Subject: [PATCH] revert to 1.5
---
diffmodels.cc | 1150 +++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 1150 insertions(+)
create mode 100644 diffmodels.cc
diff --git a/diffmodels.cc b/diffmodels.cc
new file mode 100644
index 0000000..5fb0b6c
--- /dev/null
+++ b/diffmodels.cc
@@ -0,0 +1,1150 @@
+/* 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) = dti.get_fa()<1?f2x(dti.get_fa()):f2x(0.95); // 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;
+ }
+ if (m_include_f0)
+ start(nparams)=f2x(0.001);
+
+ // 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 = std::abs(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);
+ }
+ if (m_include_f0)
+ m_f0=x2f(final_par(nparams));
+ 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;
+ if (m_include_f0)
+ sumf=m_f0;
+ for(int i=1;i<=nfib;i++){
+ sumf+=m_f(i);
+ if(sumf>=1){for(int j=i;j<=nfib;j++)m_f(j)=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);
+ }
+ if (m_include_f0)
+ p(nparams)=f2x(m_f0);
+ 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 = std::abs(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());
+ }
+ if (m_include_f0){
+ float temp_f0=x2f(p(nparams));
+ pred(i) = p(1)*(temp_f0+(1-sumf-temp_f0)*isoterm(i,_d)+val);
+ }
+ else
+ 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 = std::abs(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());
+ }
+ if (m_include_f0){
+ float temp_f0=x2f(p(nparams));
+ err = (p(1)*(temp_f0+(1-sumf-temp_f0)*isoterm(i,_d)+err) - Y(i));
+ }
+ else
+ 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 = std::abs(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(2)>0?1.0:-1.0)*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));
+ }
+ if (m_include_f0){
+ float temp_f0=x2f(p(nparams));
+ //derivative with respect to f0
+ J(i,nparams)= p(1)*(1-isoterm(i,_d)) * two_pi*sign(p(nparams))*1/(1+p(nparams)*p(nparams));
+ sig=p(1)*(temp_f0+(1-sumf-temp_f0)*isoterm(i,_d)+sig);
+ J(i,2) += (p(2)>0?1.0:-1.0)*p(1)*(1-sumf-temp_f0)*isoterm_d(i,_d);
+ }
+ else{
+ sig = p(1)*((1-sumf)*isoterm(i,_d)+sig);
+ J(i,2) += (p(2)>0?1.0:-1.0)*p(1)*(1-sumf)*isoterm_d(i,_d);
+ }
+ diff(i) = sig - Y(i);
+ J(i,1) = sig/p(1);
+ }
+
+ 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 = std::abs(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);
+ 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(2)>0?1.0:-1.0)*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));
+ }
+ if (m_include_f0){
+ float temp_f0=x2f(p(nparams));
+ //derivative with respect to f0
+ J(i,nparams)= p(1)*(1-isoterm(i,_d)) * two_pi*sign(p(nparams))*1/(1+p(nparams)*p(nparams));
+ sig=p(1)*(temp_f0+(1-sumf-temp_f0)*isoterm(i,_d)+sig);
+ J(i,2) += (p(2)>0?1.0:-1.0)*p(1)*(1-sumf-temp_f0)*isoterm_d(i,_d);
+ }
+ else{
+ sig = p(1)*((1-sumf)*isoterm(i,_d)+sig);
+ J(i,2) += (p(2)>0?1.0:-1.0)*p(1)*(1-sumf)*isoterm_d(i,_d);
+ }
+ J(i,1) = sig/p(1);
+ }
+
+ 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);
+}
+
+
+
+
+////////////////////////////////////////////////
+// PARTIAL VOLUME MODEL - MULTIPLE SHELLS
+////////////////////////////////////////////////
+
+void PVM_multi::fit(){
+
+ // initialise with simple pvm
+ PVM_single pvm1(Y,bvecs,bvals,nfib);
+ pvm1.fit();
+
+ float _a,_b;
+ _a = 1.0; // start with d=d_std
+ _b = pvm1.get_d();
+
+ ColumnVector start(nparams);
+ start(1) = pvm1.get_s0();
+ start(2) = _a;
+ start(3) = _b;
+ for(int i=1,ii=4;i<=nfib;i++,ii+=3){
+ start(ii) = pvm1.get_f(i);
+ start(ii+1) = pvm1.get_th(i);
+ start(ii+2) = pvm1.get_ph(i);
+ }
+
+ // 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 = std::abs(final_par(2)*final_par(3));
+ m_d_std = std::sqrt(std::abs(final_par(2)*final_par(3)*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(j)=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; // =1/beta
+ 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 = std::abs(p(2));
+ float _b = std::abs(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 = std::abs(p(2));
+ float _b = std::abs(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 = (std::abs(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 = std::abs(p(2));
+ float _b = std::abs(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(2)>0?1.0:-1.0)*std::abs(p(1))*fs(k)*anisoterm_a(i,_a,_b,xx);
+ // beta
+ J(i,3) += (p(3)>0?1.0:-1.0)*std::abs(p(1))*fs(k)*anisoterm_b(i,_a,_b,xx);
+ // f
+ J(i,kk) = std::abs(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) = std::abs(p(1))*fs(k)*anisoterm_th(i,_a,_b,xx,p(kk+1),p(kk+2));
+ // ph
+ J(i,kk+2) = std::abs(p(1))*fs(k)*anisoterm_ph(i,_a,_b,xx,p(kk+1),p(kk+2));
+ }
+ sig = std::abs(p(1))*((1-sumf)*isoterm(i,_a,_b)+sig);
+ diff(i) = sig - Y(i);
+ // other stuff for derivatives
+ J(i,1) = (p(1)>0?1.0:-1.0)*sig/p(1);
+ J(i,2) += (p(2)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf)*isoterm_a(i,_a,_b);
+ J(i,3) += (p(3)>0?1.0:-1.0)*std::abs(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 = std::abs(p(2));
+ float _b = std::abs(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(2)>0?1.0:-1.0)*std::abs(p(1))*fs(k)*anisoterm_a(i,_a,_b,xx);
+ // beta
+ J(i,3) += (p(3)>0?1.0:-1.0)*std::abs(p(1))*fs(k)*anisoterm_b(i,_a,_b,xx);
+ // f
+ J(i,kk) = std::abs(p(1))*(anisoterm(i,_a,_b,xx)-isoterm(i,_a,_b)) * cov;
+ // th
+ J(i,kk+1) = std::abs(p(1))*fs(k)*anisoterm_th(i,_a,_b,xx,p(kk+1),p(kk+2));
+ // ph
+ J(i,kk+2) = std::abs(p(1))*fs(k)*anisoterm_ph(i,_a,_b,xx,p(kk+1),p(kk+2));
+ }
+ sig = std::abs(p(1))*((1-sumf)*isoterm(i,_a,_b)+sig);
+ diff(i) = sig - Y(i);
+ // other stuff for derivatives
+ J(i,1) = (p(1)>0?1.0:-1.0)*sig/p(1);
+ J(i,2) += (p(2)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf)*isoterm_a(i,_a,_b);
+ J(i,3) += (p(3)>0?1.0:-1.0)*std::abs(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);
+}
+
--
GitLab