影象去噪序列——BM3D影象去噪模型實現
阿新 • • 發佈：2019-01-07
1. BM3D模型簡介
BM3D模型是一個兩階段影象去噪方法，主要包含兩個步驟:
(1) 在噪聲影象上，利用區域性區域搜尋相似塊，並進行堆疊，在變換域(DCT域、FFT域）利用硬閾值去噪方法對堆疊的影象塊進行去噪，獲得堆疊相似塊的估計值，最後，根據均值權重進行聚合；
(2) 通過步驟(1) 獲取初步估計的影象，在初步估計的影象上進行相似塊的聚合； 然後，利用維納協同濾波進行影象去噪，從而，獲取最後的去噪結果
2. 模型實現(程式碼參考網路實現)：
% BM3D_Color_Demo % BM3D 在彩色影象上去噪 % Author: HSW % Date: 2018-05-06 % clc; close all; clear all; img_org = imread('timg.png'); figure(1); imshow(img_org); title('原影象'); % 加噪聲 sigma = 25; img_noise = double(img_org)+sigma * randn(size(img_org)); figure; imshow(img_noise / 255, []); title('噪聲影象'); img_denoise = BM3D_Color(img_noise, 0, sigma, 0, 1); figure; imshow(img_denoise / 255, []); title('去噪影象');
% BM3D_Gray_Demo % BM3D 在灰度影象上去噪 % Author: HSW % Date: 2018-05-06 % clc; close all; clear all; img_org = imread('timg.png'); img_gray = rgb2gray(img_org); figure(1); imshow(img_gray); title('原影象'); % 加噪聲 sigma = 25; img_noise = double(img_gray)+sigma * randn(size(img_gray)); figure; imshow(img_noise / 255, []); title('噪聲影象'); img_denoise = BM3D_Gray(img_noise, 0, sigma, 1); figure; imshow(img_denoise / 255, []); title('去噪影象');
function img_denoise = BM3D_Color(img_noise, tran_mode, sigma, color_mode, isDisplay) % BM3D實現去噪 % Inputs: % img_noise: 噪聲影象 % tran_mode: 變換方法： 預設值為0, tran_mode: = 0, fft; = 1, dct; = 2, dwt, = 3, db1 % sigma: 噪聲水平，預設值為10 % color_mode: 彩色影象去噪時採用的顏色空間, 預設值為0, color_mode: = 0, YUV; = 1, YCbCr; = 2, OPP % Ouputs: % img_out: 去噪影象 % 參考文獻：An Analysis and Implementation of the BM3D Image Denoising Method % Inputs: % img_in: 噪聲影象，必須為矩形方陣 % tran_mode: = 0, FFT; = 1, DCT; = 2, DWT, = 3, db1 % Outputs: % img_denoise: 去噪影象 % % if ~exist('isDisplay', 'var') isDisplay = 0; end if ~exist('color_mode', 'var') color_mode = 0; end if ~exist('sigma', 'var') sigma = 10; end if ~exist('tran_mode', 'var') tran_mode = 0; end [row, col, dims] = size(img_noise); img_trans = rgb2other(img_noise, color_mode); % First Step 引數 kHard = 8; % 塊大小 pHard = 4; % 塊移動間隔 lambda_distHard = 0; % 求相似的距離時，變換後，收縮的閾值 nHard = 40; % 搜尋視窗大小 NHard = 28; % 最多相似塊個數 tauHard = 5000; % 最大的相似距離for fft % kaiser視窗的引數，實際上並沒有特別大的影響 beta=2; Wwin2D = kaiser(kHard, beta) * kaiser(kHard, beta)'; % Second Step引數 kWien = kHard; pWien = pHard; lambda_distWien = lambda_distHard; nWien = nHard; NWien = NHard; tauWien = tauHard; sigma2 = sigma*sigma; if tran_mode == 0 % FFT lambda2d=400; lambda1d=500; lambda2d_wie=50; lambda1d_wie=500; elseif tran_mode == 1 % DCT lambda2d=50; lambda1d=80; lambda2d_wie=20; lambda1d_wie=60; elseif tran_mode == 2 % DWT lambda2d=50; lambda1d=80; lambda2d_wie=20; lambda1d_wie=60; end fprintf('BM3D: First Stage Start...\n'); %block為原始影象塊， tran_block為FFT變換且硬閾值截斷後的頻域係數(頻域， 計算距離的時候採用的是變換塊) [block_ch1, tran_block_ch1, block2row_idx_ch1, block2col_idx_ch1] = im2block(img_trans(:,:,1), kHard, pHard, lambda_distHard, 0); [block_ch2, tran_block_ch2, block2row_idx_ch2, block2col_idx_ch2] = im2block(img_trans(:,:,2), kHard, pHard, lambda_distHard, 0); [block_ch3, tran_block_ch3, block2row_idx_ch3, block2col_idx_ch3] = im2block(img_trans(:,:,3), kHard, pHard, lambda_distHard, 0); %bn_r和bn_c為行和列上的影象塊個數 bn_r = floor((row - kHard) / pHard) + 1; bn_c = floor((col - kHard) / pHard) + 1; %基礎估計的影象 img_basic_sum = zeros(row, col, 3); img_basic_weight = zeros(row, col, 3); %對每個塊遍歷 for i=1:bn_r for j=1:bn_c % 利用亮度通道進行相似塊搜尋 [sim_blk_ch1, sim_num, sim_blk_idx] = search_similar_block(i, j, block_ch1, tran_block_ch1, floor(nHard/pHard), bn_r, bn_c, tauHard, NHard); % 進行亮度通道處理 % 協同濾波： 公式(2) tran3d_blk_shrink_ch1 = transform_3d(sim_blk_ch1, tran_mode, lambda2d, lambda1d); tran3d_blk_shrink_ch2 = transform_3d(block_ch2(:,:,sim_blk_idx), tran_mode, lambda2d, lambda1d); tran3d_blk_shrink_ch3 = transform_3d(block_ch3(:,:,sim_blk_idx), tran_mode, lambda2d, lambda1d); % 聚合： 公式(3)中的說明 NHard_P_ch1 = nnz(tran3d_blk_shrink_ch1); NHard_P_ch2 = nnz(tran3d_blk_shrink_ch2); NHard_P_ch3 = nnz(tran3d_blk_shrink_ch3); if NHard_P_ch1 > 1 wHard_P_ch1 = 1 / NHard_P_ch1; else wHard_P_ch1 = 1; end if NHard_P_ch2 > 1 wHard_P_ch2 = 1 / NHard_P_ch2; else wHard_P_ch2 = 1; end if NHard_P_ch3 > 1 wHard_P_ch3 = 1 / NHard_P_ch3; else wHard_P_ch3 = 1; end blk_est_ch1 = inv_transform_3d(tran3d_blk_shrink_ch1,tran_mode); blk_est_ch1 = real(blk_est_ch1); blk_est_ch2 = inv_transform_3d(tran3d_blk_shrink_ch2, tran_mode); blk_est_ch2 = real(blk_est_ch2); blk_est_ch3 = inv_transform_3d(tran3d_blk_shrink_ch3, tran_mode); blk_est_ch3 = real(blk_est_ch3); % 公式(3): 對亮度通道，即第1個通道 for k=1:sim_num idx = sim_blk_idx(k); ir = block2row_idx_ch1(idx); jr = block2col_idx_ch1(idx); img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 1) = img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 1) + wHard_P_ch1 * blk_est_ch1(:, :, k); img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 1) = img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 1) + wHard_P_ch1; img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 2) = img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 2) + wHard_P_ch2 * blk_est_ch2(:, :, k); img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 2) = img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 2) + wHard_P_ch2; img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 3) = img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 3) + wHard_P_ch3 * blk_est_ch3(:, :, k); img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 3) = img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 3) + wHard_P_ch3; end end end img_basic = img_basic_sum ./ img_basic_weight; if isDisplay figure; img_rgb = other2rgb(img_basic, color_mode); imshow(img_rgb / 255.0 ,[]); title('BM3D:Fist Stage Result'); end fprintf('BM3D: First Stage End...\n'); fprintf('BM3D: Second Stage Start...\n'); [block_basic_ch1,tran_block_basic_ch1,block2row_idx_basic_ch1,block2col_idx_basic_ch1] = im2block(img_basic(:, :, 1), kWien, pWien, lambda_distWien, 0); [block_basic_ch2,tran_block_basic_ch2,block2row_idx_basic_ch3,block2col_idx_basic_ch2] = im2block(img_basic(:, :, 2), kWien, pWien, lambda_distWien, 0); [block_basic_ch3,tran_block_basic_ch3,block2row_idx_basic_ch3,block2col_idx_basic_ch3] = im2block(img_basic(:, :, 3), kWien, pWien, lambda_distWien, 0); bn_r = floor((row - kWien) / pWien) + 1; bn_c = floor((col - kWien) / pWien) + 1; img_wien_sum = zeros(row, col, 3); img_wien_weight = zeros(row, col, 3); for i=1:1:bn_r for j=1:1:bn_c % 公式(5), 利用亮度進行相似性搜尋 [sim_blk_basic_ch1, sim_num, sim_blk_basic_idx] = search_similar_block(i, j, block_basic_ch1, tran_block_basic_ch1, floor(nWien/pWien), bn_r, bn_c, tauWien, NWien); % 公式(6) tran3d_blk_basic_ch1 = transform_3d(sim_blk_basic_ch1, tran_mode, lambda2d_wie, lambda1d_wie); tran3d_blk_basic_ch2 = transform_3d(block_basic_ch2(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); tran3d_blk_basic_ch3 = transform_3d(block_basic_ch3(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); omega_P_ch1 = (tran3d_blk_basic_ch1.^2) ./ ((tran3d_blk_basic_ch1.^2) + sigma2); omega_P_ch2 = (tran3d_blk_basic_ch2.^2) ./ ((tran3d_blk_basic_ch2.^2) + sigma2); omega_P_ch3 = (tran3d_blk_basic_ch3.^2) ./ ((tran3d_blk_basic_ch3.^2) + sigma2); % 公式(7) tran3d_blk_ch1 = transform_3d(block_ch1(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); tran3d_blk_ch2 = transform_3d(block_ch2(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); tran3d_blk_ch3 = transform_3d(block_ch3(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); blk_est_ch1 = inv_transform_3d(omega_P_ch1 .* tran3d_blk_ch1, tran_mode); blk_est_ch2 = inv_transform_3d(omega_P_ch2 .* tran3d_blk_ch2, tran_mode); blk_est_ch3 = inv_transform_3d(omega_P_ch3 .* tran3d_blk_ch3, tran_mode); blk_est_ch1 = real(blk_est_ch1); blk_est_ch2 = real(blk_est_ch2); blk_est_ch3 = real(blk_est_ch3); NWien_P_ch1 = nnz(omega_P_ch1); NWien_P_ch2 = nnz(omega_P_ch2); NWien_P_ch3 = nnz(omega_P_ch3); if NWien_P_ch1 > 1 wWien_P_ch1 = 1 / (NWien_P_ch1); else wWien_P_ch1 = 1; end if NWien_P_ch2 > 1 wWien_P_ch2 = 1/(NWien_P_ch2); else wWien_P_ch2 = 1; end if NWien_P_ch3 > 1 wWien_P_ch3 = 1 / (NWien_P_ch3); else wWien_P_ch3 = 1; end % 公式(8) for k=1:sim_num idx=sim_blk_basic_idx(k); ir=block2row_idx_basic_ch1(idx); jr=block2col_idx_basic_ch1(idx); img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 1) = img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 1) + wWien_P_ch1 * blk_est_ch1(:, :, k); img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 1) = img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 1) + wWien_P_ch1; img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 2) = img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 2) + wWien_P_ch2 * blk_est_ch2(:, :, k); img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 2) = img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 2) + wWien_P_ch2; img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 3) = img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 3) + wWien_P_ch3 * blk_est_ch3(:, :, k); img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 3) = img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 3) + wWien_P_ch3; end end end img_other = img_wien_sum ./ img_wien_weight; img_denoise = other2rgb(img_other, color_mode); fprintf('BM3D: Second Stage End\n');
function img_denoise = BM3D_Gray(img_noise, tran_mode, sigma, isDisplay)
% 參考文獻：An Analysis and Implementation of the BM3D Image Denoising Method
% Inputs:
% img_noise: 灰度噪聲影象，必須為矩形方陣
% tran_mode: = 0, fft; = 1, dct; = 2, dwt, = 3, db1
% Outputs:
% img_denoise: 去噪影象
%
if ~exist('tran_mode', 'var')
tran_mode = 0;
end
if ~exist('sigma', 'var')
sigma = 10;
end
if ~exist('isDisplay', 'var')
isDisplay = 0;
end
[row,col] = size(img_noise);
% First Step 引數
kHard = 8; % 塊大小
pHard = 4; % 塊移動間隔
lambda_distHard = 0; % 求相似的距離時，變換後，收縮的閾值
nHard = 40; % 搜尋視窗大小
NHard = 28; % 最多相似塊個數
tauHard = 5000; % 最大的相似距離for fft
% kaiser視窗的引數，實際上並沒有特別大的影響
beta=2;
Wwin2D = kaiser(kHard, beta) * kaiser(kHard, beta)';
% Second Step引數
kWien = kHard;
pWien = pHard;
lambda_distWien = lambda_distHard;
nWien = nHard;
NWien = NHard;
tauWien = tauHard;
sigma2 = sigma*sigma;
if(tran_mode==0) %fft
lambda2d=400;
lambda1d=500;
lambda2d_wie=50;
lambda1d_wie=500;
elseif(tran_mode == 1) %dct
lambda2d=50;
lambda1d=80;
lambda2d_wie=20;
lambda1d_wie=60;
elseif(tran_mode == 2) %dwt
lambda2d=50;
lambda1d=80;
lambda2d_wie=20;
lambda1d_wie=60;
end
%block為原始影象塊， tran_block為FFT變換且硬閾值截斷後的頻域係數(頻域， 計算距離的時候採用的是變換塊)
[block,tran_block,block2row_idx,block2col_idx]=im2block(img_noise,kHard,pHard,lambda_distHard,0);
%bn_r和bn_c為行和列上的影象塊個數
bn_r=floor((row-kHard)/pHard)+1;
bn_c=floor((col-kHard)/pHard)+1;
%基礎估計的影象
img_basic_sum=zeros(row,col);
img_basic_weight=zeros(row,col);
%basic處理
fprintf('BM3D: First Stage Start...\n');
%對每個塊遍歷
for i=1:bn_r
for j=1:bn_c
[sim_blk,sim_num,sim_blk_idx]=search_similar_block(i,j,block,tran_block,floor(nHard/pHard),bn_r,bn_c,tauHard,NHard);
% 協同濾波： 公式(2)
tran3d_blk_shrink=transform_3d(sim_blk,tran_mode,lambda2d,lambda1d);
% 聚合： 公式(3)中的說明
NHard_P=nnz(tran3d_blk_shrink);
if(NHard_P >1)
wHard_P=1/NHard_P;
else
wHard_P=1;
end
blk_est =inv_transform_3d(tran3d_blk_shrink,tran_mode);
blk_est=real(blk_est);
% 公式(3)
for k=1:sim_num
idx=sim_blk_idx(k);
ir=block2row_idx(idx);
jr=block2col_idx(idx);
img_basic_sum(ir:ir+kHard-1,jr:jr+kHard-1) = img_basic_sum(ir:ir+kHard-1,jr:jr+kHard-1) + wHard_P*blk_est(:,:,k);
img_basic_weight(ir:ir+kHard-1,jr:jr+kHard-1) = img_basic_weight(ir:ir+kHard-1,jr:jr+kHard-1) + wHard_P;
end
end
end
fprintf('BM3D: First Stage End...\n');
img_basic=img_basic_sum./img_basic_weight;
if isDisplay
figure;
imshow(img_basic,[]);
title('BM3D:Fist Stage Result');
end
[block_basic,tran_block_basic,block2row_idx_basic,block2col_idx_basic] = im2block(img_basic,kWien,pWien,lambda_distWien,0);
bn_r=floor((row-kWien)/pWien)+1;
bn_c=floor((col-kWien)/pWien)+1;
img_wien_sum=zeros(row,col);
img_wien_weight=zeros(row,col);
fprintf('BM3D: Second Stage Start...\n');
for i=1:1:bn_r
for j=1:1:bn_c
% 公式(5)
[sim_blk_basic,sim_num,sim_blk_basic_idx] = search_similar_block(i,j,block_basic,tran_block_basic,floor(nWien/pWien),bn_r,bn_c,tauWien,NWien);
% 公式(6)
tran3d_blk_basic = transform_3d(sim_blk_basic,tran_mode,lambda2d_wie,lambda1d_wie);
omega_P=(tran3d_blk_basic.^2)./((tran3d_blk_basic.^2)+sigma2);
% 公式(7)
tran3d_blk = transform_3d(block(:,:,sim_blk_basic_idx),tran_mode,lambda2d_wie,lambda1d_wie);
blk_est=inv_transform_3d(omega_P.*tran3d_blk,tran_mode);
blk_est=real(blk_est);
NWien_P=nnz(omega_P);
if(NWien_P >1)
wWien_P=1/(NWien_P);
else
wWien_P=1;
end
% 公式(8)
for k=1:sim_num
idx=sim_blk_basic_idx(k);
ir=block2row_idx_basic(idx);
jr=block2col_idx_basic(idx);
img_wien_sum(ir:ir+kWien-1,jr:jr+kWien-1) = img_wien_sum(ir:ir+kWien-1,jr:jr+kWien-1) + wWien_P*blk_est(:,:,k);
img_wien_weight(ir:ir+kWien-1,jr:jr+kWien-1) = img_wien_weight(ir:ir+kWien-1,jr:jr+kWien-1) + wWien_P;
end
end
end
fprintf('BM3D: Second Stage End\n');
img_denoise = img_wien_sum./img_wien_weight;
function [block,transform_block,block2row_idx,block2col_idx] =im2block(img,k,p,lambda2D,delta)
% 實現影象分塊
% Inputs：
% k: 塊大小
% p: 塊移動步長
% lambda_2D: 收縮閾值
% delta: 收縮閾值
% Outputs:
% block: 返回的塊
% transform_block: 變換後的塊
% block2row_idx: 塊索引與影象塊的左上角行座標對應關係
% block2col_idx: 塊索引與影象塊的左上角列座標對應關係
%
[row,col] = size(img);
% 頻域去噪中的硬閾值，實際上原文中，對於噪聲方差小於40時thres = 0, 具體見公式(1)的說明第2點(即距離計算)
thres = lambda2D*delta*sqrt(2*log(row*col));
% r_num 和 c_num分別表示行和列上可以採集的塊的數目
r_num = floor((row-k)/p)+1;
c_num = floor((col-k)/p)+1;
block = zeros(k,k,r_num*c_num);
block2row_idx = [];
block2col_idx = [];
cnt = 1;
for i = 0:r_num-1
rs = 1+i*p;
for j = 0:c_num-1
cs = 1+j*p;
block(:,:,cnt) = img(rs:rs+k-1,cs:cs+k-1);
block2row_idx(cnt) = rs;
block2col_idx(cnt) = cs;
tr_b = fft2(block(:,:,cnt));
idx = find(abs(tr_b)<thres);
tr_b(idx) = 0;
transform_block(:,:,cnt) = tr_b;
cnt = cnt+1;
end
end
endfunction [blk_est]=inv_transform_3d(blk_tran3d,tran_mode)
% 3D 逆變換
% Inputs:
% blk_tran3d: 在頻域中，硬閾值濾波的影象塊
% tran_mode: 變換方法
% Outputs:
% blk_est:
%
global blk_tran1d_s;
global blk_2d_s;
[m,n,blk_num]=size(blk_tran3d);
blk_invtran1d=zeros(m,n,blk_num);
blk_est=zeros(m,n,blk_num);
if(tran_mode==0) %fft
for i=1:1:m
for j=1:1:n
blk_invtran1d(i,j,:)=ifft(blk_tran3d(i,j,:));
end
end
for i=1:1:blk_num
blk_est(:,:,i)=ifft2(blk_invtran1d(:,:,i));
end
elseif(tran_mode==1) %dct
for i=1:1:m
for j=1:1:n
blk_invtran1d(i,j,:)=idct(blk_tran3d(i,j,:));
end
end
for i=1:1:blk_num
blk_est(:,:,i)=idct2(blk_invtran1d(:,:,i));
end
elseif(tran_mode==2) %dwt
blk_num=length(blk_2d_s);
blk_c=waverec2(blk_tran3d,blk_tran1d_s,'haar');
blk_est=[];
for i=1:1:blk_num
blk_est(:,:,i)=waverec2(blk_c(:,i),blk_2d_s{i},'Bior1.5');
end
else
error('tran_mode error');
end
endfunction img_trans = other2rgb(img_in, color_mode)
% 將RGB顏色空間轉為其他顏色空間
% Inputs:
% img_in: RGB顏色空間影象
% color_mode: 彩色影象去噪時採用的顏色空間, 預設值為0, color_mode: = 0, YUV; = 1, YCbCr; = 2, OPP
% Outputs:
% img_trans: 其他顏色空間
%
% Author: HSW
% Date: 2018-05-06
img_trans = zeros(size(img_in));
[row, col, dims] = size(img_in);
if color_mode == 0
color_tran = [0.30, 0.59, 0.11; -0.15, -0.29, 0.44; 0.61, -0.51, -0.10];
color_tran_inv = inv(color_tran);
for i = 1:row
for j = 1:col
other = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)];
img_trans(i, j, :) = color_tran_inv * other;
end
end
elseif color_mode == 1
color_tran = [0.30, 0.59, 0.11; -0.17, -0.33, 0.50; 0.50, -0.42, -0.08];
color_tran_inv = inv(color_tran);
for i = 1:row
for j = 1:col
other = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)];
img_trans(i, j, :) = color_tran_inv * other;
end
end
elseif color_mode == 2
color_tran = [1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0; 1.0 / 2.0, 0, -1.0 / 2.0; 1.0 / 4.0, -1.0 / 2.0, 1.0 / 4.0];
color_tran_inv = inv(color_tran);
for i = 1:row
for j = 1:col
other = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)];
img_trans(i, j, :) = color_tran_inv * other;
end
end
end
endfunction img_trans = rgb2other(img_in, color_mode)
% 將RGB顏色空間轉為其他顏色空間
% Inputs:
% img_in: RGB顏色空間影象
% color_mode: 彩色影象去噪時採用的顏色空間, 預設值為0, color_mode: = 0, YUV; = 1, YCbCr; = 2, OPP
% Outputs:
% img_trans: 其他顏色空間
%
% Author: HSW
% Date: 2018-05-06
img_trans = zeros(size(img_in));
[row, col, dims] = size(img_in);
if color_mode == 0
color_tran = [0.30, 0.59, 0.11; -0.15, -0.29, 0.44; 0.61, -0.51, -0.10];
for i = 1:row
for j = 1:col
rgb = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)];
img_trans(i, j, :) = (color_tran * rgb)';
end
end
elseif color_mode == 1
color_tran = [0.30, 0.59, 0.11; -0.17, -0.33, 0.50; 0.50, -0.42, -0.08];
for i = 1:row
for j = 1:col
rgb = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)];
img_trans(i, j, :) = (color_tran * rgb)';
end
end
elseif color_mode == 2
color_tran = [1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0; 1.0 / 2.0, 0, -1.0 / 2.0; 1.0 / 4.0, -1.0 / 2.0, 1.0 / 4.0];
for i = 1:row
for j = 1:col
rgb = [img_in(i, j, 1); img_in(i, j, 2); img_in(i, j, 3)];
img_trans(i, j, :) = (color_tran * rgb)';
end
end
end
endfunction [sim_blk,sim_num,sim_blk_idx]=search_similar_block(ik,jk,block,tran_block,np,bn_r,bn_c,tau,max_sim_num)
% 搜尋相似塊
% Inputs:
% ik, jk： 待搜尋相似塊的索引
% block: 影象塊集合
% tran_block： 影象塊FFT硬閾值過濾後的FFT係數
% k: 影象塊大小
% np: floor(nHard / pHard)， 其中nHard表示影象的搜尋區域大小， pHard表示塊的移動步長
% bn_r, bn_c: 影象總的行/列可以採集影象塊的數目
% tau: 影象塊相似性判斷閾值，見公式(1)
% max_sim_num: 最多保留相似塊的數目
% Ouputs:
% sim_blk:
% sim_num:
% sim_blk_idx:
%
% 搜尋視窗的左上角，右下角的塊索引
in_s = max(ik-floor(np/2),1);
jn_s = max(jk-floor(np/2),1);
in_e = min(ik+floor(np/2),bn_r);
jn_e = min(jk+floor(np/2),bn_c);
% 當前參考塊
ref_blk = tran_block(:,:,((ik-1)*bn_c+jk));
ii = in_s:1:in_e;
jj = jn_s:1:jn_e;
[II,JJ] = meshgrid(ii,jj);
IDX = (II-1)*bn_c+JJ;
blk_idx=IDX(:);
% 收縮範圍內的全部影象塊
cur_blk=tran_block(:,:,blk_idx);
cnt=size(cur_blk,3);
ref_blk_mat=repmat(ref_blk,[1,1,cnt]);
delta_blk=cur_blk-ref_blk_mat;
dist=sum(sum(delta_blk.*delta_blk,1),2);
[dist_sort,dist_idx]=sort(dist);
% 最大相似塊是真實相似塊和目標引數相似塊的最小值
max_num=min(cnt,max_sim_num);
if(dist_sort(max_num)<tau)
sim_num=max_num;
else
sim_num=sum(dist_sort(1:max_num)<tau);
end
cnt_idx=dist_idx(1:sim_num);
sim_blk_idx=blk_idx(cnt_idx);
sim_blk=block(:,:,sim_blk_idx);
endfunction [val]=thres_shrink(data,thres)
% 進行閾值截斷： 即 data(i) < thres ? data(i) = 0 : data(i) = data(i)
% Inputs:
% data: 閾值截斷前的資料
% thres: 閾值
% Ouputs:
% val: 閾值截斷後的資料
%
val=data;
idx=find(abs(data)<thres);
val(idx)=0;
endfunction blk_tran3d = transform_3d(blk_3d,tran_mode,lambda2d,lambda1d)
% 進行3D變換，即Collaborative Filtering: 在影象塊內進行2D變換，在影象塊間進行1D變換
% 公式(2)
% Inputs:
% blk_3d:
% tran_mode:
% Ouputs:
%
global blk_tran1d_s;
global blk_2d_s;
[m,n,blk_num]=size(blk_3d);
%變換不同時，可能需要修改？？
blk_2d_shrink=zeros(m,n,blk_num);
blk_1d_shrink=zeros(m,n,blk_num);
if(tran_mode==0) %fft
for i=1:1:blk_num
blk_tran2d = fft2(blk_3d(:,:,i));
blk_2d_shrink(:,:,i) = thres_shrink(blk_tran2d,lambda2d);
end
for i=1:1:m
for j=1:1:n
blk_tran1d = fft(blk_2d_shrink(i,j,:));
blk_1d_shrink(i,j,:) = thres_shrink(blk_tran1d,lambda1d);
end
end
blk_tran3d=blk_1d_shrink;
elseif(tran_mode==1) %dct
for i=1:1:blk_num
blk_tran2d=dct2(blk_3d(:,:,i));
blk_2d_shrink(:,:,i)=thres_shrink(blk_tran2d,lambda2d);
end
for i=1:1:m
for j=1:1:n
blk_tran1d=dct(blk_2d_shrink(i,j,:));
blk_1d_shrink(i,j,:)=thres_shrink(blk_tran1d,lambda1d);
end
end
blk_tran3d=blk_1d_shrink;
elseif(tran_mode==2) %dwt
blk_2d_s={};
blk_2d_shrink=[];%zeros()
for i=1:1:blk_num
[blk_tran2d_c,blk_tran2d_s]=wavedec2(blk_3d(:,:,i),2,'Bior1.5');
blk_2d_shrink(:,i)=thres_shrink(blk_tran2d_c,lambda2d);
blk_2d_s{i}=blk_tran2d_s;
end
%這裡應該用 wavedec.因為是對1維？？
[blk_tran1d_c,blk_tran1d_s]=wavedec2(blk_2d_shrink,1,'haar');
blk_tran3d=thres_shrink(blk_tran1d_c,lambda1d);
% elseif(strcmp(tran_mode,'db1')) %還未實現
% blk_2d_s={};
% blk_2d_shrink=[];%zeros()
% for i=1:1:blk_num
% [blk_tran2d_cA,blk_tran2d_cH,blk_tran2d_cV,blk_tran2d_cD]=...
% dwt2(blk_3d(:,:,i),'db1');
% blk_2d_shrink(:,i)=thres_shrink(blk_tran2d_c,lambda2d);
% blk_2d_s{i}=blk_tran2d_s;
% end
% [blk_tran1d_c,blk_tran1d_s]=wavedec2(blk_2d_shrink,1,'haar');
% blk_tran3d=thres_shrink(blk_tran1d_c,lambda1d);
else
error('tran_mode error');
end
end3. 模型效果：
3.1 灰度影象
3.2 彩色影象
