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

七、6.2调用的dark_channel函数

首先判断暗通道范围是否为奇数，然后对输入图像I进行四周边界扩充，类似于重复边界，这样扩充后的I进行暗通道计算时无需考虑边界值情况（比如位于第一列、第一行等）。遍历获得每一个暗通道块patch，找块在三个颜色通道（如果是彩色的话）的最小值，记作J，最小值的位置记作J_index。

function [J, J_index] = dark_channel(I, patch_size)
% function J = dark_channel(I, patch_size);

% Computes the "Dark Channel" of corresponding RGB image.
% -Finds from the input image the minimum value among all 
%  pixels within the patch centered around the location of the 
%  target pixel in the newly created dark channel image 'J'
%  J is a 2-D image (grayscale).

% Example: J = dark_channel(I, 15); % computes using 15x15 patch

% Check to see that the input is a color image
% if ndims(I) == 3
%     [M N C] = size(I);
%     J = zeros(M, N); % Create empty matrix for J
%     J_index = zeros(M, N); % Create empty index matrix
% else
%     error('Sorry, dark_channel supports only RGB images');
% end
%% for grayscale image
%
[M, N, C] = size(I);
J = zeros(M, N); % Create empty matrix for J
J_index = zeros(M, N); % Create empty index matrix

% Test if patch size has odd number
if ~mod(numel(patch_size),2) % if even number
    error('Invalid Patch Size: Only odd number sized patch supported.');
end

% pad original image
%I = padarray(I, [floor(patch_size./2) floor(patch_size./2)], 'symmetric');
I = padarray(I, [floor(patch_size./2) floor(patch_size./2)], 'replicate');%重复边界,为了不考虑块的临界

% Compute the dark channel 
 for m = 1:M
        for n = 1:N
            patch = I(m:(m+patch_size-1), n:(n+patch_size-1),:);
            tmp = min(patch, [], 3);%第三维度的最小值
            [tmp_val, tmp_idx] = min(tmp(:));
            J(m,n) = tmp_val;
            J_index(m,n) = tmp_idx;
        end
 end



end

八、6.3调用的assign_dark_channel_to_pixel函数

函数格式

function [outImg] = assign_dark_channel_to_pixel(S, dark_channel_refine, dark_channel_index, patch_size)

输入为模糊图像S、改进的暗通道dark_channel_refine（S暗通道图像去掉一些比较小的值为0），暗通道映射dark_channel_index(6.3为J_idx)，以及暗通道范围patch_size，返回S调整后的图像outImg。

8.1 padsize为patchsize/2向下取整，比如patchsize为35，那么表示以某像素为中心（第19行18列）的35*35大小的区域参数暗通道计算，某像素分别向上下左右的宽度为17（padsize），同样对S进行重复边界的扩充，宽度为padsize，为了不讨论边界点问题，比如像素位于图像边界时暗通道区域会变化。

[M N C] = size(S);
%outImge = zeros(M, N); % 

% pad original image
padsize = floor(patch_size./2);
S_padd = padarray(S, [padsize padsize], 'replicate');%重复

8.2 进行和(七)相似的操作，获得S的某个区域patch的最低值，即暗通道值，判断与输入的改进的暗通道值dark_channel_refine是否相等，如果不等，将dark_channel_refine对应的暗通道映射dark_channel_index在patch的值改为dark_channel_refine的值，将修改的patch回溯（赋值回）S_padd。

m = 1:M;
        for n = 1:N%dark_channel_refine：小于阈值置为0的暗通道
            patch = S_padd(m:(m+patch_size-1), n:(n+patch_size-1),:);
            if ~isequal(min(patch(:)), dark_channel_refine(m,n))
                patch(dark_channel_index(m,n)) = dark_channel_refine(m,n);
            end
            for cc = 1:C
                S_padd(m:(m+patch_size-1), n:(n+patch_size-1),cc) = patch(:,:,cc);
            end
        end%根据dark_channel_refine回修正S_padd

8.3 将S_padd裁剪掉在8.1步添加的边界，得到outImg，并把outImg内边界值替换为对应S的值。函数assign_dark_channel_to_pixel结束，返回outImg。

outImg = S_padd(padsize + 1: end - padsize, padsize + 1: end - padsize,:);%裁剪
%% boundary processing
outImg(1:padsize,:,:) = S(1:padsize,:,:);  outImg(end-padsize+1:end,:,:) = S(end-padsize+1:end,:,:);
outImg(:,1:padsize,:) = S(:,1:padsize,:);  outImg(:,end-padsize+1:end,:) = S(:,end-padsize+1:end,:);
%边界值替换
%figure(2); imshow([S, outImg],[]);
end

九、4.3调用的estimate_psf函数

函数格式

function psf = estimate_psf(blurred_x, blurred_y, latent_x, latent_y, weight, psf_size)

输入带有边界的模糊图像的梯度blurred_x, blurred_y（4.3输入为Bx和By），4.3中得到的清晰图像筛选后的梯度lantent_x和lantent_y和weight=2，和模糊和的尺寸，输出估计的模糊核。

9.1 分别将前四个参数转为频域。

    latent_xf = fft2(latent_x);
    latent_yf = fft2(latent_y);
    blurred_xf = fft2(blurred_x);
    blurred_yf = fft2(blurred_y);

9.2 计算b_f和p.m，并把b_f转为空间域，得到b。定义部分参数后，初始化psf后，调用自定义函数conjgrad(十)进行共轭梯度优化后，得到估计的psf，将小于psf值值和0.05的psf值置为0，重新归一化并返回，函数estimate_psf结束。

    latent_xf = fft2(latent_x);
    latent_yf = fft2(latent_y);
    blurred_xf = fft2(blurred_x);
    blurred_yf = fft2(blurred_y);
    % compute b = sum_i w_i latent_i * blurred_i
    b_f = conj(latent_xf)  .* blurred_xf ...
        + conj(latent_yf)  .* blurred_yf;
    b = real(otf2psf(b_f, psf_size));

    p.m = conj(latent_xf)  .* latent_xf ...
        + conj(latent_yf)  .* latent_yf;
    %p.img_size = size(blurred);
    p.img_size = size(blurred_xf);
    p.psf_size = psf_size;
    p.lambda = weight;

    psf = ones(psf_size) / prod(psf_size);
    psf = conjgrad(psf, b, 20, 1e-5, @compute_Ax, p);
    
    psf(psf < max(psf(:))*0.05) = 0;
    psf = psf / sum(psf(:));    
end

传入的compute_Ax函数：

function y = compute_Ax(x, p)
    x_f = psf2otf(x, p.img_size);
    y = otf2psf(p.m .* x_f, p.psf_size);
    y = y + p.lambda * x;
end

返回y，y的值为p.m与x的频域的点乘的逆傅里叶变换后与x的p.lambda的和。

十  9.2调用的conjgrad函数

函数格式

function x=conjgrad(x,b,maxIt,tol,Ax_func,func_param,visfunc)

x对应4.2的psf，maxIt迭代次数为20，tol为1e-5，Ax_func为compute_Ax,func_param为结构体p。

10.1 首先计算b与Ax_func(x,func_param)的差值，得到r，p初始值设为r，rsold为r的平方和。

    r = b - Ax_func(x,func_param);
    p = r;
    rsold = sum(r(:).*r(:));

10.2 迭代，每次以p和func_param为输入调用Ax_func，获得Ap，计算alpha为rsold除以p与Ap的点乘之和。x增加alpha倍的p，r减少alpha倍的Ap，计算新的r的平方和rsnew，如果根号下rsnew小于定值tol，那么退出迭代，否则更新p为r和rsnew/rsold倍的p的和，rsold置为rsnew。迭代完后，返回x，即估计后的模糊核psf，函数conjgrad结束。

    for iter=1:maxIt
        Ap = Ax_func(p,func_param);
        alpha = rsold/sum(p(:).*Ap(:));
        x=x+alpha*p;
        if exist('visfunc', 'var')
            visfunc(x, iter, func_param);
        end
        r=r-alpha*Ap;
        rsnew=sum(r(:).*r(:));
        if sqrt(rsnew)<tol
            break;
        end
        p=r+rsnew/rsold*p;
        rsold=rsnew;
    end
end

十一、2.8调用的adjust_psf_center函数（cho_code文件夹下）

输入psf，输出调整后的psf。

11.1 创建与psf等大的网络矢量X和Y。分别与psf点乘后求和得到xc1和yc1。xc2和yc2取值为psf最中间列号和行号。也即X每行和Y每列的均值，xshift和yshift表示了xc1与均值的偏移的四舍五入。将psf调用warpimage得到调整后的psf。

function psf = adjust_psf_center(psf)

[X Y] = meshgrid(1:size(psf,2), 1:size(psf,1));
xc1 = sum2(psf .* X);
yc1 = sum2(psf .* Y);
xc2 = (size(psf,2)+1) / 2;
yc2 = (size(psf,1)+1) / 2;
xshift = round(xc2 - xc1);
yshift = round(yc2 - yc1);
psf = warpimage(psf, [1 0 -xshift; 0 1 -yshift]);
end
function val = sum2(arr)
val = sum(arr(:));
end

11.2 warpimage函数：输入待调整的psf和调整矩阵M（[1 0 -xshift; 0 1 -yshift]）

function warped = warpimage(img, M)
if size(img,3) == 3
    warped(:,:,1) = warpProjective2(img(:,:,1), M);
    warped(:,:,2) = warpProjective2(img(:,:,2), M);
    warped(:,:,3) = warpProjective2(img(:,:,3), M);
    warped(isnan(warped))=0;
else
    warped = warpProjective2(img, M);
    warped(isnan(warped))=0;
end
end

实质调用了自定义的扭曲投影函数warpProjective2，得到扭曲后的psf，记作warped，并将其中的NaN的值改为0.返回warped。

11.3 warpProjective2输入待调整的psf和调整矩阵M（[1 0 -xshift; 0 1 -yshift]）

coords为一个k*2大小的坐标，k=m*n（图像的大小），coords每一列代表一个不重复的坐标值。取值1到m或n。homogeneousCoords在coords的基础上增加了第三行全1矩阵，扭曲坐标warpedCoords为A（[1 0 -xshift; 0 1 -yshift]）与homogeneousCoords的乘积（不是点乘），结果warpedCoords的第一行为调整后的网络矢量xprime，yprime同理，使用系统函数interp2对im图像以原网络矢量x和y线性插值到xprime和yprime网络矢量，结果为result，并调整result大小和输入图像大小im相同。result即调整后的psf。

function result = warpProjective2(im,A)

if (size(A,1)>2)
  A=A(1:2,:);
end

% Compute coordinates corresponding to input 
% and transformed coordinates for result
[x,y]=meshgrid(1:size(im,2),1:size(im,1));
coords=[x(:)'; y(:)'];
homogeneousCoords=[coords; ones(1,prod(size(im)))];
warpedCoords=A*homogeneousCoords;
xprime=warpedCoords(1,:);%./warpedCoords(3,:);
yprime=warpedCoords(2,:);%./warpedCoords(3,:);

result = interp2(x,y,im,xprime,yprime, 'linear');
result = reshape(result,size(im));

return;
end

十二  1.8调用的ringing_artifacts_removal函数

1.6步完成时已经得到了估计的模糊核kernel和中间潜在图像interim_latent(黑白图像)。函数格式如下：

function [result] = ringing_artifacts_removal(y, kernel, lambda_tv, lambda_l0, weight_ring)

该函数用来根据模糊核得到清晰图像，并抑制振铃。传入0-1的初始模糊图像y、1.6得到的模糊和kernel、两个权值参数，weight_ring用来表示抑制振铃效应的程度，为0表示没有振铃不需要去除，返回最后去除的结果result。

12.1 对y同样添加liu的边界，得到y_pad，对其每个颜色通道，分别调用自定义函数deblurring_adm_aniso(十三)，传入当前颜色通道y_pad(:,:,c)、模糊核kernel、权值参数lambda_tv和alpha=1。得到tv正则化的结果Latent_tv，将Latent_tv裁剪掉之前添加的边界。

H = size(y,1);    W = size(y,2);
y_pad = wrap_boundary_liu(y, opt_fft_size([H W]+size(kernel)-1));
Latent_tv = [];
for c = 1:size(y,3)
    Latent_tv(:,:,c) = deblurring_adm_aniso(y_pad(:,:,c), kernel, lambda_tv, 1);
end
Latent_tv = Latent_tv(1:H, 1:W, :);

12.2 如果weight_ring为0，那么不需要对Latent_tv进行后续操作，直接为函数的返回值。

if weight_ring==0
    result = Latent_tv;
    return;
end

12.3 否则，对y_pad进行L0正则化，调用自定义函数L0Restoration（五），传入y_pad、kernel、权值lambda_l0和kappa=2，得到正则化结果Latent_l0，并同样裁剪掉添加的边界。

Latent_l0 = L0Restoration(y_pad, kernel, lambda_l0, 2);
Latent_l0 = Latent_l0(1:H, 1:W, :);

12.4 得到对y_padd进行tv正则化和l0正则化结果的差异diff，对差异值调用自定义函数bilateral_filter（十四）进行滤波，最终去振铃的结果为tv正则化的结果Latent_tv减去weight_ring倍的滤波后的差异bf_diff得到result，函数ringing_artifacts_removal结束。

diff = Latent_tv - Latent_l0;
bf_diff = bilateral_filter(diff, 3, 0.1);
result = Latent_tv - weight_ring*bf_diff;

十三、12.1调用的deblurring_adm_aniso函数

函数格式

function [I] = deblurring_adm_aniso(B, k, lambda, alpha)

输入参数为模糊图像（灰度）B（对应12.1的y_padd的单个颜色通道）、模糊核k，权值lambda和alpha=1.输出正则化结果I。这个方法来自于使用超拉普拉斯先验的快速图像去卷积，采用了aniso TV正则化方法。

13.1 用于迭代和计算的beta初始化为1/lambda,迭代截至下界为beta_min，I初始化为输入值B，根据B和k调用computeDenominator计算分母，得到Nomin1、 Denom1和 Denom2。

beta = 1/lambda;
beta_min = 0.001;
[m n] = size(B); 
% initialize with input or passed in initialization
I = B; 
[Nomin1, Denom1, Denom2] = computeDenominator(B, k);

computeDenominator函数如下：输入参数为y（被调用输入的变量为x）和k。Nomin1为模糊核k的频域的共轭与y的频域的点乘，Denom1为模糊核k的频域的模的平方，Denom2为将水平和垂直梯度算子[1,-1]和[1;-1]转为频域后模平方之和。

function [Nomin1, Denom1, Denom2] = computeDenominator(y, k)
sizey = size(y);
otfk  = psf2otf(k, sizey); 
Nomin1 = conj(otfk).*fft2(y);
Denom1 = abs(otfk).^2; 
% if higher-order filters are used, they must be added here too
Denom2 = abs(psf2otf([1,-1],sizey)).^2 + abs(psf2otf([1;-1],sizey)).^2;
end

13.2 计算I的差分Ix和Iy。

Ix = [diff(I, 1, 2), I(:,1) - I(:,n)]; 
Iy = [diff(I, 1, 1); I(1,:) - I(m,:)]; 

13.3 循环，gamma值为1/(2*beta)，分母Denom为刚计算的Denom1与gamma倍的Denom2之和。Wx和Wy分别为差分Ix（Iy）的绝对值的每个元素与beta*lambda的大小，如果大于，那么Wx对应元素为Ix的值，否则为0.Wxx为Wx的水平方向上的负差分与Wy在垂直方向上的负差分之和。最后频域的结果为(Nomin1 + gamma*fft2(Wxx))./Denom，将该结果转为空间域取实数部分得到更新的I，重新计算I的差分Ix和Iy，beta缩小为之前的二分之一，直至beta小于betamin，循环结束。函数结束，返回最后的I。

while beta > beta_min
    gamma = 1/(2*beta);
    Denom = Denom1 + gamma*Denom2;
    % subproblem for regularization term
    if alpha==1
        Wx = max(abs(Ix) - beta*lambda, 0).*sign(Ix);
        Wy = max(abs(Iy) - beta*lambda, 0).*sign(Iy);
        %%
    else% 恒为1 跳过
        Wx = solve_image(Ix, 1/(beta*lambda), alpha);
        Wy = solve_image(Iy, 1/(beta*lambda), alpha);
    end
      Wxx = [Wx(:,n) - Wx(:, 1), -diff(Wx,1,2)]; 
      Wxx = Wxx + [Wy(m,:) - Wy(1, :); -diff(Wy,1,1)]; 
        
      Fyout = (Nomin1 + gamma*fft2(Wxx))./Denom; 
      I = real(ifft2(Fyout));
      % update the gradient terms with new solution
      Ix = [diff(I, 1, 2), I(:,1) - I(:,n)]; 
      Iy = [diff(I, 1, 1); I(1,:) - I(m,:)]; 
    beta = beta/2;
end 

十四、12.4调用的bilateral_filter

函数格式

function r_img = bilateral_filter(img, sigma_s, sigma, boundary_method, s_size)

在12.4中，输入的参数分别为tv正则化结果和L0正则化结果的差异diff，对应第一个输入参数img、3，对应sigma_s、0.1，对应sigma。

14.1 因为只输入了参数 所以将变量boundary_method置为replicate，代表重复边界。判断输入的img是否为0~255，满足则转化为0~1，方便操作。根据12的操作，img颜色通道和最开始选定的模糊图像的颜色通道相等，假设为彩色，将rgb通道转为lab通道，sigma变为之前的一百倍；如果为灰度，无需附加操作，lab即img。判断s_size是否存在，由于不存在，所以fr设为三倍的sigma_s的上界，即fr=9。分别对img和lab扩充边界（重复边界）边界宽度为fr，得到p_img和p_lab。

if ~exist('boundary_method', 'var')
    boundary_method = 'replicate';
end

if isinteger(img) == 1, img = single(img) / 255; end %将0-255 转为0-1

[h w d] = size(img);

if d == 3
  C = makecform('srgb2lab');
  lab = single( applycform(double(img), C) );
  sigma = sigma * 100;
else
  lab = img;
  sigma = sigma * sqrt(d);
end

if exist('s_size', 'var')
    fr = s_size;
else
    fr = ceil(sigma_s*3);
end

p_img = padarray(img, [fr fr], boundary_method);
p_lab = padarray(lab, [fr fr], boundary_method);

14.2 因为是对四周进行的扩充，那么扩充后图像的真实部分（非添加的边界部分）的两个维度的开始和结束分别是u、b和u、r，代表边界扩充后的p_img的第u行到第b行、第u列到第r列，未扩充的img。初始化r_img、w_sum，spatial_weight为2*fr+1大小（19*19）的高斯模糊核，方差为3均值为0，ss为sigma（0.1）的平方（0.01）.

u = fr+1; b = u+h-1;
l = fr+1; r = l+w-1;

r_img = zeros(h, w, d, 'single');
w_sum = zeros(h, w, 'single');

spatial_weight = fspecial('gaussian', 2*fr+1, sigma_s);
ss = sigma * sigma;

遍历y、x，取值范围为-fx到fx，y+fr+1即从1到2*fr+1，u+y的范围为1到2*fr+1，b+y范围为h到2*fr+h。w_s依次获取了高斯模糊核spatial_weight的值作为权重，n_img和n_lab对p_img和p_lab进行裁剪，依次以（1，1）到（2*fr+1,2*fr+1）为左上角，裁剪一个共h行w列的图像。f_diff为原始lab与扩增后偏移裁剪（以边界区域为左上角）得到的n_lab的差值，f_dist表示三个通道的f_diff的平方和，表示了lab与n_lab在lab三个颜色通道的平方误差的和。计算其高斯分布的概率密度函数（矩阵）w_f，w_t的值为w_f与w_s的点乘。r_img增加n_img与w_t复制三层的结果，w_sum增加w_t。

遍历完后r_img点除复制三层的w_sum得到最后返回的结果r_img，函数bilateral_filter结束。

for y = -fr:fr
  for x = -fr:fr
    
    w_s = spatial_weight(y+fr+1, x+fr+1);
    
    n_img = p_img(u+y:b+y, l+x:r+x, :);
    n_lab = p_lab(u+y:b+y, l+x:r+x, :);
    
    f_diff = lab - n_lab;
    f_dist = sum(f_diff.^2, 3);
    
    w_f = exp(-0.5 * (f_dist / ss));
    
    w_t = w_s .* w_f;
    
    r_img = r_img + n_img .* repmat(w_t, [1 1 d]);
    w_sum = w_sum + w_t;
    
  end
end

r_img = r_img ./ repmat(w_sum, [1 1 d]);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

总体算法的流程：

1.选择去模糊的图像，设定部分参数（包括权值等）。

2.建立图像金字塔，从上到下，下层的模糊核由上层的模糊核重插值得到，下层的图像由原始图像下采样得到。

3.每层金字塔进行五次迭代，每次迭代中，通过正则化得到清晰图像S（（六）或（五）），然后根据模糊图像b的梯度以及清晰图像S的梯度（阈值筛选后的）计算得到模糊核k（（九））。迭代结束，得到估计的模糊核和中间潜在图像S。（这里说的清晰图像并非论文最后结果，只是用来辅助计算模糊核。得到模糊核后再返回来计算清晰图像）

4.此时可以看作非盲复原。对原始的模糊图像y和得到的估计的模糊核k，分别进行TV正则化和L0正则化，得到两种复原结果，对他们的差异进行滤波后，用tv正则化结果减去滤波后的差异得到得到去振铃的结果即最终复原的结果；或者调用whyte的方法进行复原，亦的到最后的清晰图像。