手寫nms

計算寬高的時候加1是為什么？

本文總結自互聯網的多種nms實現，供參考，非博主原創，各原文鏈接如下，也建議大家動手寫一寫。

Ref：

淺談NMS的多種實現

目標窗口檢測算法-NMS非極大值抑制

一、faster-rcnn源碼閱讀：nms的CUDA編程

c++版 nms

nms簡介

首先還是要科普一下nms算法的思想：簡單來說就是去重框。這里的重框針對的當然是某一類的框。下面實現的時候也是默認拿到某一類所有的框。

算法思路：

For a prediction bounding box B, the model calculates the predicted probability for each category. Assume the largest predicted probability is p, the category corresponding to this probability is the predicted category of B. We also refer to pas the confidence level of prediction bounding box B. On the same image, we sort the prediction bounding boxes with predicted categories other than background by confidence level from high to low, and obtain the list L. Select the prediction bounding box B1 with highest confidence level from L as a baseline and remove all non-benchmark prediction bounding boxes with an IoU with B1 greater than a certain threshold from L. The threshold here is a preset hyper-parameter. At this point,L retains the prediction bounding box with the highest confidence level and removes other prediction bounding boxes similar to it. Next, select the prediction bounding box B2 with the second highest confidence level from L as a baseline, and remove all non-benchmark prediction bounding boxes with an IoU with B2 greater than a certain threshold from L. Repeat this process until all prediction bounding boxes in L have been used as a baseline. At this time, the IoU of any pair of prediction bounding boxes in L is less than the threshold. Finally, output all prediction bounding boxes in the list L.

前面這段話基本說出了算法的具體實現思路：先對每個框的score進行排序，首先選擇第一個，也就是score最高的框，它一定是我們要保留的框。然后拿它和剩下的框進行比較，如果IOU大于一定閾值，說明兩者重合度高，應該去掉，這樣篩選出的框就是和第一個框重合度低的框，第一次迭代結束。第二次從保留的框中選出score第一的框，重復上述過程直到沒有框保留了。

Python

版本一

該版本為 Faster RCNN 實現的版本，是網絡是最常見的版本：

def nms(dets, thresh):x1 = dets[:, 0] #xminy1 = dets[:, 1] #yminx2 = dets[:, 2] #xmaxy2 = dets[:, 3] #ymaxscores = dets[:, 4] #confidenceareas = (x2 - x1 + 1) * (y2 - y1 + 1) # 每個boundingbox的面積order = scores.argsort()[::-1] # boundingbox的置信度排序keep = [] # 用來保存最后留下來的boundingboxwhile order.size > 0:     i = order[0] # 置信度最高的boundingbox的indexkeep.append(i) # 添加本次置信度最高的boundingbox的index# 當前bbox和剩下bbox之間的交叉區域# 選擇大于x1,y1和小于x2,y2的區域xx1 = np.maximum(x1[i], x1[order[1:]]) #交叉區域的左上角的橫坐標yy1 = np.maximum(y1[i], y1[order[1:]]) #交叉區域的左上角的縱坐標xx2 = np.minimum(x2[i], x2[order[1:]]) #交叉區域右下角的橫坐標yy2 = np.minimum(y2[i], y2[order[1:]]) #交叉區域右下角的縱坐標# 當前bbox和其他剩下bbox之間交叉區域的面積w = np.maximum(0.0, xx2 - xx1 + 1)h = np.maximum(0.0, yy2 - yy1 + 1)inter = w * h# 交叉區域面積 / (bbox + 某區域面積 - 交叉區域面積)ovr = inter / (areas[i] + areas[order[1:]] - inter)#保留交集小于一定閾值的boundingboxinds = np.where(ovr <= thresh)[0]order = order[inds + 1]return keep

可以看到，基本是按照上述思路去寫的，按照score進行降序排序，然后每次拿到第一個框，也就是score最大的框，然后計算該框與其他框的IOU，最后留下iou<=thresh的框留作下次循環，這里唯一值得強調的是最后這個索引為什么要+1。這是因為我們要得到的inds是排除了當前用來比較的score最大的框，所以在其原始索引基礎上+1 ,從代碼中看就是由于order[1:]這樣寫導致的。

當然，這個大眾版本思路很清楚，但感覺不夠優雅，可以再精簡一點。

版本二

def nms(dets, thresh):                   areas=np.prod(bbox[:,2:]-bbox[:,:2],axis=1)order = scores.argsort()[::-1]keep=[]while order.size>0:i=order[0]keep.append(i)tl=np.maximum(b[:2],bbox[i+1:,:2])br=np.minimum(b[2:],bbox[i+1:,2:])inter=np.prod(br-tl,axis=1)*(br>tl).all(axis=1)ovr=inter/(areas[order[1:]]+areas[i]-inter)inds=np.where(ovr<=thresh)[0]order=order[inds+1]return keep

當然這里的思路還是要排序，只不過iou部分寫的更精簡了。

好了，鋪墊了這么久，說一下不排序怎么寫。基本思路是，依次遍歷每個框，計算這個框與其他框的iou,找到iou大于一定閾值的其他框，因為這個時候不能保證它一定是score最高的框，所以要進行判斷，如果它的score小于其他框，那就把它去掉，因為它肯定不是要保留的框。如果它的score大于其他框，那應該保留它，同時可以去掉所有其他框了。最后保留的框就是結果。

def nms(bbox, scores, thresh):area=np.prod(bbox[:,2:]-bbox[:,:2],axis=1)keep=np.ones(len(bbox),dtype=bool)for i, b in enumerate(bbox):if(keep[i]==False):continuetl=np.maximum(b[:2],bbox[i+1:,:2])br=np.minimum(b[2:],bbox[i+1:,2:])inter=np.prod(br-tl,axis=1)*(br>=tl).all(axis=1)iou=ia/(area[i+1:]+area[i]-inter)r = [ k for k in np.where(iou>thresh)[0]+i+1 if keep[k]==True]if (scores[i]>scores[r]).all():keep[r]=Falseelse:keep[i]=Falsereturn np.where(keep)[0].astype(np.int32)

這是我按照上面思路寫的，用keep表示框的去留，為True的框要保留。當然這個思路的效率不見得比前面排序的要好，只是作為其他角度思考。

版本三

補充：之前排序實現是通過每次篩除框來完成的，直到最后沒有框剩下了則循環結束。但還可以這么想：同樣先將框按score排好序，這次循環條件為遍歷所有框，在某次循環，拿到一個框，將其與所有已經保留的框進行比較，如果iou大于閾值，說明它應該被刪去，直接進行下一次循環，如果小于將其加入進要保留的框中。當遍歷完所有框時結束，拿到保留的框。

嗯，文字還是有點繞，直接看code更清楚：

def _non_maximum_suppression_cpu(bbox, thresh, score):order=np.argsort(score)[::-1]bbox=bbox[order]      bbox_area=np.prod(bbox[:, 2:]-bbox[:, :2],axis=1)keep=np.zeros(bbox.shape[0], dtype=bool)  for i, b in enumerate(bbox):tl=np.maximum(b[:2],bbox[keep, :2])  br=np.minimum(b[2:],bbox[keep, 2:])	   area = np.prod(br-tl, axis=1)*(tl<br).all(axis=1)iou=area/(bbox_area[i]+bbox_area[keep]-area) if (iou>=thresh).any():   continuekeep[i]=Truekeep = np.where(keep)[0]return keep.astype(np.int32)

至此介紹了三種NMS算法的實現思路。即使作為常規的算法題考察也很好，因為這也不牽扯到目標檢測的東西。而且還順便考察了numpy的操作。這也提醒自己有時候要多深入想想，而不是簡單copy網上現有的東西。

C++

typedef struct Bbox{int x;int y;int w;int h;float score;
}Bbox;static bool sort_score(Bbox box1,Bbox box2){return box1.score > box2.score ? true : false;
}float iou(Bbox box1,Bbox box2){int x1 = max(box1.x,box2.x);int y1 = max(box1.y,box2.y);int x2 = min(box1.x+box1.w,box2.x+box2.w);int y2 = min(box1.y+box1.h,box2.y+box2.h);// int w = max(0,x2 - x1 + 1);// int h = max(0,y2 - y1 + 1);int w = max(0,x2 - x1 + 1);int h = max(0,y2 - y1 + 1);// float over_area = w*h;float over_area = (x2 - x1) * (y2 - y1);return over_area / (box1.w * box1.h + box2.w * box2.h - over_area);
}vector<Bbox> nms(std::vector<Bbox>&vec_boxs,float threshold){vector<Bbox>results;std::sort(vec_boxs.begin(),vec_boxs.end(),sort_score);while(vec_boxs.size() > 0) {results.push_back(vec_boxs[0]);int index = 1;while(index < vec_boxs.size()){float iou_value = iou(vec_boxs[0],vec_boxs[index]);cout << "iou:" << iou_value << endl;if(iou_value > threshold)vec_boxs.erase(vec_boxs.begin() + index);elseindex++;}vec_boxs.erase(vec_boxs.begin());}return results;
}

測試：

int main(){vector<Bbox> input;Bbox box1 = {1, 2, 2, 2, 0.4};Bbox box2 = {1, 3, 2, 1, 0.5};Bbox box3 = {1, 3, 3, 1, 0.72};Bbox box4 = {1, 1, 3, 3, 0.9};Bbox box5 = {1, 1, 2, 3, 0.45};input.push_back(box1);input.push_back(box2);input.push_back(box3);input.push_back(box4);input.push_back(box5);vector<Bbox> res;res = nms(input, 0.5);for(int i = 0;i < res.size();i++){printf("%d %d %d %d %f",res[i].x,res[i].y,res[i].w,res[i].h,res[i].score);cout << endl;}return 0;
}

CUDA

cpu版驗證和理解了算法，下面來看看GPU實現加速，速度大概可以提升50X。在2080it上大概是115ms比3ms。

nms各種實現的benchmark請看：

https://github.com/fmscole/benchmark

主控函數分兩部分，第一部分計算mask，還需要的標0，不需要的（重復度大的）標1，第二部分是根據mask，選出留下來的候選框。

先簡要說一下CUDA編程模型：

GPU之所以能夠加速，是因為并行計算，即每個線程負責計算一個數據，充分利用GPU計算核心超多（幾千個）的優勢。

（1）每個計算核心相互獨立，運行同一段代碼，這段代碼稱為核函數；

（2）每個核心有自己的身份id，線程的身份id是兩個三維數組：（blockIdx.x，blockIdx.y，blockIdx.z）-（threadIdx.x，threadIdx.y，threadIdx.z）。

身份id被另兩個三維數組grid（gridDim.x,gridDim.y,gridDim.z）和block(blockDim.x,blockDim.y,blockDim.z)確定范圍。

總共有 $gridDim.x\times gridDim.y\times gridDim.z$ 個 block

每個block有 $blockDim.x\times blockDim.y\times blockDim.z$ 個 thread。

有了線程的身份id，經過恰當的安排，讓身份id（核函數可以獲取）：（blockIdx.x，blockIdx.y，blockIdx.z）-（threadIdx.x，threadIdx.y，threadIdx.z）對應到一個數據，就可以實現一個線程計算一個數據，至于如何對應，開發人員得好好安排，可以 說這是CUDA開發的一個核心問題。

gridDim.x、blockIdx.x這些是核函數可以獲取的，gridDim.x等于多少，調用核函數的時候就要定一下來。看代碼：

  dim3 blocks(DIVUP(boxes_num, threadsPerBlock),DIVUP(boxes_num, threadsPerBlock));dim3 threads(threadsPerBlock);nms_kernel<<<blocks, threads>>>(boxes_num,nms_overlap_thresh,boxes_dev,mask_dev);

這里的threadsPerBlock=8*8=64，

當boxes_num=12030時，DIVUP(12030, 64)=12030/64+12030%64>0=188

在調用核函數的時候，通過<<<blocks, threads>>>（#這是cu語法，不是標準C語言）把線程數量安排傳遞進去，核函數里就有

gridDim.x=188,gridDim.y=188,gridDim.z=1；

blockDim.x=64,blockDim.y=1,blockDim.z=1；

0<=blockIdx.x<188,

0<=blockIdx.y<188,

blockIdx.z=0,

0<=threadIdx.x<64，

threadIdx.y=threadIdx.z=0，

這樣就啟動了2,262,016個（兩百多萬個線程）來計算，兩百多萬看起來嚇人，對GPU來書毫無負擔！每個線程計算不超過64個值，后面再講。

（3）這里的grid(a,b,c），block(x,y,z)值是多少，由程序設計人員根據問題來定，在調用核函數時就要確定下來，但有一個基本限制block(x,y,z)中的x×y×z<=1024（這個值隨GPU版本確定，起碼nvidia 1080，2080都是這樣）；

（4）block中的線程每32個thread為一束，絕對同步：比如if-else語句，這32個線程中有的滿足if條件，有的滿足else。滿足else的那部分線程不能直接進入，而是要等滿足if的那部分線程運行完畢才進入else部分，而滿足if的那部分線程現在也不能結束，而是要等else部分線程運行完畢，大家才能同時結束。for語句也是一樣。因此GPU計算盡可能不要有分支語句。

不是說不能用if和for，該用還得用，用的時候要知道付出的代價。否則實現了減速都不知道為了啥。

不同的線程束之間不同步，如果同步需要請__syncthreads();

如果設置block(1)，即一個block只安排一個線程呢？事實上GPU還是要啟動32個線程，另外31個陪跑。

因此block(x,y,z)中的x×y×z應該為32的倍數。不過32×32=1024了。

（5）要并行計算，前提是數據之間沒有相互依賴，有前后依賴的部分只能放在同一個核函數里計算；

先看控制部分代碼，這部分做的事情就是：

1. 在GPU上分配內存，把數據傳到GPU
2. 調用核函數，計算mask；
3. 把數據傳回來，
4. 根據mask把獲取保留下來的候選框。

nms_kernel.cu（來源https://github.com/jwyang/faster-rcnn.pytorch）：

void _nms(int* keep_out, int* num_out, const float* boxes_host, int boxes_num,int boxes_dim, float nms_overlap_thresh, int device_id) {_set_device(device_id);
//keep_out：返回保留目標框的下標
//num_out：返回保留下來的個數float* boxes_dev = NULL;unsigned long long* mask_dev = NULL;const int col_blocks = DIVUP(boxes_num, threadsPerBlock);
//比如boxes_num=12030時，col_blocks=188，后面都以此為例CUDA_CHECK(cudaMalloc(&boxes_dev,boxes_num * boxes_dim * sizeof(float)));
//在GPU上分配內存CUDA_CHECK(cudaMemcpy(boxes_dev,boxes_host,boxes_num * boxes_dim * sizeof(float),cudaMemcpyHostToDevice));
//把候選框數據復制到GPUCUDA_CHECK(cudaMalloc(&mask_dev,boxes_num * col_blocks * sizeof(unsigned long long)));
//分配mask的內存，用于返回mask的計算結果dim3 blocks(DIVUP(boxes_num, threadsPerBlock),DIVUP(boxes_num, threadsPerBlock));
//(188,188)dim3 threads(threadsPerBlock);
//(64)//調用核函數nms_kernel<<<blocks, threads>>>(boxes_num,nms_overlap_thresh,boxes_dev,mask_dev);std::vector<unsigned long long> mask_host(boxes_num * col_blocks);
//在CPU上分配內存，用于返回mask就算結果CUDA_CHECK(cudaMemcpy(&mask_host[0],mask_dev,sizeof(unsigned long long) * boxes_num * col_blocks,cudaMemcpyDeviceToHost));
//把mask從GPU復制到CPU
//下面這段代碼最后再敘述std::vector<unsigned long long> remv(col_blocks);memset(&remv[0], 0, sizeof(unsigned long long) * col_blocks);int num_to_keep = 0;for (int i = 0; i < boxes_num; i++) {int nblock = i / threadsPerBlock;int inblock = i % threadsPerBlock;if (!(remv[nblock] & (1ULL << inblock))) {keep_out[num_to_keep++] = i;unsigned long long *p = &mask_host[0] + i * col_blocks;for (int j = nblock; j < col_blocks; j++) {remv[j] |= p[j];}}}*num_out = num_to_keep;CUDA_CHECK(cudaFree(boxes_dev));CUDA_CHECK(cudaFree(mask_dev));
//釋放GPU上的內存
}

計算mask的核函數是核心代碼，什么是mask呢？從內存里講是一段連續內存，但應該把它想象成一個矩陣：

比如候選框個數為boxes_num=12030時，由于要計算任意兩個候選框之間的IOU是否大于閾值，因此，我們需要建一個12030*12030的矩陣，1表示兩個候選框的IOU大于0.7，0表示不大于。這樣任意兩個候選框之間的關系都可以用這個12030行*12030列（行表示本框，列表示其他框）的矩陣保存。

但是，這里有好幾處值得改進的地方：

1）對角線上的值不需要計算，因為意味著自己與自己計算IOU，沒意義；

2）上三角與下三角是對稱的，只需要用到上三角即可；

3）更重要的是，為了保存0或1的值，真需要建這么大的矩陣嗎？事實上每連續的64個0、1剛好構成一個無符號的64位整數（unsigned long long），我們只用一個整數表示即可，這樣內存減少64倍！這中間設計到位運算，這不是問題。原本需要一個12030***12030的一個數組來表示這個矩陣，現在只要12030***188的一個unsigned long long數組就可以映射這個12030行，188列的矩陣。

比如：這長度為1230*188數組中的第一組188個整數（相當于矩陣的第一行）記錄的是第一個框（按分數排序之后）與其他所有框之間的重疊關系，第一個整數的第一位表示與自己，由于自己與自己不參與計算，所以這個值一定為0.

__global__ void nms_kernel(int n_boxes, float nms_overlap_thresh,float *dev_boxes, unsigned long long *dev_mask) {//n_boxes候選框的個數，比如是12030//nms_overlap_thresh=0.7閾值//dev_boxes候選框//dev_mask存放mask的值const int row_start = blockIdx.y;//把blockIdx.y想象成矩陣的行號，n_boxes=12030時，共有188行，const int col_start = blockIdx.x;//把blockIdx.x想象成矩陣的列號，n_boxes=12030時，共有188列，而每個block有64個線程，188*64=12032// if (row_start > col_start) return;const int row_size =min(n_boxes - row_start * threadsPerBlock, threadsPerBlock);const int col_size =min(n_boxes - col_start * threadsPerBlock, threadsPerBlock);//在block里，row_size和col_size最多取到64，尾部不足64就取余數。先把數據復制一份到共享內存，關于共享內存，后面詳述__shared__ float block_boxes[threadsPerBlock * 5];//每個block有64個線程，所以復制64個候選框，每個候選框有4個坐標值和一個分數值，共5個值，//所以每個block分配的共享內存大小為64*5=320if (threadIdx.x < col_size) {block_boxes[threadIdx.x * 5 + 0] =dev_boxes[(threadsPerBlock * col_start + threadIdx.x) * 5 + 0];block_boxes[threadIdx.x * 5 + 1] =dev_boxes[(threadsPerBlock * col_start + threadIdx.x) * 5 + 1];block_boxes[threadIdx.x * 5 + 2] =dev_boxes[(threadsPerBlock * col_start + threadIdx.x) * 5 + 2];block_boxes[threadIdx.x * 5 + 3] =dev_boxes[(threadsPerBlock * col_start + threadIdx.x) * 5 + 3];block_boxes[threadIdx.x * 5 + 4] =dev_boxes[(threadsPerBlock * col_start + threadIdx.x) * 5 + 4];}__syncthreads();//同步if (threadIdx.x < row_size) {//每個線程雖然運行的是同一段代碼，但看到的表示身份的threadIdx.x是不一樣的，row_size值也不一樣//具體的說，有188*188個線程看到的threadIdx.x是一樣的，因為總共有188*188個block，threadIdx.x的取值范圍是0到63const int cur_box_idx = threadsPerBlock * row_start + threadIdx.x;//int cur_box_idx這個才是真正的行號，因為grid下面還有block。//CUDA編程就是要把這些表示線程的id對應到具體數據，blockIdx.y和threadIdx.x確定當前候選框const float *cur_box = dev_boxes + cur_box_idx * 5;//取出當前的候選框，別忘了候選框用5個數值表示的，所以要乘5int i = 0;unsigned long long t = 0;//t就是存放64個0、1的整數int start = 0;if (row_start == col_start) {//對角線上的blockstart = threadIdx.x + 1;//自己跟自己就不要計算IOU了}for (i = start; i < col_size; i++) {if (devIoU(cur_box, block_boxes + i * 5) > nms_overlap_thresh) {//每一個當前框都與其他框計算IOU，（其他框存放在共享內存，是一個復制品）//本線程只負責計算第blockIdx.y×64+threadIdx.x號候選框（當前框）//與blockIdx.x×64～blockIdx.x×64+63這64個候選框（其他框）之間的關系t |= 1ULL << i;}}const int col_blocks = DIVUP(n_boxes, threadsPerBlock);//同上假定下，col_blocks=188dev_mask[cur_box_idx * col_blocks + col_start] = t;//dev_mask總長為12030*188，表示12030行，188列的矩陣//即分為12030段（行），每段188個（列）int64，每個int64標記64個0、1，所以每行的這188個int64可以標記12032個0\1，//就是說第i段（行）的12032個記錄的是第i個候選框與其他所有候選框之間的重疊關系。//cur_box_idx表示第幾段，每段col_blocks（=188）個數，col_start是本段里第幾個數，所以//dev_mask[cur_box_idx * col_blocks + col_start] =t 記錄的是第cur_box_idx個候選框//與第col_start×64到（col_start+1）×64-1之間這64個矩形框之間的重疊關系}
}

再總結一下，總共規劃了（188,188,1）個block，每個block有（64,1,1）個線程，第（c，r，0）-（t,0,0)號線程負責計算的數據是：r×64+t號候選框與c×64～c×64+63號這64個候選框之間的重疊關系，并存進（r×64+t）×188+c這個整數里。當c，r遍歷完188，t遍歷完64，這188×188×64個線程就計算完了任意兩個候選框之間的關系。

由于cuda里每個線程計算一個數據，相當于cpu里的循環。CPU里的循環是一個一個的算，而cuda里是同時在算。

本質上，這里gird（188,188,1）和block（64,1,1）來代替了三重循環：

for（int i=0；i<188;i++)for（int j=0；j<188;j++)for（int k=0；k<64;k++)。。。。。。

只要把這循環中的i，j，k替換為blockIdx.x，blockIdx.y，threadIdx.x即可。

CUDA編程就是循環替代！——這是我目前的理解。

回過頭來看把候選框復制進共享內存這部分，

block_boxes[threadIdx.x * 5 + 0] =dev_boxes[(threadsPerBlock * col_start + threadIdx.x) * 5 + 0];

可以看出這里面只用到了blockIdx.x和threadIdx.x，注意到了沒，沒用到blockIdx.y，而blockIdx.y的范圍是0~187，什么意思，意味著這段代碼被重復執行了188次！畢竟，我們啟動了188×64×188個線程，而數據只有12030個數據，因此有188個線程執行的是相同的數據！不知道我理解的對不對，望大佬幫我指正。

計算好了mask數組之后，計算保留下來的目標框：

//...................std::vector<unsigned long long> remv(col_blocks);memset(&remv[0], 0, sizeof(unsigned long long) * col_blocks);
//長為188的unsigned long long數組，初始值為0
//這12032個標記位用于記錄那些候選框已經被剔除int num_to_keep = 0;
//保留下來的個數for (int i = 0; i < boxes_num; i++) {int nblock = i / threadsPerBlock;int inblock = i % threadsPerBlock;
//把12030個候選框分成188組，每64個一組，第nblock組的第inblock個候選框就是第i個候選框if (!(remv[nblock] & (1ULL << inblock))) {
//還沒有被舍棄，這兩個條件只要有一個不成立，都意味著沒有被剔除keep_out[num_to_keep++] = i;
//把下標記錄下來unsigned long long *p = &mask_host[0] + i * col_blocks;
//再回憶一遍，mask_host總長是12030×188
//分為12030段，第i段有188個unsigned long long整數，記錄第i個候選框與其他所有候選框之間的信息
//p指向的是第i段的起始位置for (int j = nblock; j < col_blocks; j++) {
//遍歷后面的188個整數，總共12032個標記位，多出來的2位不去管它remv[j] |= p[j];
//因為第i個候選框已經保留下來，所以與第i個候選框重復的候選框都標記為去除}}}
//keep_out的前num_to_keep個數就是保留下來的候選框的下標
//.................................