cuda bilinear 紋理記憶體、共享記憶體實作與解說

之前介紹過一篇

• [原始碼] C / C++ 線性插補 bilinear 與注意事項

裡面詳細介紹注意事項與如何實作，大致簡單說一下就是等比取值

如果 A=0, B=10 ，那麼 AB=1 大概就是這個原理，如果是二維平面4個點就是取3次剛剛的過程，先兩次兩個點，然後再把兩個點的結果取一次。

CPU函式

上面那個敘述就是我們的主要函式了，大致可以寫成如下程式

__host__ __device__
inline static float bilinearRead(const float* img, 
    size_t width, float y, float x) // 線性取值
{
    // 獲取鄰點(不能用 1+)
    size_t x0 = floor(x);
    size_t x1 = ceil(x);
    size_t y0 = floor(y);
    size_t y1 = ceil(y);
    // 獲取比例(只能用 1-)
    float dx1 = x - x0;
    float dx2 = 1 - dx1;
    float dy1 = y - y0;
    float dy2 = 1 - dy1;
    // 獲取點
    const float& A = img[y0*width + x0];
    const float& B = img[y0*width + x1];
    const float& C = img[y1*width + x0];
    const float& D = img[y1*width + x1];
    // 乘出比例(要交叉)
    float AB = A*dx2 + B*dx1;
    float CD = C*dx2 + D*dx1;
    float X = AB*dy2 + CD*dy1;
    return X;
}

要操作這個函式大致上就用兩個for用x, y跑過圖中所有的點就可以了

__host__ void biliner_CPU_core(vector<float>& img, const vector<float>& img_ori, 
    size_t width, size_t height, float Ratio)
{
    int newH = static_cast<int>(floor(height * Ratio));
    int newW = static_cast<int>(floor(width  * Ratio));
    img.resize(newH*newW);
    // 跑新圖座標
    for (int j = 0; j < newH; ++j) {
        for (int i = 0; i < newW; ++i) {
            // 調整對齊
            float srcY, srcX;
            if (Ratio < 1) {
                srcY = ((j+0.5f)/Ratio) - 0.5;
                srcX = ((i+0.5f)/Ratio) - 0.5;
            } else {
                srcY = j * (height-1.f) / (newH-1.f);
                srcX = i * (width -1.f) / (newW-1.f);
            }
            // 獲取插補值
            img[j*newW + i] = bilinearRead(img_ori.data(), width, srcY, srcX);
        }
    }
}

共享記憶體

共享記憶體這邊有引用到之前寫的函式庫，詳細可以參考這篇站內文
https://charlottehong.blogspot.tw/2018/01/cuda.html

// 共享記憶體線性取值核心
__global__ void biliner_share_kernel(float* dst, const float* src, int srcW, int srcH, float ratio) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j = blockIdx.y * blockDim.y + threadIdx.y;
    int newH = (int)(floor(srcH * ratio));
    int newW = (int)(floor(srcW * ratio));
    if(i < srcW*ratio && j < srcH*ratio) { // 會多跑一點點要擋掉
        // 調整對齊
        float srcY, srcX;
        if (ratio < 1) {
            srcY = ((j+0.5f)/ratio) - 0.5;
            srcX = ((i+0.5f)/ratio) - 0.5;
        } else {
            srcY = j * (srcH-1.f) / (newH-1.f);
            srcX = i * (srcW -1.f) / (newW-1.f);
        }
        // 獲取插補值
        dst[j*newW + i] = bilinearRead(src, srcW, srcY, srcX);
    }
}

差不多就是從CPU那個函式改過來的而已，使用時候要先new兩個空間

// 共享記憶體線性取值函式
__host__ void biliner_share_core(float *dst, const float* src,
    size_t srcW, size_t srcH, float ratio)
{
    Timer T; T.priSta = 1;
    // 設置GPU所需長度
    int srcSize = srcW*srcH;
    int dstSize = srcSize*ratio*ratio;

    // 要求GPU空間
    T.start();
    CudaData<float> gpu_src(srcSize);
    T.print("  GPU new 空間1");
    T.start();
    CudaData<float> gpu_dst(dstSize);
    T.print("  GPU new 空間2");
    // 複製到GPU
    T.start();
    gpu_src.memcpyIn(src, srcSize);
    T.print("  GPU 複製");

    // 設置執行緒
    dim3 block(BLOCK_DIM, BLOCK_DIM);
    dim3 grid(ceil((float)srcW*ratio / BLOCK_DIM), ceil((float)srcH*ratio / BLOCK_DIM));
    T.start();
    biliner_share_kernel <<< grid, block >> > (gpu_dst, gpu_src, srcW, srcH, ratio);
    T.print("  核心計算");

    // 取出GPU值
    T.start();
    gpu_dst.memcpyOut(dst, dstSize);
    T.print("  GPU 取出資料");

    // 釋放GPU空間
    T.start();
    gpu_src.~CudaData();
    gpu_dst.~CudaData();
    T.print("  GPU 釋放空間");
}

紋理記憶體

這個速度最快代碼也最少，可以直接調用cuda內優化過的方法

// 宣告GPU紋理變數(只能放全域)
texture<float, 2, cudaReadModeElementType> rT;
// 紋理線性取值核心
__global__ void biliner_texture_kernel(float* dst, int srcW, int srcH, float ratio) {
    int idxX = blockIdx.x * blockDim.x + threadIdx.x,
        idxY = blockIdx.y * blockDim.y + threadIdx.y;
    if(idxX < srcW*ratio && idxY < srcH*ratio) { // 會多跑一點點要擋掉
        float srcX = idxX / ratio;
        float srcY = idxY / ratio;
        size_t idx = (idxY*srcW*ratio + idxX);
        dst[idx] = tex2D(rT, srcX+0.5, srcY+0.5);
    }
}

使用

__host__ void biliner_texture_core(float *dst, const float* src,
    size_t dstW, size_t dstH, float ratio)
{
    Timer T; T.priSta = 1;
    // 設置GPU所需長度
    int srcSize = dstW*dstH;
    int dstSize = srcSize*ratio*ratio;

    // 宣告 texture2D陣列並綁定
    T.start();
    CudaMemArr<float> cuArray(src, dstW, dstH);
    cudaBindTextureToArray(rT, cuArray);
    T.print("  GPU new 紋理空間+複製");

    // 設置 插植模式and超出邊界補邊界
    rT.filterMode = cudaFilterModeLinear;
    rT.addressMode[0] = cudaAddressModeClamp;
    rT.addressMode[1] = cudaAddressModeClamp;

    // 要求GPU空間
    T.start();
    CudaData<float> gpu_dst(dstSize);
    T.print("  GPU new 一般空間");

    // 設置執行緒
    dim3 block(BLOCK_DIM, BLOCK_DIM);
    dim3 grid(ceil((float)dstW*ratio / BLOCK_DIM), ceil((float)dstH*ratio / BLOCK_DIM));
    T.start();
    biliner_texture_kernel <<< grid, block >> > (gpu_dst, dstW, dstH, ratio);
    T.print("  核心計算");

    // 取出GPU值
    T.start();
    gpu_dst.memcpyOut(dst, dstSize);
    T.print("  GPU 取出資料");
}

總結

可以看出大部分時間其實是消耗再第一次new空間的時候，至於記憶體的複製貌似也沒消耗太多時間，建議圖可以直接讀進共享記憶體就好了，不需要先用主機記憶體讀然後再複製進去，從一開始的時候就宣告共享記憶體。

texture times = 0.327s
share times = 0.4s
CPU times = 2.578s

範例用的原始碼：https://github.com/hunandy14/Cuda_ImgBilinear

​

資料壓縮 向量量化 方法解說與實作

向量量化原理解析(未訓練)

首先我們回需要三分檔案分別是

1. 原圖：source.raw
2. 編碼簿：origin.raw
3. 索引值：idx.raw

原圖

原圖是256x256的大小，並且在實做中最小單位是一個區塊，一個區塊由4x4組合，原圖可以看成是是一個64x64大小的圖。

編碼簿

編碼簿也是以區塊為單位，其中第一次是從原圖隨機抽取256個區塊。

索引值

索引值大小是區塊化的原圖大小也就是64x64，它的獲取需要輸入sou與ori獲得，獲取的過程是以原圖為主，原圖的第一個區塊 sou.bolock[0] 去比較 ori.bolock[0~255] 跟哪一個最像獲得。

假如 sou.bolock[i] 與 ori.bolock[2] 最像，那麼就在 idx[i] 填入 2，依序填完全部

最像的定義是這樣的兩個區塊內16個像素，個別的差值平方的總和 而根據這個數據，越低代表越像。

void block_diff(unsigned char a[16], unsigned char b[16]){
    size_t sum=0;
    for(unsigned i = 0; i < 16; ++i) {
        sum += pow(abs(a[i]-b[i]), 2);
    }
}

合併回原圖

我們有了編碼簿與索引值之後就可以拼回原圖了，依照索引值內的編號將編碼簿一個區塊一個區塊慢慢補回去即可。

for(unsigned i = 0; i < 4096; ++i) {
    img.block[i] = ori.block[idx[i]];
}

訓練編碼簿

隨機抽取的編碼簿還原之後的圖與原圖的差距是很大的，為了降低這種誤差，要盡可能的讓編碼簿的碼優化。優化步驟如下。

1. 計算新的編碼簿區塊
2. 計算新的索引值
3. 疊代至收斂

計算新的編碼簿區塊

看一下圖中的 sou.raw 上面的編號是從 idx.raw 抄上來的，在這上面可以發現有兩個3號，分別是 sou.block[0] 與 sou.block[2] 他們之所以獲得相同的編號是因為，他們各自與編碼簿中的256個比較後，與 ori.block[3] 很像。

理解這個原理就可以發現 ori.block[3] 存在著優化空間，它可以取 sou.block[0] 與 sou.block[2] 的平均獲得更精準的還原；平均指的是16個點個別相加/2，獲得全新的16個點。

void tra_block(int idxNum){
    long long unsigned int s[16]{};
    size_t cnt=0;
    // 找出相同索引並累加
    for(unsigned i = 0; i < 4096; ++i) {
        if(idx[i]==idxNum) {
            for(unsigned k = 0; k < 16; ++k, ++cnt) {
                s[k] += sou.block[i][k];
            }
        }
    }
    // 根據找出的數，算出平均數並填入tra
    if (cnt != 0){
        // 算出平均
        for (int i = 0; i < 16; ++i)
            s[i] = s[i]/cnt;
        // 回填
        for (unsigned k = 0; k < 16; ++k)
            ori.block[idxNum][k] = s[k];
    }
}

ori 裡面總共有 256 個區塊，把這 256 個都做過一次優化就可以了

for(unsigned j = 0; j < 256; ++j) {
    tra_block(j);
}

計算新的索引值

有了優化過的編碼簿之後就可以重新計算新的索引值，與前面所敘述的算法是一樣的，只不過在過程中要多存一個數據。

比對的過程中每個區塊所算出來的差平方和，一共是4096個要把取平均值存下來。

疊代至收斂

根據索引值計算時所產出的 avg 做計算，依據公式需要用到本次與上次的avg

計算 avg差值/當前avg 不小於 0.01 則回到[1]

可以看到優化過後的圖片，與原本的沒有優化好了很多。

實作程式碼

GitHub連結：向量量化

​

什麼是物件導向，的差別在哪裡

什麼是物件導向，差別在哪裡

用一句話來形容就是

資料與方法的群組化再群組化

如果說編程是解放 繁瑣與重複的工作，那麼物件導向就是解放編程的繁瑣與重複的工作

下面讓我用簡單的群組概念來解釋如何解放繁瑣的工作，當然~實際上物件導向不只如此，舉例只是這個海灘的中的其中一粒沙子，但是這已經有足夠強大的功能與便利性了。

讓我們看看一個單向鍊結的實例，下面是一個用C語言實做的例子

函式

// 資料結構
typedef struct node Node;
struct node {
    int data;
    Node* next;
};
// 初始化節點
Node* link_init(int num) {
    Node* n = (Node*)malloc(sizeof(Node));;
    n->data = num;
    n->next = NULL;
    return n;
}
// 尾端新增節點
Node* link_append(Node* node, int num) {
    node->next = link_init(num);
    return node->next;
}
// 插入節點
void link_insert(Node* node, size_t idx, int data){
    // 旗標(插入點的前一個)
    for(unsigned i = 0; i < idx; ++i) {
        node = node->next;
    }
    Node* temp = node->next;
    node->next = link_init(data);
    node->next->next = temp;
}

操作範例

    // 新增鏈結
    Node* head = link_init(-1);
    Node* tail = head;

    // 新增節點
    tail = link_append(tail, 1);
    tail = link_append(tail, 2);
    tail = link_append(tail, 3);
    // 插入節點
    link_insert(head, 1, 0);
    // 查看節點
    for(Node* n = head->next; n; n=n->next) {
        printf("%d, ", n->data);
    } printf("\n");

    // 刪除鏈結
    link_delete(head);
    head=NULL;
    return 0;

在這個例子裏面可以看到有三件繁瑣的事情

1. 每一個函式呼叫都需要引入 Node*
2. 為了方便識別函式是對屬於誰的，名稱帶前缀名
3. 結束之後自己要記得 delete

由於這是一個簡單的範例所有只有幾個函式看起並不多餘，如果這個鏈結陣列裝的東西不只是int可能會需要很多個函式來操作處理，這時候每一個函式參數都要帶一個 Node* 就顯得有些多餘。

那物件導向是怎麼優化這件事情的呢，讓我們來看一下物件導向的寫法。

函式

struct Node_class{
    // 資料結構
    struct Node {
        int data;
        Node* next;
    } *node, *tail;

    // 初始化
    Node_class(){
        node = init(-1);
        tail = node;
    }
    // 自動刪除鏈結
    ~Node_class(){
        del();
    }

    // 初始化節點
    static Node* init(int num) {
        Node* node = (Node*)malloc(sizeof(Node));;
        node->data = num;
        node->next = NULL;
        return node;
    }
    // 尾端新增節點
    void append(int num) {
        tail->next = init(num);
        tail = tail->next;
    }
    // 插入節點
    void insert(size_t idx, int data){
        // 旗標(插入點的前一個)
        Node* node = this->node;
        for(unsigned i = 0; i < idx; ++i) {
            node = node->next;
        }
        Node* temp = node->next;
        node->next = init(data);
        node->next->next = temp;
    }
    // 刪除鏈結
    void del(){
        for(Node* temp=node; temp; temp=node) {
            node = temp->next;
            free(temp);
        }
    }
    // 印出節點
    void print(){
        for(Node* n = node->next; n; n=n->next) {
            printf("%d, ", n->data);
        } printf("\n");
    }
};

操作方法

    // 新增鏈結
    Node_class list;

    // 新增節點
    list.append(1);
    list.append(2);
    list.append(3);
    // 插入節點
    list.insert(1, 0);
    // 查看節點
    list.print();

    // 刪除鏈結
    // list.~Node_class();

再來回頭看那三點

1. 每一個函式呼叫都需要引入 Node*

   由於方法與資料的群組結構，每一個物件都帶有各自的資料，而這個物件可以呼叫這些方法，用的就是內帶的資料，每個物件帶有不同的資料。

2. 為了方便識別函式是對屬於誰的，名稱帶前缀名

   物件的名字取代了函式的前缀名，這個物件是誰建立的就傭有哪些方法，如果不是則無法使用，我們再也不需要為了同名而困擾了。

3. 結束之後自己要記得 delete

   之所以將她註解掉是因為並沒有影響，每一個物件再生命週期結束的時候都會自動呼叫一次解構子，也就是帶有~符號的類別名稱方法。

​

cuda 自動管理記憶體 函式庫

因為還沒有實際在專案上跑過，可能有bug或是考慮不周的情況發生，再自行修改~

另外目前最新版 VS2017 還不支援，有兩個辦法可以處理

1. 使用VS2015
2. 安裝舊版VS2017並不更新
3. 在最新版2017新增2015的套件，讓方案用2015跑

3看起來是最好的，不過實際運行有一些bug，每次編譯完之後更改程式碼在按除錯並不會編譯新的程式碼會直接執行舊的exe，清除重編可以不過每次都按會很麻煩。最佳手段是選1或2。

詳情參考：https://charlottehong.blogspot.tw/2017/12/visutal-studio-2017-cuda90.html

cuda 通常流程

宣告空間

新建 共享記憶體 這一個cpu跟gpu都可以用，一般都是用這個傳入gpu

T* gpuData;
cudaMalloc((void**)&gpuData, size*sizeof(T));

輸入資料

cudaMemcpy(gpuData, dataIn, size*sizeof(T), cudaMemcpyHostToDevice);

輸出資料

cudaMemcpy(dataIn, gpuData, size*sizeof(T), cudaMemcpyDeviceToHost);

cuda 記憶體控制類別

上述流程繁雜，不過基本上就是跟C語言一樣，為了寫作方便可以讓一個class來管理記憶體，下面是我寫的類別。

稍微有點長我把它放到外部空間gist

• Cuda記憶體管理
• 時間計時函式

使用

基本用法都包含進去，需要什麼在自己包什麼。
下面是使用範例

/*****************************************************************
Name : 
Date : 2018/01/08
By   : CharlotteHonG
Final: 2018/01/08
*****************************************************************/
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <iostream>
#include <vector>
using namespace std;

#include "timer.hpp"
#include "CudaData.cuh"

__global__ void cudacopy(float* b, float* a, int size){
    // 乾式寫法
    const int idx = blockIdx.x*blockDim.x + threadIdx.x;
    if(idx<size){
        b[idx]=a[idx];
    }
    // 迴圈寫法
    //for(int i=threadIdx.x; i<size; i+=blockDim.x){
        //b[i]=a[i];
    //}
}
void cpucopoy(float* b, float* a, int size) {
    for(int i=0; i<size; ++i){
        b[i]=a[i];
    }
}
void testCuda(size_t size) {
    Timer T;
    // 配置主機記憶體
    vector<float> img_data(size), cpu_data(size), gpu_data(size);
    float* a = img_data.data(); // 原始資料
    float* b = cpu_data.data(); // CPU計算後資料
    float* c = gpu_data.data(); // GPU輸出回來資料
    // 設置初值
    float test_val=7;
    for(int i=0; i<size; i++){
        a[i]=test_val;
    }
    // 配置顯示記憶體, 載入資料.
    T.start();
    CudaData<float> gpuDataIn(a, size), gpuDataOut(size);
    gpuDataOut.memset(0, size);
    T.print("  Cuda Data malloc and copy");

    // 網格區塊設定. (與 kernel for 的次數有關)
    const size_t blkDim=16;
    int grid(size/blkDim+1);  // 網格要含蓋所有範圍, 所以除完要加 1.
    int block(blkDim);        // 區塊設定 16x16.
    // Cuda Kernel 執行運算
    T.start();
    cudacopy<<<grid,block>>>(gpuDataOut, gpuDataIn, size);
    T.print("  Cuda-copy");
    // 取出GPU資料
    gpuDataOut.memcpyOut(c, size);

    // CPU 執行運算
    T.start();
    cpucopoy(b, a, size);
    T.print("  Cpu-copy");
    // 測試
    bool f=0;
    for(size_t i = 0; i < size; i++) {
        if(c[i] != b[i]) {f=1;}
    }
    // 測試報告
    if(f==0) {
        cout << "test ok" << endl;
    } else {
        cout << "test Error" << endl;
    }
}
int main(){
    Timer T;
    testCuda(99999);
    T.print("ALL time.");
    return 0;
}

​