碼上歡樂
相關內容    簡體    繁體

最遠點采樣介紹及CUDA實現分析

本文轉載自   查看原文   2021-04-17 21:58   321 

最遠點采樣介紹及CUDA實現分析

最遠點采樣（Farthest Point sampling/FPS）是一個基本的點雲采樣算法，在各類點雲處理算法中都有使用，如PointNet++，以及各類三維物體檢測算法。

本文從以下幾個方面對FPS算法進行介紹和分析

  • FPS邏輯描述
  • FPS算法串行實現與分析
  • FPS算法並行實現與分析
  • 串行實現與並行實現的性能比較

1. FPS邏輯描述

假設有\(n\)個點，要從中按照FPS算法，采樣出\(k\)個點（\(k<=n\)）。那么，可以將所有點歸類到兩個集合\(A,B\)中，\(A\)集合表示選中的點形成的集合，\(B\)集合表示未選中的點構成的集合。FPS所做的事情就是每次從集合\(B\)中選取一個到集合\(A\)里的點距離最大的點。

初始情況：集合\(A\)為空，集合\(B\)包含所有點

選第一個點：對所有點進行shuffle后，選取第一個點，將其加入集合\(A\)中。此時\(size(A)=1,size(B)=n-1\)

選第二個點：分別計算出集合\(B\)里面每個點到集合\(A\)中一個點的距離，選取距離最大的點，加入集合\(A\)。此時，\(size(A)=2, size(B)=n-1\)。

選第三個點及后續的點：設此時集合A中有\(m\)個點（\(m>=2\)），集合\(B\)中有\(n-m\)個點。此時需要定義集合\(B\)中的點\(p_B\)到集合\(A\)中點的距離，注意到集合\(A\)中不止一個點，這是FPS的核心。對\(p_B\)到集合\(A\)的距離\(d_B\)定義如下

\[\begin{aligned} d_B &= d(p_B,A) = min(d(p_B, p^j_A)) \quad (p_B \in B, p^j_A \in A) \end{aligned} \]

  1. 分別計算出\(p_B\)到集合\(A\)中每個點的距離。當集合\(A\)中有\(m\)個點時，可以計算出\(m\)個距離值
  2. 從計算出來的\(k\)個距離值中，取最小的距離值，作為點\(p_B\)到集合\(A\)的距離

對於集合\(B\)中的\(n-m\)個點，每個點都可以計算出一個距離值：\(\{d^1_B,d^2_B,...,d^{n-m}_B\}\)。選出最大距離值對應的點，記為第\(m+1\)個點，將其加入集合\(A\)。此時集合\(A\)包含\(m+1\)個點，集合\(B\)包含\(n-(m+1)\)個點。

重復上述步驟3，直至選出\(k\)個點為止。

2. FPS算法串行實現與分析

當集合\(A\)中有\(m\)個點時，集合\(B\)中有\(n-m\)個點，此時選取下一個點的時間復雜度為\((n-m)*m\)

顯然該步驟可以優化，分析如下：

  • 選第\(m+1\)個點時，對於集合\(B\)中的一個點\(p_B\)，需要計算出到集合\(A\)中\(m\)個點的距離。假設\(o^1_B\)表示點\(p_B\)到集合\(A\)中第一個點的距離，那么需要計算的距離為：\(\{o^1_B, o^2_B, ..., o^m_B\}\)，然后取最小值\(min(\{o^1_B, o^2_B, ..., o^m_B\})\)，作為\(p_B\)到集合\(A\)的距離。
  • 考慮選第m個點的過程，對於集合中的一個點\(p_B\)，需要計算出到集合\(A\)中\(m-1\)個點的距離，需要計算的距離為\(\{o^1_B, o^2_B,...,o^{m-1}_B\}\)。

可以看到，\(\{o^1_B, o^2_B,...,o^{m-1}_B\}\)這些值被重復計算了，所以可以進行如下的優化。

假設在選第\(m\)個點時，對於點\(p_B\in B\)：

\(t^{m-1}_B\)表示點\(p_B\)到集合\(A\)中的\(m-1\)個點的距離的最小值：

\[t^{m-1}_B=min(\{o^1_B, o^2_B,...,o^{m-1}_B\}) \]

選出第\(m\)個點之后，對於點\(p_B\)，繼續選第\(m+1\)個點時，只需要計算點\(p_B\)到選出的第\(m\)個點的距離\(o^m_B\)，因為以下公式：

\[min(\{o^1_B, o^2_B, ..., o^m_B\}) = min(min\{o^1_B, o^2_B,...,o^{m-1}_B\},o^m_B) \]

上述公式還可以寫為\(t^m_B=min(t^{m-1}_B,o^m_B)\)，這是一個簡單的動態規划問題。在實現上，只需要一個長度為n的數組，分別保存每個點到集合\(A\)的距離，每當將一個新點加入集合\(A\)后，更新該數組即可。

串行實現代碼如下

/*
dataset: (b, n, 3)
temp: (b, n)
idxs: (b, m)
*/
void farthest_point_sampling_cpu(int b, int n, int m, const float *dataset, float *temp, int *idxs){
    const float * const dataset_start = dataset;
    float * const temp_start = temp;
    int * const idxs_start = idxs;
    for(int i=0; i<b; ++i){
        dataset = dataset_start + i*n*3;
        temp = temp_start + i*n;
        idxs = idxs_start + i*m;
        int old = 0;
        idxs[0] = old;
        for(int j=1; j<m; ++j){
            int besti = 0;
            float best = -1;
            float x1 = dataset[old * 3 + 0];
            float y1 = dataset[old * 3 + 1];
            float z1 = dataset[old * 3 + 2];
            for(int k=0; k<n; ++k){
                float x2, y2, z2;
                x2 = dataset[k*3+0];
                y2 = dataset[k*3+1];
                z2 = dataset[k*3+2];

                float d = (x2 - x1)*(x2 - x1)+(y2 - y1)*(y2 - y1)+(z2 - z1)*(z2 - z1);
                float d2 = min(d, temp[k]);
                temp[k] = d2;
                besti = d2 > best ? k : besti;
                best = d2 > best ? d2 : best;
            }
            old = besti;
            idxs[j] = old;
        }
    }
}
  • hints
    • temp數組即保存了\(t_B\)的值，該數組在函數調用前，需要將所有值初始化為INF/infinity
    • 該實現是串行的batch version的FPS實現
  • 時間復雜度分析
    • 每步選一個點，需要計算\(n\)個距離，總共選\(k\)個點，時間復雜度為\(\mathcal{O}(kn)\)
    • \(k\)和\(n\)一般是常數比例關系，所以可以認為是\(\mathcal{O}(n^2)\)
    • 考慮到batch size為\(b\), 最終的時間復雜度為\(\mathcal{O}(b*n^2)\)

3. FPS算法並行實現與分析

在實際使用中，面對大規模點雲數據的下采樣需求，串行實現的FPS算法在性能和效率上不能滿足要求，因此常常使用並行化技術加速FPS算法，這里給出一個CUDA上的FPS算法實現例子和分析。

__device__ void __update(float *__restrict__ dists, int *__restrict__ dists_i, int idx1, int idx2){
    const float v1 = dists[idx1], v2 = dists[idx2];
    const int i1 = dists_i[idx1], i2 = dists_i[idx2];
    dists[idx1] = max(v1, v2);
    dists_i[idx1] = v2 > v1 ? i2 : i1;
}

template <unsigned int block_size>__global__ void farthest_point_sampling_kernel(int b, int n, int m, 
    const float *__restrict__ dataset, float *__restrict__ temp, int *__restrict__ idxs) {
    // dataset: (B, N, 3)
    // temp: (B, N)
    // output:
    //      idxs : (B, M)

    if (m <= 0) return;
    __shared__ float dists[block_size];
    __shared__ int dists_i[block_size];

    int batch_index = blockIdx.x;
    dataset += batch_index * n * 3;
    temp += batch_index * n;
    idxs += batch_index * m;

    int tid = threadIdx.x;
    const int stride = block_size;

    int old = 0;
    if (threadIdx.x == 0)
        idxs[0] = old;

    __syncthreads();
    for (int j = 1; j < m; j++) {
        int besti = 0;
        float best = -1;
        float x1 = dataset[old * 3 + 0];
        float y1 = dataset[old * 3 + 1];
        float z1 = dataset[old * 3 + 2];
        for (int k = tid; k < n; k += stride) {
            float x2, y2, z2;
            x2 = dataset[k * 3 + 0];
            y2 = dataset[k * 3 + 1];
            z2 = dataset[k * 3 + 2];

            float d = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) + (z2 - z1) * (z2 - z1);
            float d2 = min(d, temp[k]);
            temp[k] = d2;
            besti = d2 > best ? k : besti;
            best = d2 > best ? d2 : best;
        }
        dists[tid] = best;
        dists_i[tid] = besti;
        __syncthreads();
        
        //Reduce programming pattern
        //Loop unrolling
        if (block_size >= 1024) {
            if (tid < 512) {
                __update(dists, dists_i, tid, tid + 512);
            }
            __syncthreads();
        }

        if (block_size >= 512) {
            if (tid < 256) {
                __update(dists, dists_i, tid, tid + 256);
            }
            __syncthreads();
        }
        if (block_size >= 256) {
            if (tid < 128) {
                __update(dists, dists_i, tid, tid + 128);
            }
            __syncthreads();
        }
        if (block_size >= 128) {
            if (tid < 64) {
                __update(dists, dists_i, tid, tid + 64);
            }
            __syncthreads();
        }
        if (block_size >= 64) {
            if (tid < 32) {
                __update(dists, dists_i, tid, tid + 32);
            }
            __syncthreads();
        }
        if (block_size >= 32) {
            if (tid < 16) {
                __update(dists, dists_i, tid, tid + 16);
            }
            __syncthreads();
        }
        if (block_size >= 16) {
            if (tid < 8) {
                __update(dists, dists_i, tid, tid + 8);
            }
            __syncthreads();
        }
        if (block_size >= 8) {
            if (tid < 4) {
                __update(dists, dists_i, tid, tid + 4);
            }
            __syncthreads();
        }
        if (block_size >= 4) {
            if (tid < 2) {
                __update(dists, dists_i, tid, tid + 2);
            }
            __syncthreads();
        }
        if (block_size >= 2) {
            if (tid < 1) {
                __update(dists, dists_i, tid, tid + 1);
            }
            __syncthreads();
        }

        old = dists_i[0];
        if (tid == 0) 
            idxs[j] = old;
    }
}

void farthest_point_sampling_kernel_launcher(int b, int n, int m,
    const float *dataset, float *temp, int *idxs) {
    // dataset: (B, N, 3)
    // temp: (B, N)
    // output:
    //      idx: (B, M)

    cudaError_t err;
    unsigned int n_threads = opt_n_threads(n);
    
    printf("n_threads:%d\n",n_threads);
    switch (n_threads) {
        case 1024:
        farthest_point_sampling_kernel<1024><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        case 512:
        farthest_point_sampling_kernel<512><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        case 256:
        farthest_point_sampling_kernel<256><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        case 128:
        farthest_point_sampling_kernel<128><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        case 64:
        farthest_point_sampling_kernel<64><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        case 32:
        farthest_point_sampling_kernel<32><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        case 16:
        farthest_point_sampling_kernel<16><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        case 8:
        farthest_point_sampling_kernel<8><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        case 4:
        farthest_point_sampling_kernel<4><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        case 2:
        farthest_point_sampling_kernel<2><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        case 1:
        farthest_point_sampling_kernel<1><<<b, n_threads>>>(b, n, m, dataset, temp, idxs); break;
        default:
        farthest_point_sampling_kernel<512><<<b, n_threads>>>(b, n, m, dataset, temp, idxs);
    }

    err = cudaGetLastError();
    if (cudaSuccess != err) {
        fprintf(stderr, "CUDA kernel failed : %s\n", cudaGetErrorString(err));
        exit(-1);
    }
}

CUDA入門材料可以參考[4]，此處不再贅述。此部分代碼分析可見參考[3]，寫的十分詳細，此處只強調幾個地方。

for (int k = tid; k < n; k += stride) {
    float x2, y2, z2;
    x2 = dataset[k * 3 + 0];
    y2 = dataset[k * 3 + 1];
    z2 = dataset[k * 3 + 2];

    float d = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) + (z2 - z1) * (z2 - z1);
    float d2 = min(d, temp[k]);
    temp[k] = d2;
    besti = d2 > best ? k : besti;
    best = d2 > best ? d2 : best;
}

這部分代碼是用來計算\(n\)個點到old點的距離的，這里主要分析線程塊內線程和數據的映射關系（假設\(blockDim.x=1024\)）：

  • 一個thread block對應的數據為dataset中的一個batch，點數為n，很顯然常常有\(n >> blockDim.x\)成立
  • 這里的循環就是用1024個線程處理\(n\)個點的pattern
    • 每個線程訪問[tid, tid+stride, ..., tid+i*stride,...]處對應的數據
    • 這里的線程warp在開始時無控制分歧出現，只在最后會出現一部分控制分歧
    • 同時，線程的數據訪問是接合的模式，即連續線程所訪問的數據可以由一次DRAM burst獲得，這樣可以提高DRAM的訪問效率
//Reduce programming pattern
//Loop unrolling
if (block_size >= 1024) {
    if (tid < 512) {
        __update(dists, dists_i, tid, tid + 512);
    }
    __syncthreads();
}

if (block_size >= 512) {
    if (tid < 256) {
        __update(dists, dists_i, tid, tid + 256);
    }
    __syncthreads();
}
if (block_size >= 256) {
    if (tid < 128) {
        __update(dists, dists_i, tid, tid + 128);
    }
    __syncthreads();
}
if (block_size >= 128) {
    if (tid < 64) {
        __update(dists, dists_i, tid, tid + 64);
    }
    __syncthreads();
}
if (block_size >= 64) {
    if (tid < 32) {
        __update(dists, dists_i, tid, tid + 32);
    }
    __syncthreads();
}
if (block_size >= 32) {
    if (tid < 16) {
        __update(dists, dists_i, tid, tid + 16);
    }
    __syncthreads();
}
if (block_size >= 16) {
    if (tid < 8) {
        __update(dists, dists_i, tid, tid + 8);
    }
    __syncthreads();
}
if (block_size >= 8) {
    if (tid < 4) {
        __update(dists, dists_i, tid, tid + 4);
    }
    __syncthreads();
}
if (block_size >= 4) {
    if (tid < 2) {
        __update(dists, dists_i, tid, tid + 2);
    }
    __syncthreads();
}
if (block_size >= 2) {
    if (tid < 1) {
        __update(dists, dists_i, tid, tid + 1);
    }
    __syncthreads();
}

上述代碼主要有兩點需要注意

  • 循環展開
    • 消除了for循環的邊界判斷分支指令以及loop計數器更新
    • 實現性能提升
  • 規約（Reduce）模式
    • 這里就是用來求dists和dists_i數組的最大值，是一個標准的max規約操作

4. 串行實現與並行實現的性能比較

在如下的硬件環境下進行了性能比較實驗

  • CPU: Intel i7-10710U
  • GPU: Nvidia MX350
  • CUDA version: 11.2
Batch Size #InputPts #SampledPts 串行實現用時 並行實現用時
1 10000 1000 76ms 8ms
2 10000 1000 157ms 8ms
3 10000 1000 234ms 8ms
4 10000 1000 311ms 15ms
5 10000 1000 381ms 20ms
6 10000 1000 459ms 26ms
6 10000 10000 4561ms 250ms

可以看到並行實現相比於串行實現在最困難的設置下加速比為4561ms/250ms=18.244，這表明並行實現相比於串行實現的效率和速度提升是十分明顯的。

參考資料

[1] Farthest Point sampling (FPs)算法核心思想解析

[2] OpenPCDet源碼中的fps-cuda-version實現

[3] Farthest Point Sampling in 3D Object Detection

[4] 大規模並行處理器編程實戰（第2版）


免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。



猜您在找 基於sinc的音頻重采樣（二）：實現 關於傅里葉分析與香農采樣定理 OpenCV兩種畸變校正模型源代碼分析以及CUDA實現 最遠曼哈頓距離 Python實現量子態采樣 pytorch 不使用轉置卷積來實現上采樣 傅立葉變換—FFT(cuda實現) 數組的並行求和-cuda實現 主成分分析(PCA)算法介紹及matlab實現案例 CUDA
 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM