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: $latex SI= A\rho\exp\left( \frac{-TE}{T_2} \right) $

where $latex \rho$ is the proton density, is arbitrary scale factor incorporating many scanner dependent things, TE are the set of echo times and $latex T_2$ is the spin-spin relaxation time. For most things it is not possible to measure $latex A$ and $latex \rho$ independently.

Mono-Exponential Data Fitting

There are many ways to fit mono-exponential data to obtain $latex T_2$.

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.

[cc lang="matlab"]
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;
[/cc]

And now display it to make sure we understand how things should look.
[cc lang="matlab"]
figure(1),clf;
(te, y_e, ‘x’);
title(‘The test data’);
xlabel(‘TE (ms)’);
ylaabel(‘Signal’);
[/cc]

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.
[cc lang="matlab"]
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)
[/cc]

Now we can reconstruct the data based on the best guess and display it overtop the initial data.
[cc lang="matlab"]
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-’);
[/cc]

The two files to run this experiment are:

[cc lang="matlab"]
% 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-’);
[/cc]

and

[cc lang="matlab"]
function [d] = fun(x, S, te)

A = x(1);
t2 = x(2);
baseline = x(3);

S_est = A * exp(-te ./ t2) + baseline;

d = sqrt( sum( (S_est – S).^2 ) );
[/cc]