透過 SIMD 加速高斯模糊運算

contributed by <CheHsuan>, <kevinbird61>, <SwimGlass> , <carolc0708>

tags: `SIMD` `sysprog`

預期目標

延續開發紀錄: 使用 SSE/AVX 對圖檔作高斯模糊處理和 SIMD + prefetching 的成果，提供 Intel 架構 SSE/AVX 和 ARM 架構 NEON 的效能最佳化實做
整理相關共筆
簡報連結
影片連結(已更新)

程式碼產出: Image-Processing

執行環境

Architecture: x86_64
CPU 作業模式： 32-bit, 64-bit
Byte Order: Little Endian
CPU(s): 8
On-line CPU(s) list: 0-7
每核心執行緒數：2
每通訊端核心數：4
Socket(s): 1
NUMA 節點： 1
供應商識別號： GenuineIntel
CPU 家族： 6
型號： 60
Model name: Intel® Core™ i7-4700HQ CPU @ 2.40GHz
製程： 3
CPU MHz： 3199.968
CPU max MHz: 3400.0000
CPU min MHz: 800.0000
BogoMIPS: 4789.25
虛擬： VT-x
L1d 快取： 32K
L1i 快取： 32K
L2 快取： 256K
L3 快取： 6144K
NUMA node0 CPU(s): 0-7

Description: Ubuntu 16.04.1 LTS
Release: 16.04
Codename: xenial

高斯模糊

高斯模糊公式

高斯模糊，又稱為高斯平滑，主要用於降低圖片細節與層次。
如上圖，我們以一點作為運算目標，距離該點愈遠的，所佔的影響值愈小。
而權重分佈則採常態分佈（又稱正態分佈、高斯分佈）方式，距離愈近 μ ，值愈大；而距離（ σ:標準差）愈小，則代表愈集中、愈大則代表愈分散
而高斯分佈公式：

$f (x) = \frac{1}{σ \sqrt{2 π}} e^{\frac{- (x - μ)^{2}}{2 σ^{2}}} .$
而高斯模糊公式（二維空間）：

$G (u, v) = \frac{1}{2 π σ^{2}} e^{\frac{- (u^{2} + v^{2})}{2 σ^{2}}} .$
這邊 r 代表模糊半徑 r² = u² + v²，而從高斯分佈公式中，可以看到當距離大於 3*σ 後，基本上已經無影響力

請用精確的表示法，避免用 r^2，而該用 r²，或者 LaTeX 數學表達式 jserv

而此，我們可以生成以距離為主的算式：
以 u = 2 , v = 2 為例，我們可以以一個 5 x 5 的矩陣來作我們的 Gaussian 權重值
藉由加總所有矩陣值並調整到為 1 後，便以他為我們的權重值，作為 Gaussian Kernel 來使用
而我們的實例中，以normalize後的 5 x 5 Gaussian kernel 為主：








Gaussian_kernel[25] = {
    1,4,7,4,1,
    4,16,26,16,4,
    7,26,41,26,7,
    4,16,26,16,4,
    1,4,7,4,1
 }
// 並再相乘後加總，除上 273，作為我們該點的 blur 值

實作 - 兩種資料結構差異

我選擇了沒經過壓縮的圖片檔格式 -bmp 作為實驗目標
- 主要我認為沒有壓縮過的圖片來操作，相對於已經做過壓縮後的（或失真），效果比較明顯
以影像處理為實驗方式，進而檢測效能
主要實作中，分為原本的資料結構和分開的結構來實作：










// Original structure
typedef struct tagRGBTRIPLE { // (3bytes)
    BYTE rgbBlue;             // (1bytes) blue channel
    BYTE rgbGreen;            // (1bytes) green channel
    BYTE rgbRed;              // (1bytes) red channel
} RGBTRIPLE;
// split structure
color_r = (unsigned char*)malloc(bmpInfo.biWidth*bmpInfo.biHeight*sizeof(unsigned char));
color_g = (unsigned char*)malloc(bmpInfo.biWidth*bmpInfo.biHeight*sizeof(unsigned char));
color_b = (unsigned char*)malloc(bmpInfo.biWidth*bmpInfo.biHeight*sizeof(unsigned char));

執行這兩種結構來看執行效率：

	naive(ori)	naive(tri)
高斯模糊	0.395602	0.322005

可以看到，由於以 original data structure 只執行了一次 blur function，比較 split data structure 來說，相對少了兩次的 blur function call ；而兩種實作方式就帶來0.073597秒的差距
而為何我要實作兩種，主要考慮到之後實作 SSE/AVX instruction set 時所需要的操作，可能會為了符合一個 register 大小（128/256）而做出調整，儘量使用到所有空間，減少回圈執行次數。

實作: SSE/AVX

考慮到高斯模糊的實作方式：3x3 跟 5x5，以及考慮 data structure 的情況：
- 單行：3=> split: 31 = 3 bytes(24bits)、 original: 33 = 9 bytes (72 bits) ; 5 => split: 51 = 5 bytes(40 bits) 、 original: 53 = 15 bytes(120 bits)
- ( [分散] 128bits 為例)以一次 load 進一行來說，3x3可以用(n)*8 = 128 => n = 16個元素；而對操作而言，一次需要用到3個，所以一次的指令中，我們可以對16-(3-1)=14次 gaussian blur ；而如果是5x5的，再一次的指令中，則可以有16-(5-1)=12次的 gaussian blur！！
- ( [原始] 128bits 為例)以一次 load 進一行來說，3x3 or 5x5 用到(n)83 = 128 => n = 16/3 => 會有額外的浪費！ => 先實作分散版本
我們先 Load 進5行(每行16 bytes)的大小來做，找出其中5x5的元素乘上高斯模糊係數，並且加總後放回原先的位置
考慮到使用 SSE 的效益，先實作看看unroll版本

當初會想選擇試試看的原因是對於圖檔處理時，必須對大面積的資料作處理，我就想到 SSE 的操作，可以一次 load 進大量的資料來做我需要的操作；但一開始時，由於沒有很有效率的使用，導致效率奇差（比原先還要慢）；後來更改後，增加使用率的版本，便可以看出效能增加不少。

而 unroll 跟 SSE 的關係我是以兩者皆有消除兩次回圈的操作再，所以有比較一下兩者差異。kevinbird61

實作 - unroll version

unrolling 應該搭配彈性高的 macro，必要時，比照 raytracing，寫出好的 code generator jserv

– 分析執行結果：經過更改後的 makefile ，可以讓使用者決定要對該圖檔執行幾次的 Gaussian blur ；而比較原本 5x5 分開的 structure , original 來做，連續執行 3 次前後效能比較：

	naive Split	naive Original	naive Unroll split	naive Unroll original
執行時間(ms)	419.833675	332.680654	188.304807	254.358430

可以看到，unroll 後的版本都比原先的實作版本快上許多（尤其以 split 的實作更為明顯）；由於原先實作中，split 版本由於必須執行三次才能完成，相對 original 版本只需要一次執行來說，相對耗費時間；而再 unroll 之後，每次執行時間減少，原先執行三次則同時代表這三次同時受到優化，可獲得3倍的執行時間減少；計算 3 次優化所減少的時間：231.528868 ms，除以 3 後大約為77.18 ms；而看到 original 前後時間差異：78.322224 ms ，差不多相同，代表著由於 unroll 對於實驗體（ 1600 x 1200 的圖片檔）來說，對一次執行實作 unroll 的效能增加差不多為 77~78 ms
而實作 original structure unroll 版本，得到的執行時間是： 0.985009 sec；竟然比原先的還要來的慢…
- PS：後來多次執行後，執行時間跟原本差不多（忽高忽低）

如何解釋這現象？又如何避免？ jserv

再後來的實驗中，所繪製的實行時間圖，original unroll 的執行時間還是小於 naive unroll 的
kevinbird61

@Champ Yen 學長提示 " 像是 thermal control(前幾次) + DVFS(throttle)造成的 "
我在散熱良好(設定為主動散熱且效能優先並且關閉了省電與熱能控制)的環境測試的結果非常平穩結果圖片將更新在最後面 shinshipower

圖為各函式對圖形執行一次 Gaussian blur 的結果發現：
SSE 的操作方式並沒有很節省空間的使用(5x5 , 3x3 的操作，每次取完能用的有限)，導致對圖檔做三次 Gaussian blur 的結果為1.842012秒
修改再同時一個 instruction 中，為了能夠一次操作更多，改成一次操作並儲存兩個位置。執行後時間反而增加:3.015885秒
由於 SSE 本身需要 load 進 16 個位置的資訊，導致每次操作時都需要花上不少指令來拿出特定的 5 個（ per line ）來做，入不敷出。每當為了運算出一個位置，就會付出 40 幾行的程式碼，以及將近 40 個 SSE function call













































...
__m128i vg1 = _mm_loadu_si128((__m128i *)sse_g1);
__m128i vg2 = _mm_loadu_si128((__m128i *)sse_g2);
__m128i vg3 = _mm_loadu_si128((__m128i *)sse_g3);
__m128i vg4 = _mm_loadu_si128((__m128i *)sse_g4);
__m128i vg5 = _mm_loadu_si128((__m128i *)sse_g5);

__m128i vsum = _mm_set1_epi8(0),vtemplow = _mm_set1_epi8(0),vtemphigh = _mm_set1_epi8(0),vempty = _mm_set1_epi8(0);
// First element src[j*w+i]
// Load in data
__m128i L0 = _mm_loadu_si128((__m128i *)(src+(j+0)*w + i));
__m128i L1 = _mm_loadu_si128((__m128i *)(src+(j+1)*w + i));
__m128i L2 = _mm_loadu_si128((__m128i *)(src+(j+2)*w + i));
__m128i L3 = _mm_loadu_si128((__m128i *)(src+(j+3)*w + i));
__m128i L4 = _mm_loadu_si128((__m128i *)(src+(j+4)*w + i));
// Get the data we need (5 element per-line) , because we only
// need 5 element from sse instruction set , so only get low part(contain 8 elements)
__m128i v0 = _mm_unpacklo_epi8(L0,vk0);
__m128i v1 = _mm_unpacklo_epi8(L1,vk0);
__m128i v2 = _mm_unpacklo_epi8(L2,vk0);
__m128i v3 = _mm_unpacklo_epi8(L3,vk0);
__m128i v4 = _mm_unpacklo_epi8(L4,vk0);
// Multiple with specific Gaussian coef.
v0 = _mm_maddubs_epi16(v0,vg1);
v1 = _mm_maddubs_epi16(v1,vg2);
v2 = _mm_maddubs_epi16(v2,vg3);
v3 = _mm_maddubs_epi16(v3,vg4);
v4 = _mm_maddubs_epi16(v4,vg5);
// Summation the 5 line
vsum = _mm_add_epi16(vsum,v0);
vsum = _mm_add_epi16(vsum,v1);
vsum = _mm_add_epi16(vsum,v2);
vsum = _mm_add_epi16(vsum,v3);
vsum = _mm_add_epi16(vsum,v4);
// Vsum summation
// Summation all - (Summation all - (Summation with shift-off 5 number))
vtemplow = _mm_unpacklo_epi16(vsum,vempty); // 1,2,3,4
vtemphigh = _mm_unpackhi_epi16(vsum,vempty); // 5

sum += _mm_cvtsi128_si32(vtemplow); // get 1
sum += _mm_cvtsi128_si32(_mm_srli_si128(vtemplow,4)); // get 2
sum += _mm_cvtsi128_si32(_mm_srli_si128(vtemplow,8)); // get 3
sum += _mm_cvtsi128_si32(_mm_srli_si128(vtemplow,12)); // get 4
sum += _mm_cvtsi128_si32(vtemphigh); // get 5
...

你的CPU時間都浪費在load了，難怪。
要顯出SSE的威力，請減少load的步驟，簡化到load->計算->移出->（重複）
例如說load 4個連續的pixel(uint32_t) -> *gaussian_blur的值 -> 移出pixel或改變gaussian_blurgnitnaw

此部份我後來有用 original structure 重新實作一次，從原先每個 instruction 使用率（ 5/16 ）到後來的（ 15/16 ）總共提升 3 倍使用率；時間上的變化（執行一次）如下：

sse split sse original

執行時間( ms ) 633.522699 243.087927

減少了2.6倍左右
kevinbird61

據我所知，SSE只有8個128bit暫存器(還是Core i7不只8個？)，不過看你sse_gaussian_blur_5_ori裏面就不只load 8個了。
應該可以再節省點用。提示：一次load數個pixel(input)進SSE=>一次load數個對應的pixel(output)=>SSE裡的input乘以已經在SSE的gaussian55(雖說是25個但是只有6個不重複)加入在SSE裡的output=>SSE裡的output輸出到memory裡的output(記得要做完25個點)=>重複load
再多個提示：_mm_mullo_epi16 or _mm_mullo_epi32
gnitnaw

實作: POSIX Thread

關於分頭進行的想法，主要是考慮到每次進行以 5x5 矩陣操作，那麼我便可以以執行緒的數量來做分配，每個執行不同的區塊，大小為： height / threadcounts；如此一來我們便可以讓這些執行續同時操作整塊圖片，加快執行

[name=kevinbird61

考慮到可以分頭進行
- 實作 pthread 來分擔 Gaussian blur ，依據不同的 thread_id * workload（由整張圖片高度除上總共thread）作為不同 thread 起始操作位置來做；而一樣對同一張圖作三次高斯模糊（unroll split structure），執行時間為：0.402466秒（ 4 條 thread）
- 重新繪製執行時間圖如下(執行一次 blur + 4 threads)：

修改SSE部份的實作機制

為什麼要修改原本實作的原因還是再題一下比較好。

太多東西省略掉以至於閱讀性待加強。我甚至不確定這邊是講改善的方式和對應手段還是文字描述修改的演算法。建議先講你的改善目的再講改善細節。

更改取值的原因？

看起來有幾次的嘗試和結果，要不要分開來說明？
因為這領域不熟，不確定這樣是否是正常，請和老師確認
wen

修改實作方式 : 原先架構由於是 split ，為分別對 r, g , b 3個陣列分開做處理，所以問題在於一次 load 進的 16 bytes 的空間，若沒有作一個完整的利用，會導致大量 load 指令的耗費，r , g , b 3個陣列皆是一個元素 1 byte，若不改寫的情況下，對 5x5 的操作非常的傷；而修改的方式，是改以 original 的方式來操作，一次也是對 5x5 做操作，但由於 original 一個元素為 3 bytes，同時包含 R,G,B 再其中，所以使用率可以從原先 5/16 到 15/16，提升為 3 倍。這樣一來雖然load的數量不變，執行的總迴圈次數從原先 split 的三次（分別對 r,g,b 操作）到一次（直接對 original 做操作）

修改取值的方式： 主要我原先認為我最後要取出結果的操作(_mm_cvtsi128_si32 配合 _mm_srli_si128)，每次呼叫就會讓他對原本內部 128 bits的資料作操作，會不會讓執行時間上升？所以我就想說如果先把值取出來，存進一個佔有 4 個 integer( = 16 bytes)的陣列，之後對該陣列操作即可；但實作後效率並沒有比較好的表現。
目前我是覺得 _mm_storeu_si128 可能再實作上，比使用(_mm_cvtsi128_si32 配合 _mm_srli_si128)來的耗費 cpu cycles，還有更改後的作法會需要額外的空間（4 個 16 bytes的陣列）來做儲存，並再存回去的時候又有對各陣列取值的動作。
這些應該都是導致修改後執行時間增加的原因kevinbird61

由於一次 load 進 16 bytes，考慮到原本資料結構為 RGBTRIPLE（大小為3 bytes ），如果一次 load 進 5 個元素，就是 5*3 = 15 bytes ，比起原先 16 個中只能 load 5 個 bytes 來說，更加利用到空間
把原先 BGRBGRBG… 排列得資料利用 _mm_unpacklo_epi8加上為 0 的 vector 來輸出： B0G0R0B0G0R0…的樣式，此時由原先的一個 __m128i 的變數，變成兩個（總共從 5 個變成 10 個）；接著利用_mm_maddubs_epi16，變形後的來跟變形的 gaussian kernel 矩陣做相乘的動作，並利用_mm_add_epi16依序兩兩相加；最後取得兩個 __m128i 的變數，代表了所有項的總和值；為了能夠使用_mm_cvtsi128_si32來取得[32 : 0]的值，我們先利用_mm_unpacklo_epi16把原先[ Color ][ 0 ][ 0 ][ Color ][ 0 ][ 0 ]…擴展成[ Color ][ 0 ][ 0 ][ 0 ][ 0 ][ Color ][ 0 ]…的方式（每個單元以32bits為單位做操作）；最後再利用_mm_cvtsi128_si32配合_mm_srli_si128作移位，一一取得交錯的 R 、 G 、 B 的值來加總，並儲存回圖片位置。
而依照記憶體位置來看，這15 bytes 分別會以 B->G->R 的順序作排列直到底為止
實作後，可以看到變化如下：
可以從紅色的曲線（ sse original ）來看到，基本上已經超越大部分的實作
由於改成這個處理方式後，把原本 split 時所需要處理三次回圈的部份改成只需要處理一次（ r,g,b 同時做）；這樣我們再處理上，就可以使用特化過後的 gaussian kernel 來幫每個移動過後的元素作相乘得動作。最後再依序取出並加總變可以得到最後該位置的值！
更改取值機制，從原本：



sum_b += _mm_cvtsi128_si32(vtemp_1) + _mm_cvtsi128_si32(_mm_srli_si128(vtemp_1,12)) + _mm_cvtsi128_si32(_mm_srli_si128(vtemp_2,8)) + _mm_cvtsi128_si32(_mm_srli_si128(vtemp_3,4)) + _mm_cvtsi128_si32(vtemp_4);
sum_g += _mm_cvtsi128_si32(_mm_srli_si128(vtemp_1,4)) + _mm_cvtsi128_si32(vtemp_2) + _mm_cvtsi128_si32(_mm_srli_si128(vtemp_2,12)) + _mm_cvtsi128_si32(_mm_srli_si128(vtemp_3,8)) + _mm_cvtsi128_si32(_mm_srli_si128(vtemp_4,4));
sum_r += _mm_cvtsi128_si32(_mm_srli_si128(vtemp_1,8)) + _mm_cvtsi128_si32(_mm_srli_si128(vtemp_2,4)) + _mm_cvtsi128_si32(vtemp_3) + _mm_cvtsi128_si32(_mm_srli_si128(vtemp_3,12)) + _mm_cvtsi128_si32(_mm_srli_si128(vtemp_4,8));

到新的嘗試方式：










int temp1[4]= {0,0,0,0},temp2[4]= {0,0,0,0} ,temp3[4]= {0,0,0,0},temp4[4]= {0,0,0,0};

_mm_storeu_si128((__m128i *)temp1,vtemp_1); // B G R B
_mm_storeu_si128((__m128i *)temp2,vtemp_2); // G R B G
_mm_storeu_si128((__m128i *)temp3,vtemp_3); // R B G R
_mm_storeu_si128((__m128i *)temp4,vtemp_4); // B G R

sum_b += temp1[0] + temp1[3] + temp2[2] + temp3[1] + temp4[0];
sum_g += temp1[1] + temp2[0] + temp2[3] + temp3[2] + temp4[1];
sum_r += temp1[2] + temp2[1] + temp3[0] + temp3[3] + temp4[2];

察看是否執行速度上升，使用該函式對圖片操作三次高斯模糊；結果執行前時間為 793.782715 ms ；更改後執行時間變為 853.281916 ms，可能_mm_storeu_si128對於整體執行時間比原先 shift 加上取出最後 32 bits 的版本要來的耗費時間
後續：加上 prefetch 後，發現沒有比較快

實作: Pthread + SSE

使用 pthread + SSE 來再次實作=>
可以看到，執行時間幾乎壓再 0.1 秒左右（紫色三角）！！整整比起原先（黃色方形）來的快了 0.2 秒左右！！

修正繪圖曲線抖動問題

先前測試時，都有開啟其他軟體使用；關閉後做測試：
- 只有開頭部份穩定，後面執行仍然有振盪

曲線推論與測試 by @shinshipower

@Champ Yen 學長提示此狀況可能與CPU TDP控制與散熱控制有關

這 execution time 的 variation pattern 看來很像是 thermal control(前幾次) + DVFS(throttle 之後) 造成的
但是實驗沒交代實驗環境軟硬體平台及參數, 這在解讀數據很重要, 跟做物理實驗時需要把變因定下來並且做陳述一樣等義… Champ Yen @ Facebook

使用平台 Dell PowerEDGE R530

Architecture: CPU op-mode(s): Byte Order: CPU(s): On-line CPU(s) list: Thread(s) per core: Core(s) per socket: Socket(s): NUMA node(s): Vendor ID: CPU family: Model: Model name: Stepping: CPU MHz: CPU max MHz: CPU min MHz: BogoMIPS: Virtualization: L1d cache: L1i cache: L2 cache: L3 cache: NUMA node0 CPU(s): Flags: x86_64
32-bit, 64-bit
Little Endian
24
0-23
2
12
1
1
GenuineIntel
6
79
Intel® Xeon® CPU E5-2650 v4 @ 2.20GHz
1
1200.890
2900.0000
1200.0000
4394.84
VT-x
32K
32K
256K
30720K
0-23
fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc aperfmperf eagerfpu pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid dca sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch epb intel_pt tpr_shadow vnmi flexpriority ept vpid fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm cqm rdseed adx smap xsaveopt cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local dtherm ida arat pln pts

gcc (GCC) 6.2.1 20160916 (Red Hat 6.2.1-2)

在bios設定為效能優先且主動散熱的情況下使用

make perf_time

設定執行一次高斯運算與 24 Threads的情況下
得到以下perf結果

Enter the times you want to execute Gaussian blur on the input picture:1
Enter the thread number: 24

 Performance counter stats for './bmpreader img/input.bmp output.bmp 1 24' (100 runs):

            10,985      cache-misses:u            #    0.393 % of all cache refs      ( +-  2.11% )
         2,796,773      cache-references:u                                            ( +-  0.61% )

       3.178636533 seconds time elapsed                                          ( +-  0.25% )

以及 plot 出的對應結果圖
image alt

而kevinbird61的執行結果：

Enter the times you want to execute Gaussian blur on the input picture:1
Enter the thread number: 4

 Performance counter stats for './bmpreader img/input.bmp output.bmp 1 4' (100 runs):

         1,027,242      cache-misses              #   43.085 % of all cache refs      ( +-  2.07% )
         2,384,239      cache-references                                              ( +-  2.05% )

       3.284249783 seconds time elapsed                                          ( +-  1.36% )

相關探討與建議:
- 經詢問與查詢瞿旭民同學的 cpu 應為 I7 4700HQ

- 在核心數、快取頻寬、快取大小、指令架構、記憶體頻寬都大幅度提升的情況下
- 兩者跑出來的結果居然非常接近，甚至某些指標還略微下降
- 經比較兩 CPU 的其他更多 benchmark 後我認為此程式與測試方式只吃重單核心
- 就連 Pthread 的版本也是，可能思考方向如下:
- 1. 使用 pthread 優化的地方是否根本並非熱區
- 2. 使用 pthread 的方法是否不佳?
- 3. cache 的使用?

實作1D Gaussian Filter

原本我們使用二維的 gaussian filter來去做運算，這樣一來，導致每個pixel都要經過 n * n次乘法運算，因此應該藉由一些方法來降低運算量。由於 gaussian filter是一個 separable matrix，因此原本2D的 filter 可以轉換成先對 X 方向做一次1D的gaussian filter，然後再對 Y 方向做一次1D的gaussian filter，就可以達成同樣的效果。而這樣的方法所帶來的好處是降低乘法運算次數，原本每個pixel的複雜度是O(n²)，而新的方法則是O(n)。

Code更新了嗎？gnitnaw

一維陣列的實作已經更新 kevinbird61

高斯模糊1D公式

執行環境

Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Byte Order: Little Endian
CPU MHz: 3300.000
L1d cache: 32K
L1i cache: 32K
L2 cache: 256K
L3 cache: 6144K

實驗結果

修改 `naive split` 版本 - 降低每個pixel讀取次數的版本

聽取了gnitnaw的建議後，我計算了如果再一個迴圈對於 25 個 pixel 作讀取一次後，加到周圍所有可能的 pixel 上的實作方法（圖片以1600x1200為基準）：
- 原本所需要的讀取次數 (以下以 Original 做代號)：
  - 1600x1200個 pixel ，每次運算會有 25 次 read 以及 1 次 store；
- 修改後的讀取次數 (以下以 Expand 做代號)：
  - (1600/5)*(1200/5)=320*240（因為每次以 25 個為一個運算單位），每次運算共有 25 次 read 加上 81 次 store；
- 計算總成本（不考慮其他運算子的情況下，單就 store 跟 load做討論）：

	Original	Expand
迴圈總執行次數	1920000	76800
存取成本	48000000(read),1920000(store)	1920000(read),6220800(store)

由上表格中可以看出，Expand比起Original版本要來的節省 read 次數，總共減少了 46080000 次的 read ，並且增加 4300800 次的 store
以MIPS角度上來分析：
- 每個 load(read) ，都會需要先對 array 指標作移動，到正確位置後取值
- 每個 store，也需要對 array 的位置做移動，到正確位置後把值作存入

待查詢： store 是否會影響 pipeline 問題kevinbird61

如此看來，對於整體上看來，若再 read 跟 store 沒有太多差別情況下成本視為相同，這麼一來總共省下41779200 次的（read / store)
並且降低 1843200 次的 branch / jump （ for 迴圈）
總和評估過後，有實作價值存在（大幅降低執行所需要的指令數），剩下的的問題就是展開後的程式碼對於整體影響，以及新的演算法是否可以這樣使用
更改後的計算方式：
- 原本由讀取 25 的 pixel 並計算後，加入到一個 pixel 並成為他的新值
- 更改後的版本，則是每次作25個，並且把每個讀取到的 pixel ，依據位置乘上相對應的 gaussian kernel 存到周圍24個位置，並且包括自己；由原本的 5x5 擴展到 9x9大小的運算：
初步實作結果：
不過看到結果，圖片會有點馬賽克的感覺。
- 目前原因正在確認，很大的可能是我擴展部份，貼到其他 pixel 上的時候，假如不再原本 5x5 當中，也就是再加總時不會有自己出現的情形，那麼我就會用 src[…] += (其他 expand 過來的 pixel) 來做操作。若是再 5x5 內，則是以原本 src[…] = (各點做運算後除以273，並判斷是否大於255)
- 直接加上去的考量是，由於該點再現階段不是主要計算區域，所以我沒有必要現在對該點作操作
初步執行時間比較：
- 把繪圖分成兩張圖片，現在我們來看針對 naive 作繪製的圖片：
- 看到最下面unroll expand split的版本，是目前 naive 中速度最快的

實作: NEON

執行環境: Raspberry pi 3 Model B v1.2

Architecture: armv7l
Byte Order: Little Endian
CPU(s): 4
On-line CPU(s) list: 0-3
Thread(s) per core: 1
Core(s) per socket: 4
Socket(s): 1
Model name: ARMv7 Processor rev 4 (v7l)
CPU max MHz: 1200.0000
CPU min MHz: 600.0000

neon_gaussian_blur












































































...
    // const data

    uint8x16_t vk0 = vdupq_n_u8(0);

    const unsigned char sse_g1[16] = { 1,0,4,0,7,0,4,0,1,0,0,0,0,0,0,0 };
    const unsigned char sse_g2[16] = { 4,0,16,0,26,0,16,0,4,0,0,0,0,0,0,0 };
    const unsigned char sse_g3[16] = { 7,0,26,0,41,0,26,0,7,0,0,0,0,0,0,0 };
    const unsigned char sse_g4[16] = { 4,0,16,0,26,0,16,0,4,0,0,0,0,0,0,0 };
    const unsigned char sse_g5[16] = { 1,0,4,0,7,0,4,0,1,0,0,0,0,0,0,0 };

    // Operation to image
    for(int i=0; i<w-2; i++) {
        for(int j=0; j<h-2; j++) {
            int sum = 0;
            //
            uint8x16_t vg1 = vld1q_u8((uint8_t *)sse_g1);
            uint8x16_t vg2 = vld1q_u8((uint8_t *)sse_g2);
            uint8x16_t vg3 = vld1q_u8((uint8_t *)sse_g3);
            uint8x16_t vg4 = vld1q_u8((uint8_t *)sse_g4);
            uint8x16_t vg5 = vld1q_u8((uint8_t *)sse_g5);

            //
            uint8x16_t vsum = vdupq_n_u8(0),vtemplow = vdupq_n_u8(0),vtemphigh = vdupq_n_u8(0),vempty = vdupq_n_u8(0);

            // First element src[j*w+i]
            // Load /in data
            //
            uint8x16_t L0 = vld1q_u8((uint8_t *)(src+(j+0)*w + i));
            uint8x16_t L1 = vld1q_u8((uint8_t *)(src+(j+1)*w + i));
            uint8x16_t L2 = vld1q_u8((uint8_t *)(src+(j+2)*w + i));
            uint8x16_t L3 = vld1q_u8((uint8_t *)(src+(j+3)*w + i));
            uint8x16_t L4 = vld1q_u8((uint8_t *)(src+(j+4)*w + i));
            // Get the data we need (5 element per-line) , because we only
            // need 5 element from sse instruction set , so only get low part(contain 8 elements)
 //
            vzipq_u8(L0,vk0);//L0: v0
            vzipq_u8(L1,vk0);//L1: v1
            vzipq_u8(L2,vk0);//L2: v2
            vzipq_u8(L3,vk0);//L3: v3
            vzipq_u8(L4,vk0);//L4: v4

            // Multiple with specific Gaussian coef.
            //
            L0 = vmulq_u8(L0, vg1);
            L1 = vmulq_u8(L1, vg2);
            L2 = vmulq_u8(L2, vg3);
            L3 = vmulq_u8(L3, vg4);
            L4 = vmulq_u8(L4, vg5);

            // Summation the 5 line
            //
            vsum = vaddq_u8(vsum, L0);
            vsum = vaddq_u8(vsum, L1);
            vsum = vaddq_u8(vsum, L2);
            vsum = vaddq_u8(vsum, L3);
            vsum = vaddq_u8(vsum, L4);

            // Vsum summation
            // Summation all - (Summation all - (Summation with shift-off 5 number))
            //

            uint8x8x2_t result = vzip_u8(vget_low_u8(vsum),vget_low_u8(vempty));
            vtemplow = vcombine_u8(result.val[0],result.val[1]);

            result = vzip_u8(vget_high_u8(vsum),vget_high_u8(vempty));
            vtemphigh = vcombine_u8(result.val[0],result.val[1]);


            //
            sum += (uint32_t)vgetq_lane_u8(vtemplow,0);
            sum += (uint32_t)vgetq_lane_u8(vtemplow,1);
            sum += (uint32_t)vgetq_lane_u8(vtemplow,2);
            sum += (uint32_t)vgetq_lane_u8(vtemplow,3);
            sum += (uint32_t)vgetq_lane_u8(vtemphigh,0);
...

Source code裡沒這個function？ gnitnaw
這部分是新加進去的，尚未git上去，這幾天測試完後會補上 SwimGlass

neon_gaussian_blur_5_tri 改

neon_flip_horizontal_tri

這邊為了實作出 _mm_shuffle_epi8，在neon中沒有這個function，要使用vtbl2_u8來作轉換

但是vtbl2_u8沒有辦法做出那麼長的對應
uint8x8_t vtbl2_u8(uint8x8x2_t a, uint8x8_t b); // VTBL.8 d0, {d0, d1}, d0
要先將mask拆成兩部分為了符合vtbl的格式



    //const __m128i mask = _mm_set_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15);
    const uint8x8_t down_mask = {7 ,6 ,5 , 4, 3, 2, 1};
    const uint8x8_t up_mask = {15,14,13,12,11,10,9,8};







    uint8x8_t v1_up = vtbl2_u8(v1_882, up_mask);
    uint8x8_t v1_down = vtbl2_u8(v1_882, down_mask);
    uint8x8_t v2_up = vtbl2_u8(v2_882, up_mask);
    uint8x8_t v2_down = vtbl2_u8(v2_882, down_mask);

    v1 = vcombine_u8(v1_up, v1_down);
    v2 = vcombine_u8(v2_up, v2_down);

測試neon

先補上 neon 的編譯參數，將下面加入到 Makefile 中，方便做測試，之後再跟 SSE 統一規格






...
gau_blur_neon_ori: $(GIT_HOOKS) format main.c
        arm-linux-gnueabihf-gcc-5 -c -g -Wall -Wextra -Ofast -mfpu=neon -DARM -o gaussian.o gaussian_arm.c
        arm-linux-gnueabihf-gcc-5 --std gnu99 -DARM -DGAUSSIAN=1024 -Wall -Wcomment -O3 -g -Wextra -Ofast gaussian.o -o main_arm main.c
        bash execute.sh main_arm img/input.bmp output.bmp;
...

測試的時候發現，vmulq_u8這邊做錯了

// Multiple with specific Gaussian coef.
L0 = vmulq_u8(L0, vg1);
L1 = vmulq_u8(L1, vg2);
L2 = vmulq_u8(L2, vg3);
L3 = vmulq_u8(L3, vg4);
L4 = vmulq_u8(L4, vg5);

我們想要做到的是_mm_maddubs_epi16，示意圖如下

原本的 v0 跟 vg1 是 8bit * 16 組存放在register裡面 , 但是 8-bit 最大值為 256 所以如果今天相乘起來的值大於 256 後就會 overflow 所以使用 mm_maddubs_epi16 就能避免這問題
但在 neon 裡面沒有這個涵式的定義 , 所以要自己手工轉

uint16x8_t L0_l = vmovl_u8(vget_low_u8(L0)), L0_h = vmovl_u8(vget_high_u8(L0));
uint16x8_t vg1_l = vmovl_u8(vget_low_u8(vg1)), vg1_h = vmovl_u8(vget_high_u8(vg1));// or pairwise add?
uint16x8_t M0_l = vmulq_u16(L0_l,vg1_l), M0_h = vmulq_u16(L0_h,vg1_h);
//unpack M0_l and insert the first lane of M0_h to the 5th
vuzpq_u16(M0_l,M0_h);

測試結果得到的圖 , 不知道為啥會暗暗的 , 還在尋找原因

以工具 - valgrind 修正 memory leak 問題

由gnitnaw大大建議，由於再執行過程中，有許多 malloc 後並沒有 free ，導致大量的 memory leak 的問題。
使用valgrind這套軟體來作為檢查方式
先在編譯時加上-g，可以讓產生的 error log 顯示再程式碼中的第幾行有問題（讓 valgrind 更加詳細）
之後使用valgrind --leak-check=full ./$(TARGET) ....來讓 valgrind 找出問題所在
- --leak-check=[no/summary/full]：用來列出要顯示多少的訊息使用
- --log-file=[filename]：把產生的訊息導入到指定檔案中
由第一次 valgrind ，產生了50000多次的 memory error 問題，最主要在於當初處理圖片時，邊界部份沒有做好，導致 sse版本的程式會再執行期間，頻繁的到不合法的位置做取值的動作；
剩下的便是部份再 gaussian.c 中部份函式中零星的 malloc 沒有做相對應的free
而再修正這 50000 多筆錯誤後，重新執行繪製執行時間的 makefile target 後，可以看到

Enter the times you want to execute Gaussian blur on the input picture:1
Enter the thread number: 4
Enter the times on perf analysis: 100

 Performance counter stats for './bmpreader img/input.bmp output.bmp 1 4' (100 runs):

           698,735      cache-misses              #   40.627 % of all cache refs      ( +-  2.31% )
         1,719,872      cache-references                                              ( +-  1.39% )

       3.101102475 seconds time elapsed                                          ( +-  1.06% )

跟上方執行做比較，cache-misses 的量有所下降

改善目標：

鏡像 Mirror

想法

分為水平翻轉和垂直翻轉，水平翻轉也就是常見的字串反轉

SIMD 實作水平翻轉












unsigned char *s1 = src + (i*w) + j;
unsigned char *s2 = src + (i*w) + (w-j-1) - 15;
const __m128i mask = _mm_set_epi8(0, 1, 2, 3,
				    4, 5, 6, 7,
				    8, 9, 10, 11,
				    12, 13, 14, 15);
__m128i v1 = _mm_loadu_si128((__m128i *)s1);
__m128i v2 = _mm_loadu_si128((__m128i *)s2);
v1 = _mm_shuffle_epi8(v1, mask);
v2 = _mm_shuffle_epi8(v2, mask);
_mm_storeu_si128((__m128i *)s2, v1);
_mm_storeu_si128((__m128i *)s1, v2);

這邊比較特殊的是 _mm_shuffle_epi8() 的用法
下面是函式的原型

__m128i _mm_shuffle_epi8( 
   __m128i a,
   __m128i mask
);

而內部的操作為

r0 = (mask0 & 0x80) ? 0 : SELECT(a, mask0 & 0x0f)
r1 = (mask1 & 0x80) ? 0 : SELECT(a, mask1 & 0x0f)
...
r15 = (mask15 & 0x80) ? 0 : SELECT(a, mask15 & 0x0f)

SELECT()

SELECT(a, n) extracts the nth 8-bit parameter from a. 
The 0th 8-bit parameter is the least significant 8-bits.

因此當我們將 mask 的值放入 0, 1 , …, 15，在 little endian 的系統上，最後可以取得以 8-bit 為單位進行反轉的結果，簡單的來說就是可以用來做 16 bytes 的字串反轉

實驗

翻轉	naive(ori)	naive(tri)	openmp(tri)	sse(tri)
水平翻轉	0.006459	0.016355	0.016314	0.004591
垂直翻轉	0.006776	0.022558	0.022308	0.001510

垂直翻轉

水平翻轉

HSV 飽和度調整、明亮度調整

想法

先將 RGB 轉到 HSV 表示，將 S 值和 V 值上調或下調後，再將 HSV 轉回 GRB

實驗

naive : 0.107381 sec
OpenMP : 0.049145 sec

透過 SIMD 加速高斯模糊運算

tags: SIMD sysprog

預期目標

程式碼產出: Image-Processing

執行環境

高斯模糊

高斯模糊公式

實作 - 兩種資料結構差異

實作: SSE/AVX

實作 - unroll version

實作: POSIX Thread

修改SSE部份的實作機制

實作: Pthread + SSE

修正繪圖曲線抖動問題

曲線推論與測試 by @shinshipower

實作1D Gaussian Filter

執行環境

實驗結果

修改 naive split 版本 - 降低每個pixel讀取次數的版本

實作: NEON

執行環境: Raspberry pi 3 Model B v1.2

neon_gaussian_blur

neon_gaussian_blur_5_tri 改

neon_flip_horizontal_tri

測試neon

以工具 - valgrind 修正 memory leak 問題

改善目標：

鏡像 Mirror

想法

SIMD 實作水平翻轉

實驗

HSV 飽和度調整、明亮度調整

想法

實驗

Reference

tags: `SIMD` `sysprog`

修改 `naive split` 版本 - 降低每個pixel讀取次數的版本