全局拉普拉斯平滑之（1）Strucutre extraction from texture via relative total variation及稀疏矩陣求解

最近在研究圖像增強處理過程中，閱讀了關於全局拉普拉斯平滑（global laplacian smoothing），加權最小二乘平滑（weighted least squares --wls）等技術文章，深感此類方法的精妙，並且這種優化思想可以用在許多地方：例如紋理去除，這也是本篇需要重點講的paper:Structure Extraction from Texture via Relative Total Variation; 圖像融合 ,對應文章為：《Poisson Image Editing》；圖像降噪，對應文章為：《Edge-Preserving Decompositions for Multi-Scale Tone and Detail Manipulation》。

上面鏈接為項目主頁地址，利用此類思想的圖像處理文章不止上述方向，後面有時間會針對上面幾篇paper總結說明，開始本篇paper：《Strucutre extraction from texture via relative total variation》的總結，關於這篇文章的說明網上也找到挺多，有些我借用說明，加上自己的一些理解。

文章介紹

很多自然場景和人工藝術品都包含紋理，見下圖，它們有一個共同的特征就是：圖像中有意義的結構圖和紋理融合在一起，可以稱這類圖片為“結構+紋理”圖片。在不去除紋理的前提下，人可以理解這些圖像，其中圖像的結構信息是人類視覺感知的主要數據，而不是紋理。文章的目的是提取出圖像中有意義的結構。

方法介紹：

文章提出了一種基於總變差模型，該模型可以將圖像中的結構和紋理區分開，模型如下：

(1)

I為輸入圖像，p為2D圖像像素的索引，S表示要輸出的結構圖像，下三角表示一階差分運算，其中第二項可以分為x方向和y方向兩個方向上的形式：

(2)

作者的目的是提取各種不同類型的結構，所以並不針對某種特定類型的紋理結構圖像，為了求解上面的模型，作者進行了下面的轉換形式：

(3)

其中：

（4）

q是以p點為中心的一個方形區域內所有像素點的索引，其中g為高斯函數:

（5）

其中圖(a)是一幅包含紋理的圖像。(b)看出紋理和結構像素點都會產生比較大的差分值，對應像素點的亮度比鄰域高；(c)可以看出結構信息對應的值大於紋理值，一種直覺上的解釋為：在一個局部小窗口中，主要結構往往比復雜紋理具有更多相似方向的梯度。

因為自己並沒有去復現文章的求解過程，具體求解過程參加文章，這裏緊要復述一下：

(6)

上述公式的第二行是一個近似計算，結果是二次項和非線性部分的乘積：

（7）

上式中為標準差為的高斯核函數，*為卷積符號。

上面是針對X方向，Y方向同樣的道理。

數學求解時，公式（3）可以轉換成如下矩陣形式：

（8）

其中vs和vi是S和I的兩個列向量，Cx和Cy是向前差分算子，為對角矩陣。因為是求最小值，對矩陣（8）求導，得到下面的線性方程：

（9）

此時可以通過求矩陣的逆運算，或者用預處理共軛梯度法等求解。原文的代碼是公開的，如下所示：

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;  
    end  
    if (~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.5  
            sigma_iter = 0.5;  
        end  
    end  
    S = x;        
end  
  
function [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;  
     
end  
  
function 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‘);    
end  
  
function FBImg = lpfilter(FImg, sigma)       
    FBImg = FImg;  
    for ic = 1:size(FBImg,3)  
        FBImg(:,:,ic) = conv2_sep(FImg(:,:,ic), sigma);  
    end     
end  
  
function 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:ch  
            tin = IN(:,:,ii);  
            [tout, flag] = pcg(A, tin(:),0.1,100, L, L‘);   
            OUT(:,:,ii) = reshape(tout, r, c);  
        end      
    else  
        OUT = IN;  
        for ii=1:ch  
            tin = IN(:,:,ii);  
            tout = A\tin(:);  
            OUT(:,:,ii) = reshape(tout, r, c);  
        end      
    end           
end  

上面的代碼使用matlab實現的，其中求解線性方程時涉及到了稀疏矩陣的方程求解，所以如果要使用c++實現，可能要借助eigen或armadillo等數學庫工具，這篇博客對全變分模型進行了解釋，並給出了c++實現，但是並沒有使用稀疏矩陣求解，所以在處理較大尺寸的圖像時，可能會存在內存過大問題，不過可以作為參考。其代碼如下：

void CImageObj::Total_Variation(int iter, double dt, double epsilon, double lambda)  
{  
    int i, j;  
    int nx = m_width, ny = m_height;  
    double ep2 = epsilon * epsilon;  
  
    double** I_t = NewDoubleMatrix(nx, ny);  
    double** I_tmp = NewDoubleMatrix(nx, ny);  
    for (i = 0; i < ny; i++)  
        for (j = 0; j < nx; j++)  
            I_t[i][j] = I_tmp[i][j] = (double)m_imgData[i][j];  
  
    for (int t = 0; t < iter; t++)  
    {  
        for (i = 0; i < ny; i++)  
        {  
            for (j = 0; j < nx; j++)  
            {  
                int iUp = i - 1, iDown = i + 1;  
                int jLeft = j - 1, jRight = j + 1;    // 邊界處理  
                if (0 == i) iUp = i; if (ny - 1 == i) iDown = i;  
                if (0 == j) jLeft = j; if (nx - 1 == j) jRight = j;  
  
                double tmp_x = (I_t[i][jRight] - I_t[i][jLeft]) / 2.0;  
                double tmp_y = (I_t[iDown][j] - I_t[iUp][j]) / 2.0;  
                double tmp_xx = I_t[i][jRight] + I_t[i][jLeft] - 2 * I_t[i][j];  
                double tmp_yy = I_t[iDown][j] + I_t[iUp][j] - 2 * I_t[i][j];  
                double tmp_xy = (I_t[iDown][jRight] + I_t[iUp][jLeft] - I_t[iUp][jRight] - I_t[iDown][jLeft]) / 4.0;  
                double tmp_num = tmp_yy * (tmp_x * tmp_x + ep2) + tmp_xx * (tmp_y * tmp_y + ep2) - 2 * tmp_x * tmp_y * tmp_xy;  
                double tmp_den = pow(tmp_x * tmp_x + tmp_y * tmp_y + ep2, 1.5);  
  
                I_tmp[i][j] += dt*(tmp_num / tmp_den + lambda*(m_imgData[i][j] - I_t[i][j]));  
            }  
        }  // 一次叠代  
  
        for (i = 0; i < ny; i++)  
            for (j = 0; j < nx; j++)  
            {  
                I_t[i][j] = I_tmp[i][j];  
            }  
  
    } // 叠代結束  
  
    // 給圖像賦值  
    for (i = 0; i < ny; i++)  
        for (j = 0; j < nx; j++)  
        {  
            double tmp = I_t[i][j];  
            tmp = max(0, min(tmp, 255));  
            m_imgData[i][j] = (unsigned char)tmp;  
        }  
  
    DeleteDoubleMatrix(I_t, nx, ny);  
    DeleteDoubleMatrix(I_tmp, nx, ny);  
}  

後面我會嘗試將代碼轉換為c++實現，主要是稀疏矩陣構造及求解那一部分，附一段我使用eigen復現的matlab中spdiags函數：

Eigen::SparseMatrix<double> spdiags(const MatrixXd& B,const Eigen::Matrix<int, 1, 1>& d, int m, int n)
{
	Eigen::SparseMatrix<double> A(m, n);
	std::vector<Triplet < double >> triplets;
 
	for (int k = 0; k < d.size(); k++)
	{
		int i_min = std::max(0, -d(k));
		int i_max = std::min(m - 1, n - d(k) - 1);
		int B_idx_start = m >= n ? d(k) : 0;
		for (int i = i_min; i <= i_max; i++){
			if (d(k)>0)
				triplets.emplace_back(i, i+d(k), B(B_idx_start + i, k));
			else
				triplets.emplace_back(i, i-i_min, B(B_idx_start + i, k));
		}
			
	}
	A.setFromTriplets(triplets.begin(), triplets.end());
	return A;
}

int main()
{
 
//---------------------------------------
 
	Matrix<int, 1, 1> d1;
	MatrixXd d0(5,1);
 
	d0(0, 0) = 10; d0(1, 0) = 20;
	d0(2, 0) = 30; d0(3, 0) = 40;
	d0(4, 0) = 50;
        d1(0)=0;
 
	Eigen::SparseMatrix<double> Diag= spdiags(d02, d1, 5, 5);
	std::cout << Diag<< std::endl;
}

參考：

https://wenku.baidu.com/view/9caf94767375a417866f8ff1.html

http://blog.sina.com.cn/s/blog_4bdb170b0101ovi8.html

https://www.cnblogs.com/Imageshop/p/3365517.html

http://www.cnblogs.com/CCBB/archive/2010/12/29/1920884.html

https://github.com/cran/tvR/blob/7e8f900b99cdae8e65774b166423dbd39a07f6a5/src/RcppCollection_Image.cpp

轉載：https://blog.csdn.net/piaoxuezhong/article/details/79050046

轉載：全局拉普拉斯平滑之（1）Strucutre extraction from texture via relative total variation及稀疏矩陣求解