【演算法隨記六】一段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 =https://www.cnblogs.com/Imageshop/p/= 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下，

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/17742.html

標籤：其他

上一篇：Unity3D 卡通描邊之控制線條粗細

下一篇：MPV原始碼探究：背景及準備作業