数字图像处理作业3

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)=\sum_p\Bigg((u_p-f_p)^2+\lambda\sum_{q\in \mathcal{N}(p)} w_{p,q}(g)(u_p-u_p)^2 \Bigg)$$

$$w_{p,q}=exp(-\frac {\|g_p-g_q\|}{\sigma_c})$$

此方法涉及到所有点，需要求解一个很大的线性方程组，假设有$H$行(高)，$W$列(宽)，$S=W\times H$，则方程组大小是$S\times S$的。改进方法是每一行每一列分别求解这样的最小二乘问题。对于第$h$行，处理的能量问题是

$$\sum_x\Bigg((u_x^h-f_x^h)^2+\lambda_t\sum_{i\in \mathcal{N}_h(x)} w_{x,i}(g)(u_x^h-u_i^h)^2 \Bigg)$$
对于第$w$列也是同样道理。这个最小二乘问题等价于求解一个方程组
$$(I_h+\lambda_tA_h)u_h=f_h$$
也就是

$$\begin{gather} \left( \begin{matrix} b_0 & c_0 &0&\dots &0\\ \ddots &\ddots & \ddots& 0 &0 \\ 0 &a_x & b_x& c_x &0\\ 0 &0 &\ddots &\ddots & \ddots \\ 0 &\ddots&0 & a_{W-1}& b_{W-1} \end{matrix} \right) \end{gather} \times \left( \begin{array}{ll} u_0^h \\ \vdots \\ u_x^h \\ \vdots \\ u_{W-1}^h\\ \end{array} \right) = \left(\begin{array}{ll}f_0^h \\\vdots \\f_x^h \\\vdots \\f_{W-1}^h\\\end{array} \right)$$
其中
$$\begin{array}{} a_x &=-\lambda_tw_{x,x-1} \\ bx &=1+\lambda_t (w_{x,x-1}+w_{x,x+1}) \\ c_x &=-\lambda_tw_{x,x+1} \end{array}$$
取参数值 $\sigma=25$， $\lambda =900$， $\lambda_t=\frac3 2 \frac {4^{T-t}}{4^{T-1}}\lambda$ 迭代次数 $T=4$ ，由于矩阵的三对角的特点，解方程组使用追赶法，解起来非常快。

程序及实验结果

迭代smooth的主要函数

void smooth::DoSmooth(int t)
{
    MatrixXd fb;
    smooth_parameter = CalcParameter(t);
    MatrixXd matrixW(W,W),matrixH(H,H);
    for (int h = 0; 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 = 0; 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(0);
    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,0, h );
    gp1 = GetPixelValue(image_f,1 , h);
    w1 = CalcWeight(gp, gp1);
    b = 1 + lambda_t*(w1);
    c = -lambda_t*w1;
    T(0, 0) = b;
    T(0, 1) = c;

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

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

右端项赋值

MatrixXd smooth::GetValueRow(Mat img, int h)
{
    MatrixXd img_f(W, 3);
    for (int i = 0; i < W; i++)
    {
        for (int k = 0; k < 3; 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 = 1; i < size; i++)
    {
        A(i, i - 1) = A(i, i - 1) / A(i - 1, i - 1);
        A(i, i) = A(i, i) - A(i - 1, i) * A(i, i - 1);
    }
    //Chasing
    for (size_t i = 1; i < size; i++)
    {
        b.row(i) = b.row(i) - b.row(i - 1) * A(i, i - 1);
    }
    //Returning
    MatrixXd x(size, 3); x.fill(0);
    x.row(size - 1) = b.row(size - 1) / A(size - 1, size - 1);
    for (int i = size - 2; i > -1; i--)
    {
        x.row(i) = (b.row(i) - A(i, i + 1) * x.row(i + 1)) / A(i, i);
    }
    return x;
}

解出的结果对图像进行修正

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

效果

原图：

迭代四次之后：

“guided image filtering” smooth的结果：

可以发现“guided image filtering”边界保留的不如”Fast Global Image Smoothing Based on Weighted Least Squares”

细节增强的结果：

原图：

$\sigma =10$ ， $T=4$，smooth效果

细节增强4倍效果(其中把超过255的取值为255，低于0的取值为0)

程序：

void smooth::Multi_scale_detail()
{
    for (int j = 0; j < H; j++)
    {//先行后列遍历
        for (int i = 0; i < W; i++)
        {
            for (int k = 0; k < 3; k++)
            {
                if (int ( image_f.at<cv::Vec3b>(j, i)[k]) + 5*(image_g.at<cv::Vec3b>(j, i)[k] - image_f.at<cv::Vec3b>(j, i)[k])> 255)image_f.at<cv::Vec3b>(j, i)[k] = 255;
                if (int(image_f.at<cv::Vec3b>(j, i)[k]) + 5*(image_g.at<cv::Vec3b>(j, i)[k] - image_f.at<cv::Vec3b>(j, i)[k]) <0 )image_f.at<cv::Vec3b>(j, i)[k] = 0;
                else
                {
                    image_f.at<cv::Vec3b>(j, i)[k] = image_f.at<cv::Vec3b>(j, i)[k] + 5*(image_g.at<cv::Vec3b>(j, i)[k] - image_f.at<cv::Vec3b>(j, i)[k]);
                }
            }
        }
    }
}

详细程序