PCD.zip_CS yasuoganzhi_PCD_pcd算法_冗余字典_；图像修复；冗余字典；

% 此函数用PCD算法解决压缩感知（构造冗余字典用于图像修复） function PCD_Inpainting; clc;clear % 读入图像和掩模 I=imread('lena256.bmp'); % lena图像，可以尝试其它的图像 M=imread('mask_bmp.bmp'); % 掩模图像，可以尝试其它的图像 X=double(I)/255; % 原始矩阵归一化 mask=double(M>0); % 掩模(黑色的损伤区域0) N=length(X); % 原始图像大小 load ww1 % 小波变换矩阵1 load ww2 % 小波变换矩阵2 load ww3 % 小波变换矩阵3 load ww4 % 小波变换矩阵4 ww5=dct(eye(N,N)); % 离散余弦变换5 ww=[ww1;ww2;ww3;ww5]; % 构造冗余字典 N_d=4; % 字典数目 W_i=WL_for(ww,X); % 小波变换 max_e=max(max(full(abs(W_i)))); % 最大系数 threshold=0.1; % 稀疏化阈值(小波变换) th_v=(max_e*threshold); % 阈值 Y=X.*(mask); % 没有损伤的像素作为测量 X_r=zeros(N_d*N,N_d*N); % 恢复矩阵 N_iter=200; % 迭代次数 % 协方差估计谱半径和自适应阈值 diag(A'*A)^{-1} mm=10; sum1=0; for m=1:mm r_test=randn(N,N); % 高斯分布随机矩阵，协方差1 r_test=r_test.*(mask); % 没有损伤的像素作为测量 r_test=WL_for(ww,r_test); % 冗余字典变换，得到变换域图像 sum1=sum1+cov(r_test); % 协方差求和 end sum1=sum1/mm; % 平均值 sum2=1./diag(abs(sum1)); % 得到字典里面每个原子范数的逆 W=diag(sum2); % 得到下降方向的矩阵形式 gamma=W*ones(N_d*N,N_d*N); % 得到阈值的矩阵形式 mu=0.1; % 线搜索迭代参数(***需要修改保证稳定) for m=1:N_iter N_iter-m X_b=X_r; % 保留初始变换域系数 X_t=WL_back(ww,X_r); % 冗余字典反变换，得到时域图像 X_m=mask.*X_t; % 没有损伤的像素作为测量 R_m=mask.*(Y-X_m); % 测量残差 X_s=X_r+W*WL_for(ww,R_m); % 冗余字典正变换，更新变换域系数 X_r=Threshold_F(X_s,th_v*gamma); % 变换域系数软阈值 X_r=X_b+mu*(X_r-X_b); % 更新（需要线搜索自适应步长mu） th_v=th_v*0.90; % 软阈值松弛 if (th_v<1e-3) % 阈值很小时候截断迭代 break; end % 和原始图像相比的峰值信噪比PSNR errorx=sum(sum(abs(WL_back(ww,X_r)-X).^2)); % MSE误差 psnr=10*log10(1*1/(errorx/N/N)); % PSNR disp('峰值信噪比：') disp(psnr) end figure(1) subplot(2,2,1) imshow(uint8(X*255)) title('原始图像') subplot(2,2,2) imshow(uint8((Y)*255)) title('空域测量') subplot(2,2,3) imshow(uint8(full(WL_back(ww,X_r))*255)) ss=sprintf('恢复图像（PSNR: %0.2f）',psnr); title(ss) subplot(2,2,4) imshow(uint8(full(WL_back(ww,X_r)-X)*255)) title('误差图像') % 冗余字典变换和冗余字典反变换（可以不显示构造） function MM=WL_for(ww,M) MM=(ww*(M)*ww'); % 2d image function MM=WL_back(ww,M) MM=(ww'*(M)*ww); % 2d image % 阈值函数 function X_th=Threshold_F(X,threshold) % 软阈值 flag1=X>(threshold); flag2=X<-(threshold); flag3=flag1 | flag2; X_th=X-double(flag1).*threshold; X_th=X_th+double(flag2).*threshold; X_th=X_th.*flag3;