L0 gradient minimization算法笔记

L0 gradient minimization算法

论文及源代码link
[侵删]

数学基础

L0 gradient minimization算法的的目标函数：
$$\underset{S}{\min }\left\{\sum _p(S_p-I_p{)}^2+\lambda \cdot C(S)\right\}$$

$$C(S)=\#\left\{\begin{array}{ccc}p& |& |{\partial }_x{S}_p|+|{\partial }_y{S}_p|\ne 0\end{array}\right\}$$
由于原始函数一项为像素差异，另一项为全局不连续的函数，所以无论是梯度法还是离散优化方法都很难计算。所以采用了半二次分裂的方法，所包含的的子问题都是具有封闭形式的解。作者通过引入辅助变量$h_p,v_p$,对原始目标函数进行重塑，得到：
$$\underset{S.h,v}{\min }\{\sum _p(S_p-I_p{)}^2+\lambda C(h,v)+\beta (({\partial }_x{S}_p-h_p{)}^2+({\partial }_y{S}_p-v_p{)}^2)\}(1)$$
其中的$C(h,v)=\#\left\{p\right||h_p|+|v_p|\ne 0\}$用来约束梯度。

转为两个子问题：

子问题一：

解决子问题1为：
$$\left\{\sum _p(S_p-I_p{)}^2+\beta (({\partial }_x{S}_p-h_p{)}^2+({\partial }_y{S}_p-v_p{)}^2)\right\}(2)$$
公式2为二次函数，可以找到全局最优点。使用FFT进行加速，则转为计算：
$$S={\mathcal{F}}^{-1}\left(\begin{array}{c}\mathcal{F}(I)+\beta (\mathcal{F}({\partial }_x{)}^{\ast }\mathcal{F}(h)+\mathcal{F}({\partial }_y{)}^{\ast }\mathcal{F}(v))\\ \mathcal{F}(1)+\beta (\mathcal{F}({\partial }_x{)}^{\ast }\mathcal{F}({\partial }_x)+\mathcal{F}({\partial }_y{)}^{\ast }\mathcal{F}({\partial }_y))\end{array}\right)(3)$$

子问题二：

目标函数为：
$$\underset{h,v}{\min }\left\{\sum _p({\partial }_x{S}_p-h_p{)}^2+({\partial }_y{S}_p-v_p{)}^2)+\frac{\lambda }{\beta }C(h,v)\right\}(4)$$
其中的$C(h,v)$返回的是$|h|+|v|$的非零元素。
这个明显复杂的子问题实际上快速解决，因为能(4)可以进行空间分解，其中每个元素 $h_p$和$v_p$是单独估计的。拆分后变为：
$$\sum _p\underset{h_p,v_p}{\min }\left\{{\left(h_p-{\partial }_x{S}_p\right)}^2+{\left(v_p-{\partial }_y{S}_p\right)}^2+\frac{\lambda }{\beta }H(\left|h_p\right|+\left|v_p\right|)\right\}(5)$$
其中的$H(|h|+|v|)$是一个二值化函数，如果$|h_p|+|v_p|\ne 0$则返回1，否则为0。对于公式5，构建
$$E_p=\left\{(h_p-{\partial }_x{S}_p{)}^2+(v_p-{\partial }_y{S}_p{)}^2+\frac{\lambda }{\beta }H(|h_p|+|v_p|)\right\}(6)$$
那么取得最小值的$E_p^{\ast }$需要满足条件：
$$(h_p,v_p)=\{\begin{array}{ll}(0,0)& ({\partial }_x{S}_p{)}^2+({\partial }_y{S}_p{)}^2\le \lambda /\beta \\ ({\partial }_x{S}_p,{\partial }_y{S}_p)& \text{otherwise}\end{array}(7)$$

算法流程图

通过求解h,v的值，不断更新S以及$\beta$，最终求得原始公式的最优解。

代码解析

官方matlab代码

function S = L0Smoothing(Im, lambda, kappa)
%L0Smooth - Image Smoothing via L0 Gradient Minimization
%   S = L0Smooth(Im, lambda, kappa) performs L0 graidient smoothing of input
%   image Im, with smoothness weight lambda and rate kappa.
%
%   Paras: 
%   @Im    : Input UINT8 image, both grayscale and color images are acceptable.
%   @lambda: Smoothing parameter controlling the degree of smooth. (See [1]) 
%            Typically it is within the range [1e-3, 1e-1], 2e-2 by default.
%   @kappa : Parameter that controls the rate. (See [1])
%            Small kappa results in more iteratioins and with sharper edges.   
%            We select kappa in (1, 2].    
%            kappa = 2 is suggested for natural images.  
%
%   Example
%   ==========
%   Im  = imread('pflower.jpg');
%   S  = L0Smooth(Im); % Default Parameters (lambda = 2e-2, kappa = 2)
%   figure, imshow(Im), figure, imshow(S);


if ~exist('kappa','var')
    kappa = 2.0;
end
if ~exist('lambda','var')
    lambda = 2e-2;
end
S = im2double(Im);
betamax = 1e5;
fx = [1, -1];
fy = [1; -1];
[N,M,D] = size(Im);
sizeI2D = [N,M];

% psf2otf：点扩散函数转调制传递函数，这里认为是微分的FFT
otfFx = psf2otf(fx,sizeI2D);
otfFy = psf2otf(fy,sizeI2D);
Normin1 = fft2(S);

% 公式3中分母部分：微分的FFT的平方
Denormin2 = abs(otfFx).^2 + abs(otfFy ).^2;
if D>1
    Denormin2 = repmat(Denormin2,[1,1,D]);
end
beta = 2*lambda;
while beta < betamax
   
    % 公式3分母
    Denormin   = 1 + beta*Denormin2;

    % h-v subproblem
    h = [diff(S,1,2), S(:,1,:) - S(:,end,:)];
    v = [diff(S,1,1); S(1,:,:) - S(end,:,:)];
    % 公式7,求h-v问题的最小化
    if D==1
        t = (h.^2+v.^2)<lambda/beta;
    else
        t = sum((h.^2+v.^2),3)<lambda/beta;
        t = repmat(t,[1,1,D]);
    end
    h(t)=0; v(t)=0;
    
    % S subproblem
    % 计算公式3的分子：分别对h和v进行dx和dy计算
    % 等价于FFT域中dx,dy乘h,v
    Normin2 = [h(:,end,:) - h(:, 1,:), -diff(h,1,2)];
    Normin2 = Normin2 + [v(end,:,:) - v(1, :,:); -diff(v,1,1)];
    FS = (Normin1 + beta*fft2(Normin2))./Denormin;
    %计算得到当前的解
    S = real(ifft2(FS));
    % 不断增大beta，使得整体公式逼近原始公式。
    beta = beta*kappa;
    fprintf('.');
end
fprintf('\n');
end