【算法随记六】一段Matlab版本的Total Variation(TV)去噪算法的C语言翻译。

　　最近看到一篇文章讲IMAGE DECOMPOSITION，里面提到了将图像分为Texture layer和Structure layer，测试了不少方法，对于那些具备很是强烈纹理的图像，总以为用TV去燥的方法分离的结果都比其余的方法都要好（好比导向、双边），好比下图：

   

 　　再好比：

 

 　　可见TV能够把纹理很好的提取出来。

　　如今应该能找到不少的TV代码，好比IPOL上就有，详见 http://www.ipol.im/pub/art/2013/61/。

　　我在其余地方也见过一些，好比这里： http://yu-li.github.io/paper/li_eccv14_jpeg.zip，他是借助于FFT实现的，固然少不了屡次迭代，速度也是比较慢的。

　　我还收藏了好久前一位朋友写的M代码，可是如今我不知道把他QQ或者微信弄到哪里去了，也不知道他会不会介意我把他的代码分享出来。

function dualROF() clc f0=imread('rr.bmp'); f0=f0(:,:,1); [m,n]=size(f0); f0=double(f0); lamda=30; % smoothness paramter, the larger the smoother tao=.125; % fixed do not change it. p1=zeros(m,n); p2=zeros(m,n); tic for step=1:100
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% div_p=div(p1,p2); cx=Fx(div_p-f0/lamda); cy=Fy(div_p-f0/lamda);

　　　　abs_c=sqrt(cx.^2+cy.^2);
　　　　p1=(p1+tao*cx)./(1+tao*abs_c);
　　　　p2=(p2+tao*cy)./(1+tao*abs_c);

end u=f0-lamda*div_p; toc figure; imagesc(f0); colormap(gray); axis off; axis equal; figure; imagesc(u); colormap(gray); axis off; axis equal; % Compute divergence using backward derivative function f = div(a,b) f = Bx(a)+By(b); % Forward derivative operator on x with boundary condition u(:,:,1)=u(:,:,1) function Fxu = Fx(u) [m,n] = size(u); Fxu = circshift(u,[0 -1])-u; Fxu(:,n) = zeros(m,1); % Forward derivative operator on y with boundary condition u(1,:,:)=u(m,:,:) function Fyu = Fy(u) [m,n] = size(u); Fyu = circshift(u,[-1 0])-u; Fyu(m,:) = zeros(1,n); % Backward derivative operator on x with boundary condition Bxu(:,1)=u(:,1) function Bxu = Bx(u) [~,n] = size(u); Bxu = u - circshift(u,[0 1]); Bxu(:,1) = u(:,1); Bxu(:,n) = -u(:,n-1); % Backward derivative operator on y with boundary condition Bxu(1,:)=u(1,:) function Byu = By(u) [m,~] = size(u); Byu = u - circshift(u,[1 0]); Byu(1,:) = u(1,:); Byu(m,:) = -u(m-1,:);

　　M的代码，代码量不大，那是由于Matlab的向量化确实很厉害，可是这个代码仍是很慢的，256*256的灰度图迭代100次都要700ms了。

　　这里抛开一些优化不说，用这个circshift会形成很大的性能损失，咱们稍微分析下就能看到用这个地方其实就是简单的水平或者垂直方向的差分，彻底没有必要这样写。

　　直接按照代码的意思用C语言把他们展开并不作其余的优化可获得大概下面这种不怎么好的代码：

int IM_DualTVDenoising(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride,  float Lamda = 20 , int Iter = 20) { int Channel = Stride / Width; if ((Src == NULL) || (Dest == NULL))                        return IM_STATUS_NULLREFRENCE; if ((Width <= 0) || (Height <= 0))                            return IM_STATUS_INVALIDPARAMETER; if ((Channel != 1) && (Channel != 3) && (Channel != 4))        return IM_STATUS_INVALIDPARAMETER; if (Channel == 1) { float tao = 0.125; // fixed do not change it.
        float InvLamda = 1.0 / Lamda; float *p1 = (float *)malloc(Width * Height * sizeof(float)); float *p2 = (float *)malloc(Width * Height * sizeof(float)); float *div_p = (float *)malloc(Width * Height * sizeof(float)); float *cx = (float *)malloc(Width * Height * sizeof(float)); float *cy = (float *)malloc(Width * Height * sizeof(float)); float *Temp = (float *)malloc(Width * Height * sizeof(float)); int X, Y; float q1, q2, q, abs_c; float *LineP1, *LineP2, *LineP3, *LineP4; unsigned char *LinePS, *LinePD; for (int Z = 0; Z < Iter; Z++) { //Div(p1, p2, div_p);

            for (Y = 0; Y < Height; Y++) { LineP1 = p1 + Y * Width;                        //Fx(Temp, cx);
                LineP2 = cx + Y * Width; LineP2[0] = LineP1[0]; for (X = 1; X < Width; X++) { LineP2[X] = LineP1[X] - LineP1[X - 1]; } LineP2[Width - 1] = -LineP1[Width - 2]; } memcpy(cy, p2, Width * sizeof(float)); for (Y = 1; Y < Height; Y++) { LineP1 = (float *)(p2 + (Y - 1)* Width); LineP2 = (float *)(p2 + Y * Width);            //Fy(Temp, cy);
                LineP3 = (float *)(cy + Y * Width); for (X = 0; X < Width; X++) { LineP3[X] = LineP2[X] - LineP1[X]; } } LineP1 = (float *)(p2 + (Height - 2) * Width); LineP2 = (float *)(cy + (Height - 1) * Width); for (X = 0; X < Width; X++) { LineP2[X] = -LineP1[X]; } for (Y = 0; Y < Height; Y++) { LineP1 = (float *)(cx + Y * Width); LineP2 = (float *)(cy + Y * Width);            //Fy(Temp, cy);
                LineP3 = (float *)(div_p + Y * Width); for (X = 0; X < Width; X++) { LineP3[X] = LineP1[X] + LineP2[X]; } } for (Y = 0; Y < Height; Y++) { LineP1 = (float *)(div_p + Y * Width); LineP2 = (float *)(Temp + Y * Width); LinePS = Src + Y * Stride; for (X = 0; X < Width; X++) { LineP2[X] = LineP1[X] - LinePS[X] * InvLamda; } } for (Y = 0; Y < Height; Y++) { LineP1 = (float *)(Temp + Y * Width);                        //Fx(Temp, cx);
                LineP2 = (float *)(cx + Y * Width); for (X = 0; X < Width - 1; X++) { LineP2[X] = LineP1[X + 1] - LineP1[X]; } LineP2[Width - 1] = 0; } for (Y = 0; Y < Height - 1; Y++) { LineP1 = (float *)(Temp + Y * Width); LineP2 = (float *)(Temp + (Y + 1) * Width);            //Fy(Temp, cy);
                LineP3 = (float *)(cy + Y * Width); for (X = 0; X < Width; X++) { LineP3[X] = LineP2[X] - LineP1[X]; } } memset(Temp + (Height - 1) * Width, 0, Width * sizeof(float)); for (Y = 0; Y < Height; Y++) { LineP1 = (float *)(p1 + Y * Width); LineP2 = (float *)(p2 + Y * Width); LineP3 = (float *)(cx + Y * Width); LineP4 = (float *)(cy + Y * Width); for (X = 0; X < Width; X++) { abs_c = sqrt(LineP3[X] * LineP3[X] + LineP4[X] * LineP4[X]); abs_c = 1 / (1 + tao * abs_c); LineP1[X] = (LineP1[X] + tao * LineP3[X]) * abs_c; LineP2[X] = (LineP2[X] + tao * LineP4[X]) * abs_c; } } } for (Y = 0; Y < Height; Y++) { LineP1 = (float *)(div_p + Y * Width); LinePS = Src + Y * Stride; LinePD = Dest + Y * Stride; for (X = 0; X < Width; X++) { LinePD[X] = IM_ClampToByte((int)(LinePS[X] - Lamda * LineP1[X])); } } free(p1); free(p2); free(div_p); free(cx); free(cy); free(Temp); } else { } }

　　算法明显占用很大的内存，并且看起来别扭，不过速度仍是杠杠的，256*256的灰度图迭代100次都要30ms了。反编译看了下代码，编译器对代码作了很好的SIMD指令优化。

　　上面的C语言仍是能够继续优化的，这就须要你们本身的认真的去研读代码深层次的逻辑关系了，实际上能够只要上面的一半的临时内存的，并且不少计算能够集中在一个循环里完成，能够手动内嵌SIMD指令，或者直接使用编译器的优化能力，基本上这样的简单的算法逻辑编译器编译后的速度不会比咱们手写的SIMD指令慢，有的时候仍是会快一些，不得不佩服那些写编译器的大牛。优化后的速度大概在14ms左右。

　　研究TV算法须要很好的数学功底，之前朋友曾经给我寄过一本书，里面都是微分方面的数学公式，看的我吓死了，不过TV算法彷佛有不少很好的应用，也曾经流行过一段时间，惋惜如今深度学习一出来，不少人都喜欢这种直接从海量数据中建造黑盒模型，而对那些有着很明显的数学逻辑的算法嗤之以鼻了，真有点惋惜。

　　之前在基于总变差模型的纹理图像中图像主结构的提取方法 一文中曾提到那个论文附带的Matlab代码没有什么意义，由于他很难转换成C的代码，即时转换成功了，也处理不了大图，可是本文这里的TV算法总的来讲在内存占用或者速度方面都还使人满意。

　　在去噪效果上，这个算法还算能够：

          

　　本文Demo下载地址：  http://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar， 算法位于Denoise --> TV Denoising下。