数字图像处理作业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)=∑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 = 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”

细节增强的结果：

原图：

σ=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]);                }            }        }    }}

详细程序