Image deblurring via extreme channels prior 代码解析(一)

文章《 Image deblurring via extreme channels prior》源自2017年CVPR，代码地址：​​​​​​​

typeisgod/ECP (github.com)

为了和代码对照，把论文中的重要知识点整理了一下：

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