;
;  Files in this directory contain the IDL routines used in the fitting the 3-exp model to MGRE data
;  and the fitting the ML3 model to the resulting amplitudes 
;
; cplxlib, type_of, trim, regress are general support routines 
; plin_filter, pfit, polyphase, polyphase3d, fix_inv_phase are used in pre-processing of MGE data prior to fitting
; fit_3exp1_cplx contains the fitting routines for the 3exp model to the refence data
; fit_3exp1_amp fits the exp model amplitudes with fixed non-linear parameters, used for all data with preparation pulse (and reference again, for consistency)
; sim_3pool_rf contains the simulation and fitting of the ML3 model
; exchange3n contains the exchange model used in the simulation
;
; to run, define !i:
defsysv,'!i',complex(0,1)
; and compile required routines:
.r exchange3n
;  etc.
;
; Example usage: 
; 
; MGRE fitting part. im= complex data for one subject, mr a mask for selection of SCC and CSF roi's
;
im= [x,y,z,te,reps]  ; reps= contatination of ref,inv,dbl & jh1 measurements, including Ti's and repetitions
mr= [x,y,z,SCC|CSF]
;
imf= plin_filter(im,indgen(41),te,/no_sort)   ; filter all images together: the spatial fit is averaged over all and constant
imi=imf[*,*,*,*,1:5]
mrf= fix_inv_phase(im[*,*,*,*,1:5],imi,mr)
imf[*,*,*,*,1:5]= imi
sig= roi_sum(imf,mrf,/no_pfit)
fref= fit_3exp1_cplx_steps(te,sig[*,0,0]*exp(-!i*arg(sig[*,0,1])),f_init,/res)
;
;     ;fix phase of sig. The correction with CSF ROI also takes inversion phase off, restore where CSF < 0
;
sig3= sig[*,*,0]*exp(-!i*arg(sig[*,*,1]))
for i=1,4 do sig3[*,i]*= -1.0
fit3a= fltarr(6,21)
for i=0,20 do fit[*,i]= fit_3exp1_amp(te,sig3[*,i],fref[0:9,n],/res)  ; 21= 1 ref, 2*5 Ti's inv, 5 Ti's dbl, 5 Ti's jhi
;
; data for ML3 fitting derived from fit:
 data= [[reform(fit[0,1:*])],[reform(fit[1,1:*]+fit[2,1:*])]/((fltarr(20)+1)#[fit[0,0],fit[1,0]+fit[2,0]]) ; ow= axonal+ interstitial, all scans normalized by ref
;
; average over repetitions and subjects 
; estimate variance from sum of differences from averages, divided by degrees of freedom
; SE= sqrt(variance/number_of_averages)
;
;
; ML3 fitting part:
;
; Typical values for parameters:
;
r2= [16.0, 0.13, 0.033]
w=  [ 4.398, 0.22, 0.0]
rfinv= rfinv_array()
;
par14= [  0.312 ,  1.84,   1.84,   0.364,    4.650,  16.40] = [v0, r1[0], r1[1], r1[2], x01, x12] in units of 1/s
;
;   to calculate residue for one set of parameters (and return simalated signals in 'sim'):
var= res_sim_3pool_rf(data,sde,rfinv,5,[0.32,0.12,0.88],w,r1,r2,x01,x12,num=14,sim=sim,par= par14,/te_zero)
;
;   to optimize fitting parameters:
par_fit= fit_3pool(data,sde,rfinv,5,[0.326,0.12,0.88],w,[1.1,1.],r2,4e-3,10e-3,14,/five,par=par14,/te_zero)
;
; data & sde = fltarr(5,2,3) = [delay_time, pool (MW/OW), experiment (IR,2IR,DIR) ] :
;     amplitudes (and standard errors) of myelin water and other water in the three experiments normalized by the reference 
;
; Please contact Peter van Gelderen: gelderen@nih.gov for any questions
;