图像主结构的提取方法

介绍

本节主要由学习论文Structure Extraction from Texture via Relative Total Variation一文中的内容而来。我们知道，很多自然场景和人工艺术品都包含纹理(如，墙上，火车和地铁表面上的涂鸦和图案。像地毯，毛衣，和其他一些精美的工艺品包含格式各样的几何图案)。在人类历史中，马赛克被视为一种艺术形式，它可以表示人和动物这类复杂的场景，并可以用石头，玻璃，陶瓷和其他一些材料模仿油画。当用Bing或Google收索这些图像的时候，你可以很快的找到成千上万类似图片。

下图所示展示了一些代表不同形式场景的实例。他们有一个共同的特征：图像中有意义的结构和纹理单元融合在一起。所以我们通常称这类图片为“结构+纹理”图片。一个很有意思的现象：在不去除纹理的前提下，人类的视觉感知系统完全有能力理解这些图像。从心里学角度分析，图像的整体结构特才是人类视觉感知的主要数据，而不是那些个体细节(纹理)。因此从图像中提取那些有意义的结构数据是一项具有意义的工作，同时对于计算机来说也是非常有挑战性的。

理论描述

在这篇论文一文中，提出了一种基于总变差形式新模型，该模型可以有效的分解图像中的结构信息和纹理，并且无需特别指定纹理是否规则或者对称。换言之，该方法具有一般性和随意性，它适用于非统一的或各向异性的纹理。下面介绍一下该方法。

上式为总变差模型, I代表输入图像，p代表2D图像像素的索引，S代表输出结构图像。其中可以写成如下各向异性的形式：

改进的模型如下

其中：

q为以p点为中心的一个正方形区域内所有的像素点的索引，g为高斯核函数:

我们看一幅包含纹理的图像。如下图所示：

其中图(a)所示是一幅包含纹理的图像。(b)则反映了纹理和结构像素点都会产生比较大的D(D值大反应在图像中也就是对应像素点的亮度高)；(c)可以看出结构部分中的L(L值大反应在图像中也就是对应像素点的亮度高)值大于纹理部分的L值，造成这种现象的一种直觉上的解释为：在包含在一个局部的小窗口中主要结构往往产生比包含在另一个窗口的复杂纹理具有更多相似方向的梯度。(d)为 (也就是(3)式中的规则项)的映射，而这个规则项让主要结构部分更加突出。(e)则是(a)用式(3)去除纹理后的结构图像。

现在，我们介绍怎么样求解公式(3): 先讨论X方向，Y方向的计算类似。

由于引入了小上述公式的第二行实际是一个近似计算。同时重新构造二次项和非线性部分。它们可以分别表示为如下形式：

上式中为标准差为的高斯核函数，*为卷积符号。Y方向上的求解类似：

其中

有了这些算子可以将公式（3）写成如下矩阵的形式：

其中是vs和vi代表S和I的两个列矢量。Cx和Cy是向前差分梯度算子的Toeplitz Matrices.都为对角矩阵，他们对线线上的值分别为：。然后，对矩阵（4）求导得到如下线性方程：

上公式可以至直接求矩阵的逆运算，或者用预处理共轭梯度法来求解。

编码实现

function S = tsmooth(I,lambda,sigma,sharpness,maxIter)if (~exist('lambda','var'))lambda=0.01;end   if (~exist('sigma','var'))sigma=3.0;end if (~exist('sharpness','var'))sharpness = 0.02;endif (~exist('maxIter','var'))maxIter=4;end    I = im2double(I);x = I;sigma_iter = sigma;lambda = lambda/2.0;dec=2.0;for iter = 1:maxIter[wx, wy] = computeTextureWeights(x, sigma_iter, sharpness);x = solveLinearEquation(I, wx, wy, lambda);sigma_iter = sigma_iter/dec;if sigma_iter < 0.5sigma_iter = 0.5;endendS = x;
endfunction [retx, rety] = computeTextureWeights(fin, sigma,sharpness)fx = diff(fin,1,2);fx = padarray(fx, [0 1 0], 'post');fy = diff(fin,1,1);fy = padarray(fy, [1 0 0], 'post');vareps_s = sharpness;vareps = 0.001;wto = max(sum(sqrt(fx.^2+fy.^2),3)/size(fin,3),vareps_s).^(-1); fbin = lpfilter(fin, sigma);gfx = diff(fbin,1,2);gfx = padarray(gfx, [0 1], 'post');gfy = diff(fbin,1,1);gfy = padarray(gfy, [1 0], 'post');     wtbx = max(sum(abs(gfx),3)/size(fin,3),vareps).^(-1); wtby = max(sum(abs(gfy),3)/size(fin,3),vareps).^(-1);   retx = wtbx.*wto;rety = wtby.*wto;retx(:,end) = 0;rety(end,:) = 0;endfunction ret = conv2_sep(im, sigma)ksize = bitor(round(5*sigma),1);g = fspecial('gaussian', [1,ksize], sigma); ret = conv2(im,g,'same');ret = conv2(ret,g','same');
endfunction FBImg = lpfilter(FImg, sigma)     FBImg = FImg;for ic = 1:size(FBImg,3)FBImg(:,:,ic) = conv2_sep(FImg(:,:,ic), sigma);end
endfunction OUT = solveLinearEquation(IN, wx, wy, lambda)[r,c,ch] = size(IN);k = r*c;dx = -lambda*wx(:);dy = -lambda*wy(:);B(:,1) = dx;B(:,2) = dy;d = [-r,-1];A = spdiags(B,d,k,k);e = dx;w = padarray(dx, r, 'pre'); w = w(1:end-r);s = dy;n = padarray(dy, 1, 'pre'); n = n(1:end-1);D = 1-(e+w+s+n);A = A + A' + spdiags(D, 0, k, k); if exist('ichol','builtin')L = ichol(A,struct('michol','on'));    OUT = IN;for ii=1:chtin = IN(:,:,ii);[tout, flag] = pcg(A, tin(:),0.1,100, L, L'); OUT(:,:,ii) = reshape(tout, r, c);end    elseOUT = IN;for ii=1:chtin = IN(:,:,ii);tout = A\tin(:);OUT(:,:,ii) = reshape(tout, r, c);end    endend

实验结果与分析

实验中ε和ε_s是2个小正值用来避免分母出现0的情况。其中ε固定为0.001. 而ε_s稍微大点会帮助保持光滑变化的结构部分，通常设定为0.02。公式(5)中的λ是一个不可或缺的权重它用来控制图像的光滑程度，但是仅仅调节它不会使纹理分离太多。而增加λ也会造成图像的模糊并且纹理反而保留下来。一般λ选取为0.01到0.03之间。图1表示不同迭代步数显示的结果，实验发现该算法3-5步就可以达到收敛状态。

空间尺度参数σ控制了公式(4)中窗口的大小，它的选取取决于纹理的尺度大小并且在结构纹理分离过程中至关重要，经验的选取σ为0到8之间，图2说明了增强σ可以很好地抑制纹理。并且实验中发现在每一次迭代时成倍的减小,可以起到锐化边缘的效果，同时不会减弱纹理去除的能力.

当一幅物体的表面包含多重纹理形式或者可以看成非正面方向，纹理单元就可以认为是不同变化尺度的。图2和图3就是这样一类图像。因为在该算法中小于相对于尺度参数那些纹理都得到了有效的惩罚，所以文章的提出模型可以很好的处理这种类型的图像。当然，如果远处的结构和近处的纹理相似，他们也都会被去除。因为该算法依赖于局部的数据，所以我们不需要认为局部的梯度是各项同性的。只要在一个局部窗口中方向相反的梯度相互抵消，该方法就能生效，而不管梯度模式是不是各项同性的或是异性的。

图像矢量化就是把一个像素图像转化为一个矢量图。矢量图可以任意的放大和缩小而不会丢失细节部分，然而大多数矢量化的方法都不能表示好的细节部分。由于复杂的纹理模式和局部像素点的不断震荡的普遍存性，矢量化“结构+纹理”这一类图像变的更加困难。观察如下图所示：

图8(c)和(e)，经典的矢量化软件Vector Magic也不能很好地完成矢量化。在本文中，我们开始先分解纹理和结构，分解的结构图为图8(b)，然后矢量化就可以很好地运用了。在矢量化的过程中，结构图像(b)直接被放大。于此同时，纹理图像可以用双线性插值作为一个位图重新被放大。最后合成这两层图像获得图8(f)。相对于传统的方法，该矢量化算法可以产生更好地效果：不丢失边缘和细节信息。

本文的算法还可以用于边缘提取。如下图所示，展示了一个例子，该幅图像中包含很明显的前景和背景的纹理，这往往导致边缘提取的失败。图9(b)和(c)使用不同参数的额Canny边缘检测提取的边缘。很明显这样的边缘是不令人满意的。该方法可以先获得好的结构图像(d)，然后再检测该结构图像的边缘得到(e)。图(6)说明了一样的道理。

由于源纹理和目标纹理的不兼容性，有时涂鸦图像，油画，和素描不能直接运用到图像融合中。图11和图12就是一个很好的例子。

直接将图11(a)和图12(a)融入目标场景中得到图11(c)和图12(e)，不难发现融合的图像很不自然。然而将纹理分离后的图11(b)和图12(b)融入相同场景得到的图11(d)和图12(d)却很自然。改组实验说明了本文算法可以很好运用到图像融合中。

补充说明

一般来说，我们知道Matlab只适合于科研。如果想做成产品，需要用其他的比如C/C++之类的实现该算法，对于这点，这个程序可能有点困难，主要的困难在于其中的解线性方程组。因为在本算法中，方程组的系数矩阵式一个很大的稀疏矩阵，有多大呢，比如如果图像时600x600的，那么这个矩阵的大小就是500000x500000，因此，直接的实现对于内存等方面肯定不现实，必须研究稀疏矩阵的存储方式，另外，解方程所用的共轭梯度法也不是一下子就能实现的，正在研究......同时，关于C/C++版本程序也还在研究中......(^_^)。但是，此算法由于占用太多内存，那Android手机上更不用说了。不过，如果有时间还是可以好好优化一下。呵呵......

参考文献

[1] Li Xu, Qiong Yan, Yang Xia, Jiaya Jia, "Structure Extraction from Texture via Relative Total Variation", ACM Transactions on Graphics (TOG), 31(6), 139, Proceedings of ACM SIGGRAPH Asia 2012.

关于Image Engineering & Computer Vision的更多讨论与交流，敬请关注本博和新浪微博songzi_tea.