TV：Total Variation
BTV：Bilateral Total Variation
Osher等在1992 年提出了总变分（TV）超分辨率重建方法，该方法能够有效
地去除噪声和消除模糊，但是在强噪声情况下，图像的平滑区域会产生阶梯效应，
同时图像的纹理信息也不能很好地保留。
Farsiu等在2004 年提出了双边总变分（BTV）正则化方法，该方法不仅考虑了周围像素与中心像素的几何距离，同时也考虑了灰度相似性，这使得BTV 算法性能相对于TV 算法有了很大的提高。所以在几种正则化方法中， BTV 算法具有较好的重建效果，该方法不仅能够去除噪声，而且又能较好的保持图像的边缘信息。
正则化方法：
TV：很好地保持图像边缘信息。
BTV：对比原图像和将其平移整数个像素后的图像，再对两者的差求加权平均和。
p为选取窗口的半径，矩阵Sx和Sy分别为将图像Z沿x和y方向分别平移l和m个像素，alpha为尺度加权系数，其取值范围为0-1。
确定lamda值科学的方法为：在图像平滑区域，值应较大；在边缘和纹理区域，值应较小。
（在实际操作中改进，可以针对lamda值设定参数，合理lamda模型化，引入已知的先验信息。）
TV去噪模型：Rudin、Osher and Fatemi提出。
TV图像去噪模型的成功之处就在于利用了自然图像内在的正则性，易于从噪声图像的解中反映真实图像的几何正则性，比如边界的平滑性。
最小化全变分来去噪：
约束条件：
等价于最小化下式：
导出的欧拉-拉格朗日方程：
数值实现：
代码：
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 -),
%% iter - num of iterations,
%% dt - time step [0.2],
%% ep - epsilon (of gradient regularization) [],
%% lam - fidelity term lambda [],
%% I0 - input (noisy) image [I0=I]
%% (default values are in [])
%% Output: evolved image if ~exist('ep')
ep=;
end
if ~exist('dt')
dt=ep/; % dt below the CFL bound
end
if ~exist('lam')
lam=;
end
if ~exist('I0')
I0=I;
end
if ~exist('C')
C=;
end
[ny,nx]=size(I); ep2=ep^; for i=:iter, %% do iterations
% estimate derivatives
I_x = (I(:,[:nx nx])-I(:,[ :nx-]))/;
I_y = (I([:ny ny],:)-I([ :ny-],:))/;
I_xx = I(:,[:nx nx])+I(:,[ :nx-])-*I;
I_yy = I([:ny ny],:)+I([ :ny-],:)-*I;
Dp = I([:ny ny],[:nx nx])+I([ :ny-],[ :nx-]);
Dm = I([ :ny-],[:nx nx])+I([:ny ny],[ :nx-]);
I_xy = (Dp-Dm)/;
% compute flow
Num = I_xx.*(ep2+I_y.^)-*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^);
Den = (ep2+I_x.^+I_y.^).^(/);
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;
C经典版：
//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 = , 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 = ; i < ny; i++)
{
for (int j = ; j < nx; j++)
{
image0[i][j] = image[i][j] = my_image->GetPixelIndex(j, ny-i-);
}
} 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 = ;
for (int t = ; t <= iter; t++)
{
//进度条
my_image->SetProgress((long)*t/iter);
if (my_image->GetEscape())
break;
//////////////////////////////////////////////////////////////////////////
//计算I(:,[2:nx nx])和I(:,[1 1:nx-1])
//公共部分2到nx-1列
for (int i = ; i < ny; i++)
{
for (int j = ; j < nx-; j++)
{
image_tmp1[i][j] = image[i][j+];
image_tmp2[i][j+] = image[i][j];
}
}
for (int i = ; i < ny; i++)
{
image_tmp1[i][nx-] = image[i][nx-];
image_tmp2[i][] = image[i][];
} //计算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 = ; i < ny; i++)
{
for (int j = ; j < nx; j++)
{
image_x[i][j] = (image_tmp1[i][j] - image_tmp2[i][j])/;
image_xx[i][j] = (image_tmp1[i][j] + image_tmp2[i][j]) - *image[i][j];
}
} //////////////////////////////////////////////////////////////////////////
//计算I([2:ny ny],:)和I([1 1:ny-1],:)
//公共部分2到ny-1行
for (int i = ; i < ny-; i++)
{
for (int j = ; j < nx; j++)
{
image_tmp1[i][j] = image[i+][j];
image_tmp2[i+][j] = image[i][j];
}
}
for (int j = ; j < nx; j++)
{
image_tmp1[ny-][j] = image[ny-][j];
image_tmp2[][j] = image[][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 = ; i < ny; i++)
{
for (int j = ; j < nx; j++)
{
image_y[i][j] = (image_tmp1[i][j] - image_tmp2[i][j])/;
image_yy[i][j] = (image_tmp1[i][j] + image_tmp2[i][j]) - *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 = ; i < ny-; i++)
{
for (int j = ; j < nx-; j++)
{
image_tmp1[i][j] = image[i+][j+];
image_tmp2[i+][j+] = image[i][j];
}
}
for (int i = ; i < ny-; i++)
{
image_tmp1[i][nx-] = image[i+][nx-];
image_tmp2[i+][] = image[i][];
}
for (int j = ; j < nx-; j++)
{
image_tmp1[ny-][j] = image[ny-][j+];
image_tmp2[][j+] = image[][j];
}
image_tmp1[ny-][nx-] = image[ny-][nx-];
image_tmp2[][] = image[][];
//计算Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]);
for (int i = ; i < ny; i++)
{
for (int j = ; 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 = ; i < ny-; i++)
{
for (int j = ; j < nx-; j++)
{
image_tmp1[i+][j] = image[i][j+];
image_tmp2[i][j+] = image[i+][j];
}
}
for (int i = ; i < ny-; i++)
{
image_tmp1[i+][nx-] = image[i][nx-];
image_tmp2[i][] = image[i+][];
}
for (int j = ; j < nx-; j++)
{
image_tmp1[][j] = image[][j+];
image_tmp2[ny-][j+] = image[ny-][j];
}
image_tmp1[][nx-] = image[][nx-];
image_tmp2[ny-][] = image[ny-][]; //计算Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
for (int i = ; i < ny; i++)
{
for (int j = ; j < nx; j++)
{
image_dm[i][j] = image_tmp1[i][j] + image_tmp2[i][j];
}
} //////////////////////////////////////////////////////////////////////////
//计算I_xy = (Dp-Dm)/4;
for (int i = ; i < ny; i++)
{
for (int j = ; j < nx; j++)
{
image_xy[i][j] = (image_dp[i][j] - image_dm[i][j])/;
}
} //////////////////////////////////////////////////////////////////////////
//计算过程： //计算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 = ; i < ny; i++)
{
for (int j = ; j < nx; j++)
{
image_num[i][j] = image_xx[i][j]*(image_y[i][j]*image_y[i][j] + ep2)
- *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 = ; i < ny; i++)
{
for (int j = ; 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 = ; i < ny; i++)
{
for (int j = ; j < nx; j++)
{
tmp = (BYTE)image[i][j];
tmp = max(, min(tmp, ));
my_image->SetPixelIndex(j, ny-i-, 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 = ; 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 = ; i < ny; i++)
{
if (matrix[i])
{
delete[] matrix[i];
}
}
delete[] matrix; return true;
}
//////////////////////////////////////////////////////////////////////////
C简洁版：
//TV去噪函数
Mat TVDenoising(Mat img, int iter)
{
int ep = ;
int nx=img.cols;
int ny = img.rows;
double dt = 0.25f;
double lam = 0.0;
int ep2 = ep*ep; double** image = newDoubleMatrix(nx, ny);
double** image0 = newDoubleMatrix(nx, ny); for(int i=;i<ny;i++){
uchar* p=img.ptr<uchar>(i);
for(int j=;j<nx;j++){
image0[i][j]=image[i][j]=(double)p[j];
}
}
//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_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 = ; t <= iter; t++){ //for (int i = 0; i < ny; i++){
// for (int j = 0; j < nx; j++){
// //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;
// int tmp_i1=(i+1)<ny ? (i+1) :(ny-1);
// int tmp_j1=(j+1)<nx ? (j+1): (nx-1);
// int tmp_i2=(i-1) > -1 ? (i-1) : 0;
// int tmp_j2=(j-1) > -1 ? (j-1) : 0;
// image_x[i][j] = (image[i][tmp_j1] - image[i][tmp_j2])/2;
// image_y[i][j]= (image[tmp_i1][j]-image[tmp_i2][j])/2;
// image_xx[i][j] = image[i][tmp_j1] + image[i][tmp_j2]- image[i][j]*2;
// image_yy[i][j]= image[tmp_i1][j]+image[tmp_i2][j] - image[i][j]*2;
// image_dp[i][j]=image[tmp_i1][tmp_j1]+image[tmp_i2][tmp_j2];
// image_dm[i][j]=image[tmp_i2][tmp_j1]+image[tmp_i1][tmp_j2];
// image_xy[i][j] = (image_dp[i][j] - image_dm[i][j])/4;
// 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);
// image[i][j] += dt*(image_num[i][j]/image_den[i][j] + lam*(image0[i][j] - image[i][j]));
// }
//}
for (int i = ; i < ny; i++){
for (int j = ; j < nx; j++){
int tmp_i1=(i+)<ny ? (i+) :(ny-);
int tmp_j1=(j+)<nx ? (j+): (nx-);
int tmp_i2=(i-) > - ? (i-) : ;
int tmp_j2=(j-) > - ? (j-) : ;
double tmp_x = (image[i][tmp_j1] - image[i][tmp_j2])/; //I_x = (I(:,[2:nx nx])-I(:,[1 1:nx-1]))/2;
double tmp_y= (image[tmp_i1][j]-image[tmp_i2][j])/; //I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2;
double tmp_xx = image[i][tmp_j1] + image[i][tmp_j2]- image[i][j]*; //I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
double tmp_yy= image[tmp_i1][j]+image[tmp_i2][j] - image[i][j]*; //I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
double tmp_dp=image[tmp_i1][tmp_j1]+image[tmp_i2][tmp_j2]; //Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]);
double tmp_dm=image[tmp_i2][tmp_j1]+image[tmp_i1][tmp_j2]; //Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
double tmp_xy = (tmp_dp - tmp_dm)/; //I_xy = (Dp-Dm)/4;
double tmp_num = tmp_xx*(tmp_y*tmp_y + ep2)
- *tmp_x*tmp_y*tmp_xy +tmp_yy*(tmp_x*tmp_x + ep2); //Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2);
double tmp_den= pow((tmp_x*tmp_x + tmp_y*tmp_y + ep2), 1.5); //Den = (ep2+I_x.^2+I_y.^2).^(3/2);
image[i][j] += dt*(tmp_num/tmp_den+ lam*(image0[i][j] - image[i][j]));
}
} } Mat new_img;
img.copyTo(new_img);
for(int i=;i<img.rows;i++){
uchar* p=img.ptr<uchar>(i);
uchar* np=new_img.ptr<uchar>(i);
for(int j=;j<img.cols;j++){
int tmp=(int)image[i][j];
tmp=max(,min(tmp,));
np[j]=(uchar)(tmp);
}
} //////////////////////////////////////////////////////////////////////////
//删除内存
//deleteDoubleMatrix(image_x, nx, ny);
//deleteDoubleMatrix(image_y, nx, ny);
//deleteDoubleMatrix(image_xx, nx, ny);
//deleteDoubleMatrix(image_yy, 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); //imshow("Image",img);
//imshow("Denosing",new_img); return new_img;
}
C简洁版修改：
void CImageObj::Total_Variation(int iter, double dt, double epsilon, double lambda)
{
int i, j;
int nx = m_width, ny = m_height;
double ep2 = epsilon * epsilon; double** I_t = NewDoubleMatrix(nx, ny);
double** I_tmp = NewDoubleMatrix(nx, ny);
for (i = ; i < ny; i++)
for (j = ; j < nx; j++)
I_t[i][j] = I_tmp[i][j] = (double)m_imgData[i][j]; for (int t = ; t < iter; t++)
{
for (i = ; i < ny; i++)
{
for (j = ; j < nx; j++)
{
int iUp = i - , iDown = i + ;
int jLeft = j - , jRight = j + ; // 边界处理
if ( == i) iUp = i; if (ny - == i) iDown = i;
if ( == j) jLeft = j; if (nx - == j) jRight = j; double tmp_x = (I_t[i][jRight] - I_t[i][jLeft]) / 2.0;
double tmp_y = (I_t[iDown][j] - I_t[iUp][j]) / 2.0;
double tmp_xx = I_t[i][jRight] + I_t[i][jLeft] - * I_t[i][j];
double tmp_yy = I_t[iDown][j] + I_t[iUp][j] - * I_t[i][j];
double tmp_xy = (I_t[iDown][jRight] + I_t[iUp][jLeft] - I_t[iUp][jRight] - I_t[iDown][jLeft]) / 4.0;
double tmp_num = tmp_yy * (tmp_x * tmp_x + ep2) + tmp_xx * (tmp_y * tmp_y + ep2) - * tmp_x * tmp_y * tmp_xy;
double tmp_den = pow(tmp_x * tmp_x + tmp_y * tmp_y + ep2, 1.5); I_tmp[i][j] += dt*(tmp_num / tmp_den + lambda*(m_imgData[i][j] - I_t[i][j]));
}
} // 一次迭代 for (i = ; i < ny; i++)
for (j = ; j < nx; j++)
{
I_t[i][j] = I_tmp[i][j];
} } // 迭代结束 // 给图像赋值
for (i = ; i < ny; i++)
for (j = ; j < nx; j++)
{
double tmp = I_t[i][j];
tmp = max(, min(tmp, ));
m_imgData[i][j] = (unsigned char)tmp;
} DeleteDoubleMatrix(I_t, nx, ny);
DeleteDoubleMatrix(I_tmp, nx, ny);
}
【转载自】
保持图像纹理特征的超分辨率重建方法研究_百度学术 http://xueshu.baidu.com/usercenter/paper/show?paperid=a89942cdaa9f99edb04b2101216541ad&site=xueshu_se
TV全变分图像去噪的研究 - 百度文库 https://wenku.baidu.com/view/00def4edb04e852458fb770bf78a6529647d3517.html
全变分（TV）模型原理与C++实现 - cyh706510441的专栏 - CSDN博客 https://blog.csdn.net/cyh706510441/article/details/45194223
VISL http://visl.technion.ac.il/~gilboa/PDE-filt/tv_denoising.html
经典的变分法图像去噪的C++实现 - InfantSorrow - 博客园 http://www.cnblogs.com/CCBB/archive/2010/12/29/1920884.html
全变分TV图像去噪 - 小魏的修行路 - CSDN博客 https://blog.csdn.net/xiaowei_cqu/article/details/18051029
全变分（TV）模型原理与C++实现 - cyh706510441的专栏 - CSDN博客 https://blog.csdn.net/cyh706510441/article/details/45194223