CUDA版Grabcut的實現

　　在上次用 CUDA實現導向濾波 後，想着導向濾波能以很小的mask還原高分辨率下的邊緣，能不能搞點事情出來，當時正好在研究Darknet框架，而後又看到grabcut算法，用opencv試了下，感受效果有點意思，後面想了下，這幾個能夠連在一塊兒，先讀取高分辨率的圖像，而後用下降分辨率先經過yolov3算出人物框(很是穩定，不跳,幾乎不會出現有人而找不到的狀況)，再用grabcut算出低mask,而後用這個mask結合原圖用導向濾波獲得高分辨率下清晰邊緣的分割圖，最後再把CUDA算出的結果直接丟給UE4/Unity3D的對應DX11中。

　　先放出cuda版grabcut的效果圖。

　　

　　 固然opencv自己提供的grabcut是用cpu算的，416*416差很少有個幾楨左右，確定不知足要求，因此在這個前提下，用Cuda來實現整個grabcut整個流程。

　　因要集成到UE4/Untiy3D對應軟件中，在window本臺下，咱們首先須要以下。

　　1 拿到darknet源碼，C風格代碼，說實話，沒想過能看到這種短小精湛的深度學習模型框架，帶cuda與cudnn,直接在VS下建一個動態連接庫，把相應代碼放進來，很少，改動幾個Linux下的函數，包含編譯符GPU與CUDNN,而後修改darknet.h這個函數，dllexport這個文件下函數，後期還能夠改下，這個文件下的編譯符移到別的位置，pthread頭文件的引用這些，以及能夠直接傳入GPU數據這幾點。opencv4雖然已經有yolo3的實現，可是咱們要GPU的，後續的操做幾乎全在GPU上。

　　2 咱們主要是參照opencv裏的grabcut實現，爲了更好的參數一些數據，咱們最好編譯本身的opencv版本,我是用的opencv-4.0.0-alpha，比較老的一個版本，須要帶opencv_contrib，包含opencv_cuda相關的模塊，主要是後期咱們實現cuda 版grabcut若是很差確認咱們是否正常實現就能夠調試進去看數據值，看源碼，以及用GpuMat/PtrStepSz，這個簡單封裝真的很適合處理圖像數據。

　　如上二點完成後，咱們能夠簡單分析下grabcut算法主要由那種算法構成，如何完成一個流程，一次大迭代主要以下，主要是先定框，框內爲可能前景區，框外爲背景區，二個區域分別爲Kmeans分類，而後二個分類的數據分別GMM創建模型，最後用GMM算出最大流算法須要的每點的source/sink,創建graph,用GPU友好的maxflow算法push-relabel獲得最小切，也就是mask，而後把mask給GMM從新創建模型而後重複這個過程。

　　如上所說，咱們主要實現三種算法，分別是kmeans,gmm,push-relabel。下面如無別的說明，本機給的是在1070下416*416分辨率下的時間。

　　首先講下kmeans的優化，kmeans的實現比較簡單，就是一個不斷從新分配中心點至到穩定的過程，其中有個過程是把數據根據先後景分別歸類(咱們假設聚合爲5塊)，最直觀最簡單的方法就是用原子操做。

__global__ void updateClusterAtomic(PtrStepSz<uchar4> source, PtrStepSz<uchar> clusterIndex, kmeansI& meansbg, kmeansI& meansfg)
{
    const int idx = blockDim.x * blockIdx.x + threadIdx.x;
    const int idy = blockDim.y * blockIdx.y + threadIdx.y;
    if (idx < source.cols && idy < source.rows)
    {
        float4 color = rgbauchar42float4(source(idy, idx));
        //背景塊,使用kcentersbg,不然使用kcentersfg
        kmeansI& kmeans = (!inRect(idx, idy)) ? meansbg : meansfg;
        int index = clusterIndex(idy, idx);
        atomicAdd(&kmeans.kcenters[index].x, color.x);
        atomicAdd(&kmeans.kcenters[index].y, color.y);
        atomicAdd(&kmeans.kcenters[index].z, color.z);
        atomicAdd(&kmeans.kcenters[index].w, color.w);
    }
}

updateClusterAtomic

　　用Nsight看了下時間，花費超過1ms了，我是想到全部線程最後歸類時都只操做10塊地方，原子操做在這有點把多線程變成有限的幾個能同時作事了。　　

　　

　　所以須要換種想法，很簡單，不用原子操做，用一部分線程如32*8，遍歷最個紋理數據，而後在一個線程裏32*8個數據再整合，使用這步後，時間下降0.5ms左右，利用的還不是很徹底，畢竟一次kmeans穩定須要遍歷十幾回左右，這個過程就要了5ms後面就不用搞了，在這基礎上，結合共享顯存，每一個點一個線程先用共享顯存處理後放入對應grid塊大小的顯存中，而後數據除32*8的grid塊執行上面的操做，這個操做包含後面統計等加起0.12ms,這樣十幾回迭代也就不到2ms，不過我沒想到的是，從cudaMemcpyDeviceToHost一個int32的四字節Nsight給的時間很低0.001ms，可是從上面的簡隔來看，應該有0.1ms了，應該還有一些別的花費，因此不少操做，後續的統計計算我就是用GPU裏一個核心來算，也不先download給CPU算而後再upload上來，可是這個操做還沒想到辦法避免，畢竟須要在CPU上斷定是否已經穩定後關閉循環。

//把source全部收集到一塊gridDim.x*gridDim.y塊數據上。
template<int blockx, int blocky>
__global__ void updateCluster(PtrStepSz<uchar4> source, PtrStepSz<uchar> clusterIndex, PtrStepSz<uchar> mask, float4* kencter, int* kindexs)
{
    __shared__ float3 centers[blockx*blocky][CUDA_GRABCUT_K2];
    __shared__ int indexs[blockx*blocky][CUDA_GRABCUT_K2];
    const int idx = blockDim.x * blockIdx.x + threadIdx.x;
    const int idy = blockDim.y * blockIdx.y + threadIdx.y;
    const int threadId = threadIdx.x + threadIdx.y * blockDim.x;
#pragma unroll CUDA_GRABCUT_K2
    for (int i = 0; i < CUDA_GRABCUT_K2; i++)
    {
        centers[threadId][i] = make_float3(0.f);
        indexs[threadId][i] = 0;
    }
    __syncthreads();
    if (idx < source.cols && idy < source.rows)
    {
        //全部值都放入共享centers
        int index = clusterIndex(idy, idx);
        bool bFg = checkFg(mask(idy, idx));
        int kindex = bFg ? index : (index + CUDA_GRABCUT_K);
        float4 color = rgbauchar42float4(source(idy, idx));
        centers[threadId][kindex] = make_float3(color);
        indexs[threadId][kindex] = 1;
        __syncthreads();
        //每一個線程塊進行二分聚合,每次操做都保存到前一半數組裏，直到最後保存在線程塊裏第一個線程上(這塊比較費時,0.1ms)
        for (uint stride = blockDim.x*blockDim.y / 2; stride > 0; stride >>= 1)
        {
            //int tid = (threadId&(stride - 1));
            if (threadId < stride)//stride 2^n
            {
#pragma unroll CUDA_GRABCUT_K2
                for (int i = 0; i < CUDA_GRABCUT_K2; i++)
                {
                    centers[threadId][i] += centers[threadId + stride][i];
                    indexs[threadId][i] += indexs[threadId + stride][i];
                }
            }
            //if (stride > 32)
            __syncthreads();
        }
        //每塊的第一個線程集合，把共享centers存入臨時顯存塊上kencter
        if (threadIdx.x == 0 && threadIdx.y == 0)
        {
            int blockId = blockIdx.x + blockIdx.y * gridDim.x;
#pragma unroll CUDA_GRABCUT_K2
            for (int i = 0; i < CUDA_GRABCUT_K2; i++)
            {
                int id = blockId * 2 * CUDA_GRABCUT_K + i;
                kencter[id] = make_float4(centers[0][i], 0.f);
                kindexs[id] = indexs[0][i];
            }
        }
    }
}

//block 1*1,threads(暫時選32*4),對如上gridDim.x*gridDim.y的數據用blockx*blocky個線程來處理
template<int blockx, int blocky>
__global__ void updateCluster(float4* kencter, int* kindexs, kmeansI& meansbg, kmeansI& meansfg, int& delta, int gridcount)
{
    __shared__ float3 centers[blockx*blocky][CUDA_GRABCUT_K2];
    __shared__ int indexs[blockx*blocky][CUDA_GRABCUT_K2];
    const int threadId = threadIdx.x + threadIdx.y * blockDim.x;
#pragma unroll CUDA_GRABCUT_K2
    for (int i = 0; i < CUDA_GRABCUT_K2; i++)
    {
        centers[threadId][i] = make_float3(0.f);
        indexs[threadId][i] = 0;
    }
    __syncthreads();
    //int gridcount = gridDim.x*gridDim.y;
    int blockcount = blockDim.x*blockDim.y;
    int count = gridcount / blockcount + 1;
    //每塊共享變量，每一個線程操縱本身對應那塊地址，不須要同步,用線程塊操做一個大內存
    for (int i = 0; i < count; i++)
    {
        int id = i*blockcount + threadId;
        if (id < gridcount)
        {
#pragma unroll CUDA_GRABCUT_K2
            for (int j = 0; j < CUDA_GRABCUT_K2; j++)
            {
                int iid = id * CUDA_GRABCUT_K2 + j;
                centers[threadId][j] += make_float3(kencter[iid]);
                indexs[threadId][j] += kindexs[iid];
            }
        }
    }
    __syncthreads();
    for (uint stride = blockDim.x*blockDim.y / 2; stride > 0; stride >>= 1)
    {
        if (threadId < stride)
        {
#pragma unroll CUDA_GRABCUT_K2
            for (int i = 0; i < CUDA_GRABCUT_K2; i++)
            {
                centers[threadId][i] += centers[threadId + stride][i];
                indexs[threadId][i] += indexs[threadId + stride][i];
            }
        }
        //if (stride > 32)
        __syncthreads();
    }
    if (threadIdx.x == 0 && threadIdx.y == 0)
    {
        int count = 0;
        //收集全部信息，並從新更新cluster，記錄變動的大小
        for (int i = 0; i < CUDA_GRABCUT_K; i++)
        {
            meansfg.kcenters[i] = make_float4(centers[0][i], 0.f);
            if (indexs[0][i] != 0)
                meansfg.cluster[i] = meansfg.kcenters[i] / indexs[0][i];
            count += abs(indexs[0][i] - meansfg.length[i]);
            meansfg.length[i] = indexs[0][i];
        }
        for (int i = CUDA_GRABCUT_K; i < CUDA_GRABCUT_K2; i++)
        {
            meansbg.kcenters[i - CUDA_GRABCUT_K] = make_float4(centers[0][i], 0.f);
            if (indexs[0][i] != 0)
                meansbg.cluster[i - CUDA_GRABCUT_K] = meansbg.kcenters[i - CUDA_GRABCUT_K] / indexs[0][i];
            count += abs(indexs[0][i] - meansbg.length[i - CUDA_GRABCUT_K]);
            meansbg.length[i - CUDA_GRABCUT_K] = indexs[0][i];
        }
        delta = count;
    }
}
void updateCluster_gpu(PtrStepSz<uchar4> source, PtrStepSz<uchar> clusterIndex, PtrStepSz<uchar> mask, float4* kencter,
    int* kindexs, kmeansI& meansbg, kmeansI& meansfg, int& delta, cudaStream_t stream = nullptr)
{
    dim3 grid(cv::divUp(source.cols, block.x), cv::divUp(source.rows, block.y));
    updateCluster<BLOCK_X, BLOCK_Y> << <grid, block, 0, stream >> > (source, clusterIndex, mask, kencter, kindexs);
    updateCluster<BLOCK_X, BLOCK_Y> << <1, block, 0, stream >> > (kencter, kindexs, meansbg, meansfg, delta, grid.x*grid.y);
}
void KmeansCuda::compute(GpuMat source, GpuMat clusterIndex, GpuMat mask, int threshold)
{
    int delta = threshold + 1;
    initKmeans_gpu(*bg, *fg);
    while (delta > threshold)
    {
        findNearestCluster_gpu(source, clusterIndex, mask, *bg, *fg);
        updateCluster_gpu(source, clusterIndex, mask, dcenters, dlenght, *bg, *fg, *d_delta);
        //能夠每隔五次調用一次這個
        cudaMemcpy(&delta, d_delta, sizeof(int), cudaMemcpyDeviceToHost);
    }    
}

updateCluster

　　

　　這裏確定有點反直覺，十多行代碼擴展到百多行，效率還增長了十倍，可是GPU運算就是這樣，每步能多利用核心就行，還有，cudaMemcpy花費真的大，這裏有個失敗的嘗試，由於在這基礎上，都是用更多的計算核心帶來的好處，就想把這個算法再擴展了下，前景與背景二個各分五類，一共十類，故當前準備用block的z爲10，當作每塊的索引，和我想象不一樣，這個要的時間在0.4ms左右，時間關係，沒有繼續測試，後續針對這塊再仔細分析比對下。

　　以下顯示的效果圖。

　　

　　而後就是GMM的實現，這個實現和k-means過程比較類似，代碼就不貼了，當時效果完成後，個人主要疑問在若是確認個人實現是對的，可是下面這張圖出來後，我就知道應該是對了。GMM對k-means每塊求的結果，而後用這結果來獲得顏色的分類，這兩者效果應該差很少是同樣，可是確定又有點區別(image是k-means的結果，seg image是GMM一次learning而後再分配的結果)，固然最後確認仍是我比較opencv裏GMM裏各個分類與我實現GMM的值的比較，通常來講，各個分類的平均值與個數相差不大，各個分類在全部點的比例相似，固然兩者的結果不會徹底同樣，opencv中k-means中用k-means++隨機選擇五個差別比較大的中心，就算用opencv本身重複算同一張圖，每次結果都會有細小區別。

　　

　　最後是maxflow的實現，這個過程我參照了https://web.iiit.ac.in/~vibhavvinet/Software/ 的實現，要說這份代碼稍稍改下仍是能運行的，就是數據來源全在一個文件上，包含formsource,tosink,edge的信息，主要查看push與relabel與bfs的實現，能夠去掉一些不少可有可無的代碼，而且與opencv對接，全部顯存全用gpumat代碼，差很少只有200行左右的代碼。記的用這裏的數據如person/sponge大約50屢次迭代push-relabel後就能獲得很好的結果(50次640*480相似的圖只須要5-8ms,若是是416*416會更低，大約0.07ms就能一次迭代)，在這給了我信心能繼續作下去，固然後續結果另說。

　　後續就是一個組合過程，組合過程參照opencv的實現。

void GrabCutCude::renderFrame(GpuMat x3source, cv::Rect rect)
{
    cudaDeviceSynchronize();

    rgb2rgba_gpu(x3source, source);
    setMask_gpu(mask, rect);
    kmeans->compute(source, clusterIndex, mask, 0);    //testShow(showSource, clusterIndex, "image");
    gmm->learning(source, clusterIndex, mask);    //testShow(showSource, clusterIndex, "seg image");
    int index = 0;
    while (index++ < iterCount)
    {
        gmm->assign(source, clusterIndex, mask);    
        gmm->learning(source, clusterIndex, mask);
        graph->addTermWeights(source, mask, *(gmm->bg), *(gmm->fg), lambda);
        graph->addEdges(source, gamma);
        graph->maxFlow(mask, maxCount);
    }
}

renderFrame

　　

　　其中有幾處要求聚合運算，分別是求maxflow邊的權重時，整張圖的beta值，k-means更新聚合，gmm的訓練，其中求beta只費了0.02ms,而k-means在0.12ms,而gmm的訓練在1.2ms左右，beta只是把全部數據相加，每一個像素對應4byte的共享存儲，而k-means把數據取成十個類別，每一個像素對應對應4byte*(3+1)*10個字節，時間花費0.02/0.12應該是一個正常範圍，而gmm與k-means運算差很少徹底同樣，只是共享顯存佔用4byte*(3+3*3+1)*10個字節，時間花費比k-means多了十倍。這個比較頗有意思，共享顯存用越多，能夠看到，GPU利用率會變低，求beta時，occupancy能夠到達100%，而在k-meansk中只有25%，updateGMM只有6.25%。從前面k-means的代碼看，每一個聚合分紅二部分，經調整GMM的先後部分發現，上部分下降線程塊block能夠下降時間，而下部分升高線程塊block能夠下降時間，其實這麼算的話GMM的計算能夠考慮分開算，讓k-means的花費來看，最多0.3ms左右，1.2ms確定還有不少優化空間。

　　以下結合咱們本身編譯的darknet動態連接庫組合的代碼。

void testYoloGrabCutCuda()
{
    cv::VideoCapture cap;
    cv::VideoWriter video;
    cv::Mat frame;
    cv::Mat result, reframe;

    cap.open(0);
    char* modelConfiguration = "yolov3.cfg";
    char* modelWeights = "yolov3.weights";

    network *net = load_network(modelConfiguration, modelWeights, 0);
    layer l = net->layers[net->n - 1];
    set_batch_network(net, 1);

    cv::Mat netFrame(net->h, net->w, CV_32FC3);

    const char* window_name = "yolo3";
    //namedWindow(window_name);
    cv::namedWindow("image");
    cv::namedWindow("seg image");
    createUI("1 image");

    int height = net->w;
    int width = net->h;
    GrabCutCude grabCut;
    grabCut.init(width, height);
    GpuMat gpuFrame;
    while (cv::waitKey(1) < 1)
    {
        updateUI("1 image", grabCut);

        cap >> frame;
        cv::resize(frame, reframe, cv::Size(net->w, net->h));
        //reframe.convertTo(netFrame, CV_32FC3, 1 / 255.0);
        image im = mat_to_image(reframe);
        float *predictions = network_predict(net, im.data);
        int nboxes = 0;
        float thresh = 0.7f;
        detection curdet = {};
        float maxprob = thresh;
        bool bFind = false;
        detection *dets = get_network_boxes(net, im.w, im.h, thresh, 0, 0, 0, &nboxes);
        for (int i = 0; i < nboxes; i++)
        {
            auto det = dets[i];    //0 person 39 bottle        
            if (det.prob[39] > maxprob)
            {
                maxprob = det.prob[39];
                curdet = det;
                bFind = true;
            }
        }

        if (!bFind)
            continue;
        int width = curdet.bbox.w + 20;
        int height = curdet.bbox.h + 20;
        cv::Rect rectangle(curdet.bbox.x - width / 2.f, curdet.bbox.y - height / 2.f, width, height);

        if (bCpu)
        {
            cv::Mat bgModel, fgModel;
            cv::grabCut(reframe, result, rectangle, bgModel, fgModel, iterCount, cv::GC_INIT_WITH_RECT);
        }
        else
        {
            grabCut.setParams(iterCount, gamma, lambda, maxCount);
            cudaDeviceSynchronize();
            gpuFrame.upload(reframe);
            grabCut.renderFrame(gpuFrame, rectangle);
            grabCut.mask.download(result);
        }

        cv::compare(result, cv::GC_PR_FGD, result, cv::CMP_EQ);
        cv::Mat foreground(reframe.size(), CV_8UC3, cv::Scalar(0, 0, 0));
        reframe.copyTo(foreground, result);
        cv::imshow("image", foreground);
        cv::rectangle(reframe, rectangle, cv::Scalar(0, 0, 255), 3);
        cv::imshow("seg image", reframe);
    }
}

testYoloGrabCutCuda

　　這段代碼就是前面顯示的結果，可是最後我發現，效果並很差，至少相對我參考那段maxflow效果來講，迭代五十次只須要5ms,再稍微優化下，k-means加上GMM能夠控制在3ms內，這樣加上yolo3也能夠控制在30楨左右，可是，就前面那個杯子來講，須要迭代超過300次纔能有個比較好的效果，若是是人，須要迭代上千次才能達到opencv的效果，那就沒啥意義了，後續準備的優化與整合進引擎也暫停。

　　至於爲啥會這樣，我分析了下，原參考push-relabel裏數據是用種子生成的方式來生成formsource，tosink，這樣formsource與tosink原本就接近結果了。

　　

　　而按照opencv裏的流程，出來的formsource與tosink圖以下。

　　

　　這樣就須要迭代更屢次的push_relabel來完成，主要由於這二次方式生成的GMM有很大區別，種子點生成的GMM各個分類與佔比原本就好，我試驗了下，擴大生成的邊框，立刻push_relabel的迭代次數加個百次，這是由於前景中混成來的背景越多，背景分類與佔比越大，不肯定區顏色經過GMM生成的formsource與tosink差值越小，雖然push-relabel每一個點能夠獨立計算，很適合GPU算，可是他的push與relabel方式致使他每次可能就少數點在流動，大多點根本不知足要求。

　　那麼要考慮在100次迭代內獲得比較好的效果，須要種子點方式，加上一些特定需求，好比背景與人變更不大，只是人在背景中的位置在不斷變化，能夠先用幾楨算mask-rcnn獲得mask,用mask獲得相對準確的GMM模型，再把這個GMM模型算對應graph的點權重，而後用maxflow算的比較清晰的邊，最後用導向濾波還原大圖，這模式對如今的流程有些改動，故此記錄下如今流程，若是後續有比較好的效果再來繼續說。

　　2019/3/20號更新：

　　爲了驗證個人想法以及評論區有同窗提出formsource與tosink用同一例子比較，深覺得然。下面是先用幾楨效果的圖算出GMM，而後用這個GMM值直接算formsource與tosink，以下是對應的formsource與tosink圖。

　　

void GrabCutCude::renderFrame(GpuMat x3source)
{
    rgb2rgba_gpu(x3source, source);
    graph->addEdges(source, gamma);
    int index = 0;
    while (index++ < iterCount)
    {
        graph->addTermWeights(source, *(gmm->bg), *(gmm->fg), lambda);
        graph->maxFlow(mask, maxCount);
    }
}

renderFrame

　　相對於原先方框裏算k-means,GMM模型,這個是在已經獲得GMM模型的狀況下，用GMM直接算formsource與tosink，能夠看到，相對於原來來講，結果已經很乾淨了，這種狀況下，50次push-relabel就能獲得比較準確的結果。