T2 Curve Fitting – Regularized NNLS
This will be similar to the other solutions. First we will create a noisy decay curve based on some reasonable parameters. The code is runme_nnls_regularized.
I haven’t checked this over very well, so there might be some subtle mistakes.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | %%======================================================== %% %% Create the decay curve with noise. %% %%======================================================== %% Create the TE and T2 vectors te = [10:10:320]'; t2 = logspace( log10(10), log10(3000), 80); %% Create the A matrix A = exp(- kron(te, 1./t2) ); A = [A ones(nrows(A), 1)]; %% Create a dummy data vector y = 1000* ( 0.4*exp(-te / 34) + 0.6*exp(-te / 111) ); %% Add some noise. noise_std = 10; e = noise_std*randn(32, 1); |
Now we’ll solve it using a non-regularized (simple) NNLS solution.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | %%======================================================== %% %% Solve using non-regularized NNLS %% %%======================================================== %% Solve using NNLS y_e = y + e + 20; x2 = lsqnonneg(A, y_e); y_calculated = A * x2; baseline = x2(end); x2 = x2(1:end-1); x2_mu0 = x2; inds = find(x2 > 0 ); baseline [x2(inds) t2(inds)'] chi2_min = sum( (y_calculated - y_e).^2 ./ noise_std.^2 ) |
Now we’ll display the solution.
1 2 3 4 5 6 7 8 9 10 11 |
Now we’ll compute the regularized version. The basic thing here is that we are going to introduce the H matrix that will be used for regularizing. It is just the identity matrix which we’ll pre-multiply by a factor mu. Pre-multiplying by different mu will smooth out the resulting T2 spectrum and thereby increadse the χ2 a little.
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 | %%======================================================== %% %% Now do the regularized version. %% %%======================================================== %% %% Now do a regularized version %% mu = [0 ]; chi2 = [chi2_min ]; H = eye(size(A)); for ii=1:nrows(H)-1, H(ii,ii+1) = -1; end; while( ~isempty(chi2) & chi2(end) < chi2_min * 1.02 ) if( mu(end) > 0 ) mu = [ mu 2*mu(end) ]; else mu = [ mu 0.001]; end mu(end) AH = [A; mu(end)*H]; x2 = lsqnonneg(AH, [y_e; zeros(length(te), 1)]); y_calculated = A * x2; chi2 = [ chi2 sum( (y_calculated - y_e).^2 ./ noise_std.^2 )]; end |
We can plot the relationship between mu and χ2.
1 2 3 4 5 6 7 8 9 10 11 12 13 |
And then we want to find a mu which will regularize the solution such that the χ2 is within the range χ2_min*1.02 ≤ χ2 ≤ χ2_min * 1.025.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | %%======================================================== %% %% Find the mu which corresponds to a chi2 in the range %% in which we are interested. %% %%======================================================== %% Find the mu that corresponds to 1.02 * chi2_min mu_spline = 0:0.001:mu(end); chi2_spline = interp1(mu, chi2, mu_spline, 'spline'); [mm,mi] = min( abs((chi2_spline - (1.02 * chi2_min))) ); fprintf('chi2_min = %5.3f 1.02*chi2_min = %5.3f found value = %5.3f\n', chi2_min, 1.02*chi2_min, chi2_spline(mi)); h=line([mu(1) mu(end)], 1.02*[chi2_min chi2_min]); set(h, 'color', 'red'); |
And now we can plot the fitted decay curve and the spectrum.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | %%======================================================== %% %% Calculate and plot the regularized solution. %% %%======================================================== AH = [A; mu_spline(mi)*eye(size(A))]; x2 = lsqnonneg(AH, [y_e; zeros(length(te), 1)]); y_calculated = A * x2; chi2_102 = sum( (y_calculated - y_e).^2 ./ noise_std.^2 ) figure(1); hold on; plot(te, y_calculated, 'r-'); figure(3),clf; hold on; plot(t2, x2(1:end-1)); plot(t2, x2_mu0, 'r-'); set(gca, 'xscale', 'log'); xlabel('T_2 (ms)'); title('T_2 Spectrum'); grid on; legend('Non Regularized', 'Regularized'); mu_spline(mi) |
There is still a little more documentation to do on this page. That will have to wait for now.
References
- C. L. Lawson and R. J. Hanson. Solving least square problems. Prentice Hall, Englewood Cliffs NJ, 1974
This is the book in which the algorithm is originally described. Theorems to show NNLS will stop in a finite num ber of steps and will arrive at the minimum L2 solution. - K.P. Whittall and A.L. MacKay. Quantitative interpretation of NMR relaxation data. J. Magn. Reson., 84:134-152, 1989.
One of the first papers in MRI to use the NNLS algorithm. Good discussion of regularization techniques as well.
Recent Comments