Integral Histogram for fast Calculation of HOG Features

The original article can be found in

http://smsoftdev-solutions.blogspot.com/2009/08/integral-histogram-for-fast-calculation.html?showComment=1300196607644#c294288376459522527

However, there are some problems when you are going to use it:
First, the input image should be a BGR color image.
Second, the variable "xsobel" and "ysobel" must be "IPL_DEPTH_32F".
Third, the author does not write out the function "doSobel".


To avoid these problems,
you can change a little bit of original codes:

IplImage* img_gray = cvCreateImage(cvGetSize(in), IPL_DEPTH_8U,1);
if(in->nChannels>1){
cvCvtColor(in, img_gray, CV_BGR2GRAY);
}
else
{
cvCopy(in,img_gray);
}
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 = cvCreateImage(cvGetSize(in), IPL_DEPTH_16S, 1);
IplImage *ysobel = cvCreateImage(cvGetSize(in), IPL_DEPTH_16S, 1);

cvSobel(img_gray,xsobel,1, 0, 3);
cvSobel(img_gray,ysobel, 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[i]);
}
/* 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 height; y++) {
/* ptr1 and ptr2 point to beginning of the current row in the xsobel
and ysobel images respectively. ptrs[i] point to the beginning of the
current rows in the bin images */
short* ptr1 = (short*) (xsobel->imageData + y * (xsobel->widthStep));
short* ptr2 = (short*) (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[i]->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 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);
}