Skip to content
Snippets Groups Projects
Commit c1ffb543 authored by Moises Fernandez's avatar Moises Fernandez
Browse files

Routines for SINGLE SHELL OLD model

parent 899904ab
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 451 additions and 0 deletions
CUDA/PVM_single.cu 0 → 100644
+ 451
− 0
View file @ c1ffb543
#include "diffmodels_utils.h"
#include "levenberg_marquardt.cu"
#include "options.h"
//#include <fstream>
/////////////////////////////////////
/////////////////////////////////////
/// PVM_single ///
/////////////////////////////////////
/////////////////////////////////////
__device__
inline double isoterm_PVM_single(const int pt,const double _d,const double *bvals){
return exp(double(-bvals[pt]*_d));
}
__device__
inline double isoterm_d_PVM_single(const int pt,const double _d,const double *bvals){
return (-bvals[pt]*exp(double(-bvals[pt]*_d)));
}
__device__
inline double anisoterm_PVM_single(const int pt,const double _d,const double3 x, const double *bvecs, const double *bvals){
double dp = bvecs[pt]*x.x+bvecs[NDIRECTIONS+pt]*x.y+bvecs[(2*NDIRECTIONS)+pt]*x.z;
return exp(double(-bvals[pt]*_d*dp*dp));
}
__device__
inline double anisoterm_d_PVM_single(const int pt,const double _d,const double3 x,const double *bvecs, const double *bvals){
double dp = bvecs[pt]*x.x+bvecs[NDIRECTIONS+pt]*x.y+bvecs[(2*NDIRECTIONS)+pt]*x.z;
return(-bvals[pt]*dp*dp*exp(double(-bvals[pt]*_d*dp*dp)));
}
__device__
inline double anisoterm_th_PVM_single(const int pt,const double _d,const double3 x, const double _th,const double _ph,const double *bvecs, const double *bvals){
double dp = bvecs[pt]*x.x+bvecs[NDIRECTIONS+pt]*x.y+bvecs[(2*NDIRECTIONS)+pt]*x.z;
double dp1 = (cos(double(_th))*(bvecs[pt]*cos(double(_ph))+bvecs[NDIRECTIONS+pt]*sin(double(_ph)))-bvecs[(2*NDIRECTIONS)+pt]*sin(double(_th)));
return(-2*bvals[pt]*_d*dp*dp1*exp(double(-bvals[pt]*_d*dp*dp)));
}
__device__
inline double anisoterm_ph_PVM_single(const int pt,const double _d,const double3 x, const double _th,const double _ph,const double *bvecs, const double *bvals){
double dp = bvecs[pt]*x.x+bvecs[NDIRECTIONS+pt]*x.y+bvecs[(2*NDIRECTIONS)+pt]*x.z;
double dp1 = sin(double(_th))*(-bvecs[pt]*sin(double(_ph))+bvecs[NDIRECTIONS+pt]*cos(double(_ph)));
return(-2*bvals[pt]*_d*dp*dp1*exp(double(-bvals[pt]*_d*dp*dp)));
}
//in diffmodel.cc
__device__ void fix_fsum_PVM_single( //INPUT
bool m_include_f0,
int nfib,
int nparams,
//INPUT - OUTPUT){
double *params)
{
double sum=0;
if (m_include_f0)
sum=params[nparams-1];
for(int i=0;i<nfib;i++){
sum += params[2+(i*3)];
if(sum>=1){
for(int j=i;j<nfib;j++)
params[2+(j*3)]=FSMALL_gpu;
break;
}
}
}
//in diffmodel.cc
__device__ void sort_PVM_single(int nfib,int nparams,double* params)
{
double temp_f, temp_th, temp_ph;
// Order vector descending using f parameters as index
for(int i=1; i<(nfib); i++){
for(int j=0; j<(nfib-i); j++){
if (params[2+j*3] < params[2+i*3]){
temp_f = params[2+j*3];
temp_th = params[2+j*3+1];
temp_ph = params[2+j*3+2];
params[2+j*3] = params[2+i*3];
params[2+j*3+1] = params[2+i*3+1];
params[2+j*3+2] = params[2+i*3+2];
params[2+i*3] = temp_f;
params[2+i*3+1] = temp_th;
params[2+i*3+2] = temp_ph;
}
}
}
}
//in diffmodels.cc -- for calculate residuals
__device__ void forwardModel_PVM_single( //INPUT
const double* p,
const double* bvecs,
const double* bvals,
const int nfib,
const int nparams,
const bool m_include_f0,
//OUTPUT
double* predicted_signal)
{
for(int i=0;i<NDIRECTIONS;i++){
predicted_signal[i]=0; //pred = 0;
}
double val;
double _d = abs(p[1]);
////////////////////////////////////
double fs[NFIBRES];
double x[NFIBRES*3];
double sumf=0;
double3 x2;
for(int k=0;k<nfib;k++){
int kk = 2+3*k;
fs[k] = x2f_gpu(p[kk]);
sumf += fs[k];
x[k*3] = sin(p[kk+1])*cos(p[kk+2]);
x[k*3+1] = sin(p[kk+1])*sin(p[kk+2]);
x[k*3+2] = cos(p[kk+1]);
}
////////////////////////////////////
for(int i=0;i<NDIRECTIONS;i++){
val = 0.0;
for(int k=0;k<nfib;k++){
x2.x=x[k*3];
x2.y=x[k*3+1];
x2.z=x[k*3+2];
val += fs[k]*anisoterm_PVM_single(i,_d,x2,bvecs,bvals);
}
if (m_include_f0){
double temp_f0=x2f_gpu(p[nparams-1]);
predicted_signal[i] = p[0]*(temp_f0+(1-sumf-temp_f0)*isoterm_PVM_single(i,_d,bvals)+val);
}
else
predicted_signal[i] = p[0]*((1-sumf)*isoterm_PVM_single(i,_d,bvals)+val);
}
}
//in diffmodels.cc -- for calculate residuals
__device__ void get_prediction_PVM_single( //INPUT
const double* params,
const double* bvecs,
const double* bvals,
const int nfib,
const int nparams,
const bool m_include_f0,
//OUTPUT
double* predicted_signal)
{
//m_s0-myparams[0] m_d-myparams[1] m_d_std-myparams[2] m_f-m_th-m_ph-myparams[3,4,5,6 etc..] m_f0-myparams[nparams-1]
double p[NPARAMS];
p[0] = params[0];
p[1] = params[1];
for(int k=0;k<nfib;k++){
int kk = 2+3*k;
p[kk] = f2x_gpu(params[kk]);
p[kk+1] = params[kk+1];
p[kk+2] = params[kk+2];
}
if (m_include_f0)
p[nparams-1]=f2x_gpu(params[nparams-1]);
forwardModel_PVM_single(p,bvecs,bvals,nfib,nparams,m_include_f0,predicted_signal);
}
//cost function PVM_single
__device__ double cf_PVM_single( //INPUT
const double* params,
const double* mdata,
const double* bvecs,
const double* bvals,
const int nparams,
const bool m_include_f0)
{
double cfv = 0.0;
double err;
double _d = abs(params[1]);
double fs[NFIBRES];
double x[NFIBRES*3];
double sumf=0;
double3 x2;
for(int k=0;k<NFIBRES;k++){
int kk = 2+3*(k);
fs[k] = x2f_gpu(params[kk]);
sumf += fs[k];
x[k*3] = sin(params[kk+1])*cos(params[kk+2]);
x[k*3+1] = sin(params[kk+1])*sin(params[kk+2]);
x[k*3+2] = cos(params[kk+1]);
}
for(int i=0;i<NDIRECTIONS;i++){
err = 0.0;
for(int k=0;k<NFIBRES;k++){
x2.x=x[k*3];
x2.y=x[k*3+1];
x2.z=x[k*3+2];
err += fs[k]*anisoterm_PVM_single(i,_d,x2,bvecs,bvals);
}
if(m_include_f0){
double temp_f0=x2f_gpu(params[nparams-1]);
err= (params[0]*((temp_f0+(1-sumf-temp_f0)*isoterm_PVM_single(i,_d,bvals))+err))-mdata[i];
}else{
err = (params[0]*((1-sumf)*isoterm_PVM_single(i,_d,bvals)+err))-mdata[i];
}
cfv += err*err;
}
return(cfv);
}
//gradient function PVM_single
__device__ void grad_PVM_single( //INPUT
const double* params,
const double* mdata,
const double* bvecs,
const double* bvals,
const int nparams,
const bool m_include_f0,
//OUTPUT
double* grad)
{
double _d = abs(params[1]);
double fs[NFIBRES];
double x[NFIBRES*3];
double3 xx;
double sumf=0;
for(int k=0;k<NFIBRES;k++){
int kk = 2+3*(k);
fs[k] = x2f_gpu(params[kk]);
sumf += fs[k];
x[k*3] = sin(params[kk+1])*cos(params[kk+2]);
x[k*3+1] = sin(params[kk+1])*sin(params[kk+2]);
x[k*3+2] = cos(params[kk+1]);
}
double J[NPARAMS];
double diff;
double sig;
for (int p=0;p<nparams;p++) grad[p]=0;
for(int i=0;i<NDIRECTIONS;i++){
sig = 0;
for(int a=0;a<nparams;a++) J[a]=0;
for(int k=0;k<NFIBRES;k++){
int kk = 2+3*(k);
xx.x=x[k*3];
xx.y=x[k*3+1];
xx.z=x[k*3+2];
sig += fs[k]*anisoterm_PVM_single(i,_d,xx,bvecs,bvals);
J[1] += (params[1]>0?1.0:-1.0)*params[0]*fs[k]*anisoterm_d_PVM_single(i,_d,xx,bvecs,bvals);
J[kk] = params[0]*(anisoterm_PVM_single(i,_d,xx,bvecs,bvals)-isoterm_PVM_single(i,_d,bvals)) * two_pi_gpu*sign_gpu(params[kk])*1/(1+params[kk]*params[kk]);
J[kk+1] = params[0]*fs[k]*anisoterm_th_PVM_single(i,_d,xx,params[kk+1],params[kk+2],bvecs,bvals);
J[kk+2] = params[0]*fs[k]*anisoterm_ph_PVM_single(i,_d,xx,params[kk+1],params[kk+2],bvecs,bvals);
}
if(m_include_f0){
double temp_f0=x2f_gpu(params[nparams-1]);
J[nparams-1]= params[0]*(1-isoterm_PVM_single(i,_d,bvals))* two_pi_gpu*sign_gpu(params[nparams-1])*1/(1+params[nparams-1]*params[nparams-1]);
sig= params[0]*((temp_f0+(1-sumf-temp_f0)*isoterm_PVM_single(i,_d,bvals))+sig);
J[1] += (params[1]>0?1.0:-1.0)*params[0]*(1-sumf-temp_f0)*isoterm_d_PVM_single(i,_d,bvals);
}else{
sig = params[0]*((1-sumf)*isoterm_PVM_single(i,_d,bvals)+sig);
J[1] += (params[1]>0?1.0:-1.0)*params[0]*(1-sumf)*isoterm_d_PVM_single(i,_d,bvals);
}
diff = sig - mdata[i];
J[0] = sig/params[0];
for (int p=0;p<nparams;p++) grad[p] += 2*J[p]*diff;
}
}
//hessian function PVM_single
__device__ void hess_PVM_single( //INPUT
const double* params,
const double* bvecs,
const double* bvals,
const int nparams,
const bool m_include_f0,
double* hess)
{
double _d = abs(params[1]);
double fs[NFIBRES];
double x[NFIBRES*3];
double3 xx;
double sumf=0;
for(int k=0;k<NFIBRES;k++){
int kk = 2+3*(k);
fs[k] = x2f_gpu(params[kk]);
sumf += fs[k];
x[k*3] = sin(params[kk+1])*cos(params[kk+2]);
x[k*3+1] = sin(params[kk+1])*sin(params[kk+2]);
x[k*3+2] = cos(params[kk+1]);
}
double J[NPARAMS];
double sig;
for (int p=0;p<nparams;p++){
for (int p2=0;p2<nparams;p2++){
hess[p*nparams+p2] = 0;
}
}
for(int i=0;i<NDIRECTIONS;i++){
sig = 0;
for(int a=0;a<nparams;a++) J[a]=0;
for(int k=0;k<NFIBRES;k++){
int kk = 2+3*(k);
xx.x=x[k*3];
xx.y=x[k*3+1];
xx.z=x[k*3+2];
sig += fs[k]*anisoterm_PVM_single(i,_d,xx,bvecs,bvals);
J[1] += (params[1]>0?1.0:-1.0)*params[0]*fs[k]*anisoterm_d_PVM_single(i,_d,xx,bvecs,bvals);
J[kk] = params[0]*(anisoterm_PVM_single(i,_d,xx,bvecs,bvals)-isoterm_PVM_single(i,_d,bvals)) * two_pi_gpu*sign_gpu(params[kk])*1/(1+params[kk]*params[kk]);
J[kk+1] = params[0]*fs[k]*anisoterm_th_PVM_single(i,_d,xx,params[kk+1],params[kk+2],bvecs,bvals);
J[kk+2] = params[0]*fs[k]*anisoterm_ph_PVM_single(i,_d,xx,params[kk+1],params[kk+2],bvecs,bvals);
}
if(m_include_f0){
double temp_f0=x2f_gpu(params[nparams-1]);
J[nparams-1]= params[0]*(1-isoterm_PVM_single(i,_d,bvals))* two_pi_gpu*sign_gpu(params[nparams-1])*1/(1+params[nparams-1]*params[nparams-1]);
sig=params[0]*((temp_f0+(1-sumf-temp_f0)*isoterm_PVM_single(i,_d,bvals))+sig);
J[1] += (params[1]>0?1.0:-1.0)*params[0]*(1-sumf-temp_f0)*isoterm_d_PVM_single(i,_d,bvals);
}else{
sig = params[0]*((1-sumf)*isoterm_PVM_single(i,_d,bvals)+sig);
J[1] += (params[1]>0?1.0:-1.0)*params[0]*(1-sumf)*isoterm_d_PVM_single(i,_d,bvals);
}
J[0] = sig/params[0];
for (int p=0;p<nparams;p++){
for (int p2=p;p2<nparams;p2++){
hess[p*nparams+p2] += 2*(J[p]*J[p2]);
}
}
}
for (int j=0; j<nparams; j++) {
for (int i=j+1; i<nparams; i++) {
hess[i*nparams+j]=hess[j*nparams+i];
}
}
}
//in diffmodel.cc
extern "C" __global__ void fit_PVM_single_kernel( //INPUT
const double* data,
const double* bvecs,
const double* bvals,
const int nvox,
const int nfib,
const bool m_include_f0,
//INPUT-OUTPUT
double* params)
{
int id = blockIdx.x * blockDim.x + threadIdx.x;
if (id >=nvox) { return; }
int nparams;
if (m_include_f0)
nparams = nfib*3 + 3;
else
nparams = nfib*3 + 2;
double myparams[NPARAMS];
double mydata[NDIRECTIONS];
for(int i=0;i<nparams;i++){
myparams[i]=params[(id*nparams)+i];
}
for(int i=0;i<NDIRECTIONS;i++){
mydata[i]=data[(id*NDIRECTIONS)+i];
}
// do the fit
levenberg_marquardt_PVM_single_gpu(mydata, &bvecs[id*3*NDIRECTIONS], &bvals[id*NDIRECTIONS], nparams, m_include_f0, myparams);
// finalise parameters
//m_s0 in myparams[0] m_d in myparams[1] m_f-m_th-m_ph in myparams[2,3,4,5, etc..] m_f0 in myparams[nparams-1]
myparams[1] = abs(myparams[1]);
for(int k=1;k<=nfib;k++){
int kk = 2 + 3*(k-1);
myparams[kk] = x2f_gpu(myparams[kk]);
}
if (m_include_f0)
myparams[nparams-1]=x2f_gpu(myparams[nparams-1]);
sort_PVM_single(nfib,nparams,myparams);
fix_fsum_PVM_single(m_include_f0,nfib,nparams,myparams);
for(int i=0;i<nparams;i++){
params[id*nparams+i]=myparams[i];
//printf("PARAM[%i]: %.20f\n",i,myparams[i]);
}
}
//in diffmodel.cc
extern "C" __global__ void get_residuals_PVM_single_kernel( //INPUT
const double* data,
const double* params,
const double* bvecs,
const double* bvals,
const int nvox,
const int nfib,
const bool m_include_f0,
const bool* includes_f0,
//OUTPUT
double* residuals)
{
int id = blockIdx.x * blockDim.x + threadIdx.x;
if (id >=nvox) { return; }
int nparams;
if (m_include_f0)
nparams = nfib*3 + 3;
else
nparams = nfib*3 + 2;
bool my_include_f0 = includes_f0[id];
double myparams[NPARAMS];
double mydata[NDIRECTIONS];
for(int i=0;i<nparams;i++){
myparams[i]=params[(id*nparams)+i];
}
for(int i=0;i<NDIRECTIONS;i++){
mydata[i]=data[(id*NDIRECTIONS)+i];
}
double predicted_signal[NDIRECTIONS];
get_prediction_PVM_single(myparams, &bvecs[id*3*NDIRECTIONS], &bvals[id*NDIRECTIONS], nfib, nparams, my_include_f0, predicted_signal);
for(int i=0;i<NDIRECTIONS;i++){ //residuals=m_data-predicted_signal;
residuals[id*NDIRECTIONS+i]= mydata[i] - predicted_signal[i];
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment