图像重建基于matlab SIDER算法图像压缩重建含Matlab源码 2170期

Posted 2022-11-06 海神之光

一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【图像重建】基于matlab SIDER算法图像压缩重建【含Matlab源码 2170期】
点击上面蓝色字体，直接付费下载，即可。

获取代码方式2：
付费专栏图像处理（Matlab）

备注：
点击上面蓝色字体付费专栏图像处理（Matlab），扫描上面二维码，付费299.9元订阅海神之光博客付费专栏，凭支付凭证，私信博主，可免费获得5份本博客上传CSDN资源代码（有效期为订阅日起，三天内有效）；
点击CSDN资源下载链接：5份本博客上传CSDN资源代码

二、简介

1 PCA
PCA （Principle Component Analysis）是统计学中的主成分分析方法。主成分分析方法从矩阵角度讲也称K-L变换。使用PCA方法对图像进行压缩和重建的大致过程：
2 PCA图像压缩
图像压缩：首先将图像训练库里的每个二维图像拉伸成向量。然后对其进行主成分分析得到主成分的变换矩阵以及图像均值向量。图像压缩过程就是把待压缩的图像减去训练得到的图像均值向量并通过变换矩阵变换成维数很小的一个向量的过程。
3 PCA图像重建
图像重建：就是将压缩的图像通过变换矩阵的逆变换后再加上图像均值向量得到的压缩前向量的近似向量。因为是主成分分析，所以图像会有较少的信息损失，并且不能完全复原，但是这种信息损失一般是非常小的。

三、部分源代码

function [u,errAll] = SIDER(Dest,Alphaest,b,f,R,N,mu,lambda,gamma,alpha,beta,nBreg,mask,varargin)
% function [u,errAll] =
% SIDER(Dest,Alphaest,b,f,R,N,mu,lambda,gamma,alpha,beta,nBreg,mask,var
% argin)
%
% 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.
%
% 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.
%
% Code downloaded from repository:
% https://github.com/HGGM-LIM/compressed-sensing-diffusion-lung-MRI
%
% If you use this code, please, cite the following paper:
% JFPJ Abascal, M Desco, J Parra-Robles. Incorporation of prior knowledge
% of the signal behavior into the reconstruction to accelerate the
% acquisition of MR diffusion data,(submitted for publication) 2016.
%
% Departamento de Bioingenier韆 e Ingenier韆 Aeroespacial
% Universidad Carlos III de Madrid, Madrid, Spain
% [e-mail], [e-mail], [e-mail]

switch nargin
case 14
uTarget = varargin1;
end % nargin

tolKrylov = 1e-2;% 1e-2 % Krylov convergence criterion

rows = N(1);
cols = N(2);
numTime = N(3);
Nf = size(f);

R = reshape(R,Nf);
f = f.*R;

% Normalize data
normFactor = getNormalizationFactor(R,f);
f = normFactorf;
uTarget = uTargetnormFactor;

% Calculate the norm of the target on the given mask
% for ip = 1:numTime
% uTargetNorm(ip) = norm(reshape(uTarget(:,:,ip),[],1));
% end
for ip = 1:numTime
tmp = uTarget(:,:,ip);
uTargetNorm(ip) = norm(tmp(mask));
end

% Signal decay operator
dm = -[exp(-( (Destb(2:end)).^Alphaest - (Destb(1:end-1)).^Alphaest )) 0]';
L = spdiags([dm ones(N(3),1)],[-1 0],N(3),N(3));

LT = L’;
LTL = LT*L;

LMat = @(x)reshape((Lreshape(x,[prod(N(1:2)),N(3)])‘)’,N);
LTMat = @(x)reshape((LTreshape(x,[prod(N(1:2)),N(3)])‘)’,N);
LTLMat = @(x)reshape((LTL*reshape(x,[prod(N(1:2)),N(3)])‘)’,N);

uBest = zeros(rows,cols,numTime);
errBest = inf;

% Reserve memory for the auxillary variables
f0 = f;
v = zeros(rows,cols,numTime);
u = zeros(rows,cols,numTime);
x = zeros(rows,cols,numTime);
y = zeros(rows,cols,numTime);
bx = zeros(rows,cols,numTime);
by = zeros(rows,cols,numTime);
dx = zeros(rows,cols,numTime);
dy = zeros(rows,cols,numTime);

p = zeros(rows,cols,numTime);
bp = zeros(rows,cols,numTime);

murf = mu*ifft2(f);

h = waitbar(0,‘SIDER’);
h2 = figure;

% Do the reconstruction
for outer = 1:nBreg
rhs_wt = gammau;
rhs_p = lambdaLTMat(p-bp);
rhs_tv = lambda*(Dxt(x-bx)+Dyt(y-by));
rhs = murf+rhs_wt+rhs_p+rhs_tv;
u = reshape(krylov(rhs(😃),N);

% Derivatives
dx      = Dx(u);
dy      = Dy(u);
dp      = LMat(u);

% update x and y
[x,y]   = shrink2(dx+bx, dy+by,alpha/lambda);
p       = shrink1(dp+bp, beta/lambda);

% update bregman parameters
bx      = bx+dx-x;
by      = by+dy-y;
bp      = bp+dp-p;

% Bregman iteration for the data constraint
fForw   = fft2(u).*R;
f       = f + f0-fForw;
murf    = mu*(ifft2(f));

% for ip = 1:numTime
% % Solution error norm
% tmp = u(:,:,ip)-uTarget(:,:,ip);
% errAll(outer,ip) = norm(tmp(😃)/uTargetNorm(ip);
% end
% Solution error norm in a mask
for ip = 1:numTime
tmp = u(:,:,ip);
tmp2 = uTarget(:,:,ip);
errAll(outer,ip) = norm(tmp(mask)-tmp2(mask))/uTargetNorm(ip);
end

if nargin >= 14
    if any([outer ==1, outer == 100, outer == 500, rem(outer,250)==0])
        figure(h); waitbar(outer/nBreg,h);        
        ip = 1;
        close(h2);
        h2=figure;
        subplot(2,2,1);
        imagesc(abs(u(:,:,ip))); colorbar; title(['Iter ' num2str(outer)]);
        subplot(2,2,2);
        imagesc(abs(x(:,:,1))); colorbar; title(['x, nnz(x) % ' num2str(ceil(100*nnz(x(:))/prod(N)))]); colorbar;
        subplot(2,2,3);
        imagesc(abs(p(:,:,ip))); colorbar; title(['p, nnz(p) % ' num2str(ceil(100*nnz(p(:))/prod(N)))]); colorbar;
        subplot(2,2,4); plot(errAll(1:outer,:)); axis tight; title(['Sol. error' ]);
        colormap hot;
        drawnow;            
    end % rem
end % plots

if (mean(errAll(outer,:)) <= errBest)
    uBest       = u;
    errBest     = mean(errAll(outer,:));
end %errThis

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.