Vulkan移植GpuImage(三)從A到C的濾鏡

前面移植了幾個比較複雜的效果後,算是確認了複雜濾鏡不會對框架造成比較大的改動,開始從頭移植,現已把A到C的所有濾鏡用vulkan的ComputeShader實現了,講一些其中實現的過程.

AverageLuminanceThreshold 畫素亮度平均閾值比較

從名字來看,就是算整圖的平均高度,然後比較這個亮度值.

GPUImage的實現,先平均縮少3*3倍,然後讀到CPU中計算平均亮度,然後再給下一層計算.

這步回讀會浪費大量時間,我之前在CUDA測試過,1080P的回讀大約在2ms左右,就算少了9倍資料量,也可能需要0.2ms,再加上回讀CPU需要同步vulkan的cmd執行執行緒,早早提交執行command,中斷流程所導致的同步時間根據執行層的複雜度可能會比上面更長.

因此在這裡,不考慮GPUImage的這種實現方式,全GPU流程處理,使用Reduce方式算影像的聚合資料(min/max/sum)等,然後儲存結果到1x1的紋理中,現在實現效果在2070下1080P下需要0.08ms,比一般的普通計算層更短.

主要實現類似opencv cuda裡的reduce相關實現,但是他在最後使用原子操作所有區域性共享視訊記憶體的值,而在glsl中,原子操作限定太多,因此我分二步來操作.

第一步中,每個點取4x4個值,這步在我想象中,應該就需要0.2ms左右,根據opencv cuda相關的reduce的計算方法,使用周邊PATCH_SIZE_XxPATCH_SIZE_Y個執行緒組的執行緒,互相混合取PATCH_SIZE_XxPATCH_SIZE_Y個值,其中每一個執行緒會賦值給周邊PATCH_SIZE執行緒邊上地址,編號x0/y0就給每個PATCH_SIZE所有執行緒第x號塊賦值,這裡程式碼看著是有些奇怪,如果讓我來實現,我肯定就取執行緒周邊PATCH_SIZE_XxPATCH_SIZE_Y來進行操作,對於GPU來說,應該沒區別才是,在單個執行緒中PATCH_SIZE_XxPATCH_SIZE_Y個資料都是序列操作啊.

在第一步上面操作後,把sizeX(1920)/sizeX(1080)變成只有(480*270)個執行緒,其中每個執行緒組有16x16個,也就一共有510塊執行緒組,每塊執行緒組使用並行2次分操作,其中16x16=2^8,8次後就能得到所有聚合資料.

在opencv中,就是針對其中的510個執行緒進行原子操作,我看了下,glsl裡的原子操作只能針對型別int/uint,侷限太大,因此我應用了我以前CUDA版Grabcut的實現中的kmeans優化方式,把上面計算後的510個執行緒組中的資料放入510*1的臨時texture.

注意在GPGPU運算層中,不要針對BUFFER又讀又寫,以前寫CUDA相關演算法時,嘗試過幾次,不管你怎麼呼叫同步API,結果全不對,在第一步中最後把結果寫入臨時BUFFER,就需要在第二步,讀取這個臨時BUFFER.

然後開始第二步,讀取這個臨時BUFFER,因為要聚合所有資料,所以我們執行緒組就只分一個256個執行緒的組,在這個組裡,使用for步長執行緒組大小來訪問這個臨時BUFFER的所有資料,然後聚合二分,最後把結果給到一個1x1的紋理中,如下擷取第一部分的程式碼出來,有興趣可以自己根據連結檢視全部程式碼.

#version 450

layout (local_size_x = 16, local_size_y = 16) in;// gl_WorkGroupSize

#if CHANNEL_RGBA
layout (binding = 0, rgba8) uniform readonly image2D inTex;
#elif CHANNEL_R8
layout (binding = 0, r8) uniform readonly image2D inTex;
#elif CHANNEL_RGBA32F
layout (binding = 0, rgba32f) uniform readonly image2D inTex;
#elif CHANNEL_R32F
layout (binding = 0, r32f) uniform image2D inTex;
#endif

#if CHANNEL_RGBA || CHANNEL_RGBA32F
layout (binding = 1, rgba32f) uniform image2D outTex;
#elif CHANNEL_R8 || CHANNEL_R32F
layout (binding = 1, r32f) uniform image2D outTex;
#endif

shared vec4 data_shared[256];

// 一個執行緒處理每行PATCH_SIZE_X個元素
const int PATCH_SIZE_X = 4;
// 一個執行緒處理每列PATCH_SIZE_Y個元素
const int PATCH_SIZE_Y = 4;
// 每個執行緒組處理元素個數為:block size(16*4)*(16*4)

// min/max/sum 等
#if REDUCE_MIN
    #define OPERATE min
    #define ATOMIC_OPERATE atomicMin
    #define INIT_VEC4 vec4(1.0f)
#endif

#if REDUCE_MAX
    #define OPERATE max
    #define ATOMIC_OPERATE atomicMax
    #define INIT_VEC4 vec4(0.0f)
#endif

#if REDUCE_SUM
    #define OPERATE add
    #define ATOMIC_OPERATE atomicAdd
    #define INIT_VEC4 vec4(0.0f)
#endif

vec4 add(vec4 a,vec4 b){
    return a+b;
}

// 前面一個執行緒取多點的邏輯參照opencv cuda模組裡的reduce思路
void main(){
    ivec2 size = imageSize(inTex);  
    ivec2 uv = ivec2(gl_GlobalInvocationID.xy);  
    if(uv.x >= size.x || uv.y >= size.y){
        return;
    }  
    // 組內執行緒一維索引    
    int tid = int(gl_LocalInvocationIndex);  
    data_shared[tid] = INIT_VEC4;
#if REDUCE_AVERAGE
    float avg = 1000.0/(gl_WorkGroupSize.x * gl_WorkGroupSize.y * PATCH_SIZE_Y * PATCH_SIZE_X);
#endif
    memoryBarrierShared();
    barrier();
    // 執行緒塊對應記憶體塊索引
    uint x0 = gl_WorkGroupID.x * (gl_WorkGroupSize.x*PATCH_SIZE_X) + gl_LocalInvocationID.x;
    uint y0 = gl_WorkGroupID.y * (gl_WorkGroupSize.y*PATCH_SIZE_Y) + gl_LocalInvocationID.y;
    // 周邊PATCH_SIZE_X*PATCH_SIZE_Y個執行緒組的執行緒,互相混合取PATCH_SIZE_X*PATCH_SIZE_Y個值.
    // 每一個執行緒會賦值給周邊PATCH_SIZE執行緒邊上地址,編號x0/y0就給每個PATCH_SIZE所有執行緒第x號塊賦值.
    // 相比直接取本身執行緒組周邊PATCH_SIZE_X*PATCH_SIZE_Y地址進行比較來說,有什麼區別嗎?
    // 由sizex/sizey的範圍縮小到數執行緒大小(sizex/PATCH_SIZE_X,sizey/PATCH_SIZE_Y)範圍
    for (uint i = 0, y = y0; i < PATCH_SIZE_Y && y < size.y; ++i, y += gl_WorkGroupSize.y){
        for (uint j = 0, x = x0; j < PATCH_SIZE_X && x < size.x; ++j, x += gl_WorkGroupSize.x){
            vec4 rgba = imageLoad(inTex,ivec2(x,y));            
            data_shared[tid] = OPERATE(rgba,data_shared[tid]);
        }
    }
    memoryBarrierShared();
    barrier();
    // 然後執行緒組內二分比較,把值儲存在data_shared[0]中
    for (uint stride = gl_WorkGroupSize.x*gl_WorkGroupSize.y / 2; stride > 0; stride >>= 1) {
       if (tid < stride){
            data_shared[tid] = OPERATE(data_shared[tid], data_shared[tid+stride]);                     
        }
        memoryBarrierShared();
        barrier();
    }
    memoryBarrierShared();
    barrier();
    // 原子操作所有data_shared[0]
    // if(tid == 0){
    //     ATOMIC_OPERATE()
    // }
    // 原子操作限定太多,放棄
    if(tid == 0){ 
        int wid = int(gl_WorkGroupID.x + gl_WorkGroupID.y * gl_NumWorkGroups.x);         
        imageStore(outTex, ivec2(wid,0), data_shared[0]);
    }
}

pre reduce glsl程式碼

reduce glsl程式碼

AverageLuminanceThreshold C++ 實現

我在做之前根據3x3卷積需要0.2ms左右粗略估算下需要的時間應該在0.3ms左右,但是實際只有(0.07+0.01)ms,後面想了下,這其實是有個很大區別,模糊那種核是圖中每個點需要取周邊多少個點,一共取點是畫素x核長x核長,而Reduce運算最開始一個點取多個畫素,但是總值還是隻有影像畫素大小.

對於影像流(視訊,拉流)來說,可以考慮在當前graph執行完vkQueueSubmit,然後VkReduceLayer層輸出的1x1的結果,類似輸出層,然後在需要用這結果的層,在下一楨之前把這個結果寫入UNIFORM_BUFFER中,可能相比取紋理值更快,這樣不會打斷執行流程,也不需要同步,唯一問題是當前楨的引數是前面的楨的執行結果.

然後下面的AddBlend/AlphaBlend算是比較常規的影像處理,GPGPU最容易處理的型別,也就不說了.

AmatorkaFilter(顏色查詢表對映)

理解下原理就行,如GPUImage裡,需要的是一張512*512畫素,其中有8x8x(格子(64x64)),把8x8=64格子垂直平放堆疊,可以理解成一個每面64畫素的正方體,其中每個格子的水平方是紅色從0-1.0,垂直方向是綠色從0-1.0,而對於正方體的垂直方向是藍色從0-1.0.

然後就是顏色對映,把對應的顏色r/g/b,當做如上三個軸的座標,原圖的藍色確定是那二個格子(浮點數需要二個格子平均),紅色找到對應格子的水平方向,綠色是垂直方向,查詢到的顏色就是對映的顏色.

在這,如果自己linear插值程式碼就有點多了,要取周邊八個點,效能可能還受影響,利用linear sampler2D簡化下,只需要自己linear插值一次,在這我也找到以前傳入取樣器不正確的BUG,生成Texture的時間沒有加入VK_IMAGE_USAGE_SAMPLED_BIT標記,不過大部分影像處理並不需要,所以新增二個方法用來表達是否需要sampled,以及是linear/nearest插值方式.

GPUImage中,那個texPos1什麼0.125 - 1.0/512.0啥的,我一開始有點暈,後面自己推理下,其實應該是這樣.

#version 450

layout (local_size_x = 16, local_size_y = 16) in;// gl_WorkGroupSize
layout (binding = 0, rgba8) uniform readonly image2D inTex;
// 8*8*(格子(64*64)),把8*8向上平放堆,可以理解成一個每面64畫素的正方體.
// 其中每個格子的水平方是紅色從0-1.0,垂直方向是綠色從0-1.0,而對於正方體的垂直方向是藍色從0-1.0.
layout (binding = 1) uniform sampler2D inSampler;
layout (binding = 2, rgba8) uniform image2D outTex;

void main(){
    ivec2 uv = ivec2(gl_GlobalInvocationID.xy);
    ivec2 size = imageSize(inTex);
    if(uv.x >= size.x || uv.y >= size.y){
        return;
    }     
    vec4 color = imageLoad(inTex,uv);       
    // 如果b是0.175,對應b所在位置11.2,那就查詢二張圖(11,12)
    float b = color.b * 63.0;
    // 查詢11這張圖在8*8格子所在位置(1,3)
    vec2 q1 ;
    q1.y = floor(floor(b)/8.0f);
    q1.x = floor(b) - (q1.y * 8.0f);
    // 查詢12這張圖在8*8格子所在位置(1,4)
    vec2 q2 ;
    q2.y = floor(ceil(b)/8.0f);
    q2.x = ceil(b) - (q2.y * 8.0f); 
    // 格子UV+在每個格子中的UV,q1/q2[0,8)整數,rg[0.0,1.0]浮點位置轉到[0,512)
    // 整點[0,512)轉(0.0,1.0f),需要(整點位置+0.5)/512.0
    vec2 pos1 = (q1*64.0f + color.rg*63.0f + 0.5f)/512.f;
    vec2 pos2 = (q2*64.0f + color.rg*63.0f + 0.5f)/512.f;
    // 取pos1/pos2上值
    vec4 c1 = textureLod(inSampler, pos1,0.0f);  
    vec4 c2 = textureLod(inSampler, pos2,0); 
    // linear混合11,12這二張圖的值
    vec4 result = mix(c1,c2,fract(b));
    imageStore(outTex,uv,result);
}

然後再比較看了一下,其實是一樣的.

這個Lookup的實現,對應的他的Amatorka/MissEtikate/SoftElegance,本框架為了儘量少的引用第三方庫,沒有讀取圖片的類庫,需要自己把Lookup表傳入inputLayer,然後連線這層,後續會新增demo接入相應外部庫給出相關對應如上所有層的實現.

原則上主體框架不引入第三方庫,框架上的模組不引入非必要的第三方庫,而demo不限制.

BilateralFilter 雙邊濾波

網上很多講解,這簡單說下,為什麼叫雙邊,一個是高斯濾波,這個濾波只考慮了周邊畫素距離,加一個濾波對應周邊畫素顏色差值,這樣可以實現在邊緣(顏色差值大)減少模糊效果.

實現也比較簡單,拿出來說是因為暫時還沒找到優化方法,怎麼說了,這個卷積核因為不只和距離有關了,還有周邊畫素顏色差有關,所以只有計算時才能得到每個畫素與眾不同的卷積核,優化方法就需要特別一些,現在的情況的核大的話,如10核長的需要3ms,對應一般的處理層0.2ms來說,簡直誇張了,現暫時這樣實現,後續查詢相應資料後經行優化.

前面有講過優化的BoxBlur(同高斯優化)相關實現,餘下B字母全是普通的常規的影像處理,就不拿來說了.

CannyEdgeDetection Canny邊緣檢測

Canny Edge Detection Canny邊緣檢測

邏輯有點同HarrisCornerDetection,由多層構成,按上面連結來看,其中第四層Hysteresis Thresholding正確實現方法邏輯應該類似:opencv_cudaimgproc canny.cpp/canny.cu裡edgesHysteresis,不過邏輯現有些複雜,後面有時間修改成這種邏輯,現暫時使用GPUImage裡的這種邏輯.

下面的CGAColorspace是常規處理,ChromaKey前面移植過UE4 Matting的扣像處理,這個裡面還考慮到整合當前環境光的處理,暫時還沒看到ChromaKeyBlend/ChromaKey有更高明的實現邏輯,就不移植了.

Closing(膨脹(Dilation)和腐蝕(Erosion))

膨脹用於擴張放大影像中的明亮白色區域,侵蝕恰恰相反,而Closing操作就是先膨脹後侵蝕.

參考 CUDA-dilation-and-erosion-filters

很有意思的是,這個專案有個優化比較統計,抄錄如下.

I have performed some tests on a Nvidia GTX 760.

With an image of 1280x1024 and a radio ranging from 2 to 15:

Radio / Implementation	Speed-up	CPU	Naïve	Separable	Shared mem.	Radio templatized	Filter op. templatized
2	34x	0.07057s	0.00263s	0.00213s	0.00209s	0.00207s	0.00207s
3	42x	0.08821s	0.00357s	0.00229s	0.00213s	0.00211s	0.00210s
4	48x	0.10283s	0.00465s	0.00240s	0.00213s	0.00221s	0.00213s
5	56x	0.12405s	0.00604s	0.00258s	0.00219s	0.00219s	0.00221s
10	85x	0.20183s	0.01663s	0.00335s	0.00234s	0.00237s	0.00237s
15	95x	0.26114s	0.03373s	0.00433s	0.00287s	0.00273s	0.00274s

其中cpu的邏輯同cuda Separable邏輯一樣,都做了行列分離的優化.

這個GPU程式碼也有點意思,我在初看Shared mem部分程式碼時還以為邏輯有問題了,都沒有填充完需要的shared部分,看了grid/block的分配才發現,原來這樣也行,記錄下這個思路,grid還是正常分配,但是block根據核長變大,這樣執行緒組多少沒變,但是執行緒組的大小變大,總的執行緒變多.根據threadIdx確定執行緒組內索引,然後根據blockId與傳入實際執行緒組大小確定正確的輸入BUFFER索引,注意,直接使用全域性的id是對應不上相應的輸入BUFFER索引.

我在這還是使用前面的GaussianBlur row裡的優化,但是相關過程有些複雜,我簡化了下,沒有PATCH_PER_BLOCK的概念,在固定16x16的執行緒組裡,一個執行緒取三個點,如Morph row塊中,分別是左16x16,自身16x16,右16x16塊,也沒什麼判斷,核長在32以內都能滿足,程式碼抄錄第一部分如下.

#version 450

// #define IS_SHARED 1

layout (local_size_x = 16, local_size_y = 16) in;// gl_WorkGroupSize
layout (binding = 0, r8) uniform readonly image2D inTex;
layout (binding = 1, r8) uniform image2D outTex;

layout (binding = 2) uniform UBO {
	int ksize;	    
} ubo;

#if EROSION 
    #define OPERATE min
    #define INIT_VUL 1.0f
#endif

#if DILATION 
    #define OPERATE max
    #define INIT_VUL 0.0f
#endif

#if IS_SHARED
// 限定最大核為32
shared float row_shared[16][16*3];
void main(){
    ivec2 uv = ivec2(gl_GlobalInvocationID.xy);
    ivec2 size = imageSize(inTex);
    if(uv.x >= size.x || uv.y >= size.y){
        return;
    } 
    ivec2 locId = ivec2(gl_LocalInvocationID.xy);
    for(int i = 0; i < 3; i++){
        uint gIdx = max(0,min(uv.x+(i-1)*16,size.x-1));
        row_shared[locId.y][locId.x + i*16] = imageLoad(inTex,ivec2(gIdx,uv.y)).r;      
    }
    memoryBarrierShared();
	barrier();
    float result = INIT_VUL;
    for(int i =0; i < ubo.ksize; i++){
        int ix = locId.x - ubo.ksize/2 + i;
        float fr = row_shared[locId.y][16 + ix];
        result = OPERATE(fr,result);
    }
    imageStore(outTex, uv, vec4(result)); 
}

#else

void main(){
    ivec2 uv = ivec2(gl_GlobalInvocationID.xy);
    ivec2 size = imageSize(inTex);
    if(uv.x >= size.x || uv.y >= size.y){
        return;
    } 
    float result = INIT_VUL;
    for(int i = 0; i< ubo.ksize; ++i){
        int x = uv.x-ubo.ksize/2+i;
        x = max(0,min(x,size.x-1));
        float r = imageLoad(inTex,ivec2(x,uv.y)).r;
        result = OPERATE(result,r);
    }
    imageStore(outTex, uv, vec4(result));     
}

#endif

有二個版本,分別是使用共享區域性視訊記憶體和不使用共享區域性視訊記憶體的,也比較下時間.

核長	時間	區域性共享視訊記憶體
3	0.15ms*2	是
10	0.17msx2	是
20	0.21msx2	是
3	0.14ms*2	否
10	0.27msx2	否
20	0.50msx2	否

對比使用/不使用共享區域性視訊記憶體,在核不大的情況下,沒什麼區別,大核優勢開始增加.

其在大核的情況下,表現比原來的GaussianBlur裡效果更好,效能比對可以看PC平臺Vulkan運算層時間記錄裡的結果,其在20核長,原GaussianBlur要比Morph多2到3倍時間,不過GaussianBlur row這種寫法擴充套件性會更好,對執行緒組大小也沒要求.

然後是一大堆常規C開頭字母的濾鏡處理,其大部分在VkColorBlendLayer,因邏輯簡單類似,就偷懶放這一個檔案裡處理,略過說明.

LBP(Local Binary Patterns)特徵檢測

剛看GPUImage程式碼實現有點奇怪,看了OpenCV——LBP(Local Binary Patterns)特徵檢測才明白,就是用RGB每個通道8位儲存周邊比較的值而已.

ColorPackingFilter,我看了下是HarrisCornerDetection的一個子元件,然後被棄用了,所以不移植了.

FAST特徵檢測器

初看GPUImage裡的實現,以為是用了二層,後面發現根本沒有用BoxBlue層,其實現也和網上說明有區別,只查詢了周邊八個頂點,暫時用GPUImage裡的實現,後續看有時間移植opencv裡的實現不.

ContrastFilter常規影像處理,略過說明.

CropFilter

我看說明是擷取一部分,實現一個類似不復雜,不過定義的引數可能有區別,我是定義中心的UV,以及擷取長寬四個引數,其中特別加了個處理,如果長寬為0,只取一個點儲存到1x1的紋理中,我是想到後續實現如我在影像上取一點,然後用於下一層影像去處理,比如扣像中,我取影像中一點用於ChromaKeyColor用於扣像層.

CrosshairGenerator這個我看了下GPUImage實現,發現盡然有用到頂點著色器的邏輯操作,我在網上也沒找到這個類的實現效果是啥,暫時就先不移植了,後續看到實現效果再補上.

CrosshatchFilter常規影像處理,略過說明.

結尾說明

差不多到這裡,從A到C的所有層都用vulkan的ComputeShader實現完畢,加上前面移植的幾個層,感覺移植所有濾鏡的進度應該在20%左右了,不過我準備先移植所有層,然後再測試所有層,所以現在一些層可能邏輯並不正確,後續會給每層加上測試.

上面一直提到一些層的花費時間,全在PC平臺Vulkan運算層時間記錄裡,這個用於計錄一些層的效率,也用來後續優化方向,因為會一直更新,再加上對應記錄圖片全在github上,所以先只發github上的連結,大家有興趣可以看看.