# 使用SSE/AVX做圖檔處理
contributed by <`CheHsuan`>, <`kevinbird61`>, <`SwimGlass`>, <`carolc0708`>
###### tags: `sysprog`
>>中英文字間請記得以空白區隔
>>[color=red][name=課程助教]
>>>OK!
>>>[color=blue][name=kevinbird61]
## 選擇目標
- 我選擇了沒經過壓縮的圖片檔格式 -bmp 作為實驗目標
- 以影像處理為實驗方式，進而檢測效能
## Source code
[github](https://github.com/kevinbird61/Image-Processing)
## 實作項目
1. 高斯模糊
- 公式
2. 水平、垂直翻轉
3. HSV飽和度、明亮度調整
`以下多個 function 實作時間評估的圖片，皆為對圖片執行一次高斯模糊，並且使用 4 個執行序來實作結果`
### 高斯模糊
- 先使用原本的資料結構和分開的結構實作：
```C=34
// 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));
```
- 執行：
```shell
$~: make gau_blur_tri
gcc main.c -DSPLIT=7 -DGAUSSIAN=1 -o bmpreader
$~: make run
# img/wf.bmp => has the 4 element(alpha value)
./bmpreader img/input.bmp output.bmp
Picture size of picture is width: 1600 , height 1200
Read file successfully
Gaussian blur[5x5][split structure], execution time : 0.395602 sec
Save file successfully
eog output.bmp
...
$~: make gau_blur_ori
gcc main.c -DGAUSSIAN=2 -o bmpreader
$~: make run
# img/wf.bmp => has the 4 element(alpha value)
./bmpreader img/input.bmp output.bmp
Picture size of picture is width: 1600 , height 1200
Read file successfully
Gaussian blur[5x5][original structure], execution time : 0.322005 sec
Save file successfully
eog output.bmp
```
- 可以看到，由於以 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: 3*1 = 3 bytes(24bits)、 original: 3*3 = 9 bytes (72 bits) ; 5 => split: 5*1 = 5 bytes(40 bits) 、 original: 5*3 = 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)*8*3 = 128 => n = 16/3 => 會有額外的浪費！ => 先實作分散版本
- 我們先 Load 進5行(每行16 bytes)的大小來做，找出其中5x5的元素乘上高斯模糊係數，並且加總後放回原先的位置
- 考慮到使用 SSE 的效益，先實作看看unroll版本
#### 實作 - unroll version
```C=41
...
sum = src[(j-2)*w+(i-2)]*gaussian55[index++] + src[(j-2)*w+(i-1)]*gaussian55[index++]
+ src[(j-2)*w+(i)]*gaussian55[index++] + src[(j-2)*w+(i+1)]*gaussian55[index++]
+ src[(j-2)*w+(i+2)]*gaussian55[index++] + src[(j-1)*w+(i-2)]*gaussian55[index++]
+ src[(j-1)*w+(i-1)]*gaussian55[index++] + src[(j-1)*w+(i)]*gaussian55[index++]
+ src[(j-1)*w+(i+1)]*gaussian55[index++] + src[(j-1)*w+(i+2)]*gaussian55[index++]
+ src[(j)*w+(i-2)]*gaussian55[index++] + src[(j)*w+(i-1)]*gaussian55[index++]
+ src[(j)*w+(i)]*gaussian55[index++] + src[(j)*w+(i+1)]*gaussian55[index++]
+ src[(j)*w+(i+2)]*gaussian55[index++] + src[(j+1)*w+(i-2)]*gaussian55[index++]
+ src[(j+1)*w+(i-1)]*gaussian55[index++] + src[(j+1)*w+(i)]*gaussian55[index++]
+ src[(j+1)*w+(i+1)]*gaussian55[index++] + src[(j+1)*w+(i+2)]*gaussian55[index++]
+ src[(j+2)*w+(i-2)]*gaussian55[index++] + src[(j+2)*w+(i-1)]*gaussian55[index++]
+ src[(j+2)*w+(i)]*gaussian55[index++] + src[(j+2)*w+(i+1)]*gaussian55[index++]
+ src[(j+2)*w+(i+2)]*gaussian55[index++];
```
- 分析執行結果：經過更改後的 makefile ，可以讓使用者決定要對該圖檔執行幾次的 Gaussian blur ；而比較原本 5x5 分開的 structure , original 來做，連續執行 3 次前後效能比較：
- Split: 1.274697 sec
- Original: 0.961214 sec
- Unroll: 0.873029 sec
- 可以看到，Split vs unroll 一共減少`0.401668`秒；而Original vs unroll 則是減少了`0.088185`秒
- 而實作 original structure unroll 版本，得到的執行時間是：` 0.985009 sec`；竟然比原先的還要來的慢...
- PS：後來多次執行後，執行時間跟原本差不多（忽高忽低）
- 圖為各函式對圖形執行一次 Gaussian blur 的結果發現：
- sse的操作方式並沒有很節省空間的使用(5x5 , 3x3的操作，每次取完能用的有限)，導致對圖檔做三次 Gaussian blur 的結果為`1.842012`秒
- 修改再同時一個 instruction 中，為了能夠一次操作更多，改成一次操作並儲存兩個位置。執行後時間反而增加:`3.015885`秒
- 由於SSE本身需要load進16個位置的資訊，導致每次操作時都需要花上不少指令來拿出特定的5個（ per line ）來做，入不敷出。每當為了運算出一個位置，就會付出40幾行的程式碼，以及將近40個SSE function call
```C=196
...
__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
...
```
- 所以我認為對於高斯模糊，用loop unrolling是最好的解法！
- 每個元素必須乘上一定係數
- 並且必須以周圍的數值加總後除以一個常數成為自己的模糊值
- 操作過於瑣碎，並且相依，無法同時來進行
- 為了得到個別值，付出的代價太高
#### 實作 - Pthread
- 考慮到可以分頭進行
- 實作 pthread 來分擔 Gaussian blur ，依據不同的 thread_id * workload（由整張圖片高度除上總共thread） 作為不同 thread 起始操作位置來做；而一樣對同一張圖作三次高斯模糊（unroll split structure），執行時間為：`0.402466`秒（ 4 條 thread）！！
- 重新繪製執行時間圖(執行一次 blur + 4 threads)：
#### 修改SSE部份的實作機制
- 由於一次load進16bytes，考慮到原本資料結構為 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 來幫每個移動過後的元素作相乘得動作。最後再依序取出並加總變可以得到最後該位置的值！
- 後續：加上prefetch後，發現沒有比較快
#### 實作 - Pthread + SSE
- 使用 pthread + SSE 來再次實作=>
- 可以看到，執行時間幾乎壓再 0.1 秒左右（紫色三角）！！整整比起原先（黃色方形）來的**快了 0.2 秒左右**！！
> 實驗：pthread 交錯工作，並檢視效能
### Mirror
#### 想法
分為水平翻轉和垂直翻轉，水平翻轉也就是常見的字串反轉
<table>
<tr>
<td>翻轉</td>
<td>naive(ori)</td>
<td>naive(tri)</td>
<td>openmp(tri)</td>
<td>sse(tri)</td>
</tr>
<tr>
<td>水平翻轉</td>
<td>0.006459</td>
<td>0.016355</td>
<td>0.016314</td>
<td>0.004591</td>
</tr>
<tr>
<td>垂直翻轉</td>
<td>0.006776</td>
<td>0.022558</td>
<td>0.022308</td>
<td>0.001510</td>
</tr>
</table>
### HSV 飽和度調整、明亮度調整
#### 想法
先將RGB轉到HSV表示，將S值和V值上調或下調後，再將HSV轉回GRB
#### 實驗
- naive : 0.107381 sec
- openmp : 0.049145 sec
# Reference
- [Intel Intrinsics Guide](https://software.intel.com/sites/landingpage/IntrinsicsGuide/#expand=3124,3125,5604,61,3392,5208,1784,1784,5586,5577)
- [Gaussian blur](https://read01.com/2go7z4.html)
-
Sign in with Wallet
Connect another wallet