canny算法实现
关于canny算法的原理和实现 这篇博客有详细的介绍 http://blog.csdn.net/likezhaobin/article/details/6892629
下面是本人源码,和opencv的cvcanny的不同之处
1.有高斯滤波 , cvcanny 无。
2. 上下阈值自动设定,opencv手动设定。
// createGuassFilter.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include <cv.h>
#include <highgui.h>
#include <cxcore.h>
#ifdef DEBUG
#pragma comment(lib,"cxcore200d.lib")
#pragma comment(lib,"cv200d.lib")
#pragma comment(lib,"highgui200d.lib")
#else
#pragma comment(lib,"cxcore200.lib")
#pragma comment(lib,"cv200.lib")
#pragma comment(lib,"highgui200.lib")
#endif
void creatGuassFilter(float *GuassFilter, int nsize,float sigma=0);
void calGradientAndAmplitue(IplImage * src , IplImage * gradientX,IplImage * gradientY ,IplImage * amplitue );
void noMaxSuppress(IplImage * gradientX, IplImage *gradientY, IplImage *amplitue,IplImage * noMaxSupress );
void EstimateThreshold(IplImage * noMaxSuppressImg , IplImage * amplitue , double &lowthreahold ,double &highthreahold);
void hysteresis(IplImage * amplitue , IplImage * noMaxSuppressImg , float highthreshold , float lowthreshold);
void trackEdge(int x , int y , int lowthreshold, IplImage * amplitue , IplImage * noMaxSuppressImg);
void canny(IplImage * src, IplImage * dst ,float sigma , int size , float highthreshold, float lowthreshold);
int _tmain(int argc, _TCHAR* argv[])
{
IplImage * im = cvLoadImage("e.bmp",CV_LOAD_IMAGE_GRAYSCALE);
IplImage * im2 = cvCreateImage(cvSize(im->width/10,im->height/10),8,1);
cvResize(im,im2);
IplImage * dst = cvCreateImage(cvGetSize(im2),8,1);
canny(im2,dst,3,7,0,0);
cvWaitKey(-1);
return 0;
}
/*
函数名 : creatGuassFilter
功能:生成高斯滤波模板
参数:mat_GuassFilter --- out 用于存储生成的高斯滤波模板
nsize ----- in 滤波窗口尺寸
sigma ----- in 高斯半径
*/
void createGuassFilter(float *GuassFilter, int nsize,float sigma)
{
float distX,distY; //到窗口中心的距离
int ncenter = nsize/2; //窗口中心
float sigmaA = 1.0 / (sqrt(2 * CV_PI) * sigma) ;
float sum =0.0;
//如果sigma = 0 ,自动计算sigma值
if (0 == sigma)
{
sigma = (nsize /2 -1) *0.3 + 0.8;
}
for (int i = 0 ; i < nsize ; ++i)
{
for (int j = 0 ; j < nsize ; ++j)
{
distX = i - ncenter;
distY = j - ncenter;
*(GuassFilter + i * nsize + j ) = sigmaA * exp( -1.0 * ( distX * distX + distY * distY) / (2 * sigma * sigma)) ;
sum +=*(GuassFilter + i * nsize + j ) ;
}
}
//归一化
for (int i = 0 ; i <nsize ; ++i)
{
for (int j=0; j <nsize; ++j)
{
*(GuassFilter + i * nsize + j ) /=sum;
}
}
}
/*
函数:calGradientAndAmplitue
功能:计算图像梯度方向 和幅度
参数:src -- in 原图
gradientX Y -- out x,y方向梯度值 像素类型 int
amplitue -- out 幅度值 像素类型double
*/
void calGradientAndAmplitue(IplImage * src , IplImage * gradientX,IplImage * gradientY ,IplImage * amplitue )
{
//计算 x ,y 方向梯度
float filterx[3][3] = {
-1 , 0 , 1 ,
-2 , 0 , 2,
-1 , 0 , 1
};
float filtery[3][3] = {
1 , 2 , 1,
0 , 0 , 0,
-1 ,-2, -1
};
CvMat mat_filterx ;
CvMat mat_filtery;
cvInitMatHeader(&mat_filterx , 3, 3,CV_32FC1 ,filterx);
cvInitMatHeader(&mat_filtery , 3 ,3,CV_32FC1, filtery);
cvFilter2D(src,gradientX,&mat_filterx);
cvFilter2D(src,gradientY,&mat_filtery);
//计算幅度值
int nwidth = src->width ;
int nheight = src->height;
float Gx, Gy; //x y 方向梯度
float * pGx , *pGy ,*pAmp;
for (int i = 0 ; i < nheight ; ++ i)
{
pGx = (float *)(gradientX->imageData + i* gradientX->widthStep );
pGy = (float *)(gradientY->imageData + i* gradientX->widthStep );
pAmp = (float *)(amplitue->imageData + i * amplitue->widthStep);
for (int j = 0 ;j < nwidth ; ++j)
{
Gx = *(pGx + j);
Gy = *(pGy + j);
*(pAmp + j) = sqrt(1.0 * Gx * Gx + 1.0 * Gy * Gy);
}
}
}
/*
函数:noMaxSuppress
功能:非极大值抑制
*/
void noMaxSuppress(IplImage * gradientX, IplImage *gradientY, IplImage *amplitue,IplImage * noMaxSupress )
{
int nwidth = noMaxSupress->width ;
int nheight = noMaxSupress->height;
//非极大值抑制
float Gx ,Gy ; // x,y梯度值
float temp1,temp2,temp3,temp4; //用于计算插值得到梯度方向邻域幅度值的中间变量
float amp , amp1 ,amp2; // 该点幅度值 该点梯度方向上的幅度值
float angle = 0.0; //梯度方向
float rate = 0.0; // 线插比例
float *pAmp , *pnoMaxSuppress,*pgradientX,*pgradientY;
pAmp = (float *)amplitue->imageData;
pnoMaxSuppress = (float *)noMaxSupress->imageData;
pgradientX = (float *)gradientX->imageData;
pgradientY = (float *)gradientY->imageData;
for ( int i = 1 ; i <nheight -1 ; ++i)
{
for (int j = 1 ; j < nwidth -1 ; ++j)
{
amp =*( pAmp+ i *nwidth + j);
//如果幅度值为0 则为非边缘点
if (amp ==0)
{
*(pnoMaxSuppress + i * nwidth + j) = 0 ;
continue;
}
Gx = *(pgradientX + i * nwidth + j);
Gy = *(pgradientY + i *nwidth + j);
angle = atan(1.0 * Gy / Gx);
//第一种情况
// temp1
//temp3 amp temp2
//temp4
if ( (angle >=0 && angle < CV_PI/4) /* || (angle >=CV_PI && angle < 5.0/4 *CV_PI)*/ )
{
temp1 = *(pAmp + (i-1)*nwidth + j+1);
temp2 = *(pAmp + i *nwidth + j+1);
temp4 = *(pAmp + i *nwidth + j-1);
temp4 = *(pAmp + (i +1) * nwidth + j-1);
rate = fabs(1.0 *Gy)/fabs(1.0*Gx);
amp1 = temp1 * rate + temp2 * (1 -rate);
amp2 = temp3 * rate + temp4 * (1 -rate);
}
//第二种情况
// temp1 temp2
// amp
// temp3 temp4
else if ( ( angle >=CV_PI/4 && angle <= CV_PI/2)/* || (angle >=5.0/4 * CV_PI && angle <3.0/2 * CV_PI )*/)
{
temp1 = *(pAmp + (i-1) * nwidth + j );
temp2 =*(pAmp + (i-1) * nwidth + j +1);
temp3 = *(pAmp + (i +1)*nwidth + j -1);
temp4 = *(pAmp + (i +1) * nwidth + j);
rate = fabs(1.0*Gx)/fabs(1.0*Gy);
amp1 = temp2 * rate + temp1 *(1- rate);
amp2 = temp3 * rate + temp4 *(1-rate);
}
//第三种情况
// temp1 temp2
// amp
// temp3 temp4
else if( ( angle >= -1.0 * CV_PI/2 && angle < -1.0 *CV_PI /4))
{
temp1 = *(pAmp + (i-1)*nwidth + j-1);
temp2 = *(pAmp + (i-1)*nwidth + j);
temp3 =*(pAmp + (i+1)*nwidth + j);
temp4 = *(pAmp + (i +1)*nwidth +j+1);
rate = fabs(1.0*Gx) / fabs(1.0*Gy);
amp1 = temp1 * rate + temp2 *(1 - rate);
amp2 = temp4 * rate + temp3 * (1- rate);
}
//第四种情况
// temp1
// temp2 amp temp3
// temp4
else /*if ( angle >= -1.0 * CV_PI /4 && angle <0)*/
{
temp1 = *(pAmp + (i-1)*nwidth + j-1);
temp2 = *(pAmp + i * nwidth + j-1);
temp3 = *(pAmp + i * nwidth + j +1);
temp4 =*(pAmp + (i+1)*nwidth + j +1);
rate = fabs(1.0*Gy)/fabs(1.0*Gx);
amp1 = temp1 * rate + temp2 *( 1- rate);
amp2 = temp4 * rate + temp3* (1- rate);
}
//局部最大值判断 如果为最大值 128 否则0
if( amp >=amp1 && amp >= amp2)
*(pnoMaxSuppress + i * nwidth + j)=128;
else
*(pnoMaxSuppress + i * nwidth + j)=0;
}
}
}
/*
函数: EstimateThreshold
功能:估计上下阈值
参数:noMaxSuppress -- in 非极大值抑制处理结果 疑似边缘点图
amplitue -- in 幅度图
lowthreashold -- out 低阈值
highthreshold --out 高阈值
*/
void EstimateThreshold(IplImage * noMaxSuppressImg , IplImage * amplitue , float &lowthreahold ,float &highthreahold)
{
int hist[1024]; // 幅度值最大不超过 sqrt (255^2 + 255^2)
for (int i = 0 ; i < 1024 ; ++i)
{
hist[i] = 0 ;
}
int nwidth = amplitue->width;
int nheight = amplitue->height;
float * pnoMaxSuppress ,*pAmp;
pnoMaxSuppress = (float *)noMaxSuppressImg->imageData ;
pAmp = (float *)amplitue->imageData;
//直方图统计 x轴幅度值 y 在该幅度值的点数
for (int i = 0 ; i <nheight ; ++i)
{
for (int j = 0 ; j <nwidth ; ++j)
{
if(*(pnoMaxSuppress + i * nwidth + j)==128)
hist[(int)*(pAmp +i*nwidth + j)]++;
}
}
//确定最大幅度值,统计疑似边缘点个数
int nEdgeNum = 0 ;
int nmaxAmp = 0 ;
for (int i = 0 ; i< 512 ; ++i)
{
if (hist[i]!=0)
{
nmaxAmp = i ;
nEdgeNum += hist[i];
}
}
//计算上下阈值
//将疑似边缘点的幅度值排序,79%处为高阈值,显然 hist直方图是按从小到大存储的,所以省略了排序步骤
int highcount = 0.79 * nEdgeNum + 0.5 ;
int nEdgeCount = 0 ; //计数
for (int i = 0 ; i <=nmaxAmp ; ++i)
{
nEdgeCount += hist[i] ;
if (nEdgeCount >= highcount)
{
highthreahold = i ;
lowthreahold = 0.5 * highthreahold + 0.5;
break;
}
}
}
/*
函数:hysteresis
功能:滞后阈值化
参数:amplitue -- in 幅度图
nomaxsuppress -- in 疑似边缘图
highthreshold -- in 高阈值
lowthreshold -- in 低阈值
*/
void hysteresis(IplImage * amplitue , IplImage * noMaxSuppressImg , float highthreshold , float lowthreshold)
{
int nwidth = noMaxSuppressImg->width ;
int nheight = noMaxSuppressImg->height ;
float * pnoMaxSuppress , *pAmp;
pnoMaxSuppress = (float *)noMaxSuppressImg->imageData;
pAmp = (float *)amplitue->imageData ;
for (int i = 0 ; i<nheight ; ++i)
{
for (int j = 0 ; j <nwidth ; ++j)
{
if ( *(pAmp + i * nwidth + j)>=highthreshold && *(pnoMaxSuppress + i * nwidth + j)==128)
{
*(pnoMaxSuppress + i * nwidth + j) =255;
trackEdge( j , i , lowthreshold , amplitue ,noMaxSuppressImg);
}
}
}
for (int i = 0 ; i <nheight ; ++i)
{
for (int j = 0 ; j < nwidth ;++j)
{
if (*(pnoMaxSuppress + i * nwidth + j)!=255)
{
*(pnoMaxSuppress + i * nwidth + j) = 0 ;
}
}
}
}
/*
函数:trackEdge
功能:找八邻域满足条件的点
参数:x ,y -- in 中心点
lowthreshold -- in 低阈值
amplitue -- 幅度图
momaxsuppress -- 疑似边缘图
*/
void trackEdge(int x , int y , int lowthreshold, IplImage * amplitue , IplImage * noMaxSuppressImg)
{
int xnum[8] = {-1 , 0 , 1 , -1 , 1, -1 ,0 , 1};
int ynum[8] = {-1, -1 ,-1 ,0 ,0 , 1, 1, 1};
int nwidth = noMaxSuppressImg->width;
int nheight = noMaxSuppressImg->height;
float * pAmp , * pnoMaxSuppress;
pAmp = (float *)amplitue->imageData ;
pnoMaxSuppress = (float *)noMaxSuppressImg->imageData ;
for ( int i = 0 ; i < 8 ; ++i)
{
if ( *(pAmp + (y + ynum[i]) * nwidth + x + xnum[i]) >=lowthreshold &&
*(pnoMaxSuppress + (y + ynum[i]) * nwidth + x + xnum[i])==128)
{
*(pnoMaxSuppress + (y + ynum[i]) * nwidth + x + xnum[i])=255;
trackEdge(x+xnum[i] , y +ynum[i] , lowthreshold,amplitue , noMaxSuppressImg);
}
}
}
/*
函数:canny
功能:canny边缘提取
参数:src -- in 原图
dst -- in 边缘图
sigma -- in 高斯核半径
size -- in 高斯核大小
highthreshold --- 高阈值 拓展用 任意值
lowthreshold -- 低阈值 拓展用 任意值
*/
void canny(IplImage * src, IplImage * dst ,float sigma , int size , float highthreshold, float lowthreshold)
{
//高斯滤波
float *guassfilter = (float * )malloc( size * size * sizeof(float));
createGuassFilter(guassfilter,size,sigma);
CvMat mat_filter;
cvInitMatHeader(&mat_filter,size,size,CV_32FC1 , (float *)guassfilter);
IplImage * filterimg = cvCreateImage(cvGetSize(src),IPL_DEPTH_32F,1);
cvFilter2D(src,filterimg,&mat_filter);
//计算梯度方向和幅度
IplImage * gradientx = cvCreateImage(cvGetSize(src),IPL_DEPTH_32F,1);
IplImage * gradienty = cvCreateImage(cvGetSize(src),IPL_DEPTH_32F,1);
IplImage * amplitue = cvCreateImage(cvGetSize(src),IPL_DEPTH_32F,1);
calGradientAndAmplitue(filterimg , gradientx,gradienty ,amplitue);
//非极大值抑制
IplImage * noMaxSuppressImg = cvCreateImage(cvGetSize(src),IPL_DEPTH_32F,1);
cvZero(noMaxSuppressImg);
noMaxSuppress(gradientx,gradienty,amplitue,noMaxSuppressImg);
//确定上下阈值
float high , low;
EstimateThreshold(noMaxSuppressImg,amplitue,low,high);
//阈值滞后处理
hysteresis(amplitue,noMaxSuppressImg,high,low);
IplImage * canny = cvCreateImage(cvGetSize(src),8,1);
cvConvert(noMaxSuppressImg,canny);
cvShowImage("canny",canny);
memcpy(dst->imageData , noMaxSuppressImg->imageData,dst->imageSize);
cvReleaseImage(&filterimg);
cvReleaseImage(&gradientx);
cvReleaseImage(&gradienty);
cvReleaseImage(&litue);
cvReleaseImage(&noMaxSuppressImg);
free(guassfilter);
}
//
#include "stdafx.h"
#include <cv.h>
#include <highgui.h>
#include <cxcore.h>
#ifdef DEBUG
#pragma comment(lib,"cxcore200d.lib")
#pragma comment(lib,"cv200d.lib")
#pragma comment(lib,"highgui200d.lib")
#else
#pragma comment(lib,"cxcore200.lib")
#pragma comment(lib,"cv200.lib")
#pragma comment(lib,"highgui200.lib")
#endif
void creatGuassFilter(float *GuassFilter, int nsize,float sigma=0);
void calGradientAndAmplitue(IplImage * src , IplImage * gradientX,IplImage * gradientY ,IplImage * amplitue );
void noMaxSuppress(IplImage * gradientX, IplImage *gradientY, IplImage *amplitue,IplImage * noMaxSupress );
void EstimateThreshold(IplImage * noMaxSuppressImg , IplImage * amplitue , double &lowthreahold ,double &highthreahold);
void hysteresis(IplImage * amplitue , IplImage * noMaxSuppressImg , float highthreshold , float lowthreshold);
void trackEdge(int x , int y , int lowthreshold, IplImage * amplitue , IplImage * noMaxSuppressImg);
void canny(IplImage * src, IplImage * dst ,float sigma , int size , float highthreshold, float lowthreshold);
int _tmain(int argc, _TCHAR* argv[])
{
IplImage * im = cvLoadImage("e.bmp",CV_LOAD_IMAGE_GRAYSCALE);
IplImage * im2 = cvCreateImage(cvSize(im->width/10,im->height/10),8,1);
cvResize(im,im2);
IplImage * dst = cvCreateImage(cvGetSize(im2),8,1);
canny(im2,dst,3,7,0,0);
cvWaitKey(-1);
return 0;
}
/*
函数名 : creatGuassFilter
功能:生成高斯滤波模板
参数:mat_GuassFilter --- out 用于存储生成的高斯滤波模板
nsize ----- in 滤波窗口尺寸
sigma ----- in 高斯半径
*/
void createGuassFilter(float *GuassFilter, int nsize,float sigma)
{
float distX,distY; //到窗口中心的距离
int ncenter = nsize/2; //窗口中心
float sigmaA = 1.0 / (sqrt(2 * CV_PI) * sigma) ;
float sum =0.0;
//如果sigma = 0 ,自动计算sigma值
if (0 == sigma)
{
sigma = (nsize /2 -1) *0.3 + 0.8;
}
for (int i = 0 ; i < nsize ; ++i)
{
for (int j = 0 ; j < nsize ; ++j)
{
distX = i - ncenter;
distY = j - ncenter;
*(GuassFilter + i * nsize + j ) = sigmaA * exp( -1.0 * ( distX * distX + distY * distY) / (2 * sigma * sigma)) ;
sum +=*(GuassFilter + i * nsize + j ) ;
}
}
//归一化
for (int i = 0 ; i <nsize ; ++i)
{
for (int j=0; j <nsize; ++j)
{
*(GuassFilter + i * nsize + j ) /=sum;
}
}
}
/*
函数:calGradientAndAmplitue
功能:计算图像梯度方向 和幅度
参数:src -- in 原图
gradientX Y -- out x,y方向梯度值 像素类型 int
amplitue -- out 幅度值 像素类型double
*/
void calGradientAndAmplitue(IplImage * src , IplImage * gradientX,IplImage * gradientY ,IplImage * amplitue )
{
//计算 x ,y 方向梯度
float filterx[3][3] = {
-1 , 0 , 1 ,
-2 , 0 , 2,
-1 , 0 , 1
};
float filtery[3][3] = {
1 , 2 , 1,
0 , 0 , 0,
-1 ,-2, -1
};
CvMat mat_filterx ;
CvMat mat_filtery;
cvInitMatHeader(&mat_filterx , 3, 3,CV_32FC1 ,filterx);
cvInitMatHeader(&mat_filtery , 3 ,3,CV_32FC1, filtery);
cvFilter2D(src,gradientX,&mat_filterx);
cvFilter2D(src,gradientY,&mat_filtery);
//计算幅度值
int nwidth = src->width ;
int nheight = src->height;
float Gx, Gy; //x y 方向梯度
float * pGx , *pGy ,*pAmp;
for (int i = 0 ; i < nheight ; ++ i)
{
pGx = (float *)(gradientX->imageData + i* gradientX->widthStep );
pGy = (float *)(gradientY->imageData + i* gradientX->widthStep );
pAmp = (float *)(amplitue->imageData + i * amplitue->widthStep);
for (int j = 0 ;j < nwidth ; ++j)
{
Gx = *(pGx + j);
Gy = *(pGy + j);
*(pAmp + j) = sqrt(1.0 * Gx * Gx + 1.0 * Gy * Gy);
}
}
}
/*
函数:noMaxSuppress
功能:非极大值抑制
*/
void noMaxSuppress(IplImage * gradientX, IplImage *gradientY, IplImage *amplitue,IplImage * noMaxSupress )
{
int nwidth = noMaxSupress->width ;
int nheight = noMaxSupress->height;
//非极大值抑制
float Gx ,Gy ; // x,y梯度值
float temp1,temp2,temp3,temp4; //用于计算插值得到梯度方向邻域幅度值的中间变量
float amp , amp1 ,amp2; // 该点幅度值 该点梯度方向上的幅度值
float angle = 0.0; //梯度方向
float rate = 0.0; // 线插比例
float *pAmp , *pnoMaxSuppress,*pgradientX,*pgradientY;
pAmp = (float *)amplitue->imageData;
pnoMaxSuppress = (float *)noMaxSupress->imageData;
pgradientX = (float *)gradientX->imageData;
pgradientY = (float *)gradientY->imageData;
for ( int i = 1 ; i <nheight -1 ; ++i)
{
for (int j = 1 ; j < nwidth -1 ; ++j)
{
amp =*( pAmp+ i *nwidth + j);
//如果幅度值为0 则为非边缘点
if (amp ==0)
{
*(pnoMaxSuppress + i * nwidth + j) = 0 ;
continue;
}
Gx = *(pgradientX + i * nwidth + j);
Gy = *(pgradientY + i *nwidth + j);
angle = atan(1.0 * Gy / Gx);
//第一种情况
// temp1
//temp3 amp temp2
//temp4
if ( (angle >=0 && angle < CV_PI/4) /* || (angle >=CV_PI && angle < 5.0/4 *CV_PI)*/ )
{
temp1 = *(pAmp + (i-1)*nwidth + j+1);
temp2 = *(pAmp + i *nwidth + j+1);
temp4 = *(pAmp + i *nwidth + j-1);
temp4 = *(pAmp + (i +1) * nwidth + j-1);
rate = fabs(1.0 *Gy)/fabs(1.0*Gx);
amp1 = temp1 * rate + temp2 * (1 -rate);
amp2 = temp3 * rate + temp4 * (1 -rate);
}
//第二种情况
// temp1 temp2
// amp
// temp3 temp4
else if ( ( angle >=CV_PI/4 && angle <= CV_PI/2)/* || (angle >=5.0/4 * CV_PI && angle <3.0/2 * CV_PI )*/)
{
temp1 = *(pAmp + (i-1) * nwidth + j );
temp2 =*(pAmp + (i-1) * nwidth + j +1);
temp3 = *(pAmp + (i +1)*nwidth + j -1);
temp4 = *(pAmp + (i +1) * nwidth + j);
rate = fabs(1.0*Gx)/fabs(1.0*Gy);
amp1 = temp2 * rate + temp1 *(1- rate);
amp2 = temp3 * rate + temp4 *(1-rate);
}
//第三种情况
// temp1 temp2
// amp
// temp3 temp4
else if( ( angle >= -1.0 * CV_PI/2 && angle < -1.0 *CV_PI /4))
{
temp1 = *(pAmp + (i-1)*nwidth + j-1);
temp2 = *(pAmp + (i-1)*nwidth + j);
temp3 =*(pAmp + (i+1)*nwidth + j);
temp4 = *(pAmp + (i +1)*nwidth +j+1);
rate = fabs(1.0*Gx) / fabs(1.0*Gy);
amp1 = temp1 * rate + temp2 *(1 - rate);
amp2 = temp4 * rate + temp3 * (1- rate);
}
//第四种情况
// temp1
// temp2 amp temp3
// temp4
else /*if ( angle >= -1.0 * CV_PI /4 && angle <0)*/
{
temp1 = *(pAmp + (i-1)*nwidth + j-1);
temp2 = *(pAmp + i * nwidth + j-1);
temp3 = *(pAmp + i * nwidth + j +1);
temp4 =*(pAmp + (i+1)*nwidth + j +1);
rate = fabs(1.0*Gy)/fabs(1.0*Gx);
amp1 = temp1 * rate + temp2 *( 1- rate);
amp2 = temp4 * rate + temp3* (1- rate);
}
//局部最大值判断 如果为最大值 128 否则0
if( amp >=amp1 && amp >= amp2)
*(pnoMaxSuppress + i * nwidth + j)=128;
else
*(pnoMaxSuppress + i * nwidth + j)=0;
}
}
}
/*
函数: EstimateThreshold
功能:估计上下阈值
参数:noMaxSuppress -- in 非极大值抑制处理结果 疑似边缘点图
amplitue -- in 幅度图
lowthreashold -- out 低阈值
highthreshold --out 高阈值
*/
void EstimateThreshold(IplImage * noMaxSuppressImg , IplImage * amplitue , float &lowthreahold ,float &highthreahold)
{
int hist[1024]; // 幅度值最大不超过 sqrt (255^2 + 255^2)
for (int i = 0 ; i < 1024 ; ++i)
{
hist[i] = 0 ;
}
int nwidth = amplitue->width;
int nheight = amplitue->height;
float * pnoMaxSuppress ,*pAmp;
pnoMaxSuppress = (float *)noMaxSuppressImg->imageData ;
pAmp = (float *)amplitue->imageData;
//直方图统计 x轴幅度值 y 在该幅度值的点数
for (int i = 0 ; i <nheight ; ++i)
{
for (int j = 0 ; j <nwidth ; ++j)
{
if(*(pnoMaxSuppress + i * nwidth + j)==128)
hist[(int)*(pAmp +i*nwidth + j)]++;
}
}
//确定最大幅度值,统计疑似边缘点个数
int nEdgeNum = 0 ;
int nmaxAmp = 0 ;
for (int i = 0 ; i< 512 ; ++i)
{
if (hist[i]!=0)
{
nmaxAmp = i ;
nEdgeNum += hist[i];
}
}
//计算上下阈值
//将疑似边缘点的幅度值排序,79%处为高阈值,显然 hist直方图是按从小到大存储的,所以省略了排序步骤
int highcount = 0.79 * nEdgeNum + 0.5 ;
int nEdgeCount = 0 ; //计数
for (int i = 0 ; i <=nmaxAmp ; ++i)
{
nEdgeCount += hist[i] ;
if (nEdgeCount >= highcount)
{
highthreahold = i ;
lowthreahold = 0.5 * highthreahold + 0.5;
break;
}
}
}
/*
函数:hysteresis
功能:滞后阈值化
参数:amplitue -- in 幅度图
nomaxsuppress -- in 疑似边缘图
highthreshold -- in 高阈值
lowthreshold -- in 低阈值
*/
void hysteresis(IplImage * amplitue , IplImage * noMaxSuppressImg , float highthreshold , float lowthreshold)
{
int nwidth = noMaxSuppressImg->width ;
int nheight = noMaxSuppressImg->height ;
float * pnoMaxSuppress , *pAmp;
pnoMaxSuppress = (float *)noMaxSuppressImg->imageData;
pAmp = (float *)amplitue->imageData ;
for (int i = 0 ; i<nheight ; ++i)
{
for (int j = 0 ; j <nwidth ; ++j)
{
if ( *(pAmp + i * nwidth + j)>=highthreshold && *(pnoMaxSuppress + i * nwidth + j)==128)
{
*(pnoMaxSuppress + i * nwidth + j) =255;
trackEdge( j , i , lowthreshold , amplitue ,noMaxSuppressImg);
}
}
}
for (int i = 0 ; i <nheight ; ++i)
{
for (int j = 0 ; j < nwidth ;++j)
{
if (*(pnoMaxSuppress + i * nwidth + j)!=255)
{
*(pnoMaxSuppress + i * nwidth + j) = 0 ;
}
}
}
}
/*
函数:trackEdge
功能:找八邻域满足条件的点
参数:x ,y -- in 中心点
lowthreshold -- in 低阈值
amplitue -- 幅度图
momaxsuppress -- 疑似边缘图
*/
void trackEdge(int x , int y , int lowthreshold, IplImage * amplitue , IplImage * noMaxSuppressImg)
{
int xnum[8] = {-1 , 0 , 1 , -1 , 1, -1 ,0 , 1};
int ynum[8] = {-1, -1 ,-1 ,0 ,0 , 1, 1, 1};
int nwidth = noMaxSuppressImg->width;
int nheight = noMaxSuppressImg->height;
float * pAmp , * pnoMaxSuppress;
pAmp = (float *)amplitue->imageData ;
pnoMaxSuppress = (float *)noMaxSuppressImg->imageData ;
for ( int i = 0 ; i < 8 ; ++i)
{
if ( *(pAmp + (y + ynum[i]) * nwidth + x + xnum[i]) >=lowthreshold &&
*(pnoMaxSuppress + (y + ynum[i]) * nwidth + x + xnum[i])==128)
{
*(pnoMaxSuppress + (y + ynum[i]) * nwidth + x + xnum[i])=255;
trackEdge(x+xnum[i] , y +ynum[i] , lowthreshold,amplitue , noMaxSuppressImg);
}
}
}
/*
函数:canny
功能:canny边缘提取
参数:src -- in 原图
dst -- in 边缘图
sigma -- in 高斯核半径
size -- in 高斯核大小
highthreshold --- 高阈值 拓展用 任意值
lowthreshold -- 低阈值 拓展用 任意值
*/
void canny(IplImage * src, IplImage * dst ,float sigma , int size , float highthreshold, float lowthreshold)
{
//高斯滤波
float *guassfilter = (float * )malloc( size * size * sizeof(float));
createGuassFilter(guassfilter,size,sigma);
CvMat mat_filter;
cvInitMatHeader(&mat_filter,size,size,CV_32FC1 , (float *)guassfilter);
IplImage * filterimg = cvCreateImage(cvGetSize(src),IPL_DEPTH_32F,1);
cvFilter2D(src,filterimg,&mat_filter);
//计算梯度方向和幅度
IplImage * gradientx = cvCreateImage(cvGetSize(src),IPL_DEPTH_32F,1);
IplImage * gradienty = cvCreateImage(cvGetSize(src),IPL_DEPTH_32F,1);
IplImage * amplitue = cvCreateImage(cvGetSize(src),IPL_DEPTH_32F,1);
calGradientAndAmplitue(filterimg , gradientx,gradienty ,amplitue);
//非极大值抑制
IplImage * noMaxSuppressImg = cvCreateImage(cvGetSize(src),IPL_DEPTH_32F,1);
cvZero(noMaxSuppressImg);
noMaxSuppress(gradientx,gradienty,amplitue,noMaxSuppressImg);
//确定上下阈值
float high , low;
EstimateThreshold(noMaxSuppressImg,amplitue,low,high);
//阈值滞后处理
hysteresis(amplitue,noMaxSuppressImg,high,low);
IplImage * canny = cvCreateImage(cvGetSize(src),8,1);
cvConvert(noMaxSuppressImg,canny);
cvShowImage("canny",canny);
memcpy(dst->imageData , noMaxSuppressImg->imageData,dst->imageSize);
cvReleaseImage(&filterimg);
cvReleaseImage(&gradientx);
cvReleaseImage(&gradienty);
cvReleaseImage(&litue);
cvReleaseImage(&noMaxSuppressImg);
free(guassfilter);
}
作者: 小马_xiao
出处:http://www.cnblogs.com/xiaomaLV2/>
关于作者:专注halcon\opencv\机器视觉
本文版权归作者,未经作者同意必须保留此段声明,且在文章页面明显位置给出 原文链接