Preconditioned ADMM with nonlinear operator constraint
Martin Benning, Florian Knoll, Carola-Bibiane Schönlieb & Tuomo Valkonen
This is the corresponding code to [1]. In order to run this code it is necessary to either execute
setpath
first, or to add the folders manually to the path.
Author: Martin Benning - mb941@cam.ac.uk
Last updated: 07.06.2016
Load coil sensitivity maps and brain phantom
We begin by loading an artificial brain phantom ugt and measured coil sensitivities b. The coil sensitivities were generated from the measurements of a water bottle with an 8-channel head coil array. The numerical phantom is based on [2]; details can be found in [1, Section 6.1]. Courtesy of Florian Knoll, NYU.
load('Data/B1Map_190.mat')
load('Data/digital_brain_phantom_images_2d.mat')
b = b1;
ugt = u_mri_flair/max(abs(u_mri_flair(:)));
nrcoils = size(b, 3);
For comparison to later reconstructions, we visualise the modulus of the ground truth spin proton density and the coil sensitivities.
figure()
imagesc(abs(ugt), [0 1])
axis image
colormap(gray(512))
colorbar
drawnow
figure()
for i=1:nrcoils
subplot(2, 4, i), imagesc(abs(b(:, :, i)), [0 max(abs(b1(:)))]);
axis image
colorbar
colormap(gray(512))
end
drawnow
Initialise operators to create artificial k-space-data
In order to create artificial k-space data from ugt and b, we need to intialise the forward operators G and F first. Here G will realise equation (11) as described in [1], whereas F performs a two-dimensional discrete Fourier transform with circular shift from the coordinates [1, 1] to [95, 95].
G = coilsensop(size(ugt), nrcoils);
F = mridft(size(ugt), [95 95]);
We load a synthetic spiral sampling scheme samp that is sampled on a cartesisan grid.
load('Supplementary Files/spiralsampling.mat')
We further create a sub-sampling matrix S from samp and multiply it with F to create the sub-sampled Fourier operator SF.
S = speye(size(F, 1));
S = S(samp(:), :);
SF = S*F;
Create artificial k-space-data
Now we can create artificial k-space data by applying G and SF to ugt and b. First, we group ugt and b to one auxiliary variable aux1 and apply G to produce another auxiliary variable aux2.
aux1 = cat(3, ugt, b);
aux2 = G.feval(aux1(:));
Then we produce artificial, noisefree k-space data g for each coil by applying SF to each component of aux2.
g = zeros(nrcoils*size(SF, 1), 1);
for i=1:nrcoils
g(((i - 1)*size(SF, 1) + 1):(i*size(SF, 1))) = SF*aux2(((i - ...
1)*size(F, 2) + 1):(i*size(F, 2)));
end
Subsequently, we contaminate the noisefree k-space data g by adding Gaußian noise with standard deviation 
and obtain our artificial k-space data f.
sigma = 0.05;
f = g + sigma*randn(size(g));
Note that we reproduce the low noise-level results from [1] with this choice of sigma. In order to reproduce the high noise-level results, we need to set 
here.
Compute zero-filled reconstructions
For comparison purposes, we also compute a naive zero-filling reconstruction. We apply the Hermitian transpose of SF to f and store the mean as a new variable uzfmean.
uzf = zeros([size(G, 1) 1]);
for i=1:nrcoils
uzf(((i - 1)*size(SF, 2) + 1):(i*size(SF, 2))) = SF'*f(((i - ...
1)*size(SF, 1) + 1):(i*size(SF, 1)));
end
uzf = reshape(uzf, sizend(G, 1));
uzfmean = reshape(mean(uzf, 3), [size(F, 2) 1]);
We plot the computed zero-filled reconstructions to emphasise that significant artifacts are introduced when the modelling of the coil sensitivities is neglected.
figure()
for i=1:nrcoils
subplot(2, 4, i), imagesc(abs(uzf(:, :, i)), [0 max(abs(b1(:)))]);
axis image
colorbar
colormap(gray(512))
end
drawnow
Set up gradient operator for reconstruction method
For regularisation purposes we need to initialise a two-dimensional finite difference approximation of the gradient operator Grad with stepsize one in both dimensions.
Grad = buildgradient(sizend(F, 2), ones(2, 1));
The gradient has to be applied to the image and all coils, see the equation following after Equation (11) in [1]. Together with the operator G it forms the operator 
required in Algorithm 1. As we have used F for the notation of the Fourier transform already, we denote this non-linear operator as A.
Gradall = kron(speye(nrcoils + 1), Grad);
A = [G; Gradall];
Set up proximity operators for reconstruction method
The next requirement for the application of Algorithm 1 is the initialisation of the proximity operators corresponding to 
and 
. As these consist of multiple partial proximity operators, we also initialise them that way. As described in the paper, we set 
. As we have 
in our case, we initialise 9 trivial proximity operators of correct size as follows.
H0 = proxoperator(sizend(F, 2));
H1 = proxoperator(sizend(F, 2));
H2 = proxoperator(sizend(F, 2));
H3 = proxoperator(sizend(F, 2));
H4 = proxoperator(sizend(F, 2));
H5 = proxoperator(sizend(F, 2));
H6 = proxoperator(sizend(F, 2));
H7 = proxoperator(sizend(F, 2));
H8 = proxoperator(sizend(F, 2));
We concatenate them to one proxmity operator for . Note that in future releases operation like H(1:9) = proxoperator( ... ) are planned to be supported.
H = [H0, H1, H2, H3, H4, H5, H6, H7, H8];
We now do the same for the proximity operator of . According to the paper, we have 
with 
for 
, 
and 
. We initialise the proximity operators for first.
J0 = diagorthprox(S, F);
J1 = diagorthprox(S, F);
J2 = diagorthprox(S, F);
J3 = diagorthprox(S, F);
J4 = diagorthprox(S, F);
J5 = diagorthprox(S, F);
J6 = diagorthprox(S, F);
J7 = diagorthprox(S, F);
Next, we initialise 
.
J8 = onetwonormprox([size(F, 2) 2]);
And conclude with the proximity operators for 
.
J9 = twonormprox([size(F, 2)*2 1]);
J10 = twonormprox([size(F, 2)*2 1]);
J11 = twonormprox([size(F, 2)*2 1]);
J12 = twonormprox([size(F, 2)*2 1]);
J13 = twonormprox([size(F, 2)*2 1]);
J14 = twonormprox([size(F, 2)*2 1]);
J15 = twonormprox([size(F, 2)*2 1]);
J16 = twonormprox([size(F, 2)*2 1]);
As for , we concatenate all the individual proximity operators to obtain one proximity operator for .
J = [J0, J1, J2, J3, J4, J5, J6, J7, J8, J9, J10, J11, J12, J13, J14, ...
J15, J16];
Associate data to prox operators
Before initialising Algorithm 1, we equip the proximity terms that correspond to the data fidelities with the noisy k-space data f.
J0.setdata(f(1:size(SF, 1)))
J1.setdata(f((size(SF, 1) + 1):(2*size(SF, 1))))
J2.setdata(f((2*size(SF, 1) + 1):(3*size(SF, 1))))
J3.setdata(f((3*size(SF, 1) + 1):(4*size(SF, 1))))
J4.setdata(f((4*size(SF, 1) + 1):(5*size(SF, 1))))
J5.setdata(f((5*size(SF, 1) + 1):(6*size(SF, 1))))
J6.setdata(f((6*size(SF, 1) + 1):(7*size(SF, 1))))
J7.setdata(f((7*size(SF, 1) + 1):(8*size(SF, 1))))
Initialise simplified version Algorithm 1
As we now have all the ingredients necessary, we proceed to initialise an object solver of the simplified version of Algorithm 1, where (respectively A in terms of the variable name) is only non-linear in 
but not in 
.
solver = NPADMMsimple(A, H, J);
We set the parameters 
and 
fixed to 1/3, which experimentally appears to guarantee convergence - otherwise we need to compute the operator norm of the linearised A in each iteration. In addition, we fix the (maximum) number of iterations to 1500 as described in [1].
solver.param.tau = 1/3;
solver.param.delta = 1/3;
solver.param.maxiter = 1500;
Set initial values
As described in [1, Section 6.2], we initialise 
with 
being the constant one-vector, for all 
. All other initial variables (
, 
, 
) are set to zero.
solver.reinitialise;
solver.var.u = ones(9*size(F, 2), 1);
Set regularisation parameters
We fix the regularisation parameters to 
, 
and 
for 
, in order to reproduce the low noise-level results of [1]. In order to reproduce the high noise-level results, the parameters have to be adapted according to [1, Section 6.2].
lambda = 0.0621;
alpha0 = 0.0062;
alphaj = 0.9317;
Finally, we modifiy the regularisation parameters of the objects of the proximity operators.
J0.setalpha(lambda)
J1.setalpha(lambda)
J2.setalpha(lambda)
J3.setalpha(lambda)
J4.setalpha(lambda)
J5.setalpha(lambda)
J6.setalpha(lambda)
J7.setalpha(lambda)
J8.setalpha(alpha0)
J9.setalpha(alphaj)
J10.setalpha(alphaj)
J11.setalpha(alphaj)
J12.setalpha(alphaj)
J13.setalpha(alphaj)
J14.setalpha(alphaj)
J15.setalpha(alphaj)
J16.setalpha(alphaj)
Solve nonlinear coilsensitivity estimation problem
Now we run Algorithm 1 by executing
solver.solve;
We extract the reconstructed variables via
recon = solver.var.u;
and store the spin proton density as urec via
urec = reshape(recon(1:size(F, 2)), sizend(F, 2));
and the coil sensitivities as coilsens via
coilsens = zeros([sizend(F, 2) nrcoils]);
for i=1:nrcoils
coilsens(:, :, i) = reshape(recon((i*size(F, 2) + 1):(i + 1)*size(F,...
2)), sizend(F, 2));
end
Plot results
To conclude, we visualise the results.
figure()
imagesc(abs(urec), [0 1])
axis image
colorbar
colormap(gray(512))
drawnow
figure()
for i=1:nrcoils
subplot(2, 4, i), imagesc(abs(coilsens(:, :, i)), [0 max(abs(b1(:)))]);
axis image
colorbar
colormap(gray(512))
drawnow
end
References
[1] Benning, Martin, Florian Knoll, Carola-Bibiane Schönlieb, and Tuomo Valkonen. "Preconditioned ADMM with nonlinear operator constraint." Accepted for publication in IFIP TC 7 2015 conference proceedings (2016).
Preprint: arXiv preprint arXiv:1511.00425 (2015)
[2] Aubert-Broche, Berengere, Alan C. Evans, and Louis Collins. "A new improved version of the realistic digital brain phantom." NeuroImage 32, no. 1 (2006): 138-145.