典的变分法图像去噪的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++)