Harris角點算法

時間  2019-12-12
標籤 harris 算法 简体版
原文   原文鏈接

特徵點檢測普遍應用到目標匹配、目標跟蹤、三維重建等應用中，在進行目標建模時會對圖像進行目標特徵的提取，經常使用的有顏色、角點、特徵點、輪廓、紋理等特徵。如今開始講解經常使用的特徵點檢測，其中Harris角點檢測是特徵點檢測的基礎，提出了應用鄰近像素點灰度差值概念，從而進行判斷是否爲角點、邊緣、平滑區域。Harris角點檢測原理是利用移動的窗口在圖像中計算灰度變化值，其中關鍵流程包括轉化爲灰度圖像、計算差分圖像、高斯平滑、計算局部極值、確認角點。html

1、基礎知識

圖像的變化類型：算法

在特徵點檢測中常常提出尺度不變、旋轉不變、抗噪聲影響等，這些是判斷特徵點是否穩定的指標。編程

性能較好的角點：數組

  1. 檢測出圖像中「真實」的角點
  2. 準確的定位性能
  3. 很高的重複檢測率
  4. 噪聲的魯棒性
  5. 較高的計算效率

角點的類型：ide

基於圖像灰度的方法經過計算點的曲率及梯度來檢測角點，避免了第一類方法存在的缺陷，此類方法主要有Moravec算子、Forstner算子、Harris算子、SUSAN算子等。函數

image

2、Harris角點原理

一、算法思想

角點原理來源於人對角點的感性判斷，即圖像在各個方向灰度有明顯變化。算法的核心是利用局部窗口在圖像上進行移動判斷灰度發生較大的變化，因此此窗口用於計算圖像的灰度變化爲：[-1,0,1;-1,0,1;-1,0,1][-1,-1,-1;0,0,0;1,1,1]。人各個方向上移動這個特徵的小窗口，如圖3中窗口內區域的灰度發生了較大的變化，那麼就認爲在窗口內遇到了角點。如圖1中，窗口內圖像的灰度沒有發生變化，那麼窗口內就不存在角點；若是窗口在某一個方向移動時，窗口內圖像的灰度發生了較大的變化，而在另外一些方向上沒有發生變化，那麼，窗口內的圖像可能就是一條直線的線段。性能

image

二、數學模型

根據算法思想，構建數學模型，計算移動窗口的的灰度差值。學習

爲了減少計算量，利用泰勒級數進行簡化公式：ui

上圖中W函數表示窗口函數，M矩陣爲偏導數矩陣。對於矩陣能夠進行對稱矩陣的變化，假設利用兩個特徵值進行替代，其幾何含義相似下圖中的表達。在幾何url

模型中經過判斷兩個特徵值的大小，來斷定像素的屬性。

M爲梯度的協方差矩陣 ，在實際應用中爲了可以應用更好的編程，定義了角點響應函數R，經過斷定R大小來判斷像素是否爲角點。

R取決於M的特徵值，對於角點|R|很大，平坦的區域|R|很小，邊緣的R爲負值。

 

三、算法流程

1.利用水平，豎直差分算子對圖像的每一個像素進行濾波以求得Ix,Iy,進而求得M中的四個元素的值。

算法源碼：

代碼中若是array爲-1,0,1,-1,0,1,-1,0,1}則是求解X方向的，若是爲{-1,-1,-1,0,0,0,1,1,1}爲Y方向的，則Ix和Iy求解結束

//C代碼
//這裏的size/2是爲了避免把圖像邊界算進去。
//array也就是3*3的窗口函數爲：double dx[9]={-1,0,1,-1,0,1,-1,0,1};或者 //定義垂直方向差分算子並求Iy
    double dy[9]={-1,-1,-1,0,0,0,1,1,1};

CvMat *mbys(CvMat *mat,int xwidth,int ywidth,double *array,int size1,int size2)//size 
{
    int i,j;
    int i1,j1;
    int px,py;
    int m;
    CvMat *mat1;
    mat1=cvCloneMat(mat);
//爲了去除邊界，從框體一半開始遍歷
    for(i=size1/2;i<ywidth-size1/2;i++)
        for(j=size2/2;j<xwidth-size2/2;j++)                         
        {
            m=0;
            for(i1=0;i1<size1;i1++)
                for(j1=0;j1<size2;j1++)
                {
                    px=i-size1/2+i1;
                    py=j-size2/2+j1;
                    //CV_MAT_ELEM訪問矩陣元素
                  //每一個元素和窗體函數遍歷相加
                    m+=CV_MAT_ELEM(*mat,double,px,py)*array[i1*size1+j1];            
                }
               //賦值
                CV_MAT_ELEM(*mat1,double,i,j)=m;
        }
        return mat1;
}

求解IX2相對比較簡單，像素相乘便可。下面爲基於opencv的C++版本，非常簡單

//構建模板
Mat xKernel = (Mat_<double>(1,3) << -1, 0, 1);
Mat yKernel = xKernel.t();
//計算IX和IY
Mat Ix,Iy;
//可參考filter2d的定義
filter2D(gray, Ix, CV_64F, xKernel);
filter2D(gray, Iy, CV_64F, yKernel);
//計算其平方
Mat Ix2,Iy2,Ixy;
Ix2 = Ix.mul(Ix);
Iy2 = Iy.mul(Iy);
Ixy = Ix.mul(Iy);

2.對M的四個元素進行高斯平滑濾波，爲的是消除一些沒必要要的孤立點和凸起，獲得新的矩陣M。

//本例中使用5×5的高斯模板
    //計算模板參數
    //int gausswidth=5;
    //double sigma=0.8;
    double *g=new double[gausswidth*gausswidth];
    for(i=0;i<gausswidth;i++)//定義模板
        for(j=0;j<gausswidth;j++)
            g[i*gausswidth+j]=exp(-((i-int(gausswidth/2))*(i-int(gausswidth/2))+(j-int(gausswidth/2))*(j-int(gausswidth/2)))/(2*sigma));

    //歸一化：使模板參數之和爲1（其實此步能夠省略）
    double total=0;
    for(i=0;i<gausswidth*gausswidth;i++)
        total+=g[i];
    for(i=0;i<gausswidth;i++)
        for(j=0;j<gausswidth;j++)
            g[i*gausswidth+j]/=total;

    //進行高斯平滑
    mat_Ix2=mbys(mat_Ix2,cxDIB,cyDIB,g,gausswidth,gausswidth);
    mat_Iy2=mbys(mat_Iy2,cxDIB,cyDIB,g,gausswidth,gausswidth);
    mat_Ixy=mbys(mat_Ixy,cxDIB,cyDIB,g,gausswidth,gausswidth);

利用opencv接口則是很簡單

//構建高斯函數
Mat gaussKernel = getGaussianKernel(7, 1);
//卷積計算
filter2D(Ix2, Ix2, CV_64F, gaussKernel);
filter2D(Iy2, Iy2, CV_64F, gaussKernel);
filter2D(Ixy, Ixy, CV_64F, gaussKernel);

 

3.接下來利用M計算對應每一個像素的角點響應函數R，即：

也可使用改進的R：
R=[Ix^2*Iy^2-(Ix*Iy)^2]/(Ix^2+Iy^2);裏面沒有隨意給定的參數k，取值應當比第一個使人滿意。

//計算cim：即cornerness of image，咱們把它稱作‘角點量’
    CvMat *mat_cim;
    mat_cim=mbcim(mat_Ix2,mat_Iy2,mat_Ixy,cxDIB,cyDIB);
//用來求得響應度
CvMat *mbcim(CvMat *mat1,CvMat *mat2,CvMat *mat3,int xwidth,int ywidth)
{
    int i,j;
    CvMat *mat;
    mat=cvCloneMat(mat1);
    for(i = 0; i <ywidth; i++)
    {
        for(j = 0; j < xwidth; j++)
        {
            //注意：要在分母中加入一個極小量以防止除數爲零溢出
            CV_MAT_ELEM(*mat,double,i,j)=(CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j)-
                CV_MAT_ELEM(*mat3,double,i,j)*CV_MAT_ELEM(*mat3,double,i,j))/
                (CV_MAT_ELEM(*mat1,double,i,j)+CV_MAT_ELEM(*mat2,double,i,j)+0.000001);
        }
    }
    return mat;
}
//opencv代碼
Mat cornerStrength(gray.size(), gray.type());
    for (int i = 0; i < gray.rows; i++)
    {
        for (int j = 0; j < gray.cols; j++)
        {
            double det_m = Ix2.at<double>(i,j) * Iy2.at<double>(i,j) - Ixy.at<double>(i,j) * Ixy.at<double>(i,j);
            double trace_m = Ix2.at<double>(i,j) + Iy2.at<double>(i,j);
            cornerStrength.at<double>(i,j) = det_m - alpha * trace_m *trace_m;
        }
    }

四、局部極大值抑制，同時選取其極大值

//--------------------------------------------------------------------------
    //                 第四步：進行局部非極大值抑制
    //--------------------------------------------------------------------------
    CvMat *mat_locmax;
    //const int size=7;
    mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,size);
//     cout<<CV_MAT_ELEM(*mat_locmax,double,cyDIB-1,cxDIB-1)<<endl;
//用來求得局部極大值
CvMat *mblocmax(CvMat *mat1,int xwidth,int ywidth,int size)
{
    int i,j;
    double max=-1000;
    int i1,j1;
    int px,py;
    CvMat *mat;
    mat=cvCloneMat(mat1);
    for(i=size/2;i<ywidth-size/2;i++)
        for(j=size/2;j<xwidth-size/2;j++)
        {
            max=-10000;
            for(i1=0;i1<size;i1++)
                for(j1=0;j1<size;j1++)
                {
                    px=i-size/2+i1;
                    py=j-size/2+j1;
                    if(CV_MAT_ELEM(*mat1,double,px,py)>max)
                        max=CV_MAT_ELEM(*mat1,double,px,py);

                }
                if(max>0)
                    CV_MAT_ELEM(*mat,double,i,j)=max;
                else
                    CV_MAT_ELEM(*mat,double,i,j)=0;
        }
        return mat;
}

在opencv中應用maxminloc函數

// threshold
double maxStrength;
minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
Mat dilated;
Mat localMax;
//默認3*3核膨脹，膨脹以後，除了局部最大值點和原來相同，其它非局部最大值點被  
 //3*3鄰域內的最大值點取代
dilate(cornerStrength, dilated, Mat());
//與原圖相比，只剩下和原圖值相同的點，這些點都是局部最大值點，保存到localMax 
compare(cornerStrength, dilated, localMax, CMP_EQ);

5.在矩陣R中，同時知足R(i,j)大於必定閾值threshold和R(i,j)是某領域內的局部極大值，則被認爲是角點。

cv::Mat cornerMap;  
  // 根據角點響應最大值計算閾值  
  threshold= qualityLevel*maxStrength;  
  cv::threshold(cornerStrength,cornerTh,  
   threshold,255,cv::THRESH_BINARY);  
  // 轉爲8-bit圖 
 cornerTh.convertTo(cornerMap,CV_8U);// 和局部最大值圖與，剩下角點局部最大值圖，即：完成非最大值抑制  
  cv::bitwise_and(cornerMap,localMax,cornerMap);
imgDst = cornerMap.clone();
for( int y = 0; y < cornerMap.rows; y++ ) 
{  
        const uchar* rowPtr = cornerMap.ptr<uchar>(y);  
            for( int x = 0; x < cornerMap.cols; x++ ) 
    {  
               // 非零點就是角點  
              if (rowPtr[x]) 
    {  
                        points.push_back(cv::Point(x,y));  
                 }  
                 }  
 }

 

 

 

 

 

 

 

 

 

 

 

3、算法源碼

一、C版本

裏面函數最基礎的代碼解釋和註釋

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
#define B(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)*3]
#define G(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)*3+1]
#define R(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)*3+2]
#define S(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)]

//卷積計算求Ix,Iy,以及濾波
//a指向的數組是size1*size2大小的...求導
CvMat *mbys(CvMat *mat,int xwidth,int ywidth,double *a,int size1,int size2)
{
    int i,j;
    int i1,j1;
    int px,py;
    int m;
    CvMat *mat1;
    mat1=cvCloneMat(mat);
    for(i=size1/2;i<ywidth-size1/2;i++)
        for(j=size2/2;j<xwidth-size2/2;j++)
        {
            m=0;
            for(i1=0;i1<size1;i1++)
                for(j1=0;j1<size2;j1++)
                {
                    px=i-size1/2+i1;
                    py=j-size2/2+j1;
                    //CV_MAT_ELEM訪問矩陣元素
                    m+=CV_MAT_ELEM(*mat,double,px,py)*a[i1*size1+j1];            
                }
                CV_MAT_ELEM(*mat1,double,i,j)=m;
        }
        return mat1;
}
//計算Ix2,Iy2,Ixy
CvMat *mbxy(CvMat *mat1,CvMat *mat2,int xwidth,int ywidth)
{
    int i,j;
    CvMat *mat3;
    mat3=cvCloneMat(mat1);
    for(i=0;i<ywidth;i++)
        for (j=0;j<xwidth;j++)
        {
            CV_MAT_ELEM(*mat3,double,i,j)=CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j);
            
        }
        return mat3;
}

//用來求得響應度
CvMat *mbcim(CvMat *mat1,CvMat *mat2,CvMat *mat3,int xwidth,int ywidth)
{
    int i,j;
    CvMat *mat;
    mat=cvCloneMat(mat1);
    for(i = 0; i <ywidth; i++)
    {
        for(j = 0; j < xwidth; j++)
        {
            //注意：要在分母中加入一個極小量以防止除數爲零溢出
            CV_MAT_ELEM(*mat,double,i,j)=(CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j)-
                CV_MAT_ELEM(*mat3,double,i,j)*CV_MAT_ELEM(*mat3,double,i,j))/
                (CV_MAT_ELEM(*mat1,double,i,j)+CV_MAT_ELEM(*mat2,double,i,j)+0.000001);
        }
    }
    return mat;
}
//用來求得局部極大值
CvMat *mblocmax(CvMat *mat1,int xwidth,int ywidth,int size)
{
    int i,j;
    double max=-1000;
    int i1,j1;
    int px,py;
    CvMat *mat;
    mat=cvCloneMat(mat1);
    for(i=size/2;i<ywidth-size/2;i++)
        for(j=size/2;j<xwidth-size/2;j++)
        {
            max=-10000;
            for(i1=0;i1<size;i1++)
                for(j1=0;j1<size;j1++)
                {
                    px=i-size/2+i1;
                    py=j-size/2+j1;
                    if(CV_MAT_ELEM(*mat1,double,px,py)>max)
                        max=CV_MAT_ELEM(*mat1,double,px,py);

                }
                if(max>0)
                    CV_MAT_ELEM(*mat,double,i,j)=max;
                else
                    CV_MAT_ELEM(*mat,double,i,j)=0;
        }
        return mat;
}
//用來確認角點
CvMat *mbcorner(CvMat *mat1,CvMat *mat2,int xwidth,int ywidth,int size,double thresh)
{
    CvMat *mat;
    int i,j;
    mat=cvCreateMat(ywidth,xwidth,CV_32FC1);
    for(i=size/2;i<ywidth-size/2;i++)
        for(j=size/2;j<xwidth-size/2;j++)
        {
            if(CV_MAT_ELEM(*mat1,double,i,j)==CV_MAT_ELEM(*mat2,double,i,j))//首先取得局部極大值
                if(CV_MAT_ELEM(*mat1,double,i,j)>thresh)//而後大於這個閾值
                    CV_MAT_ELEM(*mat,int,i,j)=255;//知足上兩個條件，纔是角點！
            else
                CV_MAT_ELEM(*mat,int,i,j)=0;
        }
        return mat;
}

CvPoint* CHarris::harris_features(IplImage *src,int gausswidth,double sigma,int size,int threshold)
{
    CvMat *mat_I,*mat_Ix,*mat_Iy,*mat_Ixy,*mat_Ix2,*mat_Iy2;//相應的矩陣
    IplImage *pImgGray=NULL;  //灰度圖像
    IplImage *dst=NULL;    //目標圖像
    IplImage *pImgDx=NULL; //水平梯度卷積後的圖像
    IplImage *pImgDy=NULL; //豎起梯度卷積後的圖像
    IplImage *pImgDx2=NULL;//Ix2圖像
    IplImage *pImgDy2=NULL;//Iy2圖像
    IplImage *pImgDxy=NULL;//Ixy圖像

    pImgGray=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
    dst=cvCreateImage(cvGetSize(src),src->depth,3);
    pImgDx=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);//建立圖像
    pImgDy=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
    pImgDx2=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
    pImgDy2=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
    pImgDxy=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
    const int cxDIB=src->width ;         // 圖像寬度
    const int cyDIB=src->height;         // 圖像高度
    double *I=new double[cxDIB*cyDIB];
    cvCvtColor(src,pImgGray,CV_RGB2GRAY);//灰度化
    dst=cvCloneImage(src);
    int i,j;
    for(j=0;j<cyDIB;j++)
        for(i=0;i<cxDIB;i++)
        {
            I[j*cxDIB+i]=S(pImgGray,i,j);//將灰度圖像數值存入I中
        }
    mat_I=cvCreateMat(cyDIB,cxDIB,CV_64FC1);
    cvInitMatHeader(mat_I,cyDIB,cxDIB,CV_64FC1,I);//用I來初始化相應的矩陣
//    cout<<CV_MAT_ELEM(*mat_I,double,200,200)<<endl;
    //--------------------------------------------------------------------------
    //                     第一步：利用差分算子對圖像進行濾波
    //--------------------------------------------------------------------------
    //定義水平方向差分算子並求Ix
    double dx[9]={-1,0,1,-1,0,1,-1,0,1};
    mat_Ix=mbys(mat_I,cxDIB,cyDIB,dx,3,3); //求Ix矩陣
//     cout<<CV_MAT_ELEM(*mat_Ix,double,200,200)<<endl;

    //定義垂直方向差分算子並求Iy
    double dy[9]={-1,-1,-1,0,0,0,1,1,1};
    mat_Iy=mbys(mat_I,cxDIB,cyDIB,dy,3,3);//求Iy矩陣
//    cout<<CV_MAT_ELEM(*mat_Iy,double,200,200)<<endl;
    for(j=0;j<cyDIB;j++)
        for(i=0;i<cxDIB;i++)
        {
            S(pImgDx,i,j)=CV_MAT_ELEM(*mat_Ix,double,j,i);//爲相應圖像賦值
            S(pImgDy,i,j)=CV_MAT_ELEM(*mat_Iy,double,j,i);
        }

    mat_Ix2=mbxy(mat_Ix,mat_Ix,cxDIB,cyDIB);//計算Ix2,Iy2,Ixy矩陣
    mat_Iy2=mbxy(mat_Iy,mat_Iy,cxDIB,cyDIB);
    mat_Ixy=mbxy(mat_Ix,mat_Iy,cxDIB,cyDIB);
    for(j=0;j<cyDIB;j++)
        for(i=0;i<cxDIB;i++)
        {
            S(pImgDxy,i,j)=CV_MAT_ELEM(*mat_Ixy,double,j,i);//爲相應圖像賦值
            S(pImgDx2,i,j)=CV_MAT_ELEM(*mat_Ix2,double,j,i);
            S(pImgDy2,i,j)=CV_MAT_ELEM(*mat_Iy2,double,j,i);
        }
    //--------------------------------------------------------------------------
    //                  第二步：對Ix2/Iy2/Ixy進行高斯平滑，以去除噪聲
    //--------------------------------------------------------------------------
    //本例中使用5×5的高斯模板
    //計算模板參數
    //int gausswidth=5;
    //double sigma=0.8;
    double *g=new double[gausswidth*gausswidth];
    for(i=0;i<gausswidth;i++)//定義模板
        for(j=0;j<gausswidth;j++)
            g[i*gausswidth+j]=exp(-((i-int(gausswidth/2))*(i-int(gausswidth/2))+(j-int(gausswidth/2))*(j-int(gausswidth/2)))/(2*sigma));

    //歸一化：使模板參數之和爲1（其實此步能夠省略）
    double total=0;
    for(i=0;i<gausswidth*gausswidth;i++)
        total+=g[i];
    for(i=0;i<gausswidth;i++)
        for(j=0;j<gausswidth;j++)
            g[i*gausswidth+j]/=total;

    //進行高斯平滑
    mat_Ix2=mbys(mat_Ix2,cxDIB,cyDIB,g,gausswidth,gausswidth);
    mat_Iy2=mbys(mat_Iy2,cxDIB,cyDIB,g,gausswidth,gausswidth);
    mat_Ixy=mbys(mat_Ixy,cxDIB,cyDIB,g,gausswidth,gausswidth);
    
    //--------------------------------------------------------------------------
    //                        第三步：計算角點量
    //--------------------------------------------------------------------------

    //計算cim：即cornerness of image，咱們把它稱作‘角點量’
    CvMat *mat_cim;
    mat_cim=mbcim(mat_Ix2,mat_Iy2,mat_Ixy,cxDIB,cyDIB);
//     cout<<CV_MAT_ELEM(*mat_cim,double,cyDIB-1,cxDIB-1)<<endl;

    //--------------------------------------------------------------------------
    //                 第四步：進行局部非極大值抑制
    //--------------------------------------------------------------------------
    CvMat *mat_locmax;
    //const int size=7;
    mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,size);
//     cout<<CV_MAT_ELEM(*mat_locmax,double,cyDIB-1,cxDIB-1)<<endl;

    //--------------------------------------------------------------------------
    //                 第五步：得到最終角點
    //--------------------------------------------------------------------------
    CvMat *mat_corner;
    //const double threshold=4500;
    //int cornernum=0;
    mat_corner=mbcorner(mat_cim,mat_locmax,cxDIB,cyDIB,gausswidth,threshold);
    //CCommon CommonClass;
    CvPoint pt[5000];
    for(j=size/2;j<cyDIB-size/2;j++)
        for(i=size/2;i<cxDIB-size/2;i++)
        {
            if(CV_MAT_ELEM(*mat_corner,int,j,i)==255)
            {    
                pt[cornerno].x=i;
                pt[cornerno].y=j;
                cornerno++;
            //    CommonClass.DrawCross(showImg2,pt,CV_RGB(0,0,255),1,4);
            //    cvCircle(dst,pt,2,CV_RGB(255,0,0),1,8,0);
            //    cout<<i<<" "<<j<<endl;
            }
        }
    return pt;
}
C源碼

 

二、基於opencv源碼

 1 #include "imgFeat.h"
 2 
 3 void feat::detectHarrisCorners(const Mat& imgSrc, Mat& imgDst, double alpha)
 4 {
 5     Mat gray;
 6     if (imgSrc.channels() == 3)
 7     {
 8         cvtColor(imgSrc, gray, CV_BGR2GRAY);
 9     }
10     else
11     {
12         gray = imgSrc.clone();
13     }
14     gray.convertTo(gray, CV_64F);
15 
16     Mat xKernel = (Mat_<double>(1,3) << -1, 0, 1);
17     Mat yKernel = xKernel.t();
18 
19     Mat Ix,Iy;
20     filter2D(gray, Ix, CV_64F, xKernel);
21     filter2D(gray, Iy, CV_64F, yKernel);
22 
23     Mat Ix2,Iy2,Ixy;
24     Ix2 = Ix.mul(Ix);
25     Iy2 = Iy.mul(Iy);
26     Ixy = Ix.mul(Iy);
27 
28     Mat gaussKernel = getGaussianKernel(7, 1);
29     filter2D(Ix2, Ix2, CV_64F, gaussKernel);
30     filter2D(Iy2, Iy2, CV_64F, gaussKernel);
31     filter2D(Ixy, Ixy, CV_64F, gaussKernel);
32     
33 
34     Mat cornerStrength(gray.size(), gray.type());
35     for (int i = 0; i < gray.rows; i++)
36     {
37         for (int j = 0; j < gray.cols; j++)
38         {
39             double det_m = Ix2.at<double>(i,j) * Iy2.at<double>(i,j) - Ixy.at<double>(i,j) * Ixy.at<double>(i,j);
40             double trace_m = Ix2.at<double>(i,j) + Iy2.at<double>(i,j);
41             cornerStrength.at<double>(i,j) = det_m - alpha * trace_m *trace_m;
42         }
43     }
44     // threshold
45     double maxStrength;
46     minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
47     Mat dilated;
48     Mat localMax;
49     dilate(cornerStrength, dilated, Mat());
50     compare(cornerStrength, dilated, localMax, CMP_EQ);
51     
52 
53     Mat cornerMap;
54     double qualityLevel = 0.01;
55     double thresh = qualityLevel * maxStrength;
56     cornerMap = cornerStrength > thresh;
57     bitwise_and(cornerMap, localMax, cornerMap);
58     
59     imgDst = cornerMap.clone();
60     
61 }
62 
63 void feat::drawCornerOnImage(Mat& image, const Mat&binary)
64 {
65     Mat_<uchar>::const_iterator it = binary.begin<uchar>();
66     Mat_<uchar>::const_iterator itd = binary.end<uchar>();
67     for (int i = 0; it != itd; it++, i++)
68     {
69         if (*it)
70             circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1);
71     }
72 }
opencv+C++推薦

 

三、opencv中Harris接口

OpenCV的Hairrs角點檢測的函數爲cornerHairrs()，可是它的輸出是一幅浮點值圖像，浮點值越高，代表越多是特徵角點，咱們須要對圖像進行閾值化。

C++: void cornerHarris(InputArray src, OutputArray dst, int blockSize, int apertureSize, double k, int borderType = BORDER_DEFAULT);

 

  • src – 輸入的單通道8-bit或浮點圖像。
  • dst – 存儲着Harris角點響應的圖像矩陣，大小與輸入圖像大小相同，是一個浮點型矩陣。
  • blockSize – 鄰域大小。
  • apertureSize – 擴展的微分算子大。
  • k – 響應公式中的，參數$\alpha$。
  • boderType – 邊界處理的類型。
int main()
{
  Mat image = imread("../buliding.png");
  Mat gray;
  cvtColor(image, gray, CV_BGR2GRAY);

  Mat cornerStrength;
  cornerHarris(gray, cornerStrength, 3, 3, 0.01);
  threshold(cornerStrength, cornerStrength, 0.001, 255, THRESH_BINARY);
  return 0;
}

 

imageimageimage

 

 

 

 

 

 

從上面上間一幅圖像咱們能夠看到，有不少角點都是粘連在一塊兒的，咱們下面經過加入非極大值抑制來進一步去除一些粘在一塊兒的角點。

非極大值抑制原理是，在一個窗口內，若是有多個角點則用值最大的那個角點，其餘的角點都刪除，窗口大小這裏咱們用3*3，程序中經過圖像的膨脹運算來達到檢測極大值的目的，由於默認參數的膨脹運算就是用窗口內的最大值替代當前的灰度值。

int main()
{
  Mat image = imread("buliding.png");
  Mat gray;
  cvtColor(image, gray, CV_BGR2GRAY);

  Mat cornerStrength;
  cornerHarris(gray, cornerStrength, 3, 3, 0.01);

  double maxStrength;
  double minStrength;
  // 找到圖像中的最大、最小值
  minMaxLoc(cornerStrength, &minStrength, &maxStrength);

  Mat dilated;
  Mat locaMax;
  // 膨脹圖像，最找出圖像中所有的局部最大值點
  dilate(cornerStrength, dilated, Mat());
  // compare是一個邏輯比較函數，返回兩幅圖像中對應點相同的二值圖像
  compare(cornerStrength, dilated, locaMax, CMP_EQ);

  Mat cornerMap;
  double qualityLevel = 0.01;
  double th = qualityLevel*maxStrength; // 閾值計算
  threshold(cornerStrength, cornerMap, th, 255, THRESH_BINARY);
  cornerMap.convertTo(cornerMap, CV_8U);
  // 逐點的位運算黑色圖片顯示的結果
  bitwise_and(cornerMap, locaMax, cornerMap);

  drawCornerOnImage(image, cornerMap);
  namedWindow("result");
  imshow("result", image);
  waitKey();

  return 0;
}
void drawCornerOnImage(Mat& image, const Mat&binary)
{
  Mat_<uchar>::const_iterator it = binary.begin<uchar>();
  Mat_<uchar>::const_iterator itd = binary.end<uchar>();
  for (int i = 0; it != itd; it++, i++)
  {
    if (*it)
      circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1);
  }
}
View Code

 

4、Harris角點性質

 

一、閾值決定檢測點數量

增大$\alpha$的值，將減少角點響應值$R$，下降角點檢測的靈性，減小被檢測角點的數量；減少$\alpha$值，將增大角點響應值$R$，增長角點檢測的靈敏性，增長被檢測角點的數量
 

二、Harris角點檢測算子對亮度和對比度的變化不敏感

這是由於在進行Harris角點檢測時，使用了微分算子對圖像進行微分運算，而微分運算對圖像密度的拉昇或收縮和對亮度的擡高或降低不敏感。換言之，對亮度和對比度的仿射變換並不改變Harris響應的極值點出現的位置，可是，因爲閾值的選擇，可能會影響角點檢測的數量。

imageimage

 

3. Harris角點檢測算子具備旋轉不變性

Harris角點檢測算子使用的是角點附近的區域灰度二階矩矩陣。而二階矩矩陣能夠表示成一個橢圓，橢圓的長短軸正是二階矩矩陣特徵值平方根的倒數。當特徵橢圓轉動時，特徵值並不發生變化，因此判斷角點響應值也不發生變化，由此說明Harris角點檢測算子具備旋轉不變性。

4. Harris角點檢測算子不具備尺度不變性

以下圖所示，當右圖被縮小時，在檢測窗口尺寸不變的前提下，在窗口內所包含圖像的內容是徹底不一樣的。左側的圖像可能被檢測爲邊緣或曲線，而右側的圖像則可能被檢測爲一個角點。

 

image

5、函數備註

一、compare

功能：兩個數組之間或者一個數組和一個常數之間的比較

結構：

void compare(InputArray src1, InputArray src2, OutputArray dst, int cmpop)


src1 ：第一個數組或者標量，若是是數組，必須是單通道數組。
src2 ：第二個數組或者標量，若是是數組，必須是單通道數組。
dst ：輸出數組，和輸入數組有一樣的size和type=CV_8UC1
cmpop ：

標誌指明瞭元素之間的對比關係
CMP_EQ src1 相等 src2.

CMP_GT src1 大於 src2.
CMP_GE src1 大於或等於 src2.
CMP_LT src1 小於 src2.
CMP_LE src1 小於或等於 src2.
CMP_NE src1 不等於 src2.

若是對比結果爲true，那麼輸出數組對應元素的值爲255，不然爲0

//獲取角點圖
    Mat getCornerMap(double qualityLevel) {
      Mat cornerMap;
      // 根據角點響應最大值計算閾值
      thresholdvalue= qualityLevel*maxStrength;
      threshold(cornerStrength,cornerTh,
      thresholdvalue,255,cv::THRESH_BINARY);
      // 轉爲8-bit圖
      cornerTh.convertTo(cornerMap,CV_8U);
      // 和局部最大值圖與，剩下角點局部最大值圖，即：完成非最大值抑制
      bitwise_and(cornerMap,localMax,cornerMap);
      return cornerMap;
    }

二、bitwise_and（位運算函數）

功能：計算兩個數組或數組和常量之間與的關係

結構：

void bitwise_and(InputArray src1, InputArray src2, OutputArray dst, InputArray mask=noArray())


src1 ：第一個輸入的數組或常量
src2 ：第二個輸入的數組或常量
dst ：輸出數組，和輸入數組有一樣的size和type
mask ：可選擇的操做掩碼，爲8位單通道數組，指定了輸出數組哪些元素能夠被改變，哪些不能夠

操做過程爲：

\texttt{dst} (I) = \texttt{src1} (I) \wedge \texttt{src2} (I) \quad \texttt{if mask} (I) \ne0

\texttt{dst} (I) = \texttt{src1} (I) \wedge \texttt{src2} \quad \texttt{if mask} (I) \ne0

\texttt{dst} (I) = \texttt{src1} \wedge \texttt{src2} (I) \quad \texttt{if mask} (I) \ne0

若是爲多通道數組，每一個通道單獨處理

 

 

 

三、filter2D

OpenCV中的內部函數filter2D能夠直接用來作圖像卷積操做

void filter2D(InputArray src, OutputArray dst, int ddepth, InputArray kernel, Point anchor=Point(-1,-1), double delta=0, int borderType=BORDER_DEFAULT )

6、參考文章 

 

一、Opencv學習筆記（五）Harris角點檢測

二、尺度空間理論

三、圖像特徵檢測：Harris

四、圖像特徵ppt經典版

五、圖像特徵ppt經典版

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程