使用opencv作物件識別（一） —— 積分直方圖加速HOG特徵計算

方向梯度直方圖（Histograms of Oriented Gradients，簡稱HOG特徵）結合支援向量機（ support vector machine， 簡稱SVM），被廣泛應用於影象識別中，尤其在行人檢測中獲得了極大的成功。

積分直方圖可以用於快速計算原始影象矩形區域內的HOG特徵。積分直方圖的概念類似與viola和jones在臉部識別中所用的積分影象。

下面的程式碼給出了，對於一幅給定的影象，如何快速計算積分直方圖，以及如何使用其進行HOG特徵的演算（關鍵處以給出註釋）：

C程式碼

 

1. /*Function to calculate the integral histogram*/  
2. IplImage** calculateIntegralHOG(IplImage* in)   
3. {   
4.     /*Convert the input image to grayscale*/
     
5.     IplImage* img_gray = cvCreateImage(cvGetSize(in), IPL_DEPTH_8U,1);   
6.     cvCvtColor(in, img_gray, CV_BGR2GRAY);   
7.     cvEqualizeHist(img_gray,img_gray);   
8.   
9. /*Calculate the derivates of the grayscale image in the x and y directions using a sobel operator and obtain 2 gradient images for the x and y directions*/  
10.   
11.     IplImage *xsobel, *ysobel;   
12.     xsobel = doSobel(img_gray, 1, 0, 3);   
13.     ysobel = doSobel(img_gray, 0, 1, 3);   
14.     cvReleaseImage(&img_gray);   
15.   
16.   
17. /* Create an array of 9 images (9 because I assume bin size 20 degrees and unsigned gradient ( 180/20 = 9), one for each bin which will have zeroes for all pixels, except for the pixels in the original image for which the gradient values correspond to the particular bin. These will be referred to as bin images. These bin images will be then used to calculate the integral histogram, which will quicken the calculation of HOG descriptors */  
18.   
19.     IplImage** bins = (IplImage**) malloc(9 * sizeof(IplImage*));   
20.     for (int i = 0; i < 9 ; i++) {   
21.         bins[i] = cvCreateImage(cvGetSize(in), IPL_DEPTH_32F,1);   
22.         cvSetZero(bins);   
23.     }   
24.   
25.   
26. /* Create an array of 9 images ( note the dimensions of the image, the cvIntegral() function requires the size to be that), to store the integral images calculated from the above bin images. These 9 integral images together constitute the integral histogram */  
27.   
28.     IplImage** integrals = (IplImage**) malloc(9 * sizeof(IplImage*));    
29.     for (int i = 0; i < 9 ; i++) {   
30.         integrals[i] = cvCreateImage(cvSize(in->width + 1, in->height + 1),   
31.         IPL_DEPTH_64F,1);   
32.     }   
33.   
34. /* Calculate the bin images. The magnitude and orientation of the gradient at each pixel is calculated using the xsobel and ysobel images.{Magnitude = sqrt(sq(xsobel) + sq(ysobel) ), gradient = itan (ysobel/xsobel) }. Then according to the orientation of the gradient, the value of the corresponding pixel in the corresponding image is set */  
35.   
36.     int x, y;   
37.     float temp_gradient, temp_magnitude;   
38.     for (y = 0; y < in->height; y++) {   
39.   
40. /* ptr1 and ptr2 point to beginning of the current row in the xsobel and ysobel images respectively. ptrs point to the beginning of the current rows in the bin images */  
41.   
42.         float* ptr1 = (float*) (xsobel->imageData + y * (xsobel->widthStep));   
43.         float* ptr2 = (float*) (ysobel->imageData + y * (ysobel->widthStep));   
44.         float** ptrs = (float**) malloc(9 * sizeof(float*));   
45.         for (int i = 0; i < 9 ;i++){   
46.             ptrs[i] = (float*) (bins[i]->imageData + y * (bins->widthStep));   
47.         }   
48.   
49. /*For every pixel in a row gradient orientation and magnitude are calculated and corresponding values set for the bin images. */  
50.   
51.         for (x = 0; x <in->width; x++) {   
52.   
53. /* if the xsobel derivative is zero for a pixel, a small value is added to it, to avoid division by zero. atan returns values in radians, which on being converted to degrees, correspond to values between -90 and 90 degrees. 90 is added to each orientation, to shift the orientation values range from {-90-90} to {0-180}. This is just a matter of convention. {-90-90} values can also be used for the calculation. */  
54.   
55.             if (ptr1[x] == 0){   
56.                 temp_gradient = ((atan(ptr2[x] / (ptr1[x] + 0.00001))) * (180/   PI)) + 90;   
57.             }   
58.             else{   
59.                 temp_gradient = ((atan(ptr2[x] / ptr1[x])) * (180 / PI)) + 90;   
60.             }   
61.             temp_magnitude = sqrt((ptr1[x] * ptr1[x]) + (ptr2[x] * ptr2[x]));   
62.   
63. /*The bin image is selected according to the gradient values. The corresponding pixel value is made equal to the gradient magnitude at that pixel in the corresponding bin image */  
64.   
65.             if (temp_gradient <= 20) {   
66.                 ptrs[0][x] = temp_magnitude;   
67.             }   
68.             else if (temp_gradient <= 40) {   
69.                 ptrs[1][x] = temp_magnitude;   
70.             }   
71.             else if (temp_gradient <= 60) {   
72.                 ptrs[2][x] = temp_magnitude;   
73.             }   
74.             else if (temp_gradient <= 80) {   
75.                 ptrs[3][x] = temp_magnitude;   
76.             }   
77.             else if (temp_gradient <= 100) {   
78.                 ptrs[4][x] = temp_magnitude;   
79.             }   
80.             else if (temp_gradient <= 120) {   
81.                 ptrs[5][x] = temp_magnitude;   
82.             }   
83.             else if (temp_gradient <= 140) {   
84.                 ptrs[6][x] = temp_magnitude;   
85.             }   
86.             else if (temp_gradient <= 160) {   
87.                 ptrs[7][x] = temp_magnitude;   
88.             }   
89.             else {   
90.                 ptrs[8][x] = temp_magnitude;   
91.             }   
92.         }   
93.     }   
94.   
95.     cvReleaseImage(&xsobel);   
96.     cvReleaseImage(&ysobel);   
97.   
98. /*Integral images for each of the bin images are calculated*/  
99.   
100.     for (int i = 0; i <9 ; i++){   
101.         cvIntegral(bins[i], integrals[i]);   
102.     }   
103.   
104.     for (int i = 0; i <9 ; i++){   
105.         cvReleaseImage(&bins[i]);   
106.     }   
107.   
108. /*The function returns an array of 9 images which consitute the integral histogram*/  
109.   
110.     return (integrals);   
111.   
112. }  

/*Function to calculate the integral histogram*/
IplImage** calculateIntegralHOG(IplImage* in)
{
    /*Convert the input image to grayscale*/
    IplImage* img_gray = cvCreateImage(cvGetSize(in), IPL_DEPTH_8U,1);
    cvCvtColor(in, img_gray, CV_BGR2GRAY);
    cvEqualizeHist(img_gray,img_gray);

/*Calculate the derivates of the grayscale image in the x and y directions using a sobel operator and obtain 2 gradient images for the x and y directions*/

    IplImage *xsobel, *ysobel;
    xsobel = doSobel(img_gray, 1, 0, 3);
    ysobel = doSobel(img_gray, 0, 1, 3);
    cvReleaseImage(&img_gray);


/* Create an array of 9 images (9 because I assume bin size 20 degrees and unsigned gradient ( 180/20 = 9), one for each bin which will have zeroes for all pixels, except for the pixels in the original image for which the gradient values correspond to the particular bin. These will be referred to as bin images. These bin images will be then used to calculate the integral histogram, which will quicken the calculation of HOG descriptors */

    IplImage** bins = (IplImage**) malloc(9 * sizeof(IplImage*));
    for (int i = 0; i < 9 ; i++) {
        bins[i] = cvCreateImage(cvGetSize(in), IPL_DEPTH_32F,1);
        cvSetZero(bins);
    }


/* Create an array of 9 images ( note the dimensions of the image, the cvIntegral() function requires the size to be that), to store the integral images calculated from the above bin images. These 9 integral images together constitute the integral histogram */

    IplImage** integrals = (IplImage**) malloc(9 * sizeof(IplImage*)); 
    for (int i = 0; i < 9 ; i++) {
        integrals[i] = cvCreateImage(cvSize(in->width + 1, in->height + 1),
        IPL_DEPTH_64F,1);
    }

/* Calculate the bin images. The magnitude and orientation of the gradient at each pixel is calculated using the xsobel and ysobel images.{Magnitude = sqrt(sq(xsobel) + sq(ysobel) ), gradient = itan (ysobel/xsobel) }. Then according to the orientation of the gradient, the value of the corresponding pixel in the corresponding image is set */

    int x, y;
    float temp_gradient, temp_magnitude;
    for (y = 0; y < in->height; y++) {

/* ptr1 and ptr2 point to beginning of the current row in the xsobel and ysobel images respectively. ptrs point to the beginning of the current rows in the bin images */

        float* ptr1 = (float*) (xsobel->imageData + y * (xsobel->widthStep));
        float* ptr2 = (float*) (ysobel->imageData + y * (ysobel->widthStep));
        float** ptrs = (float**) malloc(9 * sizeof(float*));
        for (int i = 0; i < 9 ;i++){
            ptrs[i] = (float*) (bins[i]->imageData + y * (bins->widthStep));
        }

/*For every pixel in a row gradient orientation and magnitude are calculated and corresponding values set for the bin images. */

        for (x = 0; x <in->width; x++) {

/* if the xsobel derivative is zero for a pixel, a small value is added to it, to avoid division by zero. atan returns values in radians, which on being converted to degrees, correspond to values between -90 and 90 degrees. 90 is added to each orientation, to shift the orientation values range from {-90-90} to {0-180}. This is just a matter of convention. {-90-90} values can also be used for the calculation. */

            if (ptr1[x] == 0){
                temp_gradient = ((atan(ptr2[x] / (ptr1[x] + 0.00001))) * (180/   PI)) + 90;
            }
            else{
                temp_gradient = ((atan(ptr2[x] / ptr1[x])) * (180 / PI)) + 90;
            }
            temp_magnitude = sqrt((ptr1[x] * ptr1[x]) + (ptr2[x] * ptr2[x]));

/*The bin image is selected according to the gradient values. The corresponding pixel value is made equal to the gradient magnitude at that pixel in the corresponding bin image */

            if (temp_gradient <= 20) {
                ptrs[0][x] = temp_magnitude;
            }
            else if (temp_gradient <= 40) {
                ptrs[1][x] = temp_magnitude;
            }
            else if (temp_gradient <= 60) {
                ptrs[2][x] = temp_magnitude;
            }
            else if (temp_gradient <= 80) {
                ptrs[3][x] = temp_magnitude;
            }
            else if (temp_gradient <= 100) {
                ptrs[4][x] = temp_magnitude;
            }
            else if (temp_gradient <= 120) {
                ptrs[5][x] = temp_magnitude;
            }
            else if (temp_gradient <= 140) {
                ptrs[6][x] = temp_magnitude;
            }
            else if (temp_gradient <= 160) {
                ptrs[7][x] = temp_magnitude;
            }
            else {
                ptrs[8][x] = temp_magnitude;
            }
        }
    }

    cvReleaseImage(&xsobel);
    cvReleaseImage(&ysobel);

/*Integral images for each of the bin images are calculated*/

    for (int i = 0; i <9 ; i++){
        cvIntegral(bins[i], integrals[i]);
    }

    for (int i = 0; i <9 ; i++){
        cvReleaseImage(&bins[i]);
    }

/*The function returns an array of 9 images which consitute the integral histogram*/

    return (integrals);

}

如何使用上面的函式來計算9維的方向梯度直方圖呢？如下：

C程式碼  

1. /* The following function takes as input the rectangular cell for which the histogram of oriented gradients has to be calculated, a matrix hog_cell of dimensions 1x9 to store the bin values for the histogram, the integral histogram, and the normalization scheme to be used. No normalization is done if normalization = -1 */  
2.   
3. void calculateHOG_rect(CvRect cell, CvMat* hog_cell,   
4. IplImage** integrals, int normalization) {   
5.   
6. /* Calculate the bin values for each of the bin of the histogram one by one */  
7.   
8.     for (int i = 0; i < 9 ; i++){   
9.   
10.         float a =((double*)(integrals[i]->imageData + (cell.y) * (integrals->   
11.             widthStep)))[cell.x];   
12.   
13.         float b = ((double*) (integrals->imageData + (cell.y + cell.height) *    
14.             (integrals->widthStep)))[cell.x + cell.width];   
15.   
16.         float c = ((double*) (integrals->imageData + (cell.y) * (integrals-    
17.             >widthStep)))[cell.x + cell.width];   
18.   
19.         float d = ((double*) (integrals->imageData + (cell.y + cell.height) *    
20.             (integrals->widthStep)))[cell.x];   
21.   
22.         ((float*) hog_cell->data.fl) = (a + b) - (c + d);   
23.   
24.     }   
25.   
26.   
27.     /*Normalize the matrix*/  
28.     if (normalization != -1){   
29.         cvNormalize(hog_cell, hog_cell, 1, 0, normalization);   
30.     }   
31.   
32. }