AVX图像算法优化系列二: 使用AVX2指令集加速查表算法。

　　查表算法，无疑也是一种非常常用、有效而且快捷的算法，我们在很多算法的加速过程中都能看到他的影子，在图像处理中，尤其常用，比如我们常见的各种基于直方图的增强，可以说，在photoshop中的调整菜单里80%的算法都是用的查表，因为他最终就是用的曲线调整。

　　普通的查表就是提前建立一个表，然后在执行过程中算法计算出一个索引值，从表中查询索引对应的表值，并赋值给目标地址，比如我们常用的曲线算法如下所示：

int IM_Curve_PureC(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, unsigned char *TableB, unsigned char *TableG, unsigned char *TableR)
{int Channel = Stride / Width;if (Channel == 1){for (int Y = 0; Y < Height; Y++){unsigned char *LinePS = Src + Y * Stride;unsigned char *LinePD = Dest + Y * Stride;for (int X = 0; X < Width; X++){LinePD[X] = TableB[LinePS[X]];}}}else if (Channel == 3){for (int Y = 0; Y < Height; Y++){unsigned char *LinePS = Src + Y * Stride;unsigned char *LinePD = Dest + Y * Stride;for (int X = 0; X < Width; X++){LinePD[0] = TableB[LinePS[0]];LinePD[1] = TableG[LinePS[1]];LinePD[2] = TableR[LinePS[2]];LinePS += 3;LinePD += 3;}}}return IM_STATUS_OK;
}

　　通常我们认为这样的算法是很高效的，当然，我们其实还可以做一定的优化，比如使用下面的四路并行：

int IM_Curve_PureC(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, unsigned char *TableB, unsigned char *TableG, unsigned char *TableR)
{int Channel = Stride / Width;if ((Channel != 1) && (Channel != 3))                        return IM_STATUS_INVALIDPARAMETER;if ((Src == NULL) || (Dest == NULL))                        return IM_STATUS_NULLREFRENCE;if ((Width <= 0) || (Height <= 0))                            return IM_STATUS_INVALIDPARAMETER;int BlockSize = 4, Block = Width / BlockSize;if (Channel == 1){for (int Y = 0; Y < Height; Y++){unsigned char *LinePS = Src + Y * Stride;unsigned char *LinePD = Dest + Y * Stride;for (int X = 0; X < Block * BlockSize; X += BlockSize){LinePD[X + 0] = TableB[LinePS[X + 0]];LinePD[X + 1] = TableB[LinePS[X + 1]];LinePD[X + 2] = TableB[LinePS[X + 2]];LinePD[X + 3] = TableB[LinePS[X + 3]];}for (int X = Block * BlockSize; X < Width; X++){LinePD[X] = TableB[LinePS[X]];}}}else if (Channel == 3){for (int Y = 0; Y < Height; Y++){unsigned char *LinePS = Src + Y * Stride;unsigned char *LinePD = Dest + Y * Stride;for (int X = 0; X < Block * BlockSize; X += BlockSize){LinePD[0] = TableB[LinePS[0]];LinePD[1] = TableG[LinePS[1]];LinePD[2] = TableR[LinePS[2]];LinePD[3] = TableB[LinePS[3]];LinePD[4] = TableG[LinePS[4]];LinePD[5] = TableR[LinePS[5]];LinePD[6] = TableB[LinePS[6]];LinePD[7] = TableG[LinePS[7]];LinePD[8] = TableR[LinePS[8]];LinePD[9] = TableB[LinePS[9]];LinePD[10] = TableG[LinePS[10]];LinePD[11] = TableR[LinePS[11]];LinePS += 12;LinePD += 12;}for (int X = Block * BlockSize; X < Width; X++){LinePD[0] = TableB[LinePS[0]];LinePD[1] = TableG[LinePS[1]];LinePD[2] = TableR[LinePS[2]];LinePS += 3;LinePD += 3;}}}return IM_STATUS_OK;
}

　　这样效率能进一步的提高。

　　在早期我们的关注中，我也一直想再次提高这个算法的效率，但是一直因为他太简单了，而无法有进一步的提高，在使用SSE指令集时，我们也没有找到合适的指令，只有当查找表为16字节的表时，可以使用_mm_shuffle_epi8快速实现，详见【算法随记七】巧用SIMD指令实现急速的字节流按位反转算法。一文的描述。

　　在我们再次接触AVX指令集，正如上一篇关于AVX指令的文章所述，他增加了非常具有特色的gather系列指令，具体有哪些如下图所示：

　　有一大堆啊，其实看明白了，就只有2大类，每大类里有2个小系列，每个系列里又有4中数据类型，

　　两大类为：针对128位的类型的gather和针对256位的gather。

　　两个系列为：带mask和不带mask系列。

　　4中数据类型为： int32、int64、float、double。

　　当然，里面还有一些64为地址和32位地址的区别，因此又增加了一些列的东西，我个人认为其中最常用的函数只有4个，分别是：_mm_i32gather_epi32 、_mm256_i32gather_epi32、_mm_i32gather_ps、_mm256_i32gather_ps，我们以_mm256_i32gather_epi32为例。

　　注意，这里所以下，不要以为_mm_i32gather_ps这样的intrinsics指令以_mm开头，他就是属于SSE的指令，实际行他并不是，他是属于AVX2的，只是高级别的指令集对老指令的有效补充。

　　_mm256_i32gather_epi32的相关说明如下：

　　其作用，翻译过来就是从固定的基地址base_addr开始，燃用偏移量由 vindex提供，注意这里的vindex是一个__m256i数据类型，里面的数据要把它看成8个int32类型，即保存了8个数据的地址偏移量，最后一个scale表示地址偏移量的放大系数，容许的值只有1、2、4、8，代表了字节，双字节，四字节和把字节的意思，通常_mm256_i32gather_epi32一般都是使用的4这个数据。

　　那么注意看这些gather函数，最下的操作单位都是int32，因此，如果我们的查找表是byte或者short类型，这个就有点困难了，正如我们上面的Cure函数一样，是无法直接使用这个函数的。

　　那么我我们来看看一个正常的int型表，使用两者之间大概有什么区别呢，以及是如何使用该函数的，为了测试公平，我把正常的查找表也做了展开。

int main()
{const int Length = 4000 * 4000;int *Src = (int *)calloc(Length, sizeof(int));int *Dest = (int *)calloc(Length, sizeof(int));int *Table = (int *)calloc(65536, sizeof(int));for (int Y = 0; Y < Length; Y++)        Src[Y] = rand();    //    产生的随机数在0-65535之间，正好符号前面表的大小for (int Y = 0; Y < 65536; Y++){Table[Y] = 65535 - Y;    //    随意的分配一些数据}LARGE_INTEGER nFreq;//LARGE_INTEGER在64位系统中是LONGLONG，在32位系统中是高低两个32位的LONG，在windows.h中通过预编译宏作定义LARGE_INTEGER nBeginTime;//记录开始时的计数器的值LARGE_INTEGER nEndTime;//记录停止时的计数器的值double time;QueryPerformanceFrequency(&nFreq);//获取系统时钟频率QueryPerformanceCounter(&nBeginTime);//获取开始时刻计数值for (int Y = 0; Y < Length; Y += 4){Dest[Y + 0] = Table[Src[Y + 0]];Dest[Y + 1] = Table[Src[Y + 1]];Dest[Y + 2] = Table[Src[Y + 2]];Dest[Y + 3] = Table[Src[Y + 3]];}QueryPerformanceCounter(&nEndTime);//获取停止时刻计数值time = (double)(nEndTime.QuadPart - nBeginTime.QuadPart) * 1000 / (double)nFreq.QuadPart;//(开始-停止)/频率即为秒数，精确到小数点后6位printf("%f   \n", time);QueryPerformanceCounter(&nBeginTime);//获取开始时刻计数值for (int Y = 0; Y < Length; Y += 16){__m256i Index0 = _mm256_loadu_si256((__m256i *)(Src + Y));__m256i Index1 = _mm256_loadu_si256((__m256i *)(Src + Y + 8));__m256i Value0 = _mm256_i32gather_epi32(Table, Index0, 4);    __m256i Value1 = _mm256_i32gather_epi32(Table, Index1, 4);_mm256_storeu_si256((__m256i *)(Dest + Y), Value0);_mm256_storeu_si256((__m256i *)(Dest + Y + 8), Value1);}QueryPerformanceCounter(&nEndTime);//获取停止时刻计数值time = (double)(nEndTime.QuadPart - nBeginTime.QuadPart) * 1000 / (double)nFreq.QuadPart;//(开始-停止)/频率即为秒数，精确到小数点后6位printf("%f   \n", time);free(Src);free(Dest);free(Table);getchar();return 0;
}

　　直接使用这句即可完成查表工作：__m256i Value0 = _mm256_i32gather_epi32(Table, Index0, 4);

　　这是一个比较简单的应用场景，在我本机的测试中，普通C语言的耗时大概是27ms，AVX版本的算法那耗时大概是17ms，速度有1/3的提升。考虑到加载内存和保存数据在本代码中占用的比重明显较大，因此，提速还是相当明显的。

　　我们回到刚才的关于Curve函数的应用，因为gather相关指令最小的收集粒度都是32位，因此，对于字节版本的表是无论为力的，但是为了能借用这个函数实现查表，我们可以稍微对输入的参数做些手续，再次构造一个int类型的表格，即使用如下代码(弧度版本，Channel == 1)：

int Table[256];
for (int Y = 0; Y < 256; Y++)
{Table[Y] = TableB[Y];
}

　　这样这个表就可以用了，对于24位我们也可以用类似的方式构架一个256*3个int元素的表。

　　但是我们又面临着另外一个问题，即_mm256_i32gather_epi32这个返回的是8个int32类型的整形数，而我们需要的返回值确实字节数，所以这里就又涉及到8个int32数据转换为8个字节数并保存的问题，当然为了更为高效的利用指令集，我们这里考虑同时把2个__m256i类型里的16个int32数据同时转换为16个字节数，这个可以用如下的代码高效的实现：

for (int Y = 0; Y < Height; Y++)
{unsigned char *LinePS = Src + Y * Stride;unsigned char *LinePD = Dest + Y * Stride;for (int X = 0; X < Block * BlockSize; X += BlockSize){__m128i SrcV = _mm_loadu_si128((__m128i *)(LinePS + X));//    int32    A0    A1    A2    A3    A4    A5    A6    A7__m256i ValueL = _mm256_i32gather_epi32(Table, _mm256_cvtepu8_epi32(SrcV), 4);//    int32    B0    B1    B2    B3    B4    B5    B6    B7__m256i ValueH = _mm256_i32gather_epi32(Table, _mm256_cvtepu8_epi32(_mm_srli_si128(SrcV, 8)), 4);//    short    A0    A1    A2    A3    B0    B1    B2    B3    A4    A5    A6    A7    B4    B5    B6    B7__m256i Value = _mm256_packs_epi32(ValueL, ValueH);//    byte    A0    A1    A2    A3    B0    B1    B2    B3    0    0    0    0    0    0    0    0    A4    A5    A6    A7    B4    B5    B6    B7        0    0    0    0    0    0    0    0    Value = _mm256_packus_epi16(Value, _mm256_setzero_si256());//    byte    A0    A1    A2    A3    A4    A5    A6    A7    B0    B1    B2    B3    B4    B5    B6    B7    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    Value = _mm256_permutevar8x32_epi32(Value, _mm256_setr_epi32(0, 4, 1, 5, 2, 3, 6, 7));_mm_storeu_si128((__m128i *)(LinePD + X), _mm256_castsi256_si128(Value));}for (int X = Block * BlockSize; X < Width; X++){LinePD[X] = TableB[LinePS[X]];}

　上面的代码里涉及到了没有按常规方式出牌的_mm256_packs_epi32、_mm256_packus_epi16等等，最后我们也是需要借助于AVX2提供的_mm256_permutevar8x32_epi32才能把那些数据正确的调整为需要的格式。

　　对于彩色的图像，就要稍微复杂一些了，因为涉及到RGB格式的排布，同时考虑一些对齐问题，最友好的方式就是一次性处理8个像素，24个字节，这一部分留给有兴趣的读者自行研究。

　　在我本机的CPU中测试呢，灰度版本的查找表大概有20%的提速，彩色版本的要稍微多一些，大概有30%左右。

　　这些提速其实不太明显，因为在整个过程中处理内存耗时较多，他并不是以计算为主要过程的算法，当我们某个算法中见也有查找时，并且为了计算查找表时，需要很多的数学运算去进行隐射的坐标计算时，这个时候这些隐射计算通常都是有浮点参与，或其他各种复杂的计算参与，这个时候用SIMD指令计算这些过程是能起到很大的加速作用的，在我们没有AVX2之前，使用SSE实现时，到了进行查表时通常的做法都是把前通过SSE计算得到的坐标的_m128i元素的每个值使用_mm_extract_epi32(这个是内在的SSE指令，不是用其他伪指令拼合的)提取出每个坐标值，然后在使用_mm_set相关的函数把查找表的返回值拼接成一个行的SSE变量，以便进行后续的计算，比如下面的代码：

　　这个时候使用AVX2的这个指令就方便了，如下所示：

　　注意到上面的Texture其实是个字节类型的数组，也就是一副图像，对应的C代码如下所示：

int SampleXF = IM_ClampI(ClipXF >> 16, 0, Width - 1);            //    试着拆分VX和VY的符号情况分开写，减少ClampI的次数，结果似乎区别不是特别大，因此优化意义不大
int SampleXB = IM_ClampI(ClipXB >> 16, 0, Width - 1);
int SampleYF = IM_ClampI(ClipYF >> 16, 0, Height - 1);
int SampleYB = IM_ClampI(ClipYB >> 16, 0, Height - 1);
unsigned char *SampleF = Texture + (SampleYF * Stride + SampleXF);
unsigned char *SampleB = Texture + (SampleYB * Stride + SampleXB);
Sum += SampleF[0] + SampleB[0];

　　可见这里实际上是对字节类型进行查表，所以这里最后的那个scale参数我们取的是1，即中间的偏移是以字节为单位的，但是这里其实隐含着一个问题，即如果我们取样的是图片最右下角的那个位置的像素，因为要从那个位置开始读取四个字节的内存，除非图像原始格式是BGRA的，否则，必然会读取到超出图像内存外的内存数据，这个在普通的C语言中，已改会弹出一个系统错误框，蹦的一下说访问非法内存，但是我看用这个指令似乎目前还没有遇到这个错误，哪怕认为的输入一个会犯错误的坐标。

　　如果是这样的话，得到的一个好处就是对于那些图像扭曲滤镜、缩放图像中哪些重新计算坐标的函数来说，不用临时构建一副同样数据的int类型图了，而可以直接放心的使用这个函数了。

　　最后说明一点，经过在其他一些机器上测试，似乎有些初代即使支持AVX2的CPU，使用这些函数后相应的算法的执行速度反而有下降的可能性，不知道为什么。

　　在我提供的SIMD指令优化的DEMO中，在 Adjust-->Exposure菜单下可以看到使用C语言和使用AVX进行查表优化的功能，有兴趣的作者可以自行比较下。

　　很明显，在这里SSE优化选项是无法使用的。

本文可执行Demo下载地址： https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar，菜单中蓝色字体显示的部分为已经使用AVX加速的算法，如果您的硬件中不支持AVX2，可能这个DEMO你无法运行。

如果想时刻关注本人的最新文章，也可关注公众号：

AVX图像算法优化系列二: 使用AVX2指令集加速查表算法。相关推荐

SSE图像算法优化系列二十一：基于DCT变换图像去噪算法的进一步优化（100W像素30ms）。...
在优化IPOL网站中基于DCT(离散余弦变换)的图像去噪算法(附源代码) 一文中,我们曾经优化过基于DCT变换的图像去噪算法,在那文所提供的Demo中,处理一副1000*1000左右的灰度噪音图像耗时 ...
SSE图像算法优化系列二十二：优化龚元浩博士的曲率滤波算法，达到约1000 MPixels/Sec的单次迭代速度...
2015年龚博士的曲率滤波算法刚出来的时候,在图像处理界也曾引起不小的轰动,特别是其所说的算法的简洁性,以及算法的效果.执行效率等方面较其他算法均有一定的优势,我在该算法刚出来时也曾经有关注,不过 ...
SSE图像算法优化系列二十五:二值图像的Euclidean distance map（EDM)特征图计算及其优化。...
Euclidean distance map(EDM)这个概念可能听过的人也很少,其主要是用在二值图像中,作为一个很有效的中间处理手段存在.一般的处理都是将灰度图处理成二值图或者一个二值图处理成另外一 ...
SSE图像算法优化系列二十九：基础的拉普拉斯金字塔融合用于改善图像增强中易出现的过增强问题（一）...
拉普拉斯金字塔融合是多图融合相关算法里最简单和最容易实现的一种,我们在看网络上大部分的文章都是在拿那个苹果和橙子融合在一起,变成一个果橙的效果作为例子说明.在这方面确实融合的比较好.但是本文我们主要讲 ...
SSE图像算法优化系列十四：局部均方差及局部平方差算法的优化。
关于局部均方差有着较为广泛的应用,在我博客的基于局部均方差相关信息的图像去噪及其在实时磨皮美容算法中的应用及使用局部标准差实现图像的局部对比度增强算法中都有谈及,即可以用于去噪也可以用来增强图像,但是 ...
SSE图像算法优化系列十九：一种局部Gamma校正对比度增强算法及其SSE优化。
这是一篇2010年比较古老的文章了,是在QQ群里一位群友提到的,无聊下载看了下,其实也没有啥高深的理论,抽空实现了下,虽然不高大上,还是花了点时间和心思优化了代码,既然这样,就顺便分享下优化的思路和经 ...
SSE图像算法优化系列十六：经典USM锐化中的分支判断语句SSE实现的几种方法尝试。...
分支判断的语句一般来说是不太适合进行SSE优化的,因为他会破坏代码的并行性,但是也不是所有的都是这样的,在合适的场景中运用SSE还是能对分支预测进行一定的优化的,我们这里以某一个算法的部分代码为例进行 ...
SSE图像算法优化系列八：自然饱和度(Vibrance)算法的模拟实现及其SSE优化（附源码，可作为SSE图像入门，Vibrance算法也可用于简单的肤色调整）。...
Vibrance这个单词搜索翻译一般振动,抖动或者是响亮.活力,但是官方的词汇里还从来未出现过自然饱和度这个词,也不知道当时的Adobe中文翻译人员怎么会这样处理.但是我们看看PS对这个功能的解释: ...
SSE图像算法优化系列一：一段BGR2Y的SIMD代码解析。
一个同事在github上淘到一个基于SIMD的RGB转Y(彩色转灰度或者转明度)的代码,我抽了点时间看了下,顺便学习了一些SIMD指令,这里把学习过程中的一些理解和认识共享给大家. github上相关 ...
ElasticSearch优化系列二：机器设置（内存）
点击"阅读原文"直接打开[北京站 | GPU CUDA 进阶课程]报名链接预留一半内存给Lucence使用一个常见的问题是配置堆太大.你有一个64 GB的机器,觉得JVM内存越 ...

AVX图像算法优化系列二: 使用AVX2指令集加速查表算法。

AVX图像算法优化系列二: 使用AVX2指令集加速查表算法。相关推荐

最新文章

热门文章