经典的变分法图像去噪的C++实现

由于这学期的图像处理课程的大作业需要写一个图像处理程序，不能使用古典的线性滤波，或者基于频域（小波）或者基于统计之类的方法。只能用老师讲过的一些方法，诸如变分，PDE，微分几何等。。感觉上简单的变分法稍微要好实现一些，就打算基于最早的TV图像去噪模型，做一个VC的实现。但是找遍了网上也没有TV去噪的C++源码，与之只好自己动手写了。

关于变分法和泛函分析的一些基础原理今天就先不多说了，TV图像去噪经典论文：《Nonlinear Total Variation based noise removal algorithms》Google上可以搜得到。

关于Matlab的程序实现，有一个经典的主页： http://visl.technion.ac.il/~gilboa/PDE-filt/tv_denoising.html

下面是一个Matlab代码实现：复制到记事本用matlab打开就可以运行，要注意图像的名称和路径要对应。如果只是想学学算法思路或者看看处理效果的话，只需要Matlab的代码就行了。

function J=tv(I,iter,dt,ep,lam,I0,C)
%% Private function: tv (by Guy Gilboa).
%% Total Variation denoising.
%% Example: J=tv(I,iter,dt,ep,lam,I0)
%% Input: I    - image (double array gray level 1-256),
%%        iter - num of iterations,
%%        dt   - time step [0.2],
%%        ep   - epsilon (of gradient regularization) [1],
%%        lam  - fidelity term lambda [0],
%%        I0   - input (noisy) image [I0=I]
%%       (default values are in [])
%% Output: evolved image

clc
clear
I=imread('grids.bmp'); % load image
I = double(I);

if ~exist('ep')
   ep=1;
end
if ~exist('dt')
   dt=ep/5;  % dt below the CFL bound
end
if ~exist('lam')
   lam=0;
end
if ~exist('I0')
    I0=I;
end
if ~exist('C')
    C=0;
end
[ny,nx]=size(I); ep2=ep^2;

% params
iter=80; 

for i=1:iter,  %% do iterations
   % estimate derivatives
   I_x = (I(:,[2:nx nx])-I(:,[1 1:nx-1]))/2;
    I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2;
    I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
    I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
    Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]);
    Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
    I_xy = (Dp-Dm)/4;
   % compute flow
   Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2);
   Den = (ep2+I_x.^2+I_y.^2).^(3/2);
   I_t = Num./Den + lam.*(I0-I+C);
   I=I+dt*I_t;  %% evolve image by dt
end % for i

%% return image
%J=I*Imean/mean(mean(I)); % normalize to original mean
J=I;

figure(1); imshow(uint8(I0)); title('Noisy image');
% denoise image by using tv for some iterations
figure(2); imshow(uint8(J)); title('Denoised image');

另外我在我的图像处理框架程序里实现了这个最经典版本的TV去噪算法，核心代码如下：

//TV去噪函数
bool MyCxImage::TVDenoising(int iter /* = 80 */)
{
    if(my_image == NULL) return false;
    if(!my_image->IsValid()) return false;
    //算法目前不支持彩色图像，所以对于彩图，先要转换成灰度图。
    if(!my_image->IsGrayScale())
    {
        my_image->GrayScale();
        //return false;
    }

    //基本参数，这里由于设置矩阵C为0矩阵，不参与运算，所以就忽略之
    int ep = 1, nx = width, ny = height;
    double dt = (double)ep/5.0f, lam = 0.0;
    int ep2 = ep*ep;

    double** image = newDoubleMatrix(nx, ny);
    double** image0 = newDoubleMatrix(nx, ny);
    //注意一点是CxImage里面图像存储的坐标原点是左下角，Matlab里面图像时左上角原点
    for (int i = 0; i < ny; i++)
    {
        for (int j = 0; j < nx; j++)
        {
            image0[i][j] = image[i][j]  = my_image->GetPixelIndex(j, ny-i-1);
        }
    }

    double** image_x = newDoubleMatrix(nx, ny);   //I_x = ( I(:,[2:nx nx]) - I(:,[1 1:nx-1]))/2;
    double** image_xx = newDoubleMatrix(nx, ny);   //I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
    double** image_y = newDoubleMatrix(nx, ny);   //I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2;
    double** image_yy = newDoubleMatrix(nx, ny);   //I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
    double** image_tmp1 = newDoubleMatrix(nx, ny);
    double** image_tmp2 = newDoubleMatrix(nx, ny);

    double** image_dp = newDoubleMatrix(nx, ny);   //Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1

    double** image_dm = newDoubleMatrix(nx, ny);   //Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);

    double** image_xy = newDoubleMatrix(nx, ny);   //I_xy = (Dp-Dm)/4;

    double** image_num = newDoubleMatrix(nx, ny);   //Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2);
    double** image_den = newDoubleMatrix(nx, ny);   //Den = (ep2+I_x.^2+I_y.^2).^(3/2);

    //////////////////////////////////////////////////////////////////////////
    //对image进行迭代iter次
    iter = 80;
    for (int t = 1; t <= iter; t++)
    {
        //进度条
        my_image->SetProgress((long)100*t/iter);
        if (my_image->GetEscape())
            break;
        //////////////////////////////////////////////////////////////////////////
        //计算I(:,[2:nx nx])和I(:,[1 1:nx-1])
        //公共部分2到nx-1列
        for (int i = 0; i < ny; i++)
        {
            for (int j = 0; j < nx-1; j++)
            {
                image_tmp1[i][j] = image[i][j+1];
                image_tmp2[i][j+1] = image[i][j];
            }
        }
        for (int i = 0; i < ny; i++)
        {
            image_tmp1[i][nx-1] = image[i][nx-1];
            image_tmp2[i][0] = image[i][0];
        }

        //计算I_x, I_xx
        // I_x = ( I(:,[2:nx nx]) - I(:,[1 1:nx-1]))/2
        //I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
        for (int i = 0; i < ny; i++)
        {
            for (int j = 0; j < nx; j++)
            {
                image_x[i][j] = (image_tmp1[i][j] - image_tmp2[i][j])/2;
                image_xx[i][j] = (image_tmp1[i][j] + image_tmp2[i][j]) - 2*image[i][j];
            }
        }

        //////////////////////////////////////////////////////////////////////////
        //计算I([2:ny ny],:)和I([1 1:ny-1],:)
        //公共部分2到ny-1行
        for (int i = 0; i < ny-1; i++)
        {
            for (int j = 0; j < nx; j++)
            {
                image_tmp1[i][j] = image[i+1][j];
                image_tmp2[i+1][j] = image[i][j];
            }
        }
        for (int j = 0; j < nx; j++)
        {
            image_tmp1[ny-1][j] = image[ny-1][j];
            image_tmp2[0][j] = image[0][j];
        }
        //计算I_xx, I_yy
        // I_y = I([2:ny ny],:)-I([1 1:ny-1],:)
        //I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
        for (int i = 0; i < ny; i++)
        {
            for (int j = 0; j < nx; j++)
            {
                image_y[i][j] = (image_tmp1[i][j] - image_tmp2[i][j])/2;
                image_yy[i][j] = (image_tmp1[i][j] + image_tmp2[i][j]) - 2*image[i][j];
            }
        }

        //////////////////////////////////////////////////////////////////////////
        //计算I([2:ny ny],[2:nx nx])和I([1 1:ny-1],[1 1:nx-1])
        //公共部分分别是矩阵右下角，左上角的ny-1行和nx-1列
        for (int i = 0; i < ny-1; i++)
        {
            for (int j = 0; j < nx-1; j++)
            {
                image_tmp1[i][j] = image[i+1][j+1];
                image_tmp2[i+1][j+1] = image[i][j];
            }
        }
        for (int i = 0; i < ny-1; i++)
        {
            image_tmp1[i][nx-1] = image[i+1][nx-1];
            image_tmp2[i+1][0] = image[i][0];
        }
        for (int j = 0; j < nx-1; j++)
        {
            image_tmp1[ny-1][j] = image[ny-1][j+1];
            image_tmp2[0][j+1] = image[0][j];
        }
        image_tmp1[ny-1][nx-1] = image[ny-1][nx-1];
        image_tmp2[0][0] = image[0][0];
        //计算Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]);
        for (int i = 0; i < ny; i++)
        {
            for (int j = 0; j < nx; j++)
            {
                image_dp[i][j] = image_tmp1[i][j] + image_tmp2[i][j];
            }
        }

        //////////////////////////////////////////////////////////////////////////
        //计算I([1 1:ny-1],[2:nx nx])和I([2:ny ny],[1 1:nx-1])
        //公共部分分别是矩阵左下角，右上角的ny-1行和nx-1列
        for (int i = 0; i < ny-1; i++)
        {
            for (int j = 0; j < nx-1; j++)
            {
                image_tmp1[i+1][j] = image[i][j+1];
                image_tmp2[i][j+1] = image[i+1][j];
            }
        }
        for (int i = 0; i < ny-1; i++)
        {
            image_tmp1[i+1][nx-1] = image[i][nx-1];
            image_tmp2[i][0] = image[i+1][0];
        }
        for (int j = 0; j < nx-1; j++)
        {
            image_tmp1[0][j] = image[0][j+1];
            image_tmp2[ny-1][j+1] = image[ny-1][j];
        }
        image_tmp1[0][nx-1] = image[0][nx-1];
        image_tmp2[ny-1][0] = image[ny-1][0];

        //计算Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
        for (int i = 0; i < ny; i++)
        {
            for (int j = 0; j < nx; j++)
            {
                image_dm[i][j] = image_tmp1[i][j] + image_tmp2[i][j];
            }
        }

        //////////////////////////////////////////////////////////////////////////
        //计算I_xy = (Dp-Dm)/4;
        for (int i = 0; i < ny; i++)
        {
            for (int j = 0; j < nx; j++)
            {
                image_xy[i][j] = (image_dp[i][j] - image_dm[i][j])/4;
            }
        }

        //////////////////////////////////////////////////////////////////////////
        //计算过程：

        //计算Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2) 和 Den = (ep2+I_x.^2+I_y.^2).^(3/2);
        for (int i = 0; i < ny; i++)
        {
            for (int j = 0; j < nx; j++)
            {
                image_num[i][j] = image_xx[i][j]*(image_y[i][j]*image_y[i][j] + ep2) 
                    - 2*image_x[i][j]*image_y[i][j]*image_xy[i][j] + image_yy[i][j]*(image_x[i][j]*image_x[i][j] + ep2);

                image_den[i][j] = pow((image_x[i][j]*image_x[i][j] + image_y[i][j]*image_y[i][j] + ep2), 1.5);
            }
        }

        //计算I: I_t = Num./Den + lam.*(I0-I+C); I=I+dt*I_t;  %% evolve image by dt
        for (int i = 0; i < ny; i++)
        {
            for (int j = 0; j < nx; j++)
            {
                image[i][j] += dt*(image_num[i][j]/image_den[i][j] + lam*(image0[i][j] - image[i][j]));
            }
        }
    }
    //迭代结束
    
    //////////////////////////////////////////////////////////////////////////
    //赋值图像
    BYTE tmp;
    for (int i = 0; i < ny; i++)
    {
        for (int j = 0; j < nx; j++)
        {
            tmp = (BYTE)image[i][j];
            tmp = max(0, min(tmp, 255));
            my_image->SetPixelIndex(j, ny-i-1, tmp);
        }
    }

    //////////////////////////////////////////////////////////////////////////
    //删除内存
    deleteDoubleMatrix(image_x, nx, ny);
    deleteDoubleMatrix(image_y, nx, ny);
    deleteDoubleMatrix(image_xx, nx, ny);
    deleteDoubleMatrix(image_yy, nx, ny);
    deleteDoubleMatrix(image_tmp1, nx, ny);
    deleteDoubleMatrix(image_tmp2, nx, ny);
    deleteDoubleMatrix(image_dp, nx, ny);
    deleteDoubleMatrix(image_dm, nx, ny);
    deleteDoubleMatrix(image_xy, nx, ny);
    deleteDoubleMatrix(image_num, nx, ny);
    deleteDoubleMatrix(image_den, nx, ny);
    deleteDoubleMatrix(image0, nx, ny);
    deleteDoubleMatrix(image, nx, ny);

    return true;
}
//////////////////////////////////////////////////////////////////////////
//开辟二维数组函数
double** MyCxImage::newDoubleMatrix(int nx, int ny)
{
    double** matrix = new double*[ny];

    for(int i = 0; i < ny; i++)
    {
        matrix[i] = new double[nx];
    }
    if(!matrix)
        return NULL;
    return
        matrix;
}
//清除二维数组内存函数
bool MyCxImage::deleteDoubleMatrix(double** matrix, int nx, int ny)
{
    if (!matrix)
    {
        return true;
    }
    for (int i = 0; i < ny; i++)
    {
        if (matrix[i])
        {
            delete[] matrix[i];
        }
    }
    delete[] matrix;

    return true;
}
//////////////////////////////////////////////////////////////////////////

这个代码单独显然是无法运行的，因为还要涉及底层的图像处理的类库，图像的读取显示我用了CxIamge类，而程序界面我是用的MFC的框架。不过代码基本一直都是在做矩阵运算，如果要是能有一个比较好的矩阵运算类库的话，代码会简介许多，效率也会高一些。总体上C++代码还是要比Matlab效率高许多的。

关于变分法的算法原理和基本思想，我这两天再读一些论文在做总结。。

Email：lichao@icst.pku.edu.cn