Moravec，Forstner和Harris特征点提取opencv实现

效果图：

Moravec

 

Forstner

 

Harris:

 

代码：

//选极值点作为特征点，抑制局部非最大
Mat FindMaxValue(Mat M, int win_size, int flag)

//特征提取算子//
//Moravec算子
void Moravec(Mat img1,Mat img,int win_size1,int win_size2)
{
 
    int size1 = floor(win_size1 / 2); //一半向下取整
    int r, c;

    //创建兴趣值矩阵IV
    Mat IV = Mat::zeros(Size(img.cols, img.rows), CV_32FC1);

    for (r = size1; r < img.rows - size1; r++)
        for (c = size1; c < img.cols - size1; c++)
        {
            int V[4] = { 0 };
            int iv = 0;//兴趣值

            for (int i = -size1; i < size1; i++)
            {
                //计算4个方向灰度差平方和
                V[0] += pow(img.at(r, c + i) - img.at(r, c + i + 1), 2);
                V[1] += pow(img.at(r + i, c + i) - img.at(r + i + 1, c + i + 1), 2);
                V[2] += pow(img.at(r + i, c) - img.at(r + i + 1, c), 2);
                V[3] += pow(img.at(r - i, c + i) - img.at(r - i - 1, c + i + 1), 2);
            }

            iv = *min_element(V, V + 4);//取最小值
            IV.at(r, c) = iv; //uchar会越界
        }

    //计算经验阈值
    int iv_mean = 0;
    int iv_max = 0;

    for (r = size1; r < img.rows - size1; r++)
        for (c = size1; c < img.cols - size1; c++)
        {
            int tmp = IV.at(r, c);
            iv_mean += tmp;
            if (tmp > iv_max)
                iv_max = tmp;
        }
  
    iv_mean = iv_mean / ((img.rows - 2 * size1) * (img.cols - 2 * size1));
    cout << "iv_mean: " << iv_mean << endl;
    cout << "iv_max: " << iv_max << endl;

 
    int threshold = 3000;
    //int threshold = iv_mean;
    int pnum = 0;//候选点数
    for (r = size1; r < img.rows - size1; r++)
        for (c = size1; c < img.cols - size1; c++)
        {
            if (IV.at(r, c) <= threshold)
                IV.at(r, c) = 0;
            else
                pnum++;
        }

    //选极值点抑制局部非最大
    IV = FindMaxValue(IV, win_size2, 1);
    int size2 = floor(win_size2 / 2);//w2为抑制窗口大小
    Mat showImg;
    img1.copyTo(showImg);
    int count = ShowFeaturePoints(showImg, IV, size2, 1);

}

//Forstner算子
void Forstner(Mat oriimg,Mat input,int win_size)
{
    Mat img;
    input.copyTo(img);
    int width = img.cols;
    int height = img.rows;
    cout << "图像行数：" << height << "  图像列数：" << width << endl;

    int size = floor(win_size / 2); //窗口二分之一
    //创建存兴趣值q,w的矩阵
    Mat Mat_q = Mat::zeros(Size(width, height), CV_32FC1);
    Mat Mat_w = Mat::zeros(Size(width, height), CV_32FC1);

    double w_ave = 0;
    //计算窗口中灰度的协方差矩阵
    for (int r = size; r < height - size; r++)
        for (int c = size; c < width - size; c++)
        {
            double gu_2 = 0;
            double gv_2 = 0;
            double gu_gv = 0;

            for (int i = r - size; i < r + size; i++)
                for (int j = c - size; j < c + size; j++)
                {
                    //计算各像素的Robert梯度gu,gv
                    double gu = 0;
                    double gv = 0;
                    gu = img.at(i + 1, j + 1) - img.at(i, j);
                    gv= img.at(i, j + 1) - img.at(i + 1, j);

                    gu_2 += gu * gu;
                    gv_2 += gv * gv;
                    gu_gv += gu * gv;
                }

     
            double DetN = gu_2 * gv_2 - gu_gv * gu_gv;
            double trN = gu_2 + gv_2;
            if (trN == 0)
                continue;//cout << "矩阵N的迹为0！" << endl;

            double w = DetN / trN;
            w_ave += w;
            double q = 4 * DetN / (trN * trN);

            Mat_w.at(r, c) = w;
            Mat_q.at(r, c) = q;
        }

    w_ave = w_ave / ((height - 2 * double(size)) * (width - 2 * double(size)));
    double f = 8.0;

    double Tq = 0.75;
    double Tw = f * w_ave;

    int count = 0;//统计待选点数
    for (int r = 0; r < height; r++)
        for (int c = 0; c < width; c++)
        {
            if (Mat_q.at(r, c) > Tq && Mat_w.at(r, c) > Tw)
                count++;
            else
                Mat_w.at(r, c) = 0;
        }

    int win_size2 = 5;
    int size2 = floor(win_size2 / 2);

    Mat_w = FindMaxValue(Mat_w, win_size2, 2);

    Mat showImg;
    oriimg.copyTo(showImg);
    //count = 0;
    count = ShowFeaturePoints(showImg, Mat_w, size2, 2);

}

//Harris
void Harris(Mat oriimg, Mat input, int win_size)
{
    Mat img;
    input.copyTo(img);
    int width = img.cols;
    int height = img.rows;

    Mat gx = Mat::zeros(Size(width - 1, height - 1), CV_32FC1);
    Mat gy = Mat::zeros(Size(width - 1, height - 1), CV_32FC1);

  
    for (int r = 0; r < height - 1; r++)
        for (int c = 0; c < width - 1; c++)
        {
            gx.at(r, c) = img.at(r, c + 1) - img.at(r, c);
            gy.at(r, c) = img.at(r + 1, c) - img.at(r, c);
        }

 
    Mat Gx, Gy;
    GaussianBlur(gx, Gx, Size(5, 5), 0.8, 0.8);
    GaussianBlur(gy, Gy, Size(5, 5), 0.8, 0.8);

    Mat I = Mat::zeros(Size(width - 1, height - 1), CV_32FC1);
    int size = floor(win_size / 2);
    double I_ave = 0;

    for (int r = size; r < height - size - 1; r++)
        for (int c = size; c < width - size - 1; c++)
        {
            double gx_2 = 0;
            double gy_2 = 0;
            double gx_gy = 0;

  
            for (int i = r - size; i < r + size; i++)
                for (int j = c - size; j < c + size; j++)
                {
                    gx_2 += Gx.at(i, j) * Gx.at(i, j);
                    gy_2 += Gy.at(i, j) * Gy.at(i, j);
                    gx_gy += Gx.at(i, j) * Gy.at(i, j);
                }
            double DetM = gx_2 * gy_2 - gx_gy * gx_gy;
            double trM = gx_2 + gy_2;

            double k = 0.04;
            I.at(r, c) = DetM - k * trM * trM;
            I_ave += I.at(r, c);
        }

    I_ave = I_ave / ((height - 1 - 2 * size) * (width - 1 - 2 * size));
    
    double thre = 10 * I_ave;

    for (int r = size; r < height - size - 1; r++)
        for (int c = size; c < width - size - 1; c++)
            if (I.at(r, c) < thre)
                I.at(r, c) = 0;

    //抑制局部非最大
    int win_size2 = 5;
    int size2 = floor(win_size2 / 2);
    I = FindMaxValue(I, win_size2, 2);

    Mat showImg;
    oriimg.copyTo(showImg);
    //count = 0;
    int count = ShowFeaturePoints(showImg, I, size2, 2);
}