非局部均值滤波的指令集优化和加速(针对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

相关文章

SpringBoot3实现Gzip压缩优化的技术指南

《SpringBoot3实现Gzip压缩优化的技术指南》随着Web应用的用户量和数据量增加,网络带宽和页面加载速度逐渐成为瓶颈,为了减少数据传输量,提高用户体验,我们可以使用Gzip压缩HTTP响应,... 目录1、简述2、配置2.1 添加依赖2.2 配置 Gzip 压缩3、服务端应用4、前端应用4.1 N

Spring Boot + MyBatis Plus 高效开发实战从入门到进阶优化(推荐)

《SpringBoot+MyBatisPlus高效开发实战从入门到进阶优化(推荐)》本文将详细介绍SpringBoot+MyBatisPlus的完整开发流程,并深入剖析分页查询、批量操作、动... 目录Spring Boot + MyBATis Plus 高效开发实战:从入门到进阶优化1. MyBatis

MyBatis 动态 SQL 优化之标签的实战与技巧(常见用法)

《MyBatis动态SQL优化之标签的实战与技巧(常见用法)》本文通过详细的示例和实际应用场景,介绍了如何有效利用这些标签来优化MyBatis配置,提升开发效率,确保SQL的高效执行和安全性,感... 目录动态SQL详解一、动态SQL的核心概念1.1 什么是动态SQL?1.2 动态SQL的优点1.3 动态S

Python如何使用__slots__实现节省内存和性能优化

《Python如何使用__slots__实现节省内存和性能优化》你有想过,一个小小的__slots__能让你的Python类内存消耗直接减半吗,没错,今天咱们要聊的就是这个让人眼前一亮的技巧,感兴趣的... 目录背景:内存吃得满满的类__slots__:你的内存管理小助手举个大概的例子:看看效果如何?1.

一文详解SpringBoot响应压缩功能的配置与优化

《一文详解SpringBoot响应压缩功能的配置与优化》SpringBoot的响应压缩功能基于智能协商机制,需同时满足很多条件,本文主要为大家详细介绍了SpringBoot响应压缩功能的配置与优化,需... 目录一、核心工作机制1.1 自动协商触发条件1.2 压缩处理流程二、配置方案详解2.1 基础YAML

MySQL中慢SQL优化的不同方式介绍

《MySQL中慢SQL优化的不同方式介绍》慢SQL的优化,主要从两个方面考虑,SQL语句本身的优化,以及数据库设计的优化,下面小编就来给大家介绍一下有哪些方式可以优化慢SQL吧... 目录避免不必要的列分页优化索引优化JOIN 的优化排序优化UNION 优化慢 SQL 的优化,主要从两个方面考虑,SQL 语

MySQL中慢SQL优化方法的完整指南

《MySQL中慢SQL优化方法的完整指南》当数据库响应时间超过500ms时,系统将面临三大灾难链式反应,所以本文将为大家介绍一下MySQL中慢SQL优化的常用方法,有需要的小伙伴可以了解下... 目录一、慢SQL的致命影响二、精准定位问题SQL1. 启用慢查询日志2. 诊断黄金三件套三、六大核心优化方案方案

Redis中高并发读写性能的深度解析与优化

《Redis中高并发读写性能的深度解析与优化》Redis作为一款高性能的内存数据库,广泛应用于缓存、消息队列、实时统计等场景,本文将深入探讨Redis的读写并发能力,感兴趣的小伙伴可以了解下... 目录引言一、Redis 并发能力概述1.1 Redis 的读写性能1.2 影响 Redis 并发能力的因素二、

使用国内镜像源优化pip install下载的方法步骤

《使用国内镜像源优化pipinstall下载的方法步骤》在Python开发中,pip是一个不可或缺的工具,用于安装和管理Python包,然而,由于默认的PyPI服务器位于国外,国内用户在安装依赖时可... 目录引言1. 为什么需要国内镜像源?2. 常用的国内镜像源3. 临时使用国内镜像源4. 永久配置国内镜

C#原型模式之如何通过克隆对象来优化创建过程

《C#原型模式之如何通过克隆对象来优化创建过程》原型模式是一种创建型设计模式,通过克隆现有对象来创建新对象,避免重复的创建成本和复杂的初始化过程,它适用于对象创建过程复杂、需要大量相似对象或避免重复初... 目录什么是原型模式?原型模式的工作原理C#中如何实现原型模式?1. 定义原型接口2. 实现原型接口3