【数字图像处理】Hough变换C语言实现
文章目录
(一)Hough圆检测
经过图像的预处理后,我们得到了只保留边缘信息的二值图像,统计可以发现在一个320x240的图像中最后有效的信息点只有5000个左右现在就是通过对这5000个点的统计与运算完成对圆与直线的检测。直线检测最容易受到环境干扰,如果直接对一幅图像进行直线检测,检测出最长的那根直线不一定,是你想要的直线。所以我们首先要进行的是Hough圆的检测。
Hough变换不仅适用于直线检测,还适用于任何形式的f(x,a)=0所表示的图形的检测,其中x 表示坐标向量,a表示系数向量。下边我们对Hough变换检测圆的原理做简要介绍。
对于一个半径为r,圆心为(a,b)的圆,我们将其表示为:
此时x=[x,y]T,a=[a,b,r]T,其参数空间为三维。
显然,图像空间上的一点(x,y),在参数空间中对应着一个圆锥,如下图所示。
而图像空间的一个圆就对应着这一簇圆锥相交的一个点,这个特定点在参数空间的三维参数一定,就表示一定半径一定圆心坐标的图像空间的那个圆。
上述方法是经典的Hough圆检测方法的原理,它具有精度高,抗干扰能力强等优点,但由于该方法的参数空间为三维,要在三维空间上进行证据累计的话,需要的时间和空间都是庞大的,在实际应用中不适用。
程序设计
Hough变换的具体算法步骤如下:
• 适当的量化参数空间。
• 将参数空间的每一个单元看作一个累加器。
• 初始化累加器为0。
• 对图像空间的每一点,在其所满足参数方程对应的累加器上加1。
• 累加器存储的最大值即为对应的图形的参数。
程序实际思路:
(1) 因为如果直接进行Hough圆的投票的话,我们需要一个3维的表保存行,列,角度信息。我们将这个投票过程分成两个步骤。
(2) 首先完成对于固定半径的圆的投票。然后在对每一个半径的圆的投票进行比较。
1) 建立一张与原图大小一致的投票表,用来表示每个点作为圆心会有多少个点与之匹配:
int AccuArrLength=img->width * img->height;;
pAccumulateArr=(int*)malloc(AccuArrLength*sizeof(int));
2)在根据极坐标运算公式,对每一个点可能对应的圆心进行投票,通过遍历每一个点,然后对经过这个点所有的圆的圆心位置进行一次计票
int x = i - R * cos(theta);//得到理想圆心x坐标
int y = j - R * sin(theta);//得到理想圆心y坐标
if(x>0&&x<img->width&&y>0&&y<img->height&&x<range[1]&&x>range[0]&&y<range[3]&&y>range[2])
{
pAccumulateArr[y * img->width + x]++;
}
3) 找到对应票数最高的圆心位置
for(int i = 0; i < img->height; i++)
for(int j = 0; j < img->width; j++)
{
int value_ = pAccumulateArr[i * img->width + j];
if(value_>max_value)
{
//printf("(%d,%d,%d,%d)",max_value,value_,i,j);
max_value=value_;
circles->x=j;
circles->y=i;
circles->r=R;
}
}
然后我们就可以找到想要的圆了。
(3) 对一个半径范围内的所有圆进行一次(1)中的操作,然后,从中选取票数最高的圆,那么这个圆就是我们想要的圆了。
(4) 具体程序请参考Litecv/Imgproc/li_canny.c
(二) Hough直线检测
Hough变换的方法基本思想可以从检测图像中的直线这个简单问题中看到。直线由两点A=(X1,Y1)和B=(X2,Y2)定义,所下图1(a)示。通过点A的所有直线由y1=kx1+q表示,k和q是某些值。这意味着同一个方程可以解释为参数空间k,q的方程。因此通过点A的所须直线可以表示为方程q=-x1k+y1图1.(b)。类似地通过点B的直线可以表示q=-x2*k+y2。在参数空间k和q中,两条直线的唯一公共点是在原图像空间中表示连接点A和B的唯一存在的直线。
这意味着图像中的每条直线在参数空k,q中由单独一个点表示,直线的任何一部分都变换为同一个点。直线检测的主要思想是确定图中所在的直线像素,将通过这些像素的所在直线变换到参数空间的对应点,在参数空间检测点(a,b),此点是图像中出现的直线y=ax+b的Hough变换的结果。
我们可以注意到直线的参数方程y =kx+q只适合解释Hough变换原理,在检测垂直线条和参数的非线性离散化时会遇到困难。如果直线表示成s=xcosθ+ysinθ。Hough变换就没有这些局限性。直线还是被变换为单个点,因此可用该原理进行直线检测。如下图2所示:
我们要注意到Hough变换的重要性质是对图像中直线的殘缺部分、噪声以及其它共存的非直线结构不敏感。因为从图像空间到累计空间的变换的鲁棒性引起的,直线殘缺的部分只会造成较低的局部极值。
(三)代码实现
#ifndef LI_CANNY_C
#define LI_CANNY_C
#include "cv.h"
#include "li_image_proc.h"
#include <stdio.h>
#include <math.h>
/**
* @name: Li_Canny
* @msg: 参考文章 https://blog.csdn.net/HUSTER_Gy/article/details/102942452?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522160498444419724838560446%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fall.%2522%257D&request_id=160498444419724838560446&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~first_rank_v2~rank_v28-11-102942452.first_rank_ecpm_v3_pc_rank_v2&utm_term=canny%E8%BE%B9%E7%BC%98%E6%A3%80%E6%B5%8B%E7%AE%97%E6%B3%95c%E5%AE%9E%E7%8E%B0&spm=1018.2118.3001.4449
* 图像砍尼检测
* @param {Li_Image* img 原图像
* BYTE CannyType选择算子
* BYTE min 最大阈值
* BYTE max} 最小阈值
* @return {*}
*/
LI_API
Li_Image* Li_Canny(Li_Image* img,BYTE CannyType,BYTE min,BYTE max)
{
if(img==NULL||img->imgdepth!=LI_DEP_8U)return NULL;
LILOG("CANNY");
Li_Image* out =Li_Copy_Image(img);
Li_Image* PP =Li_Copy_Image(img);
Li_Image* QQ =Li_Copy_Image(img);
Li_Kernel* SX,*SY;
double *P =(double*)malloc(sizeof(double)*img->width*img->height);/*x方向偏导*/
double *Q =(double*)malloc(sizeof(double)*img->width*img->height); /*y方向偏导*/
/*PI*-1/2* -- PI*1/2*/
double *Threa =(double*)malloc(sizeof(double)*img->width*img->height);
#ifdef DEBUG
Li_Save_Image("before_minus.bmp",out);
#endif
/*开始计算梯度与方向角*/
switch (CannyType)
{
case LI_CANNY_MYDEFINE:
{
for(int i=0;i<img->height-1;i++)
for(int j=0;j<img->width-1;j++)
{
BYTE* ptr[4];
BYTE* ptr2;
/**
* 2 3
* 0 1
*/
ptr[0]=img->at(img,j,i);
ptr[1]=img->at(img,j+1,i);
ptr[2]=img->at(img,j,i+1);
ptr[3]=img->at(img,j+1,i+1);
P[i*img->width+j]=(double)((double)*ptr[3]+*ptr[1]-*ptr[0]-*ptr[2])/2;
Q[i*img->width+j]=(double)((double)*ptr[0]+*ptr[1]-*ptr[2]-*ptr[3])/2;
Threa[i*img->width+j]=atan(Q[i*img->width+j]/P[i*img->width+j]);
ptr2=out->at(out,j,i);
*ptr2=sqrt(P[i*img->width+j]*P[i*img->width+j]+Q[i*img->width+j]*Q[i*img->width+j]);
}
}
break;
case LI_CANNY_SOBEL:
{
for(int i=1;i<img->height-1;i++)
for(int j=1;j<img->width-1;j++)
{
BYTE* ptr[9];
BYTE* ptr2;
/**6 7 8
* 3 4 5
* 0 1 2
*/
if(j-1>=0&&i-1>=0)
ptr[0]=(BYTE*)img->at(img,j-1,i-1);
if(j>=0&&i-1>=0)
ptr[1]=(BYTE*)img->at(img,j+0,i-1);
if(j+1<=img->width&&i-1>=0)
ptr[2]=(BYTE*)img->at(img,j+1,i-1);
if(j-1>=0&&i>=0)
ptr[3]=(BYTE*)img->at(img,j-1,i+0);
if(j>=0&&i>=0)
ptr[4]=(BYTE*)img->at(img,j+0,i+0);
if(j+1<=img->width&&i>=0)
ptr[5]=(BYTE*)img->at(img,j+1,i+0);
if(j-1>=0&&i+1<=img->height)
ptr[6]=(BYTE*)img->at(img,j-1,i+1);
if(j>=0&&i+1<=img->height)
ptr[7]=(BYTE*)img->at(img,j+0,i+1);
if(j+1<=img->width&&i+1<=img->height)
ptr[8]=(BYTE*)img->at(img,j+1,i+1);
P[i*img->width+j]=(double)((double)*ptr[2]+*ptr[5]*2+*ptr[8]-*ptr[0]-*ptr[3]*2-*ptr[6]);
Q[i*img->width+j]=(double)((double)*ptr[6]+*ptr[7]*2+*ptr[8]-*ptr[0]-*ptr[1]*2-*ptr[2]);
Threa[i*img->width+j]=atan(Q[i*img->width+j]/P[i*img->width+j]);
ptr2=out->at(out,j,i);
*ptr2=sqrt(P[i*img->width+j]*P[i*img->width+j]+Q[i*img->width+j]*Q[i*img->width+j]);
}
}
break;
case LI_CANNY_ROBERTS:
{
for(int i=0;i<img->height-1;i++)
for(int j=0;j<img->width-1;j++)
{
BYTE* ptr[4];
BYTE* ptr2;
/**
* 2 3
* 0 1
*/
ptr[0]=img->at(img,j,i);
ptr[1]=img->at(img,j+1,i);
ptr[2]=img->at(img,j,i+1);
ptr[3]=img->at(img,j+1,i+1);
ptr2=out->at(out,j,i);
*ptr2=abs(*ptr[2]-*ptr[1])+abs(*ptr[0]-*ptr[3]);
Threa[i*img->width+j]=atan(abs(*ptr[2]-*ptr[1])/abs(*ptr[0]-*ptr[3]));
}
}
break;
case LI_CANNY_PREWITT:
{
for(int i=1;i<img->height-1;i++)
for(int j=1;j<img->width-1;j++)
{
BYTE* ptr[9];
BYTE* ptr2;
/**6 7 8
* 3 4 5
* 0 1 2
*/
if(j-1>=0&&i-1>=0)
ptr[0]=(BYTE*)img->at(img,j-1,i-1);
if(j>=0&&i-1>=0)
ptr[1]=(BYTE*)img->at(img,j+0,i-1);
if(j+1<=img->width&&i-1>=0)
ptr[2]=(BYTE*)img->at(img,j+1,i-1);
if(j-1>=0&&i>=0)
ptr[3]=(BYTE*)img->at(img,j-1,i+0);
if(j>=0&&i>=0)
ptr[4]=(BYTE*)img->at(img,j+0,i+0);
if(j+1<=img->width&&i>=0)
ptr[5]=(BYTE*)img->at(img,j+1,i+0);
if(j-1>=0&&i+1<=img->height)
ptr[6]=(BYTE*)img->at(img,j-1,i+1);
if(j>=0&&i+1<=img->height)
ptr[7]=(BYTE*)img->at(img,j+0,i+1);
if(j+1<=img->width&&i+1<=img->height)
ptr[8]=(BYTE*)img->at(img,j+1,i+1);
P[i*img->width+j]=(double)((double)*ptr[2]+*ptr[5]+*ptr[8]-*ptr[0]-*ptr[3]-*ptr[6]);
Q[i*img->width+j]=(double)((double)*ptr[6]+*ptr[7]+*ptr[8]-*ptr[0]-*ptr[1]-*ptr[2]);
Threa[i*img->width+j]=atan(Q[i*img->width+j]/P[i*img->width+j]);
ptr2=out->at(out,j,i);
*ptr2=sqrt(P[i*img->width+j]*P[i*img->width+j]+Q[i*img->width+j]*Q[i*img->width+j]);
}
}
break;
default:
break;
}
#ifdef DEBUG
Li_Save_Image("after_minus.bmp",out);
#endif
/*非极大值抑制*/
for(int j=1;j<out->height-1;j++)
for(int i=1;i<out->width-1;i++)
{
double t=Threa[j*img->width+i];
BYTE* ptr=out->at(out,i,j);
double g=(double) *ptr;
double g0, g1;
if ((t >= -(3*M_PI/8)) && (t < -(M_PI/8)))
{
ptr=out->at(out,i-1,j-1);
g0=(double) *ptr;
ptr=out->at(out,i+1,j+1);
g1=(double) *ptr;
}
else if ((t >= -(M_PI/8)) && (t < M_PI/8))
{
ptr=out->at(out,i-1,j);
g0=(double) *ptr;
ptr=out->at(out,i+1,j);
g1=(double) *ptr;
}
else if ((t >= M_PI/8) && (t < 3*M_PI/8))
{
ptr=out->at(out,i+1,j-1);
g0=(double) *ptr;
ptr=out->at(out,i-1,j+1);
g1=(double) *ptr;
}
else
{
ptr=out->at(out,i,j-1);
g0=(double) *ptr;
ptr=out->at(out,i,j+1);
g1=(double) *ptr;
}
if (g <= g0 || g <= g1) {
ptr=out->at(out,i,j);
*ptr=0;
}
}
/*阈值化操作*/
#ifdef DEBUG
Li_Save_Image("before_thre.bmp",out);
#endif
Li_Image*out1=Li_Double_Threshold(out,min,max);
#ifdef DEBUG
Li_Save_Image("after_thre.bmp",out1);
#endif
/*边缘链接*/
for (int j = 1; j < out1->height-2; j++)
for (int i = 1; i < out1->width-2; i++) {
BYTE* ptr=out1->at(out1,i,j);
if(*ptr==255)
{
for (int m = -1; m < 1; m++) {
for (int n = -1; n < 1; n++) {
BYTE* temp=out1->at(out1,i+n,j+m);
if(*ptr!=0&&*ptr!=255)
*ptr=255;
}
}
}
}
for (int j = 0; j < out1->height-1; j++) {
for (int i = 0; i < out1->width-1; i++) {
// 如果该点依旧是弱边缘点,及此点是孤立边缘点
BYTE* ptr=out1->at(out1,i,j);
if(*ptr!=255&&*ptr!=0)
*ptr=0;
}
}
Li_Destroy_Image(PP);
Li_Destroy_Image(QQ);
return out1;
}
#endif // !LI_CANNY_C
(四)效果展示
(五)写在后面
因为LiteCV项目才刚刚写了一个开头,代码中有错误的地方还望指出。我已经将项目同步到了github,我会实时更新这个代码仓库。
项目github地址:
LiteCV
下一篇: Java实现模拟五子棋游戏