查表演算法，無疑也是一種非常常用、有效而且快捷的演算法，我們在很多演算法的加速程序中都能看到他的影子，在影像處理中，尤其常用，比如我們常見的各種基于直方圖的增強，可以說，在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 =https://www.cnblogs.com/Imageshop/p/= 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 = https://www.cnblogs.com/Imageshop/p/(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 =https://www.cnblogs.com/Imageshop/p/ _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 =https://www.cnblogs.com/Imageshop/p/ _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 = https://www.cnblogs.com/Imageshop/p/_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你無法運行，
如果想時刻關注本人的最新文章，也可關注公眾號：
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/514096.html
標籤：其他