ESPIRiT Reconstruction Demo
This is a demo on how to generate ESPIRiT maps and use them to perform ESPIRiT reconstruction for parallel imaging. It is based on the paper Uecker et. al, MRM 2013 DOI 10.1002/mrm.24751. ESPIRiT is a method that finds the subspace of multi-coil data from a calibration region in k-space using a series of eigen-value decompositions in k-space and image space. Here we also use the "soft" sense idea (Uecker et. al, "ESPIRiT Reconstruction using Soft-SENSE", Proceedings of the ISMRM 2013, pp-127) by using the eigen values to weight the eigen-vectors.
Here we perform ESPIRiT calibration on data which has strong aliasing in the phase-encode direction. SENSE often fails with this type of data.
load brain_alias_8ch DATA = DATA/max(max(max(abs(ifft2c(DATA))))) + eps; %DATA = crop(DATA,[256,256,8]); ncalib = 24; % use 24 calibration lines to compute compression ksize = [6,6]; % ESPIRiT kernel-window-size eigThresh_k = 0.02 % threshold of eigenvectors in k-space eigThresh_im = 0.9; % threshold of eigenvectors in image space [sx,sy,Nc] = size(DATA); % create a sampling mask to simulate x2 undersampling with autocalibration % lines mask = zpad(ones(sx,ncalib),[sx,sy]); mask = repmat(mask,[1,1,8]); mask(:,1:2:end,:) = 1; DATAc = DATA.*mask; calib = crop(DATAc,[sx,ncalib,Nc]);
eigThresh_k = 0.0200
Display coil images:
im = ifft2c(DATAc); figure, imshow3(abs(im),,[1,Nc]); title('magnitude of physical coil images'); colormap((gray(256))); colorbar; figure, imshow3(angle(im),,[1,Nc]); title('phase of physical coil images'); colormap('default'); colorbar;
Maps are computed in two steps.
% compute Calibration matrix, perform 1st SVD and convert singular vectors % into k-space kernels [k,S] = dat2Kernel(calib,ksize); idx = max(find(S >= S(1)*eigThresh_k));
Display the singular vectors and values of the calibration matrix
kdisp = reshape(k,[ksize(1)*ksize(2)*Nc,ksize(1)*ksize(2)*Nc]); figure, subplot(211), plot([1:ksize(1)*ksize(2)*Nc],S,'LineWidth',2); hold on, plot([1:ksize(1)*ksize(2)*Nc],S(1)*eigThresh_1,'r-','LineWidth',2); plot([idx,idx],[0,S(1)],'g--','LineWidth',2) legend('signular vector value','threshold') title('Singular Vectors') subplot(212), imagesc(abs(kdisp)), colormap(gray(256)); xlabel('Singular value #'); title('Singular vectors')
crop kernels and compute eigen-value decomposition in image space to get maps
[M,W] = kernelEig(k(:,:,:,1:idx),[sx,sy]);
show eigen-values and eigen-vectors. The last set of eigen-vectors corresponding to eigen-values 1 look like sensitivity maps
figure, imshow3(abs(W),,[1,Nc]); title('Eigen Values in Image space'); colormap((gray(256))); colorbar; figure, imshow3(abs(M),,[Nc,Nc]); title('Magnitude of Eigen Vectors'); colormap(gray(256)); colorbar; figure, imshow3(angle(M),,[Nc,Nc]); title('Magnitude of Eigen Vectors'); colormap(jet(256)); colorbar;
Warning: Image is too big to fit on screen; displaying at 33% Warning: Image is too big to fit on screen; displaying at 33%
crop sensitivity maps according to eigenvalues==1. Note that we have to use 2 sets of maps. Here we weight the 2 maps with the eigen-values
maps = M(:,:,:,end-1:end); % Weight the eigenvectors with soft-senses eigen-values weights = W(:,:,end-1:end) ; weights = (weights - eigThresh_im)./(1-eigThresh_im).* (W(:,:,end-1:end) > eigThresh_im); weights = -cos(pi*weights)/2 + 1/2; % create and ESPIRiT operator ESP = ESPIRiT(maps,weights); nIterCG = 12;
ESPIRiT CG reconstruction with soft-sense and 2 sets of maps
disp('Performing ESPIRiT reconstruction from 2 maps') tic; [reskESPIRiT, resESPIRiT] = cgESPIRiT(DATAc,ESP, nIterCG, 0.01,DATAc*0); toc % ESPIRiT CG reconstruction with 1 map disp('Performing SENSE reconstruction from 1 set of maps') SNS = ESPIRiT(maps(:,:,:,end)); tic;[reskSENSE, resSENSE] = cgESPIRiT(DATAc, SNS, nIterCG,0.01 ,DATAc*0);toc % GRAPPA reconstruction disp('Performing GRAPPA reconstruction ') tic; reskGRAPPA = GRAPPA(DATAc,calib,[5,5],0.01);toc resGRAPPA = ifft2c(reskGRAPPA);
Performing ESPIRiT reconstruction from 2 maps Elapsed time is 5.465768 seconds. Performing SENSE reconstruction from 1 set of maps Elapsed time is 4.023993 seconds. Performing GRAPPA reconstruction reconstructing coil 1 reconstructing coil 2 reconstructing coil 3 reconstructing coil 4 reconstructing coil 5 reconstructing coil 6 reconstructing coil 7 reconstructing coil 8 Elapsed time is 12.164346 seconds.
Note the typical center FOV aliasing in SENSE. Also, note that ESPIRiT has (very slightly) less error than GRAPPA
figure, imshow(cat(2,abs(resSENSE),sos(resESPIRiT), sos(resGRAPPA)),) title('SENSE reconstruction with 1 map vs ESPIRiT with 2 vs GRAPPA') figure, imshow(cat(2,sos(ifft2c(reskSENSE-DATA)),sos(ifft2c(reskESPIRiT-DATA)),sos(ifft2c(reskGRAPPA-DATA))).^(1/2),) title('SENSE reconstruction error with 1 map vs ESPIRiT with 2 vs GRAPPA')