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 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 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; 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;
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 = 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; } //////////////////////////////////////////////////////////////////////////
C简洁版：
//TV去噪函数 Mat TVDenoising(Mat img, int iter) { int ep = 1; 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=0;i<ny;i++){ uchar* p=img.ptr<uchar>(i); for(int j=0;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 = 1; 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 = 0; i < ny; i++){ for (int j = 0; j < nx; j++){ 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; double tmp_x = (image[i][tmp_j1] - image[i][tmp_j2])/2; //I_x = (I(:,[2:nx nx])-I(:,[1 1:nx-1]))/2; double tmp_y= (image[tmp_i1][j]-image[tmp_i2][j])/2; //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]*2; //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]*2; //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)/4; //I_xy = (Dp-Dm)/4; double tmp_num = tmp_xx*(tmp_y*tmp_y + ep2) - 2*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=0;i<img.rows;i++){ uchar* p=img.ptr<uchar>(i); uchar* np=new_img.ptr<uchar>(i); for(int j=0;j<img.cols;j++){ int tmp=(int)image[i][j]; tmp=max(0,min(tmp,255)); 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 = 0; i < ny; i++) for (j = 0; j < nx; j++) I_t[i][j] = I_tmp[i][j] = (double)m_imgData[i][j]; for (int t = 0; t < iter; t++) { for (i = 0; i < ny; i++) { for (j = 0; j < nx; j++) { int iUp = i - 1, iDown = i + 1; int jLeft = j - 1, jRight = j + 1; // 边界处理 if (0 == i) iUp = i; if (ny - 1 == i) iDown = i; if (0 == j) jLeft = j; if (nx - 1 == 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] - 2 * I_t[i][j]; double tmp_yy = I_t[iDown][j] + I_t[iUp][j] - 2 * 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) - 2 * 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 = 0; i < ny; i++) for (j = 0; j < nx; j++) { I_t[i][j] = I_tmp[i][j]; } } // 迭代结束 // 给图像赋值 for (i = 0; i < ny; i++) for (j = 0; j < nx; j++) { double tmp = I_t[i][j]; tmp = max(0, min(tmp, 255)); 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