专栏首页一心无二用,本人只专注于基础图像算法的实现与优化。【算法随记六】一段Matlab版本的Total Variation(TV)去噪算法的C语言翻译。

【算法随记六】一段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);

  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下。

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • SSE图像算法优化系列二:高斯模糊算法的全面优化过程分享(一)。

         这里的高斯模糊采用的是论文《Recursive implementation of the Gaussian filter》里描述的递归算法。 ? ...

    用户1138785
  • 【短道速滑一】OpenCV中cvResize函数使用双线性插值缩小图像到长宽大小一半时速度飞快(比最近邻还快)之异象解析和自我实现。

      今天,一个朋友想使用我的SSE优化Demo里的双线性插值算法,他已经在项目里使用了OpenCV,因此,我就建议他直接使用OpenCV,朋友的程序非常注意效率...

    用户1138785
  • PhotoShop算法原理解析系列 - 像素化---》碎片。

          接着上一篇文章的热度,继续讲讲一些稍微简单的算法吧。       本文来讲讲碎片算法,先贴几个效果图吧: ? ? ? ?   这是个破坏性的...

    用户1138785
  • 使用Apache Common CSV读写CSV文件

    jar包下载地址:http://commons.apache.org/proper/commons-csv/,点击Download进行下载!

    卡尔曼和玻尔兹曼谁曼
  • 6.3 VR扫描:工信部近期将发放5G商用牌照;添田武人卸任索尼互娱中国区总裁

    今日,工信部正式确认,中国将于近期发放5G商用牌照。此前,工信部曾多次表示,在5G产业链成熟的情况下,工信部将会适时发放5G牌照。此次消息的确认,意味着我国包括...

    VRPinea
  • Vs Code搭建 TypeScript 开发环境

      更多细节参考:http://www.typescriptlang.org/docs/handbook/tsconfig-json.html

    喝茶去
  • 深度学习Pytorch检测实战 - Notes - 第3章 网络骨架

    Sigmoid函数将特征压缩到了(0,1)区间,0端对应抑制状态,而1对应激活状态,中间部分梯度较大。Sigmoid函数可以用来做二分类,但其计算量较大,并且容...

    肉松
  • VUE 插件注册

    通过全局方法Vue.use()使用插件,它需要你调用new Vue()启动应用之前完成

    用户7533190
  • 机器学习为CAD插上一双翅膀(上)

    基于机器学习在无数行业中得到了充分利用,从网络上的提示性搜索到照片库存图像推荐。其核心是,推荐引擎可以在大量数据库中查询相关信息(文本、图像等),并在用户与给定...

    AiTechYun
  • 3. 经典卷积网络之GooleInceptionNet

    GooleInceptionNet首次出现是在2014年的ILSVRC的比赛中,当时是第一名,最大的特点就是控制计算量的同时获得了比较好的分类性能--top-5...

    和蔼的zhxing

扫码关注云+社区

领取腾讯云代金券