opencv学习_13 (trajkovic 角点检测)

来源:http://blog.csdn.net/songzitea/article/details/13614977

背景引言

本节主要内容来源于是由 Miroslav Trajkovic和Mark Hedley[1]在1998年提出Trajkovic算子，其论文为FastCorner Detection.和Trajkovic Operator (4-Neighbours)[3](注：本节图片主要是来源于此)。Trajkovic算子角点提取方法存在价值是在角点提取的效果上，它优于同时期的其他角点(如Moravec角点，Harris角点)提取方法，同时，从算法的运行速度角度而言，它比同是期的其他角点提取方法要快很多。

基本理论

作者对角点的定义，与Moravec角点、 Harris 角点的定义是一样的：图像灰度值在各个方向变化都比较大的点，即认为是角点。角点量的定义也类似于Moravec,即角点量的是在各个方向上灰度变化的最小值（当在某个方向上，灰度值的变化最小，并且这个最小值也大于某某个设定的阈值，那么认为这个点就是一个角点）。与Moravec 相比 ，Trajkovic 的性能更好，因为Trajkovic通过使用像素插值的方法比较每一个方向上的灰度值变化程度，而moravec只有限个方向的比较，与harris相比，两都效果差不多，但是Trajkovic的速度远快于haris.

Trajkovic在速度方面的表现于决定它能够应用对实时要求比较高的系统中。但是Trajkovic也有缺点是它不具备旋转不变性，对噪声比较敏感，对对角点的响应值比较大等。对于这些缺点也会有相应的方法，但不能完全的解决这些缺点。Trajkovic算子计算角点量是在一个小的圆形窗口内，并且考虑了所有通过圆心的直线，即计算了每一个可能的方向上的灰度值变化程度）。假设用C表示圆的圆心，那么通过圆心C的任意一条直线会与圆在两个交点，如图1所示：

Figure 1: Notation for Trajkovic Operator

那么，Trajkovic对角点量的定义是

C(x,y)表示的是图像I上，任意一点(x,y)处的角点量。这个公式能断出是否是角点。点(x,y)在图像目的位置有四种情况。如下图所示：

第一情况：点(x,y)在内部区域：此时经点(x,y)位中心的圆形窗中完全在一个比较平坦的区域，或者大部分的圆形窗口是在这个平坦区域的（如上图A所示）。这种情况下，至少存大一条直线过中心点C的直线使得Ic约等于Ip和Ic约等于Ip‘,根据角点量计算公式，计算出的C(x,y)一定是一个很小的值，那么就可以排除点(x,y)不是角点。并且，大多数的情况下，会有多个这样的点P存大，使得对噪声一定的承受能力。

第二情况：点(x,y)在边缘上：此时圆形窗口的中心点(x,y)刚好落在边缘上（如上图B所示）。仅存在一个点对P,P'使得Ic约等于Ip和Ic约等于Ip‘。在理想状况下，根据角点量的公式计算出来的角点量C(x,y)是比较小的。但是由于只存在一个条这样直线，所以容易受到噪声的干扰。

第三情况：点(x,y)在角点上：当圆窗口的中心(x, y)位于角点上时（如上图C 所示）。任意一条通过点(x,y)与圆相交的交点P、P', P、P' 中至少有一个点的灰度值与中心的点的灰度值差别比较大。所以计算出来的角点量也比较大（大于设定的阈值，所以会被认为是角点）。

第四情况：点(x,y)在一个孤立点上：如果在孤立点上（如上图D所示），对于每一个直线上的Ip, Ip’,都分与Ic的值相差比较大，此时计算出的角点量也比较大。此时这个点是噪声,但是由于角点量比较大，会被认为是角点。为了削除这种噪声的影响，可以先使用高斯平滑进行去噪。

根据上面的分析，Trajkovic 受噪声的影响比较大，所以，可以先实现高斯平滑去噪再使用Trajkovic算子进行角点检测。

角点量的计算方法

目前，关于Trajkovic计算角点量的并未具体介绍如何时计算。首先，如下图2所示，水平方向的灰度值变化量rA,垂直方向的灰度值变化量rB很容易计算。

如何计算任意一个方向上灰度值的变化程度呢？首先，我们角点量就取其最小值即可：

 

Figure 2: Interpixel positions               Figure3: Interpixel approximation for a 3x3 window using 4-neighbours

现在，我们可将问题可以转换为，任意一条通过点C 的直线与圆的交点是P，P',Q,Q'即可：

把公式(2)代入公式(1) 中即可：

那么，最小值是即可结论是

多格算法（Multigrid Algorithm）

角点可以分为两类：几何角点和纹理角点。几何角点是由图像中物体边缘的相交。而纹理角点有由物体表面的纹理产生的角点（如：草地，衣服的纹理等）。通常情况下，一幅图像中的几何角点的数量要远少于纹理角点。多格算法的目的是希望能多检测到几何角点，少检测到纹理点。因Trajkovic 认为几何角点比纹理角点更加的稳定。所以，从实践使用来看，减少纹理角点是合理的。

经过观察发现，纹理角点一般都是非常的密集，并且是在一个很小区域内灰度值发生变化，所以，采用把原图缩小以后，再来提取角点，通过缩小图像，可以消除区域内灰度值的变换。缩小图像的方法不是使用插值法，而是使用平均法。如下图所示：

Trajkovic角点检测，首先使用多格法对原如图像进行缩小，这样不仅能够减少检测到纹理角点数目，而且加快了Trajkovic角点检测的速度，在Trajkovic角点检测缩小版的图像上，使用角为简单的角点量计算公式初步判断是不是侯选角点。如果是，再到原始图像上进一步判定。

Figure 4: Corner points detected at different image resolutions

如图4所示，原始的大小是256*256,图像中有一片草地，而草地上的角点，对于我们进行角点匹配时，是没有什么意义。当把图像缩小到128*128时，草地上的角点已经消除了。当然，图像也不能缩小太多，否则，几何角点也会消失，当缩小到64*64时，部分几何角点也消失了。

算法流程

代码：

[cpp] view plain copy

 print?

1. #include <iostream>  
2. #include "cv.h"  
3. #include "highgui.h"  
4. #include "cxcore.h"  
5. using namespace std;  
6.   
7. void getTrajkovic4NCorner(IplImage* src,IplImage* srcResize , float T1 , float T2 ,CvSeq* corners,int maximumSize)  
8. {  
9.     int x,y,maxChar=255,scaleX = src->width/srcResize->width,scaleY = src->height/srcResize->height;  
10.       
11.       
12.     IplImage* srcResizeMap ,*srcMap;  
13.       
14.     srcResizeMap = cvCreateImage(cvGetSize(srcResize),8,1);     // 用来保存缩小图像的角点量  
15.     srcMap = cvCreateImage(cvGetSize(src),32,1);                // 用来保存源图像计算得到的角点量  
16.       
17.     cvZero(srcMap);  
18.     cvZero(srcResizeMap);  
19.   
20.     for( y=1;y<srcResize->height-1;y++)  
21.     {  
22.         uchar* preRow = (uchar*)(srcResize->imageData + (y-1)*srcResize->widthStep);  
23.         uchar* curRow = (uchar*)(srcResize->imageData + y*srcResize->widthStep);  
24.         uchar* nextRow = (uchar*)(srcResize->imageData + (y+1)*srcResize->widthStep);  
25.           
26.         uchar* MapData = (uchar*)(srcResizeMap->imageData + y*srcResizeMap->widthStep);  
27.       
28.         for(x=1;x<srcResize->width-1;x++)  
29.         {  
30.               
31.             int IC,IA,IB,IAA,IBB;  
32.             int rA,rB,C_Simple;  
33.               
34.             // curRow[x]是uchar型,而IC是int型,如果直接赋值int类型的其它三个字节可能会产生随机值  
35.             IC = curRow[x]&maxChar;  
36.               
37.             IA = curRow[x+1]&maxChar;  
38.             IAA = curRow[x-1]&maxChar;  
39.   
40.             IB = preRow[x]&maxChar;  
41.             IBB = nextRow[x]&maxChar;  
42.               
43.             rA = (IA-IC)*(IA-IC) + (IAA-IC)*(IAA-IC);  
44.             rB = (IB-IC)*(IB-IC) + (IBB-IC)*(IBB-IC);  
45.               
46.             C_Simple = rA < rB ? rA : rB;  
47.           
48.             if(C_Simple > T1)  
49.             {  
50.                 MapData[x]=1;  
51.             }  
52.         }  
53.       
54.     }  
55.   
56.     for( y=1;y<srcResize->height-1;y++)  
57.     {  
58.           
59.         uchar* srcResizeMapData = (uchar*)(srcResizeMap->imageData + y*srcResizeMap->widthStep);  
60.         //uchar* srcMapData = (uchar*)(srcMap->imageData + y*srcMap->widthStep);  
61.   
62.         for(x=1;x<srcResize->width-1;x++)  
63.         {  
64.               
65.             if(srcResizeMapData[x]==0)  
66.             {  
67.                 continue;  
68.             }  
69.               
70.             int originX,originY,IC,IA,IB,IAA,IBB;  
71.             int rA,rB,C_Simple;  
72.               
73.             originX = x*scaleX;  
74.             originY = y*scaleY;  
75.               
76.             uchar* srcPreRow = (uchar*)(src->imageData + (originY-1)*src->widthStep);  
77.             uchar* srcCurRow = (uchar*)(src->imageData + originY*src->widthStep);  
78.             uchar* srcNextRow = (uchar*)(src->imageData + (originY+1)*src->widthStep);  
79.               
80.             IC = srcCurRow[originX]&maxChar;  
81.               
82.             IA = srcCurRow[originX+1]&maxChar;  
83.             IAA = srcCurRow[originX-1]&maxChar;  
84.   
85.             IB = srcPreRow[x]&maxChar;  
86.             IBB = srcNextRow[x]&maxChar;  
87.   
88.             rA = (IA-IC)*(IA-IC) + (IAA-IC)*(IAA-IC);  
89.             rB = (IB-IC)*(IB-IC) + (IBB-IC)*(IBB-IC);  
90.               
91.             C_Simple = rA < rB ? rA : rB;  
92.               
93.             if(C_Simple>T2)  
94.             {  
95.                 float B1,B2,C,A,B,C_InterPixel;  
96.                   
97.                 B1 = (IB-IA)*(IA-IC)+(IBB-IAA)*(IAA-IC);  
98.                 B2 = (IB-IAA)*(IAA-IC)+(IBB-IA)*(IA-IC);  
99.                 C = rA;  
100.                 B = B1<B2 ? B1 : B2;  
101.                 A = rB-rA-2*B;  
102.                   
103.                 if(B<0 &&(A+B)>0)  
104.                 {  
105.                     C_InterPixel = C-(B*B)/A;  
106.                 }  
107.                 else  
108.                 {  
109.                     C_InterPixel=C_Simple;  
110.                 }  
111.                   
112.                 if(C_InterPixel>T2)  
113.                 {  
114.   
115.                     float* srcMapData = (float*)(srcMap->imageData + originY*srcMap->widthStep);  
116.                     srcMapData[originX]=C_InterPixel;  
117.                     //cvSetReal2D(srcMap,originY,originX,C_InterPixel);  
118.                 }  
119.             }  
120.               
121.         }  
122.   
123.     }  
124.   
125.   
126.     //计算局部极大值 及 极大值是否大于阈值  
127.     int beginY,endY,beginX,endX;  
128.     int halfWinSize = maximumSize/2;  
129.   
130.     beginY = halfWinSize;  
131.     endY = srcMap->height - halfWinSize;  
132.   
133.     beginX = halfWinSize;  
134.     endX = srcMap->width - halfWinSize;  
135.       
136.     for(y=beginY;y<endY;)  
137.     {  
138.         for(x=beginX;x<endX;)  
139.         {  
140.             //寻找局部极大值 及其位置信息  
141.             float maxValue=0;  
142.             int flag = 0 ;  
143.             CvPoint maxLoc;  
144.             maxLoc.x = -1;  
145.             maxLoc.y = -1;  
146.   
147.   
148.             //首先计算以点(x,y)位中心的maximumSize*maximumSize的窗口内部的局部极大值  
149.             for(int winy=-halfWinSize;winy<=halfWinSize;winy++)  
150.             {  
151.                 for(int winx=-halfWinSize;winx<=halfWinSize;winx++)  
152.                 {  
153.                     float value ;  
154.                     value = cvGetReal2D(srcMap,y+winy,x+winx);  
155.                       
156.                     //计算该窗口内 最大值 保存到max 并保存其坐标到maxLoc  
157.                     if(value>maxValue)  
158.                     {  
159.                         maxValue = value;  
160.                         maxLoc.x = x+winx;  
161.                         maxLoc.y = y+winy;  
162.                         flag = 1;  
163.                     }  
164.                 }  
165.             }  
166.               
167.   
168.   
169.               
170.             //如果找到局部极大值 并且该值大于预先设定的阈值 则认为是角点  
171.             if(flag==1 && maxValue>T2)  
172.             {  
173.                 cvSeqPush(corners,&maxLoc);  
174.               
175.             }  
176.               
177.           
178.             x = x+halfWinSize;  
179.   
180.         }  
181.           
182.         y = y + halfWinSize;  
183.     }  
184.   
185.     cvReleaseImage(&srcResizeMap);  
186.     cvReleaseImage(&srcMap);  
187. }  
188.   
189. int main(int argc,char* argv[])  
190. {  
191.     //相关变量  
192.     int scale = 2;  
193.     IplImage* src,*srcGray,*srcGrayResize;  
194.     CvMemStorage* mem = cvCreateMemStorage(0);  
195.     CvSeq* TrajkovicPoints;  
196.   
197.     src = cvLoadImage("E:\\study_opencv_video\\lesson17_2\\2.jpg");//源图像  
198.       
199.     srcGray = cvCreateImage(cvGetSize(src),8,1);  
200.       
201.     if(!src)  
202.     {  
203.         cout << "src is null" << endl;  
204.         return 0;  
205.     }  
206.       
207.     cvCvtColor(src,srcGray,CV_BGR2GRAY);  
208.     srcGrayResize = cvCreateImage(cvSize(srcGray->width/scale,srcGray->height/scale),8,1);  
209.     cvResize(srcGray,srcGrayResize);//请将resize修改为多格算法  
210.   
211.     //Trajkovic 4 邻域角点角点保存的空间  角点坐标保存在一个序列中  
212.     TrajkovicPoints = cvCreateSeq(0,sizeof(CvSeq),sizeof(CvPoint),mem);  
213.       
214.   
215.     //进行Trajkovic 4 邻域角点检测  
216.     int T1=120,T2=150,localArea = 8;   // localArea 是局部极大值抑制时窗口的大小  
217.     getTrajkovic4NCorner(srcGray,srcGrayResize,T1,T2,TrajkovicPoints,localArea);  
218.       
219.     //获取每一个角点的坐标  
220.     for(int x=0;x<TrajkovicPoints->total;x++)  
221.     {  
222.         //获取第x个角点的坐标  
223.         CvPoint* pt = (CvPoint*)cvGetSeqElem(TrajkovicPoints,x);  
224.   
225.         //以角点坐标为中心  绘制一个半径为5的圆  
226.         cvCircle(src,*pt,2,cvScalar(255,0,255,0));  
227.     }  
228.   
229.   
230.     //显示图像  
231.     cvNamedWindow("dst");  
232.     cvShowImage("dst",src);  
233.     cvWaitKey(0);  
234.   
235.     //释放资源  
236.     cvReleaseImage(&src);  
237.     cvReleaseImage(&srcGray);  
238.     cvReleaseImage(&srcGrayResize);  
239.     cvReleaseMemStorage(&mem);  
240.     return 0;  
241.   
242. }  

英文链接：http://kiwi.cs.dal.ca/~dparks/CornerDetection/trajkovic.htm