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

典的变分法图像去噪的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*(im

age0[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效率高许多的。

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

本文来源：https://www.bwwdw.com/article/pgc3.html