局部边沿方向特征的快速算法

    这几天是元旦假期，正好有时间写个实验程序验证一下新想法。原来写的这个 “用局部边沿方向特征做图像识别目标检测” 检测效果还好，但是速度比较慢。这次的新算法可以快速计算图像的局部边沿方向。新算法首先计算图像的积分图，然后根据积分图计算图像的近似梯度场，在图像的边沿上梯度方向是法线方向，切线方向就是图像的局部边沿方向，这个算法比原来基于卷积的边沿方向算法快很多。

    图像的近似梯度场图：

求近似梯度场的算法：

inline  HaarGradient  img_haar_gradient_l4(IplImage  *img_gray, CvMat* mat_Integral, CvPoint p)
{
        int  l = 4;

        double   fix01,fix02,fix03,fix04;

        double   sum;

        int x,y;

        sum=0.0;
        x = p.x - l;
        y = p.y - l;
        float  ve = CV_IMAGE_ELEM(img_gray, uchar, y, x);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y, x+1);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y+1, x);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y+1, x+1);
        sum+=ve;
        fix01 = sum;

        sum=0.0;
        x = p.x - l;
        y = p.y + l;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y, x);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y, x+1);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y-1, x);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y-1, x-1);
        sum+=ve;
        fix02 = sum;

        sum=0.0;
        x = p.x + l;
        y = p.y + l;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y, x);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y, x-1);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y-1, x);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y-1, x-1);
        sum+=ve;
        fix03 = sum;

        sum=0.0;
        x = p.x + l;
        y = p.y - l;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y, x);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y, x-1);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y+1, x);
        sum+=ve;
        ve = CV_IMAGE_ELEM(img_gray, uchar, y+1, x-1);
        sum+=ve;
        fix04 = sum;


        HaarGradient  grad;

        int  x1 = p.x - l, x2 = p.x;
        int  y1 = p.y - l, y2 = p.y + l;

        float   *p_x1y1,*p_x1y2,*p_x2y1,*p_x2y2;
        double  inte_x1y1,inte_x1y2,inte_x2y1,inte_x2y2;

        p_x1y1 = (float*)cvPtr2D(mat_Integral, y1, x1);
        p_x1y2 = (float*)cvPtr2D(mat_Integral, y2, x1);
        p_x2y1 = (float*)cvPtr2D(mat_Integral, y1, x2);
        p_x2y2 = (float*)cvPtr2D(mat_Integral, y2, x2);

        double   inte_xy01 = *p_x2y2 - *p_x1y2 - *p_x2y1 + *p_x1y1 - fix01 - fix02;

        x1 = p.x,        x2 = p.x + l;
        y1 = p.y - l,    y2 = p.y + l;


        p_x1y1 = (float*)cvPtr2D(mat_Integral, y1, x1);
        p_x1y2 = (float*)cvPtr2D(mat_Integral, y2, x1);
        p_x2y1 = (float*)cvPtr2D(mat_Integral, y1, x2);
        p_x2y2 = (float*)cvPtr2D(mat_Integral, y2, x2);

        double   inte_xy02 = *p_x2y2 - *p_x1y2 - *p_x2y1 + *p_x1y1 - fix03 - fix04;

        double   haar_horizontal = inte_xy01 -  inte_xy02;

        x1 = p.x - l, x2 = p.x + l;
        y1 = p.y - l, y2 = p.y;

        *p_x1y1,*p_x1y2,*p_x2y1,*p_x2y2;
        inte_x1y1,inte_x1y2,inte_x2y1,inte_x2y2;

        p_x1y1 = (float*)cvPtr2D(mat_Integral, y1, x1);
        p_x1y2 = (float*)cvPtr2D(mat_Integral, y2, x1);
        p_x2y1 = (float*)cvPtr2D(mat_Integral, y1, x2);
        p_x2y2 = (float*)cvPtr2D(mat_Integral, y2, x2);

        inte_xy01 = *p_x2y2 - *p_x1y2 - *p_x2y1 + *p_x1y1 - fix01 - fix04;

        x1 = p.x - l, x2 = p.x + l;
        y1 = p.y,     y2 = p.y + l;


        p_x1y1 = (float*)cvPtr2D(mat_Integral, y1, x1);
        p_x1y2 = (float*)cvPtr2D(mat_Integral, y2, x1);
        p_x2y1 = (float*)cvPtr2D(mat_Integral, y1, x2);
        p_x2y2 = (float*)cvPtr2D(mat_Integral, y2, x2);

        inte_xy02 = *p_x2y2 - *p_x1y2 - *p_x2y1 + *p_x1y1 - fix02 - fix03;

        double   haar_vertical = inte_xy02 -  inte_xy01;

        double   haar_vh = fabs(haar_vertical/haar_horizontal);

        //grad.angle = atan2(haar_vertical, haar_horizontal);
        grad.angle = (CV_PI/2.0 * haar_vh) / (1.0 + haar_vh); //(pi/2 * y) / (1 + y)

        if (haar_vertical < 0.0 && haar_horizontal < 0.0) grad.angle += CV_PI;

        if (haar_vertical < 0.0 && haar_horizontal > 0.0) grad.angle = 2*CV_PI - grad.angle;

        if (haar_vertical > 0.0 && haar_horizontal < 0.0) grad.angle = CV_PI - grad.angle;

        //if (haar_vertical < 0.0)  grad.angle += CV_PI;

        grad.mod = sqrt(haar_vertical*haar_vertical + haar_horizontal*haar_horizontal);

        if (p.x == 50 && p.y == 50)
        {
                std::cout << "haar_vertical,  haar_horizontal  " << haar_vertical << ", " << haar_horizontal  << std::endl;
                std::cout << " haar_vh = " << haar_vh;
                std::cout << "     ang = " << grad.angle/CV_PI * 180.0 << std::endl;
        }

        return  grad;
}

根据梯度可以得到局部边沿方向图，和原来的基于卷积的算法结果一致。

  这样特征提取问题就解决了，下一步还要改进匹配算法，提高匹配效率，加强匹配算法的弹性适应一定的视角变化，试试分组投票等算法，不过元旦假期马上就结束了，只能以后放假再说了，估计春节假期能写完新程序。

                                                                  2013.01.03

                                                                   zhujiang