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 end
扫描二维码关注公众号,回复:
814528 查看本文章
function [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 end
function 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 end
function 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 end
function [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); end
function [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; end
function 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 end
3. 模型效果:
3.1 灰度图像
3.2 彩色图像