Digital Image Processing(3)

實作文章(Fast Global Image Smoothing Based on Weighted Least Squares)

與“Guided image filtering”作對比

文章的實作

原理及思路

文章的方法是基于前人Weighted Least Squares的方法上進行改進的，達到加速的效果。

假設 g 表示指導圖檔，f為需要smooth的圖檔， u 為smooth之後的圖檔。那麼smooth過程可以描述成求解最小二乘問題

J(u)=∑p((up−fp)2+λ∑q∈N(p)wp,q(g)(up−up)2)

wp,q=exp(−∥gp−gq∥σc)

此方法涉及到所有點，需要求解一個很大的線性方程組，假設有 H 行(高)，W列(寬)， S=W×H ，則方程組大小是 S×S 的。改進方法是每一行每一列分别求解這樣的最小二乘問題。對于第 h 行，處理的能量問題是

∑x((uhx−fhx)2+λt∑i∈Nh(x)wx,i(g)(uhx−uhi)2)

對于第 w 列也是同樣道理。這個最小二乘問題等價于求解一個方程組

(Ih+λtAh)uh=fh

也就是

⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜b0⋱000c0⋱ax0⋱0⋱bx⋱0…0cx⋱aW−1000⋱bW−1⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟×⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜uh0⋮uhx⋮uhW−1⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜fh0⋮fhx⋮fhW−1⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟

其中

axbxcx=−λtwx,x−1=1+λt(wx,x−1+wx,x+1)=−λtwx,x+1

取參數值 σ=25 ， λ=900 ， λt=324T−t4T−1λ 疊代次數 T=4 ，由于矩陣的三對角的特點，解方程組使用追趕法，解起來非常快。

程式及實驗結果

疊代smooth的主要函數

void smooth::DoSmooth(int t)
{
    MatrixXd fb;
    smooth_parameter = CalcParameter(t);
    MatrixXd matrixW(W,W),matrixH(H,H);
    for (int h = ; h < H; h++)
    {
        FitMatrix_x(h, matrixW);
        fb = GetValueRow(image_f, h);
        fb = recursive_solver(matrixW, fb);
        GiveValueRow(image_f, fb, h);
    }
    for (int w = ; w < W; w++)
    {
        FitMatrix_y(w, matrixH);
        fb = GetValueCol(image_f, w);
        fb = recursive_solver(matrixH, fb);
        GiveValueCol(image_f, fb, w);
    }
}

按行裝配矩陣(按列裝配的類似)

void smooth::FitMatrix_x(int h, MatrixXd &T)
{//固定高度，選擇某一行，x進行變化
    T.fill();
    double lambda_t = smooth_parameter;
    double w_1, w1;//w(x,x-1),w(x,x+1)
    VectorXd gp, gp_1, gp1;
    double a, b, c;
    gp = GetPixelValue(image_f,, h );
    gp1 = GetPixelValue(image_f, , h);
    w1 = CalcWeight(gp, gp1);
    b =  + lambda_t*(w1);
    c = -lambda_t*w1;
    T(, ) = b;
    T(, ) = c;

    for (int i=;i<W-1;i++)
    {
        gp = GetPixelValue(image_f,  i,h);
        gp1 = GetPixelValue(image_f,i+, h);
        gp_1 = GetPixelValue(image_f, i-,h );
        w_1 = CalcWeight(gp, gp_1); w1 = CalcWeight(gp, gp1);
        a = -lambda_t*w_1;
        b =  + lambda_t*(w_1 + w1);
        c = -lambda_t*w1;

        T(i, i) = b;
        T(i, i - )= a;
        T(i, i + )= c;
    }
    gp = GetPixelValue(image_f, W-, h);
    gp_1 = GetPixelValue(image_f, W-, h);
    w_1 = CalcWeight(gp, gp_1);
    a = -lambda_t*w_1;
    b =  + lambda_t*(w_1);
    T(W-, W-) = b;
    T(W-, W - ) = a;
}

右端項指派

MatrixXd smooth::GetValueRow(Mat img, int h)
{
    MatrixXd img_f(W, );
    for (int i = ; i < W; i++)
    {
        for (int k = ; k < ; k++)
        {
            img_f(i, k) = double(img.at<cv::Vec3b>(h, i)[k]);
        }
    }
    return img_f;
}

追趕法解方程組

MatrixXd smooth::recursive_solver(MatrixXd A, MatrixXd b)
{
    int size = b.rows();
    //LU decomposition
    for (size_t i = ; i < size; i++)
    {
        A(i, i - ) = A(i, i - ) / A(i - , i - );
        A(i, i) = A(i, i) - A(i - , i) * A(i, i - );
    }
    //Chasing
    for (size_t i = ; i < size; i++)
    {
        b.row(i) = b.row(i) - b.row(i - ) * A(i, i - );
    }
    //Returning
    MatrixXd x(size, ); x.fill();
    x.row(size - ) = b.row(size - ) / A(size - , size - );
    for (int i = size - ; i > -; i--)
    {
        x.row(i) = (b.row(i) - A(i, i + ) * x.row(i + )) / A(i, i);
    }
    return x;
}

解出的結果對圖像進行修正

void smooth::GiveValueRow(Mat &img, MatrixXd v,int h)
{
    for (int i = ; i < W; i++)
    {
        GiveValueToPixel(img, v.row(i).transpose(), i, h);
    }
}

效果

原圖：