Curve Fitting – Mono-Exponential
Given a set of data from an MRI scanner there are several ways to fit the data. Given the assumption that the data is exponential (mono or multi).
The standard signal equation is:
where 
Mono-Exponential Data Fitting
There are many ways to fit mono-exponential data to obtain
First we will create some test data. It will be based on 4 acquired echoes, this can change to any number you wish to use. We will also add on a little noise just to make it interesting. The noise in this case is Gaussian, which is not what is normal for standard magnitude MRI images. That we can deal with later.
1 2 3 4 5 6 7 |
And now display it to make sure we understand how things should look.
Now for the curve fitting. This is the important part. I use Matlab’s function fmincon which is a nonlinear function minimizer. The function we are going to minimize is called “fun” and is the first parameter. There are many linear inequality (A and B) and equality (Aeq and Beq) constraints we can add. For this I am not going to use any. Bounds on our desired variables are added. The lower bound will make sure it stays positive, and the upper bound is to make sure things don’t blow up. The function fmincon can also be used to pass in information to the function “fun” (in our case y_e and te) as that is needed to calculate the L2 norm of the signal intensities at a current iteration. When the minimization is done, the output vector will be the minimum L2 best guess.
1 2 3 4 5 6 7 8 9 10 | x0 = [ 10^5 30 0 ]; % A T2 baseline A = []; B = []; Aeq = []; Beq = []; LB = [ 0 0 0 ]; UB = [ 10^6 200 10^4]; ops = optimset('fmincon'); x = fmincon('fun', x0, A, B, Aeq, Beq, LB, UB, [], ops, y_e, te); A_est = x(1) t2_est = x(2) baseline_est = x(3) |
Now we can reconstruct the data based on the best guess and display it overtop the initial data.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | y_recon = A_est * exp(-te ./ t2_est) + baseline_est; figure(1), hold on; (te, y_recon, 'ro-'); dd a baseline component x0 = [ 10^5 30 0 ]; % A T2 baseline A = []; B = []; Aeq = []; Beq = []; LB = [ 0 0 0 ]; UB = [ 10^6 200 10^4]; ops = optimset('fmincon'); x = fmincon('fun', x0, A, B, Aeq, Beq, LB, UB, [], ops, y_e, te); A_est = x(1) t2_est = x(2) baseline_est = x(3) %% Reconstruct the data. y_recon = A_est * exp(-te ./ t2_est) + baseline_est; figure(1), hold on; plot(te, y_recon, 'r-'); |
The two files to run this experiment are:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 | % Create fake data te = [24 34 44 54]'; y = 2*10^5 * exp(-te ./ 22 ); noise_factor = 5; e = noise_factor * randn(size(y)); y_e = y + e; figure(1),clf; plot(te, y_e, 'x'); title('The fake data'); xlabel('TE (ms)'); ylabel('Signal'); % Add a baseline component x0 = [ 10^5 30 0 ]; % A T2 baseline A = []; B = []; Aeq = []; Beq = []; LB = [ 0 0 0 ]; UB = [ 10^6 200 10^4]; ops = optimset('fmincon'); x = fmincon('fun', x0, A, B, Aeq, Beq, LB, UB, [], ops, y_e, te); A_est = x(1) t2_est = x(2) baseline_est = x(3) %% Reconstruct the data. y_recon = A_est * exp(-te ./ t2_est) + baseline_est; figure(1), hold on; plot(te, y_recon, 'r-'); |
and
Recent Comments