【图像重建】基于SIDER算法实现图像的压缩重建附matlab代码

✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信。

🍎个人主页:Matlab科研工作室

🍊个人信条:格物致知。

更多Matlab仿真内容点击👇

智能优化算法  神经网络预测 雷达通信  无线传感器

信号处理 图像处理 路径规划 元胞自动机 无人机

⛄ 内容介绍

Diffusion MRI measurements using hyperpolarized gases are generally acquired during patient breath hold, which yields a compromise between achievable image resolution, lung coverage and number of b-values. In this work, we propose a novel method that accelerates the acquisition of MR diffusion data by undersampling in both spatial and b-value dimensions, thanks to incorporating knowledge about the signal decay into the reconstruction (SIDER). SIDER is compared to total variation (TV) reconstruction by assessing their effect on both the recovery of ventilation images and estimated mean alveolar dimensions (MAD). Both methods are assessed by retrospectively undersampling diffusion datasets of normal volunteers and COPD patients (n=8) for acceleration factors between x2 and x10. TV led to large errors and artefacts for acceleration factors equal or larger than x5. SIDER improved TV, presenting lower errors and histograms of MAD closer to those obtained from fully sampled data for accelerations factors up to x10. SIDER preserved image quality at all acceleration factors but images were slightly smoothed and some details were lost at x10. In conclusion, we have developed and validated a novel compressed sensing method for lung MRI imaging and achieved high acceleration factors, which can be used to increase the amount of data acquired during a breath-hold. This methodology is expected to improve the accuracy of estimated lung microstructure dimensions and widen the possibilities of studying lung diseases with MRI.

⛄ 部分代码

function Demo_SIDER

% function Demo_SIDER.m

%

% The proposed method incorporates the knowledge of the signal decay into

% the reconstruction (SIDER) to accelerate the acquisition of MR diffusion

% data by undersampling in both spatial and b-value dimensions.

%

% Demo for reconstructing diffusion lung MR images with spatial total

% variation (S-TV) and spatiotemporal total variation (ST-TV) using the

% Split Bregman formulation. SIDER combines total variation (TV) with a

% penalty function that promotes sparsity across the b-direction as follows

% min_u beta|grad_x,y u|_1 + gamma|M u|_1 st. ||Fu-f||^2 < delta, where the

% first term corresponds to spatial TV and M is is an operator that encodes

% the relationship between ventilation images for consecutives values of b,

% based on a stretched exponential model.

%

% The demo loads fully sampled data, undersampled data and reconstruct

% images using both TV and SIDER methods.

% Data is undersampled using a modified version of Lustig's variable

% density pdf, downloaded from (SparseMRI V0.2)

% http://web.stanford.edu/~mlustig/SparseMRI.html, see M. Lustig, D.L

% Donoho and J.M Pauly "Sparse MRI: The Application of Compressed Sensing

% for Rapid MR Imaging" Magnetic Resonance in Medicine, 2007 Dec;

% 58(6):1182-1195.

%

% -------------------------------------------------------------------------

% Load complex image

load('DataControl','b','uTarget');

% Choose one slice

zSlice      = 4;

uTarget     = uTarget(:,:,:,zSlice); 

b           = b(1:5);

N           = size(uTarget);

[nel_Lm,xcenters_Lm]    = hist(150:20:450,30);

% -------------------------------------------------------------------------

% MASK FOR THE LUNG to define the region where to estimate D, alpha and Lm

% maps (this estimation is a pixel-by-pixel fitting so its calculation is

% not affected by the mask). The mask is automatically computed by

% thresholding the ventilation image

se      = strel('disk',3);

im      = imclose(abs(uTarget(:,:,1)),se);   % imdilate

[val, ind] = sort(im(:),'descend');

im      = im/val(5);

mask    = im>=0.3;                  % thresholding

se      = strel('disk',1);          % remove small objects

mask    = (imopen(abs(mask),se))>0;

h=figure;

subplot(2,1,1);

imagesc(abs((uTarget(:,:,1))));colorbar; colormap hot; axis image;

title('b=0');

subplot(2,1,2);

imagesc(abs((mask(:,:,1))));colorbar; colormap hot; axis image;

title('Mask');

pause(2);

close(h);

% -------------------------------------------------------------------------

% ESTIMATE D, alpha and Lm MAPS

%

% Ventilation images are filtered before the estimation (as images are very

% noisy for patient data sets, especially for large values of b)

k       = fspecial('gaussian',3,1); 

h=figure;

for ip = 1:N(3)

    uTargetSmooth(:,:,ip) = imfilter(uTarget(:,:,ip),k);

    subplot(2,1,1); imagesc(abs(uTarget(:,:,ip)));colorbar; colormap hot; axis image

    title(['Original ventilation image, b ' num2str(ip)]);

    subplot(2,1,2); imagesc(abs(uTargetSmooth(:,:,ip)));colorbar; colormap hot; axis image;

    title('Filtered ventilation image');

    pause(1);

end

close(h);

% Estimate D and alpha using a STRETCHED EXPONENTIAL MODEL

options     = optimset('display','iter','Algorithm','levenberg-marquardt','Jacobian','on');

x0_D        = 0.2*ones(N(1:2));

x0_D        = x0_D(mask);

x0_a        = 0.8*ones(N(1:2));

x0_a        = x0_a(mask);

x0          = [x0_D; x0_a];

x           = lsqnonlin(@objVecPotExpJac,x0,[],[],options,...

    abs(uTargetSmooth(:,:,1)),abs(uTargetSmooth(:,:,2:end)),b(2:end),mask);

Drec0        = zeros(N(1:2));

Drec0(mask)  = abs(x(1:nnz(mask)));

alpharec0    = zeros(N(1:2));

alpharec0(mask)= abs(x(nnz(mask)+1:2*nnz(mask)));

% Apply a threshold to the estimated values

Drec0(abs(Drec0)>0.9)=0; Drec0(abs(Drec0)<0.05)=0;

alpharec0(abs(alpharec0)>1.3)=0; alpharec0(abs(alpharec0)<0.3)=0;

% ESTIMATE Lm map

tmp         = alpharec0;

tmp(logical(mask.*alpharec0>1)) = 0.9; % It doesnt work for values of alpha larger than 1

Lm0         = ComputeLmImage(N,Drec0,tmp);

h=figure;

subplot(2,2,1); imagesc(abs(Drec0)); colormap hot; axis image; colorbar; title('D fully sampled');caxis([0 0.9]);

subplot(2,2,2); imagesc(abs(alpharec0)); colormap hot; axis image; colorbar; title('alpha');caxis([0 1.3]);

subplot(2,2,3); imagesc(abs(Lm0)); colormap hot; axis image; colorbar; title('Lm');caxis([150 450]);

pause(2);

close(h);

% -------------------------------------------------------------------------

% -------------------------------------------------------------------------

% COMPRESSED SENSING

%

% Pseudo-random undersampling pattern for all values of b

rand('state',1);

% Select parameters for Lustig's variable density pdf: radius, sparsity and

% P. Below we provide the set of paired sparsity and P values used in the

% paper:

% sparsity [0.5 0.25 0.20 0.15 0.1 0.07]

% P        [4   4    4    6    9   13]

radius      = 0.017;

DN          = [N(1),1];

sparsity    = 0.15;      % Corresponds to acceleration factor x7

P           = 6;

RAll    = zeros(N);

h=figure;

for ip = 1:N(3)

    % Loop across of values of b

    pdf     = genPDF(DN,P,sparsity,2,radius,0);

    temp    = genSampling_LIM(pdf,20,1);

    indR    = temp(:,1)==1;

    R       = zeros(N(1:2));

    R(indR,:)= 1;

    R(1,:)= 1;

    spy(R); title(['Undersampling b ' num2str(ip)]);

    RAll(:,:,ip) = R;

    data    = fft2(uTarget(:,:,ip)).*R;

    dataAll(:,:,ip) = data;

    pause(1);

end % ip

close(h);

% -------------------------------------------------------------------------

% ZERO-FILLING RECONSTRUCTION

u_fft       = ifft2(dataAll);

% Smooth ventilation images

for ip = 1:N(3)

    u_fft_s(:,:,ip) = imfilter(u_fft(:,:,ip),k);

end

% POTENTIAL MODEL

options     = optimset('display','iter','Algorithm','levenberg-marquardt','Jacobian','on');

x0_D        = 0.2*ones(N(1:2));

x0_D        = x0_D(mask);

x0_a        = 0.8*ones(N(1:2));

x0_a        = x0_a(mask);

x0          = [x0_D; x0_a];

x           = lsqnonlin(@objVecPotExpJac,x0,[],[],options,abs(u_fft_s(:,:,1)),abs(u_fft_s(:,:,2:end)),b(2:end),mask);

% D and A

Drecfft      = zeros(N(1:2));

Drecfft(mask)= abs(x(1:nnz(mask)));

alpharecfft  = zeros(N(1:2));

alpharecfft(mask)= abs(x(nnz(mask)+1:2*nnz(mask)));

Drecfft(abs(Drecfft)>0.9)=0; Drecfft(abs(Drecfft)<0.05)=0;

alpharecfft(abs(alpharecfft)>1.3)=0; alpharecfft(abs(alpharecfft)<0.3)=0;

% Lm

tmp         = alpharecfft;

tmp(logical(mask.*alpharecfft>1)) = 0.9;%!!!!1 Othewise, Lms doesnt work!!!

Lmfft       = ComputeLmImage(N,Drecfft,tmp);

Lmfft(isnan(Lmfft)) = 0;

h=figure;

subplot(2,2,1); imagesc(abs(Drecfft)); colormap hot; axis image; colorbar; title('D IFFT');caxis([0 0.9]);

subplot(2,2,2); imagesc(abs(alpharecfft)); colormap hot; axis image; colorbar; title('alpha');caxis([0 1.3]);

subplot(2,2,3); imagesc(abs(Lmfft)); colormap hot; axis image; colorbar; title('Lm');caxis([150 450]);

pause(2);

close(h);

% -------------------------------------------------------------------------

% STATIC Spatial Total Variation reconstruction using Split Bregman

% Code download from

% http://www.ece.rice.edu/~tag7/Tom_Goldstein/Split_Bregman.html

%

% Goldstein's spatial TV using the Split Bregman formulation

% u_tv = mrics(RAll(:,:,1),data(:,:,1), mu, lambda, gamma, nInner, nBreg);

%

% SpatialTVSB.m: same as mrics.m but it computes the solution error

% norm

mu          = 1;

lambda      = 1;

gamma       = 0.01;

nInner      = 1;

nBreg       = 800;

% Loop across all values of b to reconstruct ventilation images

for ip = 1:N(3)

    % The figure displayed shows (at several iteration numbers) the

    % reconstructed image (u), the dummy variable x that outputs the

    % shrinkage (for the gradient along x-direction) indicating the

    % sparsity as the number of nonzero coefficients, and the solution

    % error norm for all ventilation images

    [uThis,errThis] = SpatialTVSB(RAll(:,:,ip),dataAll(:,:,ip), mu, lambda, gamma, nInner, nBreg,uTargetSmooth(:,:,ip),mask);

    u_tv(:,:,ip) = uThis;

    err_tv(:,ip) = errThis;

end % ip

% Smooth ventilation images

for ip = 1:N(3)

    u_tv_s(:,:,ip) = imfilter(u_tv(:,:,ip),k);

end

% POTENTIAL MODEL

options     = optimset('display','iter','Algorithm','levenberg-marquardt','Jacobian','on');

x0_D        = 0.2*ones(N(1:2));

x0_D        = x0_D(mask);

x0_a        = 0.8*ones(N(1:2));

x0_a        = x0_a(mask);

x0          = [x0_D; x0_a];

x           = lsqnonlin(@objVecPotExpJac,x0,[],[],options,abs(u_tv_s(:,:,1)),abs(u_tv_s(:,:,2:end)),b(2:end),mask);

% D and A

Drectv      = zeros(N(1:2));

Drectv(mask)= abs(x(1:nnz(mask)));

alpharectv  = zeros(N(1:2));

alpharectv(mask)= abs(x(nnz(mask)+1:2*nnz(mask)));

Drectv(abs(Drectv)>0.9)=0; Drectv(abs(Drectv)<0.05)=0;

alpharectv(abs(alpharectv)>1.3)=0; alpharectv(abs(alpharectv)<0.3)=0;

% Lm

tmp         = alpharectv;

tmp(logical(mask.*alpharectv>1)) = 0.9;

Lmtv        = ComputeLmImage(N,Drectv,tmp);

Lmtv(isnan(Lmtv)) = 0;

h=figure;

subplot(2,2,1); imagesc(abs(Drectv)); colormap hot; axis image; colorbar; title('D TV');caxis([0 0.9]);

subplot(2,2,2); imagesc(abs(alpharectv)); colormap hot; axis image; colorbar; title('alpha');caxis([0 1.3]);

subplot(2,2,3); imagesc(abs(Lmtv)); colormap hot; axis image; colorbar; title('Lm');caxis([150 450]);

pause(2);

close(h);

% -------------------------------------------------------------------------

% SIDER

%

% Get mean estimated values of D and alpha for the model

Dest        = mean(Drectv(Drectv>0));

Alphaest    = mean(alpharectv(alpharectv>0));

mu      = 1;

lambda  = 1;

gamma   = 1e-2;

nBreg   = 800;

alpha   = 0.2;

beta    = 0.2;

% The figure displayed shows (at several iteration numbers) the

% reconstructed image (u), the dummy variable x that outputs the

% shrinkage (for the gradient along x-direction) and the dummy variable p

% (for the operator along b-direction), and the solution error norm for all

% ventilation images

[u_t,err_t] = SIDER(Dest,Alphaest,b,dataAll,RAll,N(1:3),mu,lambda,gamma,alpha,beta,nBreg,mask,uTargetSmooth);

% Smooth ventilation images

for ip = 1:N(3)

    u_t_s(:,:,ip) = imfilter(u_t(:,:,ip),k);

end

% POTENTIAL MODEL

options     = optimset('display','iter','Algorithm','levenberg-marquardt','Jacobian','on','TolFun',1e-06);

x0_D        = Dest*ones(N(1:2));

x0_D        = x0_D(mask);

x0_a        = 0.8*ones(N(1:2));

x0_a        = x0_a(mask);

x0          = [x0_D; x0_a];

x           = lsqnonlin(@objVecPotExpJac,x0,[],[],options,abs(u_t_s(:,:,1)),abs(u_t_s(:,:,2:end)),b(2:end),mask);

% D and A

Drect       = zeros(N(1:2));

Drect(mask)= abs(x(1:nnz(mask)));

alpharect   = zeros(N(1:2));

alpharect(mask)= abs(x(nnz(mask)+1:2*nnz(mask)));

Drect(abs(Drect)>0.9)=0; Drect(abs(Drect)<0.05)=0;

alpharect(abs(alpharect)>1.3)=0; alpharect(abs(alpharect)<0.3)=0;

% Lm

tmp         = alpharect;

tmp(logical(mask.*alpharect>1)) = 0.9;

Lmt         = ComputeLmImage(N,Drect,tmp);

Lmt(isnan(Lmt)) = 0;

h=figure;

subplot(2,2,1); imagesc(abs(Drect)); colormap hot; axis image; colorbar; title('D SIDER'); caxis([0 0.9]);

subplot(2,2,2); imagesc(abs(alpharect)); colormap hot; axis image; colorbar; title('alpha'); caxis([0 1.3]);

subplot(2,2,3); imagesc(abs(Lmt)); colormap hot; axis image; colorbar; title('Lm'); caxis([150 450]);

pause(2);

close(h);

% --------------------------------------------------------------------

% Plot solution error norm

figure; plot(mean(err_tv,2)); hold on; plot(mean(err_t,2),'r--');

xlabel('Iteration number'); ylabel('Solution error norm'); legend('TV','SIDER');

title('Solution error vs. iteration number (mean of ventilation images)');

% Ventilation images

figure;

subplot(3,2,1);

imagesc(abs(uTargetSmooth(:,:,1))); colormap hot; axis image; axis off; colorbar; title('u(b=0) fully sampled');

subplot(3,2,2);

imagesc(abs(uTargetSmooth(:,:,5))); colormap hot; axis image; axis off; colorbar; title('u(b=6.4) fully sampled');

subplot(3,2,3);

imagesc(abs(u_tv(:,:,1))); colormap hot; axis image; axis off; colorbar; title('u(b=0) TV');

subplot(3,2,4);

imagesc(abs(u_tv(:,:,5))); colormap hot; axis image; axis off; colorbar; title('u(b=6.4) TV');

subplot(3,2,5);

imagesc(abs(u_t(:,:,1))); colormap hot; axis image; axis off; colorbar; title('u(b=0) SIDER');

subplot(3,2,6);

imagesc(abs(u_t(:,:,5))); colormap hot; axis image; axis off; colorbar; title('u(b=6.4) SIDER');

figure;

subplot(3,3,1); imagesc(abs(Drec0)); colormap hot; axis image; axis off; colorbar; title('D fully sampled'); caxis([0 0.9]);

subplot(3,3,2); imagesc(abs(alpharec0)); colormap hot; axis image; axis off; colorbar; title('\alpha fully sampled'); caxis([0 1.3]);

subplot(3,3,3); imagesc(abs(Lm0)); colormap hot; axis image; axis off; colorbar; title('Lm fully sampled'); caxis([150 450]);

subplot(3,3,4); imagesc(abs(Drectv)); colormap hot; axis image; axis off; colorbar; title('D TV'); caxis([0 0.9]);

subplot(3,3,5); imagesc(abs(alpharectv)); colormap hot; axis image; axis off; colorbar; title('\alpha TV'); caxis([0 1.3]);

subplot(3,3,6); imagesc(abs(Lmtv)); colormap hot; axis image; axis off; colorbar; title('Lm TV');caxis([150 450]);

subplot(3,3,7); imagesc(abs(Drect)); colormap hot; axis image; axis off; colorbar; title('D SIDER');caxis([0 0.9]);

subplot(3,3,8); imagesc(abs(alpharect)); colormap hot; axis image; axis off; colorbar; title('\alpha SIDER');caxis([0 1.3]);

subplot(3,3,9); imagesc(abs(Lmt)); colormap hot; axis image; axis off; colorbar; title('Lm SIDER');caxis([150 450]);

pause(1);

end

% -------------------------------------------------------------------------

% -------------------------------------------------------------------------

function [obj,JacThis] = objVecPotExpJac(x,A0,A,b,mask)

% Fit vntilation images to a stretch exponential model to get maps

% of diffusion (D) and heterogeneity index (alpha)

obj     = objVecPotExp(x,A0,A,b,mask);

JacThis = objVecPotJac(x,A0,A,b,mask);

%     end

end

function obj = objVecPotExp(x,A0,A,b,mask)

% Cost function for nonlinear LS

n           = size(A);

nm          = nnz(mask);

Dval        = x(1:nm);

alphaval    = x(nm+1:2*nm);

DThis       = zeros(n(1:2));

DThis(mask) = Dval;

alphaThis   = zeros(n(1:2));

alphaThis(mask) = alphaval;

obj         = [];

for iq = 1:size(A,3)

    tmp     = A0.*exp( -( (b(iq)*DThis).^alphaThis) ) - A(:,:,iq);

    obj     = [obj; (tmp(mask))];

end

end

function JacThis = objVecPotJac(x,A0,A,b,mask)

% Jacobian for nonlinear LS

nm          = nnz(mask);

Dval        = x(1:nm);

alphaval    = x(nm+1:2*nm);

JacThis     = [];

for ig  = 1:size(A,3)

    % Derivative wrt D

    A0e         = A0(mask).*exp(-((b(ig)*Dval).^alphaval));

    ba          = b(ig).^alphaval;

    Da          = Dval.^(alphaval-1);

    JacThisD    = -alphaval.*ba.*Da.*A0e;

    

    % Derivative wrt a

    logbD       = log(b(ig)*Dval);

    bDa         = (b(ig)*Dval).^alphaval;

    JacThisA    = -bDa.*logbD.*A0e;

    

    tmp         = [spdiags(JacThisD,0,nm,nm) spdiags(JacThisA,0,nm,nm)];

    JacThis     = [JacThis; tmp];

end

end

function [Lm,Lstd] = ComputeLmImage(N,Drec0,alpharec0)

% Compute maps of mean alveolar length (Lm) given maps of diffusion

% (D) and heterogeneity index (alpha)

%

Lm      = zeros(N(1:2));

Lstd    = zeros(N(1:2));

indThis = find( (Drec0 > 0).*(alpharec0 > 0) );

for ia = 1:length(indThis)

    [Lm_val,Lstd_val] = Lmdiststretchedexp(Drec0(indThis(ia)),alpharec0(indThis(ia)));

    Lm(indThis(ia))      = Lm_val;

    Lstd(indThis(ia))    = Lstd_val;

end % ia

end

function [Lm,Lmstd,Hk,ldd]=Lmdiststretchedexp(DDC,alpha,diftime, Do)

if nargin==3,

    Do=0.9;

end

if nargin==2,

    Do=0.88;

    diftime=1.6e-3;

end

D=0.01:0.001:Do;

beta=alpha;

b=0:0.5:10;

t0=1./DDC;

k=D;

betavect=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9];

betai=0.01:0.01:1;

B=[0.145 0.197 0.243 0.285 0.382 0.306 0.360 0.435 0.7];

C=[0.89 0.50 0.35 0.25 0 0.13 0.22 0.4015 0.33];

Bi=interp1(betavect,B,betai,'cubic');

Ci=interp1(betavect,C,betai,'cubic');

delta=beta.*(beta-0.5)./(1-beta);

fk=1+Ci(round(beta.*100)).*(k.*t0).^delta;

Hk=(t0.*Bi(round(beta.*100))./(k*t0).^((1-beta/2)/(1-beta))).*exp(-((1-beta).*beta.^(beta./(1-beta)))./(k.*t0).^(beta./(1-beta))).*fk;

Hk=Hk/sum(Hk);

ldd=1e4*sqrt(2*D*diftime);

Lm=sum(Hk.*ldd);

Lmstd=sqrt(sum(((ldd-Lm).^2).*Hk));

end

%

⛄ 运行结果

⛄ 参考文献

[1] Abascal J ,  Desco M ,  Parra-Robles J . Incorporation of prior knowledge of the signal behavior into the reconstruction to accelerate the acquisition of MR diffusion data[J]. arXiv e-prints, 2017.

❤️ 关注我领取海量matlab电子书和数学建模资料

❤️部分理论引用网络文献,若有侵权联系博主删除

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

matlab科研助手

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值