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