/*----------------------------------------------------------------------
This is for the 3-parameter DEMRI model from Bob Cox and Steve Fung,
implemented for John Butman and Steve Fung.
Author: Rick Reynolds, Daniel Glen 22 Oct 2006
*----------------------------------------------------------------------
*/
#include <stdlib.h>
#include "NLfit_model.h"
extern int AFNI_needs_dset_ijk(void) ;
#define M_D3_R1_MIN 0.0 /* these are exclusive limits */
#define M_D3_R1_MAX 1000.0
#define M_D3_R_MIN 0.001
#define M_D3_R_MAX 10000.0
#define M_D3_THETA_MIN 0.0
#define M_D3_THETA_MAX 90.0
#define M_D3_TR_MIN 0.0
#define M_D3_TR_MAX 10000.0
#define M_D3_TF_MIN 0.1
#define M_D3_TF_MAX 30.0
#define M_D3_NFIRST 5
#define EPSILON 0.0001
#define Mmmmmm_PIiiiii 3.1415926535897932385
void signal_model
(
float * params, /* parameters for signal model */
int ts_len, /* length of time series data */
float ** x_array, /* independent variable matrix */
float * ts_array /* estimated signal model time series */
);
typedef struct
{
float K, kep, ve, fvp; /* fit params (one of kep or ve given) */
float r1, RIB, RIT; /* given params (via env) */
float theta, TR, TF; /* TR & inter-frame TR (TR of input dset) */
float cos0; /* cos(theta) */
int nfirst; /* num TRs used to compute mean Mp,0 */
int ijk; /* voxel index*/
int debug;
double * comp; /* computation data, and elist */
double * elist; /* (for easy allocation, etc.) */
float * mcp; /* (for easy allocation, etc.) */
} demri_params;
static int g_use_ve = 0; /* can be modified in initialize model() */
static THD_3dim_dataset *dset_R1I = NULL; /* input 3d+time data set */
static float *R1I_data_ptr = NULL;
static int alloc_param_arrays(demri_params * P, int len);
static int compute_ts (demri_params *P, float *ts, int ts_len);
static int c_from_ct_cp (demri_params * P, int len);
static int convert_mp_to_cp (demri_params * P, int mp_len);
static int ct_from_cp (demri_params * P, int len);
static int disp_demri_params(char * mesg, demri_params * p );
static int get_env_params (demri_params * P);
static int get_Mp_array (demri_params * P, int * mp_len);
static int get_time_in_seconds(char * str, float * time);
static int model_help (void);
static int Mx_from_R1 (demri_params * P, float * ts, int len);
static int R1_from_c (demri_params * P, int len);
/*----------------------------------------------------------------------
Initialize MODEL_interface struct with default values and model name.
*/
DEFINE_MODEL_PROTOTYPE /* for C++ */
MODEL_interface * initialize_model ()
{
MODEL_interface * M;
if(AFNI_yesenv("AFNI_MODEL_HELP_DEMRI_3") ||
AFNI_yesenv("AFNI_MODEL_HELP_ALL")) model_help();
if( AFNI_yesenv("AFNI_MODEL_D3_USE_VE") )
{
fprintf(stderr,"-d will use Ve for param #2, to set k_ep\n");
g_use_ve = 1;
}
/* get space for new struct */
M = (MODEL_interface *)malloc(sizeof(MODEL_interface));
strcpy(M->label, "demri_3"); /* 3-param DEMRI */
M->model_type = MODEL_SIGNAL_TYPE; /* signal model (not noise) */
M->params = 3; /* number of user parameters */
/* set param labels */
strcpy(M->plabel[0], "K_trans");
if( g_use_ve ) strcpy(M->plabel[1], "Ve");
else strcpy(M->plabel[1], "k_ep");
strcpy(M->plabel[2], "f_pv");
/* set min/max constraints */
M->min_constr[0] = 0.0; M->max_constr[0] = 0.99;
M->min_constr[1] = 0.0; M->max_constr[1] = 0.99;
M->min_constr[2] = 0.0; M->max_constr[2] = 0.99;
M->call_func = &signal_model; /* set the signal model generator callback */
return (M);
}
/*----------------------------------------------------------------------
model the time series, called repeatedly
*/
void signal_model (
float * params, /* parameters for signal model */
int ts_len, /* length of time series data */
float ** x_array, /* independent variable matrix */
float * ts_array /* estimated signal model time series */
)
{
static demri_params P = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0, 0, 0, NULL, NULL, NULL };
static int first_call = 1;
int mp_len; /* length of mcp list */
int rv, errs = 0;
if( first_call > 1 ) return; /* fail, and be quick about it */
/* first time here, get set params, and process Mp data */
if( first_call )
{
if( get_env_params( &P ) ||
get_Mp_array(&P, &mp_len) ||
alloc_param_arrays(&P, mp_len) ||
convert_mp_to_cp(&P, mp_len) )
errs++; /* bad things man, bad things */
/* verify TR (maybe we don't need TR) */
if( !errs && P.TR != x_array[1][1] - x_array[1][0] )
fprintf(stderr, "** warning: TR (%f) != x_array time (%f)\n",
P.TR, x_array[1][1] - x_array[1][0]);
/* verify time series length */
if( !errs && mp_len != ts_len )
{
fprintf(stderr,
"** mp_len (%d) != ts_len (%d), going home in snit.\n",
mp_len, ts_len);
errs++;
}
if( P.debug ) disp_demri_params("ready to rock: ", &P);
if( errs )
{
first_call = 2; /* prevent progression, but don't exit */
memset(ts_array, 0, ts_len*sizeof(float));
fprintf(stderr,"** MD3: failing out, and clearing time series\n");
}
else
first_call = 0;
}
/* note passed parameters */
P.K = params[0];
P.fvp = params[2];
if( g_use_ve )
{
P.ve = params[1];
/* it is not clear what to do if P.ve is small */
if( P.ve < 10e-6 )
{
if( P.debug > 2 )
fprintf(stderr, "** warning, ve, K = %f, %f --> kep = %f\n",
P.ve, P.K, P.kep);
memset(ts_array, 0, ts_len*sizeof(float));
return;
}
else
P.kep = P.K / P.ve;
if( P.kep > 10.0 )
{
if( P.debug > 2 )
fprintf(stderr, "** warning, ve, K = %f, %f ---> kep = %f\n",
P.ve, P.K, P.kep);
memset(ts_array, 0, ts_len*sizeof(float));
return;
}
}
else P.kep = params[1];
if( R1I_data_ptr )
{
P.ijk = AFNI_needs_dset_ijk();
P.RIT = R1I_data_ptr[P.ijk]; /* get R1I voxelwise from dataset */
if( P.RIT < 0.02 || P.RIT > 20 )
{
if( P.debug > 1 )
fprintf(stderr,"** warning, bad RIT value %f\n", P.RIT);
memset(ts_array, 0, ts_len*sizeof(float));
return;
/* do something here? panic into error?? */
}
if(P.debug > 3) printf("Voxel index %d, R1I value %f\n", P.ijk, P.RIT);
}
(void)compute_ts( &P, ts_array, ts_len );
if( (rv = thd_floatscan(ts_len, ts_array)) > 0 )
{
static int nbad = 0;
static int nprint = 1;
nbad++;
if( (nbad % nprint) == 0 )
{
char mesg[128];
sprintf(mesg, "\n** MD3: %d bad results (occurance %d):", rv, nbad);
disp_demri_params(mesg, &P);
if( nprint < 10e+6 ) nprint <<= 1; /* slower and slower */
}
memset(ts_array, 0, ts_len*sizeof(float));
}
}
/*----------------------------------------------------------------------
* compute a time series
*
*/
static int compute_ts( demri_params * P, float * ts, int ts_len )
{
if( ct_from_cp(P, ts_len) ) return 1;
if( c_from_ct_cp(P, ts_len) ) return 1;
if( R1_from_c(P, ts_len) ) return 1;
if( Mx_from_R1(P, ts, ts_len) ) return 1;
return 0;
}
/*----------------------------------------------------------------------
compute Ct from Cp and params
Ct[n] = Ktrans * ( sum [ Cp[k] * exp( -kep(n-k)*TF )
------ * ( exp(kep*TF/2) - exp(-kep*TF/2) ) ]
kep + Cp[n] * (1 - exp(-kep*TF/2)) )
where sum is over k=m..n-1 (m is nfirst), and n is over 0..len-1
note: Cp[n] = Ct[n] = 0, for n in {0..m-1}
Let P1 = Ktrans / kep
P2 = exp(kep*TF/2) - exp(-kep*TF/2)
= exp(-P3/2) - exp(P3/2)
P3 = -kep*TF
P4 = 1 - exp(-kep*TF/2)
= 1 - exp(P3/2)
Ct[n] = P1*P2 * sum [ Cp[k] * exp(P3*(n-k)) ] + P1*P4 * Cp[n]
note: in exp_list, max power is (P3*(len-1-m)), m is nfirst
note: Ct is stored in P->comp
** This is the only place that the dataset TR (which we label as TF,
the inter-frame TR) is used in this model.
*/
static int ct_from_cp(demri_params * P, int len)
{
static int first = 1;
double * ct = P->comp;
double * elist = P->elist;
double P12, P3, P14;
double dval;
float * cp = P->mcp;
int k, n;
/* we don't need to fill with zeros every time, but be explicit */
if( first ){ for(n=0; n < P->nfirst; n++) ct[n] = 0.0; first = 0; }
/* assign P*, and then fill elist[] from P3 */
P12 = P->K / P->kep; /* P1 for now */
P3 = -P->kep * P->TF; /* P3 */
dval = exp(P3/2.0);
P14 = P12 * (1.0 - dval); /* P1 * P4 */
P12 = P12 * (1.0/dval - dval); /* P1 * P2 */
/* In a test, float accuracy was not lost until ~1 million products.
Computing powers over the range of a time series should be okay. */
/* fill elist with powers of e(P3*i), i = 1..len-1-nfirst */
elist[0] = 1.0;
dval = exp(P3); /* first val, elist[i] is a power of this */
for( k = 1; k < len - P->nfirst; k++ ) /* fill the list */
elist[k] = dval * elist[k-1];
for( n = P->nfirst; n < len; n++ ) /* ct[k] = 0, for k < nfirst */
{
dval = 0.0; /* dval is sum here */
for( k = P->nfirst; k < n; k++ )
dval += cp[k]*elist[n-k];
ct[n] = P12 * dval + P14 * cp[n];
}
/* return tissue concentration of Gd over time (in ct[] via P->comp[]) */
return 0;
}
/*----------------------------------------------------------------------
C[n] = { 0, for n < nfirst
{ Ct[n] + fpv * Cp[n], n = nfirst..len-1
note: C will replace Ct in P->comp
*/
static int c_from_ct_cp(demri_params * P, int len)
{
double * C = P->comp;
float * cp = P->mcp;
int n;
for( n = P->nfirst; n < len; n++ )
C[n] += P->fvp * cp[n]; /* C already holds Ct */
return 0;
}
/*----------------------------------------------------------------------
R1[n] = RIT + r1 * C[n]
note: R1 will replace C in P->comp
R1[i] will be constant P->RIT over i=0..nfirst-1
*/
static int R1_from_c(demri_params * P, int len)
{
double * R1 = P->comp;
int n;
for( n = P->nfirst; n < len; n++ )
R1[n] = P->RIT + P->r1 * R1[n];
return 0;
}
/*----------------------------------------------------------------------
compute the final time series, M_transpose, from the R1 array
Mx[n] = 1 * (1 - exp(-R1[n] * TR)) * (1 - exp(-RIT*TR)cos0)
---------------------------------------------
(1 - exp(-R1[n] * TR)cos0) * (1 - exp(-RIT*TR))
= (1 - e[n]) * P1c) / [(1 - cos0*e[n]) * P1]
where P1 = 1 - exp(-RIT*TR)
P1c = 1 - exp(-RIT*TR)*cos0
e[n] = exp(-R1[n] * TR)
notes: R1 is in P->comp
The '1' is a placeholder for Mx(t=0), which should have
been factored out of the input time series.
Mx[i] is identically 1, for i = 0..nfirst-1 .
I can see no speed-up for e[n]. :'(
*/
static int Mx_from_R1(demri_params * P, float * ts, int len)
{
double * R1 = P->comp;
double P1, P1c, e, cos0;
int n;
cos0 = P->cos0;
P1 = exp(-P->RIT * P->TR);
P1c = 1 - P1 * cos0;
P1 = 1 - P1;
/* the first nfirst values 1.0 */
for( n = 0; n < P->nfirst; n++ )
ts[n] = 1.0;
/* and compute the last ones */
for( n = P->nfirst; n < len; n++ )
{
e = exp(-R1[n] * P->TR);
ts[n] = (1 - e)*P1c / ( (1-cos0*e) * P1);
}
return 0;
}
/* get required environment vars, return 0 on success (number of errors) */
static int get_env_params(demri_params * P)
{
char * envp;
int errs = 0;
envp = my_getenv("AFNI_MODEL_D3_R1");
if( !envp )
{
fprintf(stderr,"\n** NLfim: need env var AFNI_MODEL_D3_R1\n");
fprintf(stderr," (in 1/(mMol * sec))\n");
errs++;
}
else
{
P->r1 = atof(envp);
if( P->r1 <= M_D3_R1_MIN || P->r1 >= M_D3_R1_MAX )
{
fprintf(stderr, "** error: r1 (%f) is not in (%f, %f)\n",
P->r1, M_D3_R1_MIN, M_D3_R1_MAX);
errs++;
}
}
envp = my_getenv("AFNI_MODEL_D3_RIB");
if( !envp )
{
fprintf(stderr,"\n** NLfim: need env var AFNI_MODEL_D3_RIB\n");
fprintf(stderr," (in seconds (reciprocal will be taken))\n");
errs++;
}
else
{
if( get_time_in_seconds(envp, &P->RIB) )
{
fprintf(stderr,"** cannot process time '%s' for RIB\n", envp);
errs++;
}
else if( P->RIB <= M_D3_R_MIN || P->RIB >= M_D3_R_MAX )
{
fprintf(stderr, "** error: RIB (%f) is not in (%f, %f)\n",
P->RIB, M_D3_R_MIN, M_D3_R_MAX);
errs++;
}
P->RIB = 1.0 / P->RIB; /* and take the reciprocal */
}
envp = my_getenv("AFNI_MODEL_D3_RIT");
if( !envp )
{
fprintf(stderr,"\n** NLfim: need env var AFNI_MODEL_D3_RIT\n");
fprintf(stderr," (in seconds (reciprocal will be taken))\n");
errs++;
}
else
{
if( get_time_in_seconds(envp, &P->RIT) )
{
fprintf(stderr,"** cannot process time '%s' for RIT\n", envp);
errs++;
}
else if( P->RIT <= M_D3_R_MIN || P->RIT >= M_D3_R_MAX )
{
fprintf(stderr, "** error: RIT (%f) is not in (%f, %f)\n",
P->RIT, M_D3_R_MIN, M_D3_R_MAX);
errs++;
}
P->RIT = 1.0 / P->RIT; /* and take the reciprocal */
}
envp = my_getenv("AFNI_MODEL_D3_THETA");
if( !envp )
{
fprintf(stderr,"\n** NLfim: need env var AFNI_MODEL_D3_THETA\n");
fprintf(stderr," (in degrees)\n");
errs++;
}
else
{
P->theta = atof(envp);
P->cos0 = cos(P->theta * Mmmmmm_PIiiiii / 180.0);
if( P->theta <= M_D3_THETA_MIN || P->theta >= M_D3_THETA_MAX )
{
fprintf(stderr, "** error: theta (%f) is not in (%f, %f)\n",
P->theta, M_D3_THETA_MIN, M_D3_THETA_MAX);
errs++;
}
}
envp = my_getenv("AFNI_MODEL_D3_TR");
if( !envp )
{
fprintf(stderr,"\n** NLfim: need env var AFNI_MODEL_D3_TR\n");
fprintf(stderr," (TR in seconds)\n");
errs++;
}
else
{
if( get_time_in_seconds(envp, &P->TR) )
{
fprintf(stderr,"** cannot process time '%s' for TR\n", envp);
errs++;
}
else if( P->TR <= M_D3_TR_MIN || P->TR >= M_D3_TR_MAX )
{
fprintf(stderr, "** error: TR (%f) is not in (%f, %f)\n",
P->TR, M_D3_TR_MIN, M_D3_TR_MAX);
errs++;
}
}
envp = my_getenv("AFNI_MODEL_D3_TF");
if( !envp )
{
fprintf(stderr,"\n** NLfim: need env var AFNI_MODEL_D3_TF\n");
fprintf(stderr," (TF in seconds)\n");
errs++;
}
else
{
if( get_time_in_seconds(envp, &P->TF) )
{
fprintf(stderr,"** cannot process time '%s' for TF\n", envp);
errs++;
}
else if( P->TF <= M_D3_TF_MIN || P->TF >= M_D3_TF_MAX )
{
fprintf(stderr, "** error: TF (%f) is not in (%f, %f)\n",
P->TF, M_D3_TF_MIN, M_D3_TF_MAX);
errs++;
}
}
envp = my_getenv("AFNI_MODEL_D3_DEBUG");
if( envp ) P->debug = atoi(envp);
/* check if non-uniform intrinsic Relaxivity map is assigned */
envp = my_getenv("AFNI_MODEL_D3_R1I_DSET");
if(envp)
{
/* verify R1I dataset existence and open dataset */
dset_R1I = THD_open_one_dataset (envp);
if (dset_R1I == NULL)
{ fprintf(stderr,"Unable to open R1I dataset: %s", envp); return 1; }
DSET_mallocize (dset_R1I);
DSET_load(dset_R1I);
if( !DSET_LOADED((dset_R1I)) )
{ fprintf(stderr,"Can't load dataset %s",envp) ; return 1; }
if( DSET_BRICK_TYPE(dset_R1I, 0) != MRI_float )
{ fprintf(stderr,"dset %s is not of type float\n",envp); return 1; }
R1I_data_ptr = DSET_ARRAY(dset_R1I, 0);
if(P->debug>0) printf("Set R1I_data_ptr\n");
}
else
{ /* should I close any open datasets? */
if(R1I_data_ptr)
{
fprintf(stderr,"** MD3: warning, we should not be here\n");
R1I_data_ptr = NULL;
dset_R1I = NULL;
}
}
if( envp && P->debug>1 && !errs ) disp_demri_params("env params set: ", P);
return errs;
}
static int get_Mp_array(demri_params * P, int * mp_len)
{
MRI_IMAGE * im;
char * envp;
envp = my_getenv("AFNI_MODEL_D3_MP_FILE");
if( !envp )
{
fprintf(stderr,"\n** NLfim: need env var AFNI_MODEL_D3_MP_FILE\n");
fprintf(stderr," (a 1D file of Mp values)\n");
return 1;
}
if( P->debug ) fprintf(stderr,"-d reading Mp file, %s...\n", envp);
im = mri_read_1D(envp);
if( !im )
{
fprintf(stderr,"** failed to open Mp file %s\n", envp);
return 1;
}
P->mcp = MRI_FLOAT_PTR(im); /* do not free this */
*mp_len = im->nx;
if(P->debug>1) fprintf(stderr,"-d Mp (len, ny) = (%d, %d)\n",im->nx,im->ny);
if( ! P->mcp )
{
fprintf(stderr,"** missing data in Mp file %s\n", envp);
return 1;
}
envp = my_getenv("AFNI_MODEL_D3_NFIRST");
if( !envp )
{
fprintf(stderr,"\nNLfim: no AFNI_MODEL_D3_NFIRST, assuming %d\n",
M_D3_NFIRST);
P->nfirst = M_D3_NFIRST;
}
else
{
P->nfirst = atoi(envp);
if( P->nfirst < 0 || P->nfirst > 1000 )
{
fprintf(stderr, "** NLfim: nfirst = %d must be in [0,1000]\n",
P->nfirst);
return 1;
}
}
if( P->debug > 1 )
{
int c;
fprintf(stderr,"-d Mp array is:\n ");
for( c = 0; c < *mp_len; c++ )
fprintf(stderr," %s", MV_format_fval(P->mcp[c]));
fputc('\n', stderr);
}
return 0;
}
/* convert Mp array to Cp(t) array
Cp(t) = 1 * ln[ 1-exp(-RIB*TR)cos0 - (1-exp(-RIB*TR))cos0*Mp(t)/Mp(0) ]
----- [ ---------------------------------------------------- ]
r1*TR [ 1-exp(-RIB*TR)cos0 - (1-exp(-RIB*TR))*Mp(t)/Mp(0) ]
- RIB/r1
subject to Cp(t) >= 0
*/
static int convert_mp_to_cp(demri_params * P, int mp_len)
{
float * mp = P->mcp;
double m0; /* mean of first nfirst values */
double rTR, R_r1; /* first and last simple term */
double ertr, ertr_c0, c0_ertr_c0; /* three 1-exp terms (1 repeat) */
double dval;
int c;
/* use local vars for those in P, for readability */
float r1 = P->r1, RIB = P->RIB, TR = P->TR, cos0 = P->cos0;
int nfirst = P->nfirst;
/* compute m0 equal to mean of first 'nfirst' values */
dval = 0.0;
for(c = 0; c < nfirst; c++)
dval += mp[c];
m0 = dval / nfirst;
if( m0 < EPSILON ) /* negative is bad, too */
{
fprintf(stderr,"** m0 == %s is too small (for my dreams of konquest)\n",
MV_format_fval(m0));
return 1;
}
/* simple terms */
rTR = 1.0/(r1*TR);
R_r1 = RIB/r1;
/* exponential terms */
ertr = 1 - exp(-RIB * TR);
ertr_c0 = 1 - exp(-RIB * TR) * cos0;
c0_ertr_c0 = (1 - exp(-RIB * TR)) * cos0;
if(P->debug > 1)
fprintf(stderr,
"+d mp_len, m0, rTR, R_r1 = %d, %f, %f, %f\n"
" ertr, ertr_c0, c0_ertr_c0 = %f, %f, %f\n"
, mp_len, m0, rTR, R_r1, ertr, ertr_c0, c0_ertr_c0);
/* we don't have to be too fast here, since this is one-time-only */
/* start with setting nfirst terms to 0 */
for( c = 0; c < nfirst; c++ )
mp[c] = 0.0;
/* and compute the remainder of the array */
for( ; c < mp_len; c++)
{
/* start with ratio */
dval = (ertr_c0 - c0_ertr_c0 * mp[c] / m0) /
(ertr_c0 - ertr * mp[c] / m0 );
if( dval < 1.0 ) mp[c] = 0.0; /* if ln < 0, then c < 0, so skip */
else mp[c] = rTR * log(dval) - R_r1;
if( mp[c] < 0.0 ) mp[c] = 0.0; /* don't allow result < 0 */
}
if( P->debug > 1 ) /* maybe print out the list */
{
fprintf(stderr,"+d Cp =");
for( c = 0; c < mp_len; c++ )
fprintf(stderr," %s", MV_format_fval(mp[c]));
fputc('\n', stderr);
}
return 0;
}
static int disp_demri_params( char * mesg, demri_params * p )
{
if( mesg ) fputs( mesg, stderr );
if( !p ) { fprintf(stderr,"demri_params: p == NULL\n"); return 1; }
fprintf(stderr, "demri_params struct at %p:\n"
" K = %f ( K trans (plasma Gd -> tissue Gd) )\n"
" kep = %f ( back-transfer rate ( Gd_t -> Gd_p ) )\n"
" ve = %f ( external cellular volume fraction)\n"
" fvp = %f ( fraction of voxel occupied by blood )\n"
" r1 = %f ( 1/(mMol*seconds) )\n"
" RIB = %f ( 1/seconds )\n"
" RIT = %f ( 1/seconds )\n"
" theta = %f ( degrees )\n"
" TR = %f ( seconds (~0.007) )\n"
" TF = %f ( seconds (~20) )\n"
" cos0 = %f ( cos(theta) )\n"
" nfirst = %d\n\n"
" debug = %d\n"
" ijk = %d\n"
" comp = %p\n"
" elist = %p\n"
" mcp = %p\n"
, p,
p->K, p->kep, p->ve, p->fvp,
p->r1, p->RIB, p->RIT, p->theta, p->TR, p->TF, p->cos0,
p->nfirst, p->debug, p->ijk, p->comp, p->elist, p->mcp);
return 0;
}
static int model_help(void)
{
printf(
"------------------------------------------------------------\n"
"DEMRI - Dynamic (contrast) Enhanced MRI\n"
"\n"
" MODEL demri_3: 3-parameter DEMRI model\n"
"\n"
" model parameters to fit:\n"
" K_trans in [0, 0.99]\n"
" k_ep in [0, 0.99]\n"
" f_pv in [0, 0.99]\n"
"\n"
" model parameters passed via environment variables:\n"
" AFNI_MODEL_D3_R1 : to set r1 in [%f, %f]\n"
" in 1/(mMol*second)\n"
" e.g. 4.8 (@ 1.5T)\n"
"\n"
" AFNI_MODEL_D3_RIB : to set blood R_1,i in [%f, %f]\n"
" ** R_1,i, given as reciprocal, in seconds\n"
" e.g. 1.3 or 1.3s or 1300ms\n"
"\n"
" AFNI_MODEL_D3_RIT : to set tissue R_1,i in [%f, %f]\n"
" ** R_1,i, given as reciprocal, in seconds\n"
" e.g. 0.8 or 0.8s or 800ms\n"
"\n"
" AFNI_MODEL_D3_THETA : to set theta in [%f, %f]\n"
" flip angle, in degrees\n"
" e.g. 30\n"
"\n"
" AFNI_MODEL_D3_TR : to set TR in [%f, %f]\n"
" time between shots, in seconds\n"
" e.g. 0.0076 or 0.0076s or 7.6ms\n"
"\n"
" AFNI_MODEL_D3_TF : to set TF in [%f, %f]\n"
" inter-frame TR (TR of input dataset)\n"
" e.g. 20 or 20sec\n"
"\n"
" ** note: default units are with respect to seconds\n"
"\n"
" environment variables to control Mp(t):\n"
" AFNI_MODEL_D3_MP_FILE : file containing Mp(t) data\n"
" AFNI_MODEL_D3_NFIRST : to set nfirst (injection TR)\n"
" index of input dataset, e.g. 5\n"
"\n"
" optional environment variables:\n"
" AFNI_MODEL_HELP_DEMRI_3 (Y/N) : to get this help\n"
" AFNI_MODEL_D3_USE_VE (Y/N) : param #2 is Ve, not k_ep\n"
" AFNI_MODEL_D3_DEBUG (0..2): to set debug level\n"
"\n"
" -----------------------------------------------\n"
"\n"
" Assumptions:\n"
"\n"
" - nfirst is the number of TRs before Gd injection\n"
" - the input M_trans is normalized:\n"
" M_trans[0] = mean( M_trans[0..nfirst-1] )\n"
" M_trans[k] = M_trans[k] / M_trans[0], k=0..N-1\n"
" - m0 is set to average of first nfirst Mp(t) values\n"
"\n"
" -----------------------------------------------\n"
" sample command-script:\n"
"\n"
" setenv AFNI_MODEL_D3_R1 4.8\n"
" setenv AFNI_MODEL_D3_RIB 1500ms\n"
" setenv AFNI_MODEL_D3_RIT 1100ms\n"
" setenv AFNI_MODEL_D3_TR 8.0ms\n"
" setenv AFNI_MODEL_D3_THETA 30\n"
" setenv AFNI_MODEL_D3_TF 20s\n"
" \n"
" setenv AFNI_MODEL_D3_MP_FILE ROI_mean.1D\n"
" setenv AFNI_MODEL_D3_NFIRST 7\n"
" \n"
" 3dNLfim -input scaled_data.nii \\\n"
" -signal demri_3 \\\n"
" -noise Zero \\\n"
" -sconstr 0 0 .95 \\\n"
" -sconstr 1 0 .95 \\\n"
" -sconstr 2 0 .99 \\\n"
" -nconstr 0 0.0 0.0 \\\n"
" -SIMPLEX \\\n"
" -nabs \\\n"
" -ignore 0 \\\n"
" -nrand 10000 \\\n"
" -mask my_mask+orig \\\n"
" -nbest 10 \\\n"
" -jobs 2 \\\n"
" -voxel_count \\\n"
" -sfit subj7.sfit \\\n"
" -bucket 0 subj7.buc\n"
"\n"
" -----------------------------------------------\n"
" R Reynolds, D Glen, Oct 2006\n"
" thanks to RW Cox\n"
"------------------------------------------------------------\n",
M_D3_R1_MIN, M_D3_R1_MAX,
M_D3_R_MIN, M_D3_R_MAX,
M_D3_R_MIN, M_D3_R_MAX,
M_D3_TR_MIN, M_D3_TR_MAX,
M_D3_THETA_MIN, M_D3_THETA_MAX,
M_D3_TF_MIN, M_D3_TF_MAX
);
return 0;
}
static int alloc_param_arrays(demri_params * P, int len)
{
P->comp = (double *)malloc(len * sizeof(double));
P->elist = (double *)malloc(len * sizeof(double));
if( !P->comp || !P->elist )
{
fprintf(stderr,"** alloc_param_arrays failure for %d doubles\n",len);
return 1;
}
return 0;
}
/* str can be of the form "7.2", "7.2s" or "7.2ms" */
static int get_time_in_seconds( char * str, float * time )
{
char * cp;
*time = strtod(str, &cp);
if( cp == str ) return 1;
if( *cp && !strncmp(cp, "ms", 2) )
*time /= 1000.0;
return 0;
}
syntax highlighted by Code2HTML, v. 0.9.1