C++ Forstner算子进行点特征提取

C++实现Forster点特征提取

---

---

前言

该算子由德国伯恩大学摄影测量系主任、教授 Forster1984年提出的。本文介绍如何使用C++实现该算子

---

一、原理

计算各像素的Robert’s梯度和像素（c，r）为中心的一个窗口的灰度协方差矩阵，在影像中寻找具有尽可能小而接近圆的误差椭圆的点作为特征点

二、使用步骤

1.实现步骤

step1: 计算各个像素的Robert梯度
step2：窗口中灰度的协方差矩阵
step3：兴趣值q和权w的计算公式
step4: 阈值的选定
step5: 确定待选点
step6: 选取极值点
（阈值范围：）

2.C++代码

// forstner.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//
#include <iostream>
#include<string>
#include<cmath>
#include<vector>
#include<gdal.h>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/opencv.hpp>
#include "cpl_conv.h"

#include<cstring>

using namespace std;
using namespace cv;


void Forstner(Mat img,Mat img2,string str2)
{
    
    /*
    step1:计算各个像素的Robert梯度
    f_x = f(x + 1, y + 1) - f(x, y)
    f_y = f(x + 1, y) - f(x, y + 1)
    */
    
    int gu=0, gv=0;//xy方向梯度
    int r, c;
    int add;
    int rows = img.rows - 1;
    int cols = img.cols - 1;
   // Mat G(rows + 1, cols + 1, CV_32FC3, Scalar::all(0));//计算各像素的Robert梯度
    Mat GU(rows + 1, cols + 1, CV_32FC1, Scalar::all(0));//用来存放gu的矩阵
    Mat GV(rows + 1, cols + 1, CV_32FC1, Scalar::all(0));//用来存放gv的矩阵
   
    for ( r = 0; r < rows +1; r++)
    {
        for ( c = 0; c < cols +1; c++)
        {
            add = 1;//des=1;
            //if (c == 0 || r == 0) des = 0;
            if (r == rows || c == cols) add = 0;
            GU.at<float>(r, c)= (img.at<uchar>(r + add, c + add) - img.at<uchar>(r, c));
            GV.at<float>(r, c)= (img.at<uchar>(r + add, c) - img.at<uchar>(r, c + add));
            //cout << "GU:"<< GU.at<float>(r, c) <<"GV:"<<GV.at<float>(r, c)<< endl;      
        }
        //cout << "lala" << endl;
    }

    float MaxValue, MaxW;//极值点的大小
    int win_r,win_c;//窗口内的移动位置
    int size = 7;//窗口的大小
    int maxR, maxC;//极值点的位置
    float gu2, gv2, guv;
    float detN = 0, trN = 0;//行列式和迹
    float w ,q;//兴趣值
    float Tq = 0.6;//阈值选取
    
    float T_w =0;
    double w_ba = 0;//w的平均值
    Mat N(rows+1,cols+1,CV_32FC3,Scalar::all(0));//协方差矩阵
    Mat W(rows+1,cols+1,CV_32FC1,Scalar::all(0));//存放w的矩阵
    Mat M_q(rows + 1, cols + 1, CV_32FC1, Scalar::all(0));//存放q的矩阵
    double num =0,number=0;
    //cout << "95" << endl;
    for ( r = 0; r <= rows - size; r+=size)
    {
        for ( c = 0; c <= cols - size; c+=size)
        {
            
            gu2 = 0,gv2=0, guv=0;
            w = 0, q = 0;
            //Tw = 0;
            //cout << "123" << endl;
            for ( win_r = 0; win_r < size; win_r++)
            {
                for (win_c = 0; win_c < size; win_c++)
                {
                    //step2:计算窗口内灰度的协方差矩阵
                    gu = GU.at<float>(r + win_r, c + win_c);
                    gv = GV.at<float>(r + win_r, c + win_c);

                    gu2 = gu2 + pow(gu,2);//求窗口范围内x方向梯度和
                    gv2 = gv2 + pow(gv, 2);
                    guv = guv + gu * gv;
                    //cout << "gu2:"<< gu2 <<"gv2"<< gv2 <<"guv"<< guv << endl;
                }
            }
            //step3:计算兴趣值q与w
            //cout << "dudu" << endl;
            detN = gu2 * gv2 - pow(guv, 2);
            trN = gu2 + gv2;
            q = (4.0 * detN) / (pow(trN, 2));
            if (q>Tq)
            {
                  w = detN / trN; 
                  num = num + w;//求w的总和
            }
            else
            {
                w = 0;
                num = num + w;
            }
          
            //cout << "w:"<< w << endl;
            
            for (win_r = 0; win_r < size; win_r++)
            {
                for (win_c = 0; win_c < size; win_c++)
                {
                    W.at<float>(r + win_r, c + win_c) = w;//将计算得到的wq存入矩阵中
                    M_q.at<float>(r + win_r, c + win_c) = q;
                    //cout << "M_q.at<float>(r + win_r, c + win_c):"<< M_q.at<float>(r + win_r, c + win_c) <<" W.at<float>(r + win_r, c + win_c)"<< W.at<float>(r + win_r, c + win_c) << endl;
                }
            }
            number = number + 1;//窗口个数
            //cout << "number:"<< number << endl;
            //cout << "num"<< num << endl;
        }
       // cout << "123" << endl;
    }
    //至此我们得到整幅图像中所有像素的wq值
    cout << "789" << endl;
    //cout << "num:"<<num<<"number:"<< number<< endl;
    w_ba = num / number;//w的均值
    //w_ba = num / ((rows + 1) * (cols +1));
    T_w = 1.0 * w_ba;
    cout << "T_w:"<< T_w << endl;

    //step4:确定待选点 Tq=0.6;
    for ( r = 0; r <= rows -size; r+=size)
    {
        for ( c = 0; c <= cols-size; c+=size)
        {
            maxR = 0, maxC = 0;
            MaxW = 0;
            for (win_r = 0; win_r < size; win_r++)//逐窗口选点
            {
                for ( win_c = 0; win_c < size ; win_c++)
                {
                    //待选点
                    if (M_q.at<float>(r+win_r, c+win_c) > Tq && W.at<float>(r+win_r, c+win_c) > T_w )
                    { 
                        //step5:局部非最大抑制,在一个适当窗口中选择最大的待选点
                        if (W.at<float>(r + win_r, c + win_c) > MaxW)
                        {
                            MaxW = W.at<float>(r + win_r, c + win_c);
                            maxR = r + win_r;
                            maxC = c + win_c;
                        }
                        //cout << "maxR"<< maxR << "maxC" << maxC<< endl;
                    }
                }
            }

            //画图
            Point pt = Point(maxC, maxR);
            int pointColor = (0, 0, 255);
            circle(img2, pt, 3, pointColor);
        }
    }

   

    namedWindow("Forster",WINDOW_NORMAL);
    imshow("Forster", img2);
    imwrite(str2,img2);
    cv::waitKey(0);
}


int main()
{

    string str = "D://学习课件//数字摄影测量//实验指导//实验指导//实验数据//实验1//DSC00438.JPG"; 
    string str2 = "D://学习课件//数字摄影测量//实验指导//实验指导//Forstner.JPG";//存放图片的位置
    Mat img,img2;
    img = imread(str,IMREAD_GRAYSCALE);//读取灰色图像
    img2 = imread(str,IMREAD_COLOR);//彩色图像
    if (img.empty() || img2.empty())
    {
        cout << "无法打开图像" << endl;
    }
    else
    {
        cout << "OKKKKKKKKKKKK" << endl;
    }
    Forstner(img,img2,str2);
    return 0;

   // std::cout << "Hello World!\n";
}

// 运行程序: Ctrl + F5 或调试 >“开始执行(不调试)”菜单
// 调试程序: F5 或调试 >“开始调试”菜单

// 入门使用技巧: 
//   1. 使用解决方案资源管理器窗口添加/管理文件
//   2. 使用团队资源管理器窗口连接到源代码管理
//   3. 使用输出窗口查看生成输出和其他消息
//   4. 使用错误列表窗口查看错误
//   5. 转到“项目”>“添加新项”以创建新的代码文件，或转到“项目”>“添加现有项”以将现有代码文件添加到项目
//   6. 将来，若要再次打开此项目，请转到“文件”>“打开”>“项目”并选择 .sln 文件

---

三、运行结果

原图：

结果：

补充：评论区有位友友指出了代码中的一处疏漏，说的很在理。因为我为了跑的快一点 一次性平移了7个像元（窗口大小），所以没法求出所有所有像元的wq值，就干脆都赋给窗口内了，按理说是不应该这样的。可能这也是结果中一些细节没有提取出来的原因，谢谢这位友友~
从效果上看，大部分点还是可以提取出来的，如果大家要求比较高，可以自行改成逐像元计算wq。