非局部均值滤波的指令集优化和加速(针对5*5的搜索特例,可达到单核1080P灰度图 28ms/帧的速度)。

本文主要是介绍非局部均值滤波的指令集优化和加速(针对5*5的搜索特例,可达到单核1080P灰度图 28ms/帧的速度)。,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

       非局部均值滤波(Non Local Means)作为三大最常提起来的去燥和滤波算法之一(双边滤波、非局部均值、BM3D),也是有着很多的论文作为研究和比较的对象,但是也是有着致命的缺点,速度慢,严重的影响了算法的应用范围。目前在已有的文献中尚未看到在不对算法的本质原理上做更改的情况下,能取得实时的效果,本文呢,也不求得到这个目的,只是对现有的开放的资源上来取得更进一步的提升。

  标准的NL-Means算法中,一般有三个参数,搜索半径SearchRadius,块半径PatchRadius,以及一个决定平滑程度的高斯函数参数Delta。在百度上能够搜索到的大部分文章所描述的提速算法都是使用积分图来提升NL-Means的速度,这也是目前来说唯一比较靠谱的优化技术,通过积分图,可以做到算法和块半径PatchRadius的大小基本无关,和Delta也无关,和SearchRadius成平方关系。

  因此,在我们很多的严重的噪音图像中,SearchRadius至少需要取到7以上(涉及15*15= 225个领域范围)才有明显的效果,因此,这就相当于要计算225次全图的某种计算,即使是每次这种计算只需要1ms(通常,任何图像处理算法无法超越同等内存大小的memcpy的大小的),也需要225ms的,因此,确实比较感尴。

  通过多线程方式可以适当对这个过程进行加速,毕竟每个像素点的处理相对来说还是独立的,但是,这个加速也收到物理核心的限制,就是8核的机器,满利用,也无法达到8倍的加速效果。

       话说回来,这么大的计算量,用GPU也都是很吃力的。

  目前使用积分图来记性NL-Means算法的比较好的文章也是大家常看的还是这一篇:

       非局部均值滤波(NL-means)算法的积分图加速原理与C++实现

  其实积分图一直有个问题,可能很多搞图像的人都没有注意到,或者说这个问题可能对某些算法的影响还不是很大,不足有对大家注意到或者关注到,即积分图存在着一下2个方面的问题和缺点:

  1、当图像较大时,积分图无法使用int类型来保存,我们必须选择能够容纳更大数据范围的数据类型来存储。

  如果我们保存的是一副字节图像的积分图,考虑极端情况,每个图像的值都是255,则积分图像最多只可保存uint.Maxvalue / 255 = 16843009个像素,大约4100*4100大小的图像,一般来说,这个规模对于实际的应用来说是足够了。

  但是我们在实际中用到的很多情况,不是直接的图像积分图,还常用的有平方积分图,即求图像平方值后的积分图,这个时候每一个元素的最大值达到了65025,极限情况下uint类型只可保存uint.Maxvalue / 255 =66051个元素的总和,只大概是指256*256个图像的大小了,已经远远的无法满足实际的应用需求。

  因此,为了实现这个要求,我们必须选择能够容纳更大数据范围的数据类型,这里有三个选择:long long(int64) / float / double

  第一个long long类型,即64位的整形数据,这个数据的表达范围已经完全够我们在图像中使用积分图了,而且保存的数据是非常准确的,特别是对于图像方面的数据来说(都是整形的),但是有个致命的问题: 速度相当相当的慢,特别是同样的计算和int类型比较的话,那真的不是一个档次上的。我一直不太理解,现在大部分都是64位系统,为什么对64位的数据的支持还是这么的弱。而且我们看大部分指令集优化的函数对64位整形的支持都比较少。因此,非常不建议使用这个类型。

  第二个float类型。如果使用这个类型,保存的数据范围是没有什么大的问题的,我们在网络上看到的文章大部分也是使用这个类型来保存结果的。但是,我在实践中多次遇到用float类型得不到正确的结果的问题,后来发现核心的原因是float的计算精度严重不足,特别是对于积分图这种连续的加法的计算,累计误差会越来越严重,当计算量足够大时,就会出现明显的误差瑕疵。因此,这个数据类型从本质上来说,对积分图是不够安全的。

  关于这一点,实际上已经有不少作者注意到了,我在博文:SSE图像算法优化系列十四:局部均方差及局部平方差算法的优化 中也有提及。

  第三个是double类型,这个类型也是64位的,因此,数据范围毫无问题,计算精度经过测试,也是没有什么问题的,而且,编译器和指令集的支持和优化做的都还很不错, 因此,个人认为这个数据类型是用来保存积分图最为合适的类型,但是有一个不友好的特点,计算速度慢,而且指令集对其能加速的空间有限。

  2、积分图虽然能做到某些算法和参数无关,但是其并不是最佳的最速度

  使用积分图技术,首先是要分配积分图占用的那部分额外的内存(而且是相当客观的内存),其次,积分图本身的计算也是需要一定时间的。还有,积分图的计算必须对边界部分做特别的判断,这个当参数较大时,计算量有的时候还是相当可观的。

  还是回到我们的非局部均值滤波上吧。上面说了这么多,意思就是虽然非局部均值可以用积分图去优化,但是还是不是很好,那有没有更好的更快的实现呢,其实在我的下面两篇博客里就已经有了相关的技术。

  一是 : SSE图像算法优化系列十三:超高速BoxBlur算法的实现和优化(Opencv的速度的五倍) 

       二是: SSE图像算法优化系列十四:局部均方差及局部平方差算法的优化 

  因为非局部均值里使用的积分图是差的平方的积分图,因此我们只要对上述参考博文一里面的参与积分的数据稍微换一下即可,一个简单的C++代码如下所示:

//    SearchRadius    搜索的领域半径
//    PatchRadius      计算相似度的块的半径
//    Delta     高斯平滑的参数
int IM_NLM_Denoising(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int SearchRadius, int PatchRadius, float Delta)
{int Status = IM_STATUS_OK;int Channel = Stride / Width;if (Src == NULL)                        return IM_STATUS_NULLREFRENCE;if ((Width <= 0) || (Height <= 0))        return IM_STATUS_INVALIDPARAMETER;if ((Channel != 1))                        return IM_STATUS_INVALIDPARAMETER;    int ExpandW = Width + SearchRadius + PatchRadius + SearchRadius + PatchRadius;int ExpandH = Height + SearchRadius + PatchRadius + SearchRadius + PatchRadius;int DiffW = Width + PatchRadius + PatchRadius;int DiffH = Height + PatchRadius + PatchRadius;float *Weight = (float *)calloc(Width * Height , sizeof(float));float *Sum = (float *)calloc(Width * Height , sizeof(float));int *ColValue = (int *)malloc((Width + PatchRadius + PatchRadius) * sizeof(int));unsigned char *Expand = (unsigned char *)malloc(ExpandH * ExpandW * sizeof(unsigned char));if ((Weight == NULL) || (Sum == NULL) || (ColValue == NULL) || (Expand == NULL)){Status = IM_STATUS_OUTOFMEMORY;goto FreeMemory;}Status = IM_GetExpandImage(Src, Expand, Width, Height, Stride, ExpandW, SearchRadius + PatchRadius, SearchRadius + PatchRadius, SearchRadius + PatchRadius, SearchRadius + PatchRadius, IM_EDGE_MIRROR);if (Status != IM_STATUS_OK)    goto FreeMemory;int Area = (2 * PatchRadius + 1) * (2 * PatchRadius + 1);float Inv = -1.0f / Area / Delta / Delta;for (int YY = -SearchRadius; YY <= SearchRadius; YY++){for (int XX = -SearchRadius; XX <= SearchRadius; XX++){for (int Y = 0; Y < Height; Y++){float *LinePS = Sum + Y * Width;float *LinePW = Weight + Y * Width;unsigned char *LinePE = Expand + (Y + YY + SearchRadius + PatchRadius) * ExpandW + XX + SearchRadius + PatchRadius;if (Y == 0){memset(ColValue, 0, DiffW * sizeof(int));for (int Z = -PatchRadius; Z <= PatchRadius; Z++){unsigned char *LineP1 = Expand + (Z + PatchRadius + YY + SearchRadius) * ExpandW + XX + SearchRadius;unsigned char *LineP2 = Expand + (Z + PatchRadius + SearchRadius) * ExpandW + SearchRadius;for (int X = 0; X < DiffW; X++){int Value = LineP2[X] - LineP1[X];ColValue[X] += Value * Value;}}}else{unsigned char *LineOut1 = Expand + (Y - 1 + YY + SearchRadius) * ExpandW + XX + SearchRadius;unsigned char *LineOut2 = Expand + (Y - 1 + SearchRadius) * ExpandW + SearchRadius;unsigned char *LineIn1 = Expand + (Y + PatchRadius + PatchRadius + YY + SearchRadius) * ExpandW + XX + SearchRadius;unsigned char *LineIn2 = Expand + (Y + PatchRadius + PatchRadius + SearchRadius) * ExpandW + SearchRadius;for (int X = 0; X < DiffW; X++){int Out = LineOut2[X] - LineOut1[X];int In = LineIn2[X] - LineIn1[X];ColValue[X] -= Out * Out - In * In;                                    //    更新列数据}}int SumA = IM_SumofArray(ColValue, PatchRadius * 2 + 1);                //    处理每行第一个数据    float W = IM_Fexp(SumA * Inv);LinePW[0] += W;LinePS[0] += W * LinePE[0];for (int X = 1; X < Width; X++){SumA = SumA - ColValue[X - 1] + ColValue[X + PatchRadius + PatchRadius];float W = IM_Fexp(SumA * Inv);LinePW[X] += W;LinePS[X] += W * LinePE[X];}}}}for (int Y = 0; Y < Height; Y++){int Index = Y * Width;unsigned char *LinePD = Dest + Y * Stride;for (int X = 0; X < Width; X++){LinePD[X] = Sum[Index + X] / (Weight[Index + X]);}}
FreeMemory:if (Expand != NULL)            free(Expand);if (Weight != NULL)            free(Weight);if (Sum != NULL)               free(Sum);if (ColValue != NULL)         free(ColValue);return Status;
}

  对于常用的搜索半径为10,块半径为3, 500 * 500的灰度图耗时大概是500ms, 相当的慢啊。 

  通过我一系列博文里的资料,可以知道上面的循环内的代码其实是很容易进行指令集优化的,基本上我那个方框模糊的优化是同一个技巧,经过指令集优化后,500*500的灰度图的耗时大概在200ms左右,如果加上线程技术,可以优化到75ms左右,这个时间和 非局部均值滤波(NL-means)算法的CUDA优化加速 文章里提的CUDA的优化速度基本已经差不多了。

  

  另外,在搜索半径较小时,一种可行的优化方式是进行行列分离的卷积,即先计算中心行的结果,然后已这个结果为原始数据,在计算列方向的卷积,注意不是同时计算中心列和中心行的累加,而是有顺序的,这样就可以利用到周边领域所有的信息,这个时候计算量就会大为下降,计算的耗时也可以明显提高。这种优化方式必须测试是否对去燥的效果有很大的影响,因为他有可能会产生较为明显的水平或垂直线条效果。

  另外,如果搜索半径较大,还可以尝试在上述基础上再进行45度和135度两个方向的卷积,以便抵消这种线条效果,同样提速也还是很明显的。

  作为一个特例,有些情况下我们可能需要搜索半径为2,块大小也为2的非局部均值滤波,这种尺寸的滤波对于高强度的高斯噪音是没有什么去燥效果的,但是对于小规模的噪音,比如一般视频会议中的噪音还是有较强的抑制作用,因此,也还是有应用场景的,但是这种场景对算法的速度提出了极高的要求,如果考虑流畅性,给这个算法的处理不易超过20ms, 而现在的视频流越来越高清了,因此,对这个算法的优化处理就必须做的更到位。

  针对这个特例,我又做了一些优化,首先,因为是小半径,所以可以使用行列分离的算法,即先计算如下区域的合成结果:

                      

   计算得到中间的结果后再通过中间结果计算下述区域的合成结果作为最终值:

                      

  这里面的优化技巧有很多很多。 

  我这里提几个想法供参考:  

  1、中心点C的处的权重不需要计算,因为必然为1。

  2、在进行水平计算时,去除掉中线后,只有4个点了,PP/P/N/NN,我们可以通过一些组合手段把他们的权重以及权重乘以值的量一次性计算出来,这样就可以直接把中间值给搞定了,如此做的好处是,不用多次保存累加的权重值和乘以值,也就是说舍弃掉了前面代码里的Weight和Sum变量,这对算法速度的提高也是有极大的贡献的,垂直的计算也是一样的。 

       3、那个exp函数是个比较耗时的函数,在http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/中一个比较快速的函数,虽然精度一般,但是经过测试,可以满足本算法的需求,但是其计算量小很多。

  4、去掉中心点后,无论在垂直或者水平方向都只有4个数据,这为SIMD指令的某些读取和优化提供了无尚的便利和方便,而有些数据真的是巧合加天工设计,比如水平方向处理时,我们需要一次性处理4个点(出去中间那个),而我们用_mm_loadl_epi64恰好可以读取8个字节,这个时候,读取的8个字节如下:

        A0   A1   A2   A3   A4   A5   A6   A7

  正好可以组成4对5个组合 :

          A0  A1  A2  A3  A4

          A1   A2   A3   A4   A5 

          A2   A3  A4   A5   A6

          A3   A4  A5   A6       A7

  如果多一个字节,都不好处理了(主要是不好读取),多么完美的事情啊。

  通过多种优化方式后,对于常见的1080P视频(1920*1080),处理其Y分量(灰度)耗时能做到28ms了,如果使用2个线程或4个线程,则可以完全满足实时的需求,如果是720P的视频,则单核也能到20ms。

   提供两个测试DEMO做速度比较吧:

         5X5的特例非局部均值去燥Demo:  https://files.cnblogs.com/files/Imageshop/NL-Means5x5.rar?t=1696752228&download=true

      任意尺寸的非局部均值去燥Demo:  https://files.cnblogs.com/files/Imageshop/NL-Means.rar?t=1696752228&download=true

  如果想时刻关注本人的最新文章,也可关注公众号或者添加本人微信:  laviewpbt

                             

这篇关于非局部均值滤波的指令集优化和加速(针对5*5的搜索特例,可达到单核1080P灰度图 28ms/帧的速度)。的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/174090

相关文章

Vue3 的 shallowRef 和 shallowReactive:优化性能

大家对 Vue3 的 ref 和 reactive 都很熟悉,那么对 shallowRef 和 shallowReactive 是否了解呢? 在编程和数据结构中,“shallow”(浅层)通常指对数据结构的最外层进行操作,而不递归地处理其内部或嵌套的数据。这种处理方式关注的是数据结构的第一层属性或元素,而忽略更深层次的嵌套内容。 1. 浅层与深层的对比 1.1 浅层(Shallow) 定义

HDFS—存储优化(纠删码)

纠删码原理 HDFS 默认情况下,一个文件有3个副本,这样提高了数据的可靠性,但也带来了2倍的冗余开销。 Hadoop3.x 引入了纠删码,采用计算的方式,可以节省约50%左右的存储空间。 此种方式节约了空间,但是会增加 cpu 的计算。 纠删码策略是给具体一个路径设置。所有往此路径下存储的文件,都会执行此策略。 默认只开启对 RS-6-3-1024k

认识、理解、分类——acm之搜索

普通搜索方法有两种:1、广度优先搜索;2、深度优先搜索; 更多搜索方法: 3、双向广度优先搜索; 4、启发式搜索(包括A*算法等); 搜索通常会用到的知识点:状态压缩(位压缩,利用hash思想压缩)。

hdu1240、hdu1253(三维搜索题)

1、从后往前输入,(x,y,z); 2、从下往上输入,(y , z, x); 3、从左往右输入,(z,x,y); hdu1240代码如下: #include<iostream>#include<algorithm>#include<string>#include<stack>#include<queue>#include<map>#include<stdio.h>#inc

使用opencv优化图片(画面变清晰)

文章目录 需求影响照片清晰度的因素 实现降噪测试代码 锐化空间锐化Unsharp Masking频率域锐化对比测试 对比度增强常用算法对比测试 需求 对图像进行优化,使其看起来更清晰,同时保持尺寸不变,通常涉及到图像处理技术如锐化、降噪、对比度增强等 影响照片清晰度的因素 影响照片清晰度的因素有很多,主要可以从以下几个方面来分析 1. 拍摄设备 相机传感器:相机传

Open3D 基于法线的双边滤波

目录 一、概述 1.1原理 1.2实现步骤 1.3应用场景 二、代码实现 2.1关键函数 输入参数: 输出参数: 参数影响: 2.2完整代码 三、实现效果 3.1原始点云 3.2滤波后点云 Open3D点云算法汇总及实战案例汇总的目录地址: Open3D点云算法与点云深度学习案例汇总(长期更新)-CSDN博客 一、概述         基于法线的双边

MySQL高性能优化规范

前言:      笔者最近上班途中突然想丰富下自己的数据库优化技能。于是在查阅了多篇文章后,总结出了这篇! 数据库命令规范 所有数据库对象名称必须使用小写字母并用下划线分割 所有数据库对象名称禁止使用mysql保留关键字(如果表名中包含关键字查询时,需要将其用单引号括起来) 数据库对象的命名要能做到见名识意,并且最后不要超过32个字符 临时库表必须以tmp_为前缀并以日期为后缀,备份

SWAP作物生长模型安装教程、数据制备、敏感性分析、气候变化影响、R模型敏感性分析与贝叶斯优化、Fortran源代码分析、气候数据降尺度与变化影响分析

查看原文>>>全流程SWAP农业模型数据制备、敏感性分析及气候变化影响实践技术应用 SWAP模型是由荷兰瓦赫宁根大学开发的先进农作物模型,它综合考虑了土壤-水分-大气以及植被间的相互作用;是一种描述作物生长过程的一种机理性作物生长模型。它不但运用Richard方程,使其能够精确的模拟土壤中水分的运动,而且耦合了WOFOST作物模型使作物的生长描述更为科学。 本文让更多的科研人员和农业工作者

hdu 4517 floyd+记忆化搜索

题意: 有n(100)个景点,m(1000)条路,时间限制为t(300),起点s,终点e。 访问每个景点需要时间cost_i,每个景点的访问价值为value_i。 点与点之间行走需要花费的时间为g[ i ] [ j ] 。注意点间可能有多条边。 走到一个点时可以选择访问或者不访问,并且当前点的访问价值应该严格大于前一个访问的点。 现在求,从起点出发,到达终点,在时间限制内,能得到的最大

AI基础 L9 Local Search II 局部搜索

Local Beam search 对于当前的所有k个状态,生成它们的所有可能后继状态。 检查生成的后继状态中是否有任何状态是解决方案。 如果所有后继状态都不是解决方案,则从所有后继状态中选择k个最佳状态。 当达到预设的迭代次数或满足某个终止条件时,算法停止。 — Choose k successors randomly, biased towards good ones — Close