文章《 Image deblurring via extreme channels prior》源自2017年CVPR,代码地址:
为了和代码对照,把论文中的重要知识点整理了一下:
1.作者提出的模型:
(10)
D(l)和B(l)分别代表潜在图像l的暗通道和亮通道。
2.亮通道:
(2)
Ω(x)表示以x为中心的区域。
暗通道:把上式的两个max更换为min。
3.熟悉的场景,分别优化求解潜在图像l和模糊核k。
(11)
(12)
4.用半二次分裂求解l的L0正则化:
(13)
引入了变量p、q和g,分别近似于D(l)、1-B(l)和▽l。其中
然后交替求解l、p、q、g,同时固定其他变量为定值,比如求解l时认为p q 和g为一个定值:
(14)
像比于上一个公式,μ||g||0、λ||p||0和η||q||0为定值,与l无关,所以在求解l时被去掉。为了省事,将1-B(l)替换为D(1-l)。因为D(·)是非线性的,将他转为一个线性操作M,代表像素点到他的暗通道值的映射(该位置像素点映射到以该位置为中心一定范围内的最低像素值)
(15)
在去模糊过程中,我们使用中间潜在图像来计算M,当潜在图像接近清晰图像时,Ml与D(l)十分接近。使用映射分别代表D(l)和D(1-l),所以,(8)式转化为
(16)
(16)式的解为:
(17)
其中
▽h和▽v代表水平和垂直梯度。
再分别求解g、p和q:
(18-20)
(18)的解为
(21)
(19-20)的解同理
5.模糊核k的估计:为了快速求解,将||l*k-b||替换为||▽l*k-▽b||
(22)
最后解为:
(23)
6.预设值:μ=λ=η=0.004,γ=2 暗通道块的尺寸为35(以当前像素为中心的,35*35大小的块,即往上下左右各17行/列)
7.总体流程
该代码的主入口为demo_deblurring.m,下载的压缩包中附带了两张测试图像,real_blur_img3.png和summerhouse.jpg 参数设置分别再115-117行、22-24行,取消原有的注释即可,两张图像的运行结果如下:
分别对应输入图像(模糊图像)、中间潜在图像、模糊核、输出图像(清晰图像)
一、demo_deblurring 主函数
1.1 首先添加了三个文件夹的地址,里面存放了之后被调用的他人的方法以及输入图像:
addpath(genpath('image'));
addpath(genpath('whyte_code'));
addpath(genpath('cho_code'));
1.2 之后定义了结构体opts的四个变量:
opts.prescale = 1; %%downsampling 没用到
opts.xk_iter = 5; %% the iterations 4.2中用到,每次计算k和潜在图像的迭代次数
opts.gamma_correct = 1.0;%用于2.1
opts.k_thresh = 20;%用于2.9
分别代表(?)、迭代次数为五次、灰度拉伸值为1(不进行灰度拉伸)、 阈值为20.
1.3 接下来,定义了输入图像名(filename)、模糊核尺寸(kernel_size)、是否过饱和(saturation)、暗通道权值(lambda_dark)、梯度权值(lambda_grad)、颜色校正值(gamma_correct)、tv正则化参数(lambda_tv)、l0正则化参数(lambda_l0)、是否要去振铃(weight_ring)。
filename = 'image\real_blur_img3.png'; opts.kernel_size = 35; saturation = 0;
lambda_dark = 4e-3; lambda_grad = 4e-3;opts.gamma_correct = 1.0;
lambda_tv = 0.001; lambda_l0 = 5e-4; weight_ring = 1;
1.4 接下来,读入图像,记作y。
y = imread(filename);
1.5 将图像转为double格式(之前为uint8格式)并且值缩放到0~1,如果是彩色图像,则转为灰度,储存为yg。
if size(y,3)==3
yg = im2double(rgb2gray(y));
else
yg = im2double(y);
end
1.6 调用自定义的blind_deconv(二)函数,传入yg、lambda_dark和lambda_grad、opts参数,返回得到估计的模糊核kernel和中间潜在图像interim_latent。
[kernel, interim_latent] = blind_deconv(yg, lambda_dark, lambda_grad, opts);
1.7 将y转为0-1的double类型:
y = im2double(y);
1.8 如果输入时非饱和的,则调用自定义的ringing_artifacts_removal(十二)函数,传入y、之前得到的模糊核kernel、lambda_tv、lambda_l0和weight_ring等参数,返回得到最后清晰图像Latent;否则调用自定义函数whyte_deconv(由于篇幅过长,这一块之后再整理)(在whyte_code文件夹下),传入y和模糊核kernel,得到清晰图像Latent。
if ~saturation
%% 1. TV-L2 denoising method
Latent = ringing_artifacts_removal(y, kernel, lambda_tv, lambda_l0, weight_ring);
else
%% 2. Whyte's deconvolution method (For saturated images)
Latent = whyte_deconv(y, kernel);
end
1.9 对模糊核k进行归一化,并分别将生成的模糊核k、清晰图像Latent和中间潜在图像interim_latent(1.6返回结果)写入到指定位置。
k = kernel - min(kernel(:));
k = k./max(k(:));
imwrite(k,['results\' filename(7:end-4) '_kernel.png']);
imwrite(Latent,['results\' filename(7:end-4) '_result.png']);
imwrite(interim_latent,['results\' filename(7:end-4) '_interim_result.png']);
二 、 1.6步调用的自定义函数blind_deconv
函数格式
function [kernel, interim_latent] = blind_deconv(y, lambda_dark, lambda_grad, opts)
y为灰度模糊图像,double类型,取值0-1,对应于(一)的yg变量。该函数用来多尺度地盲卷积。
2.1 如果输入的opts.gamma_correct不为1,则进行灰度拉伸:
if opts.gamma_correct~=1
y = y.^opts.gamma_correct;%灰度拉伸
end
如果gamma_correct小于1,那么会向高对比度拉伸(更亮);大于1则会向低对比度拉伸(更暗)。
2.2 计算金字塔层数num_scales
ret = sqrt(0.5);
%%
maxitr=max(floor(log(5/min(opts.kernel_size))/log(ret)),0); %迭代次数 6次
num_scales = maxitr + 1;
2.3 计算每一层的模糊核k1*k2,分别存在k1list和k2list中:
retv表示缩放倍数,第一个值为1,第二个值为1/sqrt(2)倍...
retv=ret.^[0:maxitr];
k1list=ceil(opts.kernel_size*retv);%每次迭代估计的模糊核大小
k1list=k1list+(mod(k1list,2)==0);%保证大小为奇数
k2list=ceil(opts.kernel_size*retv);
k2list=k2list+(mod(k2list,2)==0);%k1*k2大小
例如模糊核大小为35时,k1和k2都为35 25 19 13 9 7
2.4 迭代,从粗到细(模糊核由大向小,金字塔从上到下),每一步的模糊核初始为上一步模糊核的值进行重插值,如果是第一次迭代,则初始化为一个方形矩阵ks,第(minisize-1)/2行的最中间两个值为0.5,其他为0,比如一个5*5的:
(为啥不是最中间那行呢)
如果不是第一次,那么根据上一步的结果ks调整当前大小。
for s = num_scales:-1:1
if (s == num_scales)
%%
% at coarsest level, initialize kernel
ks = init_kernel(k1list(s));%第一层 初始化
k1 = k1list(s);
k2 = k1; % always square kernel assumed
else
% upsample kernel from previous level to next finer level
k1 = k1list(s);
k2 = k1; % always square kernel assumed
% resize kernel from previous level
ks = resizeKer(ks,1/ret,k1list(s),k2list(s));
end
%for循环还没有结束 直至2.9
初始化模糊核函数
function [k] = init_kernel(minsize)
k = zeros(minsize, minsize);
k((minsize - 1)/2, (minsize - 1)/2:(minsize - 1)/2+1) = 1/2;%最中间两个为0.5 其他为0
end
调整大小函数,部分解释加在注释中
调用fixsize的目的是,直接resize的话大小会根据倍数向上取整,可能和预期的大小(k1 k2)不等,所以通过fixsize调整到k1*k2。(要是我直接resize到[k1 k2])
function k=resizeKer(k,ret,k1,k2)
%%
% levin's code
k=imresize(k,ret);%放大根2倍,大小向上取整
k=max(k,0);%去掉小于0的值
k=fixsize(k,k1,k2);
if max(k(:))>0
k=k/sum(k(:));%归一化
end
end
%%
function nf=fixsize(f,nk1,nk2)
[k1,k2]=size(f);
while((k1~=nk1)|(k2~=nk2))%循环增改模糊核的大小直至k1=nk1 k2=nk2
if (k1>nk1)%如果resize的k的行数比预期多
s=sum(f,2);%每行的sum求和
if (s(1)<s(end))
f=f(2:end,:);%如果第一行的和比最后一行小,那么模糊核删除第一行
else
f=f(1:end-1,:);%否则删除最后一行
end
end
if (k1<nk1)%如果小于的话则需要增补一行
s=sum(f,2);
if (s(1)<s(end))
tf=zeros(k1+1,size(f,2));
tf(1:k1,:)=f;%在最后一行后补一行0
f=tf;
else
tf=zeros(k1+1,size(f,2));
tf(2:k1+1,:)=f;%在第一行前补一行0
f=tf;
end
end
%列和行同理
if (k2>nk2)
s=sum(f,1);
if (s(1)<s(end))
f=f(:,2:end);
else
f=f(:,1:end-1);
end
end
if (k2<nk2)
s=sum(f,1);
if (s(1)<s(end))
tf=zeros(size(f,1),k2+1);
tf(:,1:k2)=f;
f=tf;
else
tf=zeros(size(f,1),k2+1);
tf(:,2:k2+1)=f;
f=tf;
end
end
[k1,k2]=size(f);
end
nf=f;
end
2.5 对图像进行下采样得到ys:
cret=retv(s);
ys=downSmpImC(y,cret);
cret表示当前下采样倍数。
下采样函数(原理并没有很明白,就没有过多的解释):
function sI=downSmpImC(I,ret)
%% refer to Levin's code
if (ret==1) %如果ret为1,结果为本身
sI=I;
return
end
%%%%%%%%%%%%%%%%%%%
sig=1/pi*ret;
g0=[-50:50]*2*pi;
sf=exp(-0.5*g0.^2*sig^2);
sf=sf/sum(sf);
csf=cumsum(sf);%csf是sf的累积
csf=min(csf,csf(end:-1:1));%csf和csf反转的最小
ii=find(csf>0.05);
sf=sf(ii);
sum(sf);
I=conv2(sf,sf',I,'valid');%首先求I的各列与向量 sf 的卷积,然后求每行结果与向量 sf' 的卷积。
%valid:舍去需要补零的部分:裁剪为(ma-mb+1)*(na-nb+1) mb*nb为卷积核大小
[gx,gy]=meshgrid([1:1/ret:size(I,2)],[1:1/ret:size(I,1)]);
sI=interp2(I,gx,gy,'bilinear');%重插值
end
2.6 如果是第一次迭代,则调用自定义函数threshold_pxpy_v1(三),输入当前层图像ys和模糊核大小,输入为阈值threshold。
if (s == num_scales)
[~, ~, threshold]= threshold_pxpy_v1(ys,max(size(ks)));
end
2.7 调用自定义函数blind_deconv_mainBDF(四),输入当前层的图像ys、模糊核ks、lambda_dark和lambda_grad、阈值threshold和opts,得到估计的模糊核ks、新的权值lambda_dark和lambda_grad以及潜在图像interim_latent。
[ks, lambda_dark, lambda_grad, interim_latent] = blind_deconv_mainBDF(ys, ks, lambda_dark,...
lambda_grad, threshold, opts);
2.8 中心化模糊核,调用自定义函数adjust_psf_center(十一),输入模糊核ks,输出调整后的模糊核ks,将小于0的值置为0,并归一化。
ks = adjust_psf_center(ks);
ks(ks(:)<0) = 0;
sumk = sum(ks(:));
ks = ks./sumk;
2.9 如果是最后一层,那么ks为最终的模糊核kernel。如果opts.k_thresh(1.2最开始的预设)大于0,将模糊核kernel中,值小于最大值除以阈值opts.k_thresh的值置为0,比如kernel最大值为0.8,opts.k_thresh为20,那么置kernel小于0.04的值为0;否则置小于0的值为0.最后归一化。
if (s == 1)
kernel = ks;
if opts.k_thresh>0
kernel(kernel(:) < max(kernel(:))/opts.k_thresh) = 0;
else
kernel(kernel(:) < 0) = 0;
end
kernel = kernel / sum(kernel(:));
end
end%for迭代结束
此时,2.4迭代中的第一个for循环结束
2.10 函数blind_deconv结束
end