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

原创

佛怒火莲 2022-10-24 09:12:27 ©著作权

文章标签 matlab 算法开发语言 ide 主成分分析 文章分类 数据结构与算法人工智能

©著作权归作者所有：来自51CTO博客作者佛怒火莲的原创作品，请联系作者获取转载授权，否则将追究法律责任

二、简介

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
 % juanabascal78@gmail.com, jmparra@hggm.es, desco@hggm.esswitch nargin
 case 14
 uTarget = varargin{1};
 end % nargintolKrylov = 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