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( &amplitue);
    cvReleaseImage( &noMaxSuppressImg);

    free(guassfilter);
}

通过 Wiz 发布

转载于:https://www.cnblogs.com/xiaomaLV2/archive/2012/01/06/2314850.html