Sensitivity encoded MRI reconstruction¶
Nan-kuei Chen
2021-04-15
BME 639: Magnetic Resonance Imaging
The University of Arizona
The goal is to demonstrate the use of SENSE reconstruction to recover full-FOV images (i.e., free-from aliasing artifact) from under-sampled k-space data
The matlab version of this program is described here; and k-space data in matlab format is available in google drive
Section 1¶
We load two sets of fully-sampled k-space data: 4-shot spin-echo echo-planar imaging obtained at different echo times.
The matrix size for each echo image is \(128 \times 128 \times 8 \) (in-plane matrix size: 128x128; number of RF channels: 8)
using PyPlot
using FFTW
using Statistics
fid = open("data/kdata_4s_seepi","r");
kdata_4s_seepi = zeros(ComplexF64, 128,128,8,2);
read!(fid,kdata_4s_seepi);
close(fid);
@show size(kdata_4s_seepi);
size(kdata_4s_seepi) = (128, 128, 8, 2)
kdata_echo1 = kdata_4s_seepi[:,:,:,1]; # k-space data for the first echo
kdata_echo2 = kdata_4s_seepi[:,:,:,2]; # k-space data for the second echo
imgdata_echo1 = fftshift(fft(fftshift(kdata_echo1,[1,2]),[1,2]),[1,2]); # Image reconstruction with Fourier transform
imgdata_echo2 = fftshift(fft(fftshift(kdata_echo2,[1,2]),[1,2]),[1,2]);
# showing the magnitude and phase information of both echo images
function convertImage(input)
showImage = reverse(permutedims(input,[2,1,3]),dims=1);
showImage = reshape(showImage, size(showImage,1), size(showImage,2)*size(showImage,3));
end
figure(1,figsize=(9,3));imshow(abs.(convertImage(imgdata_echo1)), cmap="gray"); axis("off"); title("first echo: magnitude")
figure(2,figsize=(9,3));imshow(abs.(convertImage(imgdata_echo2)), cmap="gray"); axis("off") ; title("second echo: magnitude")
figure(3,figsize=(9,3));imshow(angle.(convertImage(imgdata_echo1)), cmap="gray"); axis("off"); title("first echo: phase")
figure(4,figsize=(9,3));imshow(angle.(convertImage(imgdata_echo2)), cmap="gray"); axis("off"); title("second echo: phase");
Section 2¶
Here we compute coil sensitivity profiles from the first echo images (fully-sampled in k-space). The information could be used in SENSE reconstruction of under-sampled k-space data (see sections 3 and 4)
We could combine the information from 8 RF channels with “root-mean-square”
function rmsCombineDataCoils(inputdata)
szl = ndims(inputdata)
if szl == 2
outputdata = deepcopy(inputdata);
end
if szl > 2
outputdata = sqrt.(mean(abs.(inputdata).^2,dims=szl));
end
return outputdata
end;
# Showing the RMS images from both echoes (fully-sampled in k-space)
combinedData_echo1 = rmsCombineDataCoils(imgdata_echo1);
combinedData_echo2 = rmsCombineDataCoils(imgdata_echo2);
figure(1,figsize=(2,2));imshow(convertImage(combinedData_echo1),cmap="gray"); axis("off"); title("first echo: RMS image ")
figure(2,figsize=(2,2));imshow(convertImage(combinedData_echo2),cmap="gray"); axis("off"); title("second echo: RMS image ");
Calculation of coil sensitivity profiles from the first echo image
coilSensitivityProfile1 = imgdata_echo1./combinedData_echo1;
figure(1,figsize=(9,3));imshow(abs.(convertImage(coilSensitivityProfile1)), cmap="gray"); axis("off"); title("magnitude of coil sensitivity profiles (measured from the first echo)")
figure(2,figsize=(9,3));imshow(angle.(convertImage(coilSensitivityProfile1)), cmap="gray"); axis("off"); title("phase information of coil sensitivity profiles (from the first echo)");
Section 3¶
Although we have fully-sampled k-space data for the second echo image, here we only use 50% of its k-space data (i.e., simulating under-sampled scan)
After Fourier transform, images reconstructed from under-sampled k-space show aliasing artifact
Assuming that we have acquied only 50% of the k-space data for the second echo image (with only odd ky lines acquired), what would the reconstructed images look like?
kData_acquired = kdata_echo2[:,1:2:end,:];
imgdata_acquired = fftshift(fft(fftshift(kData_acquired,[1,2]),[1,2]),[1,2]);
figure(1,figsize=(9,3));imshow(abs.(convertImage(imgdata_acquired)), cmap="gray"); axis("off"); title("second echo: magnitude of undersampled image")
figure(2,figsize=(9,3));imshow(angle.(convertImage(imgdata_acquired)), cmap="gray"); axis("off"); title("second echo: phase of undersampled image")
figure(3,figsize=(2,2));imshow(convertImage(rmsCombineDataCoils(imgdata_acquired)),cmap="gray"); axis("off"); title("second echo: RMS of undersampled data ");
Would zero-filling k-space data change anything?
kData_acquired_zeroFilled = zeros(ComplexF64,size(kData_acquired)[1],size(kData_acquired)[2]*2,size(kData_acquired)[3]);
kData_acquired_zeroFilled[:,1:2:end,:] = kData_acquired;
imgdata_acquired_zeroFilled = fftshift(fft(fftshift(kData_acquired_zeroFilled,[1,2]),[1,2]),[1,2]);
figure(1,figsize=(9,3));imshow(abs.(convertImage(kData_acquired_zeroFilled)).^0.3, cmap="gray"); axis("off"); title("undersampled and zero-filled k-space data")
figure(2,figsize=(2,2));imshow(convertImage(rmsCombineDataCoils(imgdata_acquired_zeroFilled)),cmap="gray"); axis("off"); title("Image-domain representation of undersampled and zero-filled k-space data");
figure(3,figsize=(9,3));imshow(abs.(convertImage(kData_acquired)).^0.3, cmap="gray"); axis("off"); title("undersampled k-space data (without zero-filling)")
figure(4,figsize=(2,2));imshow(convertImage(rmsCombineDataCoils(imgdata_acquired)),cmap="gray"); axis("off"); title("Image-domain representation of undersampled k-space data (without zero-filling)");
Section 4¶
Here we use coil sensitivity information estimated from the first echo image (fully-sampled: section 2) to remove aliasing artifact in under-sampled k-space data of the second echo image (section 3)
SENSE reconstruction for recovering full-FOV images from under-sampled data using the coil sensitivity profiles
sense_reconstructed_image = zeros(ComplexF64,size(kData_acquired)[1],size(kData_acquired)[2]*2);
for cntx = 1:size(kData_acquired)[1]
for cnty = 1: size(kData_acquired)[2]
mat1 = imgdata_acquired_zeroFilled[cntx,cnty,:];
mat2 = hcat(coilSensitivityProfile1[cntx,cnty,:],coilSensitivityProfile1[cntx,cnty+size(kData_acquired)[2],:])
mat3 = mat2\mat1;
sense_reconstructed_image[cntx,cnty]=mat3[1];
sense_reconstructed_image[cntx,cnty+size(kData_acquired)[2]]=mat3[2];
end
end
figure(1,figsize=(3,3));imshow(abs.(sense_reconstructed_image)',cmap="gray", origin="lower"); axis("off"); title("magnitude of SENSE produced image")
figure(2,figsize=(3,3));imshow(angle.(sense_reconstructed_image)',cmap="gray", origin="lower") ; axis("off"); title("phase of SENSE produced image")
figure(3,figsize=(3,3));imshow((abs.(sense_reconstructed_image)*2-combinedData_echo2[:,:,1])',cmap="gray", origin="lower"); axis("off"); title("Residual errors of SENSE produced image; compared with the ground truth");
Questions and exercise¶
Can we simply average information from 8 RF channels instead using room-mean-square?
What if we have only 4 RF channels? will the SENSE reconstruction work?
What if we have only 2 RF channels? will the SENSE reconstruction work?
Should we smooth coil sensitivity profiles? Will that improve the SNR of SENSE-reconstructed images?
How can we use regularization to reduce noise amplification?
Can we reconstruct full FOV images from only 25% of the k-space data (e.g., \(k_y\) line # 1, 5, 9 …)?