2014-10-05 16:17:46 gongem 阅读数 12126

1. 边缘检测的概念

边缘检测是图像处理与计算机视觉中极为重要的一种分析图像的方法,至少在我做图像分析与识别时,边缘是我最喜欢的图像特征。边缘检测的目的就是找到图像中亮度变化剧烈的像素点构成的集合,表现出来往往是轮廓。在对现实世界的图像采集中,有下面4种情况会表现在图像中时形成一个边缘。

  1. 深度的不连续(物体处在不同的物平面上);
  2. 表面方向的不连续(如正方体的不同的两个面);
  3. 物体材料不同(这样会导致光的反射系数不同);
  4. 场景中光照不同(如被树萌投向的地面);


    image上面的图像是图像中水平方向7个像素点的灰度值显示效果,我们很容易地判断在第4和第5个像素之间有一个边缘,因为它俩之间发生了强烈的灰度跳变。在实际的边缘检测中,边缘远没有上图这样简单明显,我们需要取对应的阈值来区分出它们。

    原图 Roberts边缘 Prewitt边缘
    image image image
    Sobel边缘 Log边缘 Canny边缘
    image image image
     


    2. 边缘检测的基本方法

    2.1 一阶微分边缘算子

    一阶微分边缘算子也称为梯度边缘算子,它是利用图像在边缘处的阶跃性,即图像梯度在边缘取得极大值的特性进行边缘检测。梯度是

    一个矢量,它具有方向θ和模|ΔI|

    梯度的方向提供了边缘的趋势信息,因为梯度方向始终是垂直于边缘方向,梯度的模值大小提供了边缘的强度信息。

    在实际使用中,通常利用有限差分进行梯度近似。对于上面的公式,我们有如下的近似:

    2.2 Roberts边缘检测算子

    1963年,Roberts提出了这种寻找边缘的算子。Roberts边缘算子是一个2x2的模板,采用的是对角方向相邻的两个像素之差。从图像处理的实际效果来看,边缘定位较准,对噪声敏感。在Roberts检测算子中:

    可以导出Roberts在点(i+1/2,j+1/2)处的水平与竖直边缘检测卷积核为:

    2.3 Prewitt边缘检测算子

    Prewitt利用周围邻域8个点的灰度值来估计中心的梯度,它的梯度计算公式如下:

    所以,Prewitt的卷积核为:

    2.4 Sobel边缘检测算子

    比起Prewitt算子,Sobel也是用周围8个像素来估计中心像素的梯度,但是Sobel算子认为靠近中心像素的点应该给予更高的权重,所以Sobel算子把与中心像素4邻接的像素的权重设置为2或-2。

    Sobel边缘检测算子的卷积核为:

    Sobel进行边缘检测的实现可以参考我原来写的一篇博文:图像特征检测:sobel边缘检测,重要的是梯度图像计算后的阈值的确定与边缘的非极大值抑制算法,Roberts与Prewitt原理与sobel一致。

    3. 二阶微分算子

    学过微积分我们都知道,边缘即是图像的一阶导数局部最大值的地方,那么也意味着该点的二阶导数为零。二阶微分边缘检测算子就是利用图像在边缘处的阶跃性导致图像二阶微分在边缘处出现零值这一特性进行边缘检测的。

    对于图像的二阶微分可以用拉普拉斯算子来表示:

    我们在像素点(i,j)3×3的邻域内,可以有如下的近似:

    对应的二阶微分卷积核为:

    所以二阶微分检测边缘的方法就分两步:1)用上面的Laplace核与图像进行卷积;2)对卷积后的图像,取得那些卷积结果为0的点。

    虽然上述使用二阶微分检测边缘的方法简单,但它的缺点是对噪声十分敏感,同时也没有能够提供边缘的方向信息。为了实现对噪声的抑制,Marr等提出了LOG的方法。

    为了减少噪声对边缘的影响,首先图像要进行低通滤波,LOG采用了高斯函数作为低通滤波器。高斯函数为:

    上面的公式中σ决定了对图像的平滑程度。高斯函数生成的滤波模板尺寸一般设定为6σ+1(加1是会了使滤波器的尺寸为奇数)。使用高斯函数对图像进行滤波并对图像滤波结果进行二阶微分运算的过程,可以转换为先对高斯函数进行二阶微分,再利用高斯函数的二阶微分结果对图像进行卷积运算:

    关于LOG算子的计算在上一篇文章:图像特征提取:斑点检测中有实现的代码。

    4. Canny边缘检测

    canny边缘检测实际上是一种一阶微分算子检测算法,但为什么这里拿出来说呢,因为它几乎是边缘检测算子中最为常用的一种,也是个人认为现在最优秀的边缘检测算子。Canny提出了边缘检测算子优劣评判的三条标准:

    Canny边缘检测之所以优秀是因为它在一阶微分算子的基础上,增加了非最大值抑制和双阈值两项改进。利用非极大值抑制不仅可以有效地抑制多响应边缘,而且还可以提高边缘的定位精度;利用双阈值可以有效减少边缘的漏检率。

    Canny边缘检测主要分四步进行:

    1. 去噪声;
    2. 计算梯度与方向角;
    3. 非最大值抑制;
    4. 滞后阈值化;

    其中前两步很简单,先用一个高斯滤波器对图像进行滤波,然后用Sobel水平和竖直检测子与图像卷积,来计算梯度和方向角。

    非极大值抑制

    图像梯度幅值矩阵中的元素值越大,说明图像中该点的梯度值越大,但这不不能说明该点就是边缘(这仅仅是属于图像增强的过程)。在Canny算法中,非极大值抑制是进行边缘检测的重要步骤,通俗意义上是指寻找像素点局部最大值,将非极大值点所对应的灰度值置为0,这样可以剔除掉一大部分非边缘的点。

    image

    根据上图可知,要进行非极大值抑制,就首先要确定像素点C的灰度值在其8值邻域内是否为最大。图中蓝色的线条方向为C点的梯度方向,这样就可以确定其局部的最大值肯定分布在这条线上,也即出了C点外,梯度方向的交点dTmp1和dTmp2这两个点的值也可能会是局部最大值。因此,判断C点灰度与这两个点灰度大小即可判断C点是否为其邻域内的局部最大灰度点。如果经过判断,C点灰度值小于这两个点中的任一个,那就说明C点不是局部极大值,那么则可以排除C点为边缘。这就是非极大值抑制的工作原理。

    在理解的过程中需要注意以下两点:

    1. 中非最大抑制是回答这样一个问题:“当前的梯度值在梯度方向上是一个局部最大值吗?” 所以,要把当前位置的梯度值与梯度方向上两侧的梯度值进行比较;
    2. 梯度方向垂直于边缘方向。

    但实际上,我们只能得到C点邻域的8个点的值,而dTmp1和dTmp2并不在其中,要得到这两个值就需要对该两个点两端的已知灰度进行线性插值,也即根据图中的g1和g2对dTmp1进行插值,根据g3和g4对dTmp2进行插值,这要用到其梯度方向。

    滞后阈值化

    由于噪声的影响,经常会在本应该连续的边缘出现断裂的问题。滞后阈值化设定两个阈值:一个为高阈值Th,一个为低阈值Tl。如果任何像素边缘算子的影响超过高阈值,将这些像素标记为边缘;响应超过低阈值(高低阈值之间)的像素,如果与已经标记为边缘的像素4-邻接或8-邻接,则将这些像素也标记为边缘。所以不整个过程描述如下:

    1. 如果该像素的梯度值小于Tl,则该像素为非边缘像素;
    2. 如果该像素的梯度值大于Th,则该像素为边缘像素;
    3. 如果该像素的梯度值介于TlTh之间,需要进一步检测该像素的3×3邻域内的8个点,如果这8个点内有一个或以上的点梯度超过了Th,则该像素为边缘像素,否则不是边缘像素。
     


    Canny边缘检测的实现


    #include <iostream>


    #include "opencv2/core/core.hpp"
    #include "opencv2/highgui/highgui.hpp"
    #include "opencv2/imgproc/imgproc.hpp"




    using namespace std;
    using namespace cv;




    void getCannyEdge(const Mat& imgSrc, Mat& imgDst, double lowThresh = -1, double highThresh = -1, double sigma = 1);
    double getGaussianThresh(const Mat& inputArray, double percentage);


    int main(int argc, char** argv)
    {
        Mat image = imread(argv[1]);
        Mat imgCanny;
        getCannyEdge(image,imgCanny);
        t = ((double)getTickCount() - t)/getTickFrequency();
        imwrite("imgCanny.png", imgCanny);


        return 0;
    }


    void getCannyEdge(const Mat& imgSrc, Mat& imgDst, double lowThresh, double highThresh, double sigma)
    {
        Mat gray;
        if (imgSrc.channels() == 3)
        {
            cvtColor(imgSrc, gray, CV_BGR2GRAY);
        }
        else
        {
            gray = imgSrc.clone();
        }
        gray.convertTo(gray, CV_64F);
        gray = gray / 255;
        
        double gaussianDieOff = .0001;
        double percentOfPixelsNotEdges = .7; // Used for selecting thresholds
        double thresholdRatio = .4;   // Low thresh is this fraction of the high


        int possibleWidth = 30;
        double ssq = sigma * sigma;
        for (int i = 1; i <= possibleWidth; i++)
        {
            if (exp(-(i * i) / (2* ssq)) < gaussianDieOff)
            {
                possibleWidth = i - 1;
                break;
            }
        }


        if (possibleWidth == 30)
        {
            possibleWidth = 1; // the user entered a reallly small sigma
        }


        // get the 1D gaussian filter
        int winSz = 2 * possibleWidth + 1;
        Mat gaussKernel1D(1, winSz, CV_64F);
        double* kernelPtr = gaussKernel1D.ptr<double>(0);
        for (int i = 0; i < gaussKernel1D.cols; i++)
        {
            kernelPtr[i] = exp(-(i - possibleWidth) * (i - possibleWidth) / (2 * ssq)) / (2 * CV_PI * ssq);
        }


        
        // get the derectional derivatives of gaussian kernel
        Mat dGaussKernel(winSz, winSz, CV_64F);
        for (int i = 0; i < dGaussKernel.rows; i++)
        {
            double* linePtr = dGaussKernel.ptr<double>(i);
            for (int j = 0; j< dGaussKernel.cols; j++)
            {
                linePtr[j] = - (j - possibleWidth) * exp(-((i - possibleWidth) * (i - possibleWidth) + (j - possibleWidth) * (j - possibleWidth)) / (2 * ssq)) / (CV_PI * ssq);
            }
        }




        /* smooth the image out*/
        Mat imgSmooth;
        filter2D(gray, imgSmooth, -1, gaussKernel1D);
        filter2D(imgSmooth, imgSmooth, -1, gaussKernel1D.t());
        /*apply directional derivatives*/


        Mat imgX, imgY;
        filter2D(imgSmooth, imgX, -1, dGaussKernel);
        filter2D(imgSmooth, imgY, -1, dGaussKernel.t());


        Mat imgMag;
        sqrt(imgX.mul(imgX) + imgY.mul(imgY), imgMag);
        double magMax;
        minMaxLoc(imgMag, NULL, &magMax, NULL, NULL);


        if (magMax > 0 )
        {
            imgMag = imgMag / magMax;
        }


        
        if (lowThresh == -1 || highThresh == -1)
        {
            highThresh = getGaussianThresh(imgMag, percentOfPixelsNotEdges);
            lowThresh = thresholdRatio * highThresh;
        }








        Mat imgStrong = Mat::zeros(imgMag.size(), CV_8U);
        Mat imgWeak = Mat::zeros(imgMag.size(), CV_8U);
        
        
        for (int dir = 1; dir <= 4; dir++)
        {
            Mat gradMag1(imgMag.size(), imgMag.type());
            Mat gradMag2(imgMag.size(), imgMag.type());
            Mat idx = Mat::zeros(imgMag.size(), CV_8U);
            if (dir == 1)
            {
                Mat dCof = abs(imgY / imgX);
                idx = (imgY <= 0 & imgX > -imgY) | (imgY >= 0 & imgX < -imgY);
                idx.row(0).setTo(Scalar(0));
                idx.row(idx.rows - 1).setTo(Scalar(0));
                idx.col(0).setTo(Scalar(0));
                idx.col(idx.cols - 1).setTo(Scalar(0));
                for (int i = 1; i < imgMag.rows - 1; i++)
                {
                    for (int j = 1; j < imgMag.cols - 1; j++)
                    {
                        gradMag1.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i,j + 1) + dCof.at<double>(i,j) * imgMag.at<double>(i - 1,j + 1);
                        gradMag2.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i,j - 1) + dCof.at<double>(i,j) * imgMag.at<double>(i + 1,j - 1);
                    }
                }
            }
            else if(dir == 2)
            {
                Mat dCof = abs(imgX / imgY);
                idx = (imgX > 0 & -imgY >= imgX) | (imgX < 0 & -imgY <= imgX);
                for (int i = 1; i < imgMag.rows - 1; i++)
                {
                    for (int j = 1; j < imgMag.cols - 1; j++)
                    {
                        gradMag1.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i - 1,j) + dCof.at<double>(i,j) * imgMag.at<double>(i - 1,j + 1);
                        gradMag2.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i + 1,j) + dCof.at<double>(i,j) * imgMag.at<double>(i + 1,j - 1);
                    }
                }
            }
            else if(dir == 3)
            {
                Mat dCof = abs(imgX / imgY);
                idx = (imgX <= 0 & imgX > imgY) | (imgX >= 0 & imgX < imgY);
                for (int i = 1; i < imgMag.rows - 1; i++)
                {
                    for (int j = 1; j < imgMag.cols - 1; j++)
                    {
                        gradMag1.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i - 1,j) + dCof.at<double>(i,j) * imgMag.at<double>(i - 1,j - 1);
                        gradMag2.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i + 1,j) + dCof.at<double>(i,j) * imgMag.at<double>(i + 1,j + 1);
                    }
                }
            
            }
            else
            {
                Mat dCof = abs(imgY / imgX);
                idx = (imgY <0 & imgX <= imgY) | (imgY > 0 & imgX >= imgY);
                for (int i = 1; i < imgMag.rows - 1; i++)
                {
                    for (int j = 1; j < imgMag.cols - 1; j++)
                    {
                        gradMag1.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i,j - 1) + dCof.at<double>(i,j) * imgMag.at<double>(i - 1,j - 1);
                        gradMag2.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i,j + 1) + dCof.at<double>(i,j) * imgMag.at<double>(i + 1,j + 1);
                    }
                }
            }


            Mat idxLocalMax = idx & ((imgMag >= gradMag1) & (imgMag >= gradMag2));




            imgWeak = imgWeak | ((imgMag > lowThresh) & idxLocalMax);
            imgStrong= imgStrong| ((imgMag > highThresh) & imgWeak);


        }


        imgDst = Mat::zeros(imgWeak.size(),imgWeak.type());
        for (int i = 1; i < imgWeak.rows - 1; i++)
        {
            uchar* pWeak = imgWeak.ptr<uchar>(i);
            uchar* pDst = imgDst.ptr<uchar>(i);
            uchar* pStrPre = imgStrong.ptr<uchar>(i - 1);
            uchar* pStrMid = imgStrong.ptr<uchar>(i);
            uchar* pStrAft = imgStrong.ptr<uchar>(i + 1);
            for (int j = 1; j < imgWeak.cols - 1; j++)
            {
                if (!pWeak[j])
                {
                    continue;
                }
                if (pStrMid[j])
                {
                    pDst[j] = 255;
                }
                if (pStrMid[j-1] || pStrMid[j+1] || pStrPre[j-1] || pStrPre[j] || pStrPre[j+1] || pStrAft[j-1] || pStrAft[j] ||pStrAft[j+1])
                {
                    pDst[j] = 255;
                }
            }
        }
    }


    double getGaussianThresh(const Mat& inputArray, double percentage)
    {
        double thresh = -1.0;
        // compute the 64-hist of inputArray
        int nBins = 64;
        double minValue, maxValue;
        minMaxLoc(inputArray, &minValue, &maxValue, NULL, NULL);
        double step = (maxValue - minValue) / nBins;


        vector<unsigned> histBin(nBins,0);
        for (int i = 0; i < inputArray.rows; i++)
        {
            const double* pData = inputArray.ptr<double>(i);
            for(int j = 0; j < inputArray.cols; j++)
            {


                int index = (pData[j] - minValue) / step;
                histBin[index]++;
            }
        }
        unsigned cumSum = 0; 
        for (int i = 0; i < nBins; i++)
        {
            cumSum += histBin[i];


            if (cumSum > percentage * inputArray.rows * inputArray.cols)
            {
                thresh = (i + 1) / 64.0;
                break;
            }
        }
        return thresh;
        
    }
  5. 文章引子:

         

    http://www.cnblogs.com/ronny/p/4001910.html

原图 Roberts边缘 Prewitt边缘
image image image
Sobel边缘 Log边缘 Canny边缘
image image image
2014-11-16 16:50:47 jia20003 阅读数 60251

图像处理之Canny 边缘检测

一:历史

Canny边缘检测算法是1986年有John F. Canny开发出来一种基于图像梯度计算的边缘

检测算法,同时Canny本人对计算图像边缘提取学科的发展也是做出了很多的贡献。尽

管至今已经许多年过去,但是该算法仍然是图像边缘检测方法经典算法之一。

二:Canny边缘检测算法

经典的Canny边缘检测算法通常都是从高斯模糊开始,到基于双阈值实现边缘连接结束

。但是在实际工程应用中,考虑到输入图像都是彩色图像,最终边缘连接之后的图像要

二值化输出显示,所以完整的Canny边缘检测算法实现步骤如下:

1.      彩色图像转换为灰度图像

2.      对图像进行高斯模糊

3.      计算图像梯度,根据梯度计算图像边缘幅值与角度

4.      非最大信号压制处理(边缘细化)

5.      双阈值边缘连接处理

6.      二值化图像输出结果

三:各步详解与代码实现

1.      彩色图像转灰度图像

根据彩色图像RGB转灰度公式:gray  =  R * 0.299 + G * 0.587 + B * 0.114

将彩色图像中每个RGB像素转为灰度值的代码如下:

int gray = (int) (0.299 * tr + 0.587 * tg + 0.114 * tb);

2.      对图像进行高斯模糊

图像高斯模糊时,首先要根据输入参数确定高斯方差与窗口大小,这里我设置默认方

差值窗口大小为16x16,根据这两个参数生成高斯卷积核算子的代码如下:

		float kernel[][] = new float[gaussianKernelWidth][gaussianKernelWidth];
		for(int x=0; x<gaussianKernelWidth; x++)
		{
			for(int y=0; y<gaussianKernelWidth; y++)
			{
				kernel[x][y] = gaussian(x, y, gaussianKernelRadius);
			}
		}

获取了高斯卷积算子之后,我们就可以对图像高斯卷积模糊,关于高斯图像模糊更详

细的解释可以参见这里:http://blog.csdn.net/jia20003/article/details/7234741实现

图像高斯卷积模糊的代码如下:

// 高斯模糊 -灰度图像
int krr = (int)gaussianKernelRadius;
for (int row = 0; row < height; row++) {
	for (int col = 0; col < width; col++) {
		index = row * width + col;
		double weightSum = 0.0;
		double redSum = 0;
		for(int subRow=-krr; subRow<=krr; subRow++)
		{
			int nrow = row + subRow;
			if(nrow >= height || nrow < 0)
			{
				nrow = 0;
			}
			for(int subCol=-krr; subCol<=krr; subCol++)
			{
				int ncol = col + subCol;
				if(ncol >= width || ncol <=0)
				{
					ncol = 0;
				}
				int index2 = nrow * width + ncol;
				int tr1 = (inPixels[index2] >> 16) & 0xff;
				redSum += tr1*kernel[subRow+krr][subCol+krr];
				weightSum += kernel[subRow+krr][subCol+krr];
			}
		}
		int gray = (int)(redSum / weightSum);
		outPixels[index] = gray;
	}
}

3.      计算图像X方向与Y方向梯度,根据梯度计算图像边缘幅值与角度大小

高斯模糊的目的主要为了整体降低图像噪声,目的是为了更准确计算图像梯度及边缘

幅值。计算图像梯度可以选择算子有Robot算子、Sobel算子、Prewitt算子等。关于

图像梯度计算更多的解释可以看这里:

http://blog.csdn.net/jia20003/article/details/7664777。

这里采用更加简单明了的2x2的算子,其数学表达如下:


// 计算梯度-gradient, X放与Y方向
data = new float[width * height];
magnitudes = new float[width * height];
for (int row = 0; row < height; row++) {
	for (int col = 0; col < width; col++) {
		index = row * width + col;
		// 计算X方向梯度
		float xg = (getPixel(outPixels, width, height, col, row+1) - 
				getPixel(outPixels, width, height, col, row) + 
				getPixel(outPixels, width, height, col+1, row+1) -
				getPixel(outPixels, width, height, col+1, row))/2.0f;
		float yg = (getPixel(outPixels, width, height, col, row)-
				getPixel(outPixels, width, height, col+1, row) +
				getPixel(outPixels, width, height, col, row+1) -
				getPixel(outPixels, width, height, col+1, row+1))/2.0f;
		// 计算振幅与角度
		data[index] = hypot(xg, yg);
		if(xg == 0)
		{
			if(yg > 0)
			{
				magnitudes[index]=90;						
			}
			if(yg < 0)
			{
				magnitudes[index]=-90;
			}
		}
		else if(yg == 0)
		{
			magnitudes[index]=0;
		}
		else
		{
			magnitudes[index] = (float)((Math.atan(yg/xg) * 180)/Math.PI);					
		}
		// make it 0 ~ 180
		magnitudes[index] += 90;
	}
}

在获取了图像每个像素的边缘幅值与角度之后

4.      非最大信号压制

信号压制本来是数字信号处理中经常用的,这里的非最大信号压制主要目的是实现边

缘细化,通过该步处理边缘像素进一步减少。非最大信号压制主要思想是假设3x3的

像素区域,中心像素P(x,y) 根据上一步中计算得到边缘角度值angle,可以将角度分

为四个离散值0、45、90、135分类依据如下:

其中黄色区域取值范围为0~22.5 与157.5~180

绿色区域取值范围为22.5 ~ 67.5

蓝色区域取值范围为67.5~112.5

红色区域取值范围为112.5~157.5

分别表示上述四个离散角度的取值范围。得到角度之后,比较中心像素角度上相邻

两个像素,如果中心像素小于其中任意一个,则舍弃该边缘像素点,否则保留。一

个简单的例子如下:


// 非最大信号压制算法 3x3
Arrays.fill(magnitudes, 0);
for (int row = 0; row < height; row++) {
	for (int col = 0; col < width; col++) {
		index = row * width + col;
		float angle = magnitudes[index];
		float m0 = data[index];
		magnitudes[index] = m0;
		if(angle >=0 && angle < 22.5) // angle 0
		{
			float m1 = getPixel(data, width, height, col-1, row);
			float m2 = getPixel(data, width, height, col+1, row);
			if(m0 < m1 || m0 < m2)
			{
				magnitudes[index] = 0;
			}
		}
		else if(angle >= 22.5 && angle < 67.5) // angle +45
		{
			float m1 = getPixel(data, width, height, col+1, row-1);
			float m2 = getPixel(data, width, height, col-1, row+1);
			if(m0 < m1 || m0 < m2)
			{
				magnitudes[index] = 0;
			}
		}
		else if(angle >= 67.5 && angle < 112.5) // angle 90
		{
			float m1 = getPixel(data, width, height, col, row+1);
			float m2 = getPixel(data, width, height, col, row-1);
			if(m0 < m1 || m0 < m2)
			{
				magnitudes[index] = 0;
			}
		}
		else if(angle >=112.5 && angle < 157.5) // angle 135 / -45
		{
			float m1 = getPixel(data, width, height, col-1, row-1);
			float m2 = getPixel(data, width, height, col+1, row+1);
			if(m0 < m1 || m0 < m2)
			{
				magnitudes[index] = 0;
			}
		}
		else if(angle >=157.5) // 跟零度是一致的,感谢一位网友发现了这个问题
		{
			float m1 = getPixel(data, width, height, col+1, row);
			float m2 = getPixel(data, width, height, col-1, row);
			if(m0 < m1 || m0 < m2)
			{
				magnitudes[index] = 0;
			}
		}
	}
}

1.      双阈值边缘连接

非最大信号压制以后,输出的幅值如果直接显示结果可能会少量的非边缘像素被包

含到结果中,所以要通过选取阈值进行取舍,传统的基于一个阈值的方法如果选择

的阈值较小起不到过滤非边缘的作用,如果选择的阈值过大容易丢失真正的图像边

缘,Canny提出基于双阈值(Fuzzy threshold)方法很好的实现了边缘选取,在实际

应用中双阈值还有边缘连接的作用。双阈值选择与边缘连接方法通过假设两个阈值

其中一个为高阈值TH另外一个为低阈值TL则有

a.      对于任意边缘像素低于TL的则丢弃

b.      对于任意边缘像素高于TH的则保留

c.      对于任意边缘像素值在TL与TH之间的,如果能通过边缘连接到一个像素大于

TH而且边缘所有像素大于最小阈值TL的则保留,否则丢弃。代码实现如下:

Arrays.fill(data, 0);
int offset = 0;
for (int row = 0; row < height; row++) {
	for (int col = 0; col < width; col++) {
		if(magnitudes[offset] >= highThreshold && data[offset] == 0)
		{
			edgeLink(col, row, offset, lowThreshold);
		}
		offset++;
	}
}
基于递归的边缘寻找方法edgeLink的代码如下:

private void edgeLink(int x1, int y1, int index, float threshold) {
	int x0 = (x1 == 0) ? x1 : x1 - 1;
	int x2 = (x1 == width - 1) ? x1 : x1 + 1;
	int y0 = y1 == 0 ? y1 : y1 - 1;
	int y2 = y1 == height -1 ? y1 : y1 + 1;
	
	data[index] = magnitudes[index];
	for (int x = x0; x <= x2; x++) {
		for (int y = y0; y <= y2; y++) {
			int i2 = x + y * width;
			if ((y != y1 || x != x1)
				&& data[i2] == 0 
				&& magnitudes[i2] >= threshold) {
				edgeLink(x, y, i2, threshold);
				return;
			}
		}
	}
}

6.      结果二值化显示 - 不说啦,直接点,自己看吧,太简单啦

// 二值化显示
for(int i=0; i<inPixels.length; i++)
{
	int gray = clamp((int)data[i]);
	outPixels[i] = gray > 0 ? -1 : 0xff000000;     
}
最终运行结果:


四:完整的Canny算法源代码

package com.gloomyfish.filter.study;

import java.awt.image.BufferedImage;
import java.util.Arrays;

public class CannyEdgeFilter extends AbstractBufferedImageOp {
	private float gaussianKernelRadius = 2f;
	private int gaussianKernelWidth = 16;
	private float lowThreshold;
	private float highThreshold;
	// image width, height
	private int width;
	private int height;
	private float[] data;
	private float[] magnitudes;

	public CannyEdgeFilter() {
		lowThreshold = 2.5f;
		highThreshold = 7.5f;
		gaussianKernelRadius = 2f;
		gaussianKernelWidth = 16;
	}

	public float getGaussianKernelRadius() {
		return gaussianKernelRadius;
	}

	public void setGaussianKernelRadius(float gaussianKernelRadius) {
		this.gaussianKernelRadius = gaussianKernelRadius;
	}

	public int getGaussianKernelWidth() {
		return gaussianKernelWidth;
	}

	public void setGaussianKernelWidth(int gaussianKernelWidth) {
		this.gaussianKernelWidth = gaussianKernelWidth;
	}

	public float getLowThreshold() {
		return lowThreshold;
	}

	public void setLowThreshold(float lowThreshold) {
		this.lowThreshold = lowThreshold;
	}

	public float getHighThreshold() {
		return highThreshold;
	}

	public void setHighThreshold(float highThreshold) {
		this.highThreshold = highThreshold;
	}

	@Override
	public BufferedImage filter(BufferedImage src, BufferedImage dest) {
		width = src.getWidth();
		height = src.getHeight();
		if (dest == null)
			dest = createCompatibleDestImage(src, null);
		// 图像灰度化
		int[] inPixels = new int[width * height];
		int[] outPixels = new int[width * height];
		getRGB(src, 0, 0, width, height, inPixels);
		int index = 0;
		for (int row = 0; row < height; row++) {
			int ta = 0, tr = 0, tg = 0, tb = 0;
			for (int col = 0; col < width; col++) {
				index = row * width + col;
				ta = (inPixels[index] >> 24) & 0xff;
				tr = (inPixels[index] >> 16) & 0xff;
				tg = (inPixels[index] >> 8) & 0xff;
				tb = inPixels[index] & 0xff;
				int gray = (int) (0.299 * tr + 0.587 * tg + 0.114 * tb);
				inPixels[index] = (ta << 24) | (gray << 16) | (gray << 8)
						| gray;
			}
		}
		
		// 计算高斯卷积核
		float kernel[][] = new float[gaussianKernelWidth][gaussianKernelWidth];
		for(int x=0; x<gaussianKernelWidth; x++)
		{
			for(int y=0; y<gaussianKernelWidth; y++)
			{
				kernel[x][y] = gaussian(x, y, gaussianKernelRadius);
			}
		}
		// 高斯模糊 -灰度图像
		int krr = (int)gaussianKernelRadius;
		for (int row = 0; row < height; row++) {
			for (int col = 0; col < width; col++) {
				index = row * width + col;
				double weightSum = 0.0;
				double redSum = 0;
				for(int subRow=-krr; subRow<=krr; subRow++)
				{
					int nrow = row + subRow;
					if(nrow >= height || nrow < 0)
					{
						nrow = 0;
					}
					for(int subCol=-krr; subCol<=krr; subCol++)
					{
						int ncol = col + subCol;
						if(ncol >= width || ncol <=0)
						{
							ncol = 0;
						}
						int index2 = nrow * width + ncol;
						int tr1 = (inPixels[index2] >> 16) & 0xff;
						redSum += tr1*kernel[subRow+krr][subCol+krr];
						weightSum += kernel[subRow+krr][subCol+krr];
					}
				}
				int gray = (int)(redSum / weightSum);
				outPixels[index] = gray;
			}
		}
		
		// 计算梯度-gradient, X放与Y方向
		data = new float[width * height];
		magnitudes = new float[width * height];
		for (int row = 0; row < height; row++) {
			for (int col = 0; col < width; col++) {
				index = row * width + col;
				// 计算X方向梯度
				float xg = (getPixel(outPixels, width, height, col, row+1) - 
						getPixel(outPixels, width, height, col, row) + 
						getPixel(outPixels, width, height, col+1, row+1) -
						getPixel(outPixels, width, height, col+1, row))/2.0f;
				float yg = (getPixel(outPixels, width, height, col, row)-
						getPixel(outPixels, width, height, col+1, row) +
						getPixel(outPixels, width, height, col, row+1) -
						getPixel(outPixels, width, height, col+1, row+1))/2.0f;
				// 计算振幅与角度
				data[index] = hypot(xg, yg);
				if(xg == 0)
				{
					if(yg > 0)
					{
						magnitudes[index]=90;						
					}
					if(yg < 0)
					{
						magnitudes[index]=-90;
					}
				}
				else if(yg == 0)
				{
					magnitudes[index]=0;
				}
				else
				{
					magnitudes[index] = (float)((Math.atan(yg/xg) * 180)/Math.PI);					
				}
				// make it 0 ~ 180
				magnitudes[index] += 90;
			}
		}
		
		// 非最大信号压制算法 3x3
		Arrays.fill(magnitudes, 0);
		for (int row = 0; row < height; row++) {
			for (int col = 0; col < width; col++) {
				index = row * width + col;
				float angle = magnitudes[index];
				float m0 = data[index];
				magnitudes[index] = m0;
				if(angle >=0 && angle < 22.5) // angle 0
				{
					float m1 = getPixel(data, width, height, col-1, row);
					float m2 = getPixel(data, width, height, col+1, row);
					if(m0 < m1 || m0 < m2)
					{
						magnitudes[index] = 0;
					}
				}
				else if(angle >= 22.5 && angle < 67.5) // angle +45
				{
					float m1 = getPixel(data, width, height, col+1, row-1);
					float m2 = getPixel(data, width, height, col-1, row+1);
					if(m0 < m1 || m0 < m2)
					{
						magnitudes[index] = 0;
					}
				}
				else if(angle >= 67.5 && angle < 112.5) // angle 90
				{
					float m1 = getPixel(data, width, height, col, row+1);
					float m2 = getPixel(data, width, height, col, row-1);
					if(m0 < m1 || m0 < m2)
					{
						magnitudes[index] = 0;
					}
				}
				else if(angle >=112.5 && angle < 157.5) // angle 135 / -45
				{
					float m1 = getPixel(data, width, height, col-1, row-1);
					float m2 = getPixel(data, width, height, col+1, row+1);
					if(m0 < m1 || m0 < m2)
					{
						magnitudes[index] = 0;
					}
				}
				else if(angle >=157.5) // angle 0
				{
					float m1 = getPixel(data, width, height, col, row+1);
					float m2 = getPixel(data, width, height, col, row-1);
					if(m0 < m1 || m0 < m2)
					{
						magnitudes[index] = 0;
					}
				}
			}
		}
		// 寻找最大与最小值
		float min = 255;
		float max = 0;
		for(int i=0; i<magnitudes.length; i++)
		{
			if(magnitudes[i] == 0) continue;
			min = Math.min(min, magnitudes[i]);
			max = Math.max(max, magnitudes[i]);
		}
		System.out.println("Image Max Gradient = " + max + " Mix Gradient = " + min);

		// 通常比值为 TL : TH = 1 : 3, 根据两个阈值完成二值化边缘连接
		// 边缘连接-link edges
		Arrays.fill(data, 0);
		int offset = 0;
		for (int row = 0; row < height; row++) {
			for (int col = 0; col < width; col++) {
				if(magnitudes[offset] >= highThreshold && data[offset] == 0)
				{
					edgeLink(col, row, offset, lowThreshold);
				}
				offset++;
			}
		}
		
		// 二值化显示
		for(int i=0; i<inPixels.length; i++)
		{
			int gray = clamp((int)data[i]);
			outPixels[i] = gray > 0 ? -1 : 0xff000000;     
		}
		setRGB(dest, 0, 0, width, height, outPixels );
		return dest;
	}
	
	public int clamp(int value) {
		return value > 255 ? 255 :
			(value < 0 ? 0 : value);
	}
	
	private void edgeLink(int x1, int y1, int index, float threshold) {
		int x0 = (x1 == 0) ? x1 : x1 - 1;
		int x2 = (x1 == width - 1) ? x1 : x1 + 1;
		int y0 = y1 == 0 ? y1 : y1 - 1;
		int y2 = y1 == height -1 ? y1 : y1 + 1;
		
		data[index] = magnitudes[index];
		for (int x = x0; x <= x2; x++) {
			for (int y = y0; y <= y2; y++) {
				int i2 = x + y * width;
				if ((y != y1 || x != x1)
					&& data[i2] == 0 
					&& magnitudes[i2] >= threshold) {
					edgeLink(x, y, i2, threshold);
					return;
				}
			}
		}
	}
	
	private float getPixel(float[] input, int width, int height, int col,
			int row) {
		if(col < 0 || col >= width)
			col = 0;
		if(row < 0 || row >= height)
			row = 0;
		int index = row * width + col;
		return input[index];
	}
	
	private float hypot(float x, float y) {
		return (float) Math.hypot(x, y);
	}
	
	private int getPixel(int[] inPixels, int width, int height, int col,
			int row) {
		if(col < 0 || col >= width)
			col = 0;
		if(row < 0 || row >= height)
			row = 0;
		int index = row * width + col;
		return inPixels[index];
	}
	
	private float gaussian(float x, float y, float sigma) {
		float xDistance = x*x;
		float yDistance = y*y;
		float sigma22 = 2*sigma*sigma;
		float sigma22PI = (float)Math.PI * sigma22;
		return (float)Math.exp(-(xDistance + yDistance)/sigma22)/sigma22PI;
	}

}
转载请务必注明出自本博客-gloomyfish

2017-04-19 19:54:00 The_star_is_at 阅读数 8472

实验三   图像轮廓提取与边缘检测


一、实验目的:

理解并掌握对二值图像进行轮廓提取的相关算法(比如,掏空内部点法),以及用于图像边缘检测和提取的典型微分算子(梯度算子和拉普拉斯算子)。

二、实验环境:

计算机、Windows XP操作系统,Matlab7.0

二、实验内容:

1、根据掏空内部点算法,运用Matlab编程实现二值图像的轮廓提取。

%以下适用于黑色背景白色前景的二值图像轮廓提取(以二值图像circles为例)

BW=imread('circles.png');     %二值图像circlesuint80黑,255

subplot(1,2,1);  imshow(BW);  title('二值图像');

[M, N]=size(BW);     %M行,N

BW_Buffer=BW;

for i=2: M-1

for j=2: N-1

if (BW(i, j)==255 & BW(i-1, j)==255 & BW(i+1, j)==255 & BW(i, j-1)==255 & BW(i, j+1)==255 & BW(i-1, j-1)==255 & BW(i-1, j+1)==255 & BW(i+1, j-1)==255 & BW(i+1, j+1)==255)    %说明BW(i, j)是前景中的一个内部点

       BW_Buffer(i, j)=0;    %掏空该内部点,即将该内部点置成与背景相同灰度

end

end

end

subplot(1,2,2);  imshow(BW_Buffer);  title('提取轮廓');


%以下适用于白色背景黑色前景的二值图像轮廓提取(以二值图像source为例)

BW=imread('source.bmp');      %二值图像sourceuint80黑,255

subplot(1,2,1);  imshow(BW);  title('二值图像');

[M, N]=size(BW);   %M行,N

BW_Buffer=BW;

for i=2: M-1

for j=2: N-1

if (BW(i, j)==0 & BW(i-1, j)==0 & BW(i+1, j)==0 & BW(i, j-1)==0 & BW(i, j+1)==0 & BW(i-1, j-1)==0 & BW(i-1, j+1)==0 & BW(i+1, j-1)==0 & BW(i+1, j+1)==0)      %说明BW(i, j)是前景中的一个内部点

       BW_Buffer(i, j)=255;   %掏空该内部点,即将该内部点置成与背景相同灰度

end

end

end

subplot(1,2,2);  imshow(BW_Buffer);  title('提取轮廓');

注意:使用掏空内部点的方法来提取二值图像的轮廓时,不能直接在原始二值图像矩阵上判断一个点掏空一个点,否则对前面像素的掏空操作会影响到对后面像素的判断结果。

解决方法:创建原始二值图像矩阵的副本(即图像矩阵BW_Buffer),在原始二值图像矩阵上执行判断操作,即依次判断每个像素点是否为前景中的内部点,如果是,则在图像矩阵BW_Buffer上执行掏空内部点的操作。

2、以灰度图像ricecameraman为例,利用Matlab图像处理工具箱中的edge函数,分别使用Roberts 算子、Sobel算子、Prewitt 算子对其进行边缘检测。

(1)函数格式: BW = edge(I, 'method', thresh)

(2)格式说明:edge函数输入灰度图像矩阵I,输出二值图像矩阵BW;参数'method'用于指定所使用的边缘检测算子,可以是'roberts''sobel''prewitt''log''canny';参数thresh用于指定梯度门限值(也称梯度阈值),图像中梯度值大于等于门限值thresh的像素用白色(1)表示,说明这些地方对应边缘,梯度值小于门限值thresh的像素用黑色(0)表示(edge function will ignore all edges that are not stronger than thresh)。若不指定参数thresh,则edge函数会自动选择阈值。所以edge函数最终将原始灰度图像中的边缘和背景用二值图像的形式展现出来,以突出边缘的位置,达到边缘检测的目的。

(3)程序如下:

I=imread('rice.png');

subplot(2,2,1);  imshow(I);  title('原始图像');

[BW1,thresh1]=edge(I,'roberts');  %进行Roberts算子边缘检测并返回门限值

[BW2,thresh2]=edge(I,'sobel');    %进行Sobel算子边缘检测并返回门限值

[BW3,thresh3]=edge(I,'prewitt');  %进行Prewitt算子边缘检测并返回门限值

subplot(2,2,2);  imshow(BW1);  title('Roberts算子边缘检测结果');

subplot(2,2,3);  imshow(BW2);  title('Sobel算子边缘检测结果');

subplot(2,2,4);  imshow(BW3);  title('Prewitt算子边缘检测结果');

若向原始图像中加入随机噪声(比如高斯噪声),之后再对噪声图像分别运用Roberts 算子、Sobel算子、Prewitt 算子、Log算子(高斯-拉普拉斯算子)进行边缘检测,观察检测结果,试比较4种边缘检测算子的抗噪声干扰能力。

I=imread('rice.png');

subplot(2,3,1);  imshow(I);  title('原始图像');

G=imnoise(I, 'gaussian');  %向原始图像中加入高斯噪声

subplot(2,3,2);  imshow(G);  title('噪声图像');

BW1=edge(G, 'roberts');  %进行Roberts算子边缘检测

BW2=edge(G, 'sobel');    %进行Sobel算子边缘检测

BW3=edge(G, 'prewitt');  %进行Prewitt算子边缘检测

BW4=edge(G, 'log');      %进行Log算子边缘检测

subplot(2,3,3);  imshow(BW1);  title('Roberts算子边缘检测结果');

subplot(2,3,4);  imshow(BW2);  title('Sobel算子边缘检测结果');

subplot(2,3,5);  imshow(BW3);  title('Prewitt算子边缘检测结果');

subplot(2,3,6);  imshow(BW4);  title('Log算子边缘检测结果');

2018-10-02 23:06:53 qq_36554582 阅读数 10006

MATLAB中有几种算法可以对图像进行边缘提取,其中一种就是edge算法,这个edge算法中有好几个算子,每一个算子分别对应着一种边缘提取的原理,接下来就来看一下几种方法的异同

%读取一张图片,并显示
original_picture=imread('D:\SoftWare\matlab2016a\Project\Picture\cat.jpg');
Pic2=im2bw(original_picture,thresh);
figure(1)
subplot(2,2,1);
imshow(original_picture);
title('原始RGB图像')
subplot(222)
imshow(Pic2)
title('二值化图像')

%用edge算法对二值化图像进行边缘提取
PicEdge1=edge(Pic2,'log');
subplot(223);
imshow(PicEdge1);
title('log算子')

PicEdge2 = edge(Pic2,'canny');
subplot(224);
imshow(PicEdge2);
title('canny算子');

PicEdge3=edge(Pic2,'sobel');
figure(2)
subplot(221)
imshow(PicEdge3);
title('sobel算子')

PicEdge4=edge(Pic2,'prewitt');
subplot(222)
imshow(PicEdge4);
title('sprewitt算子')

PicEdge5=edge(Pic2,'zerocross');
subplot(223)
imshow(PicEdge5);
title('zerocross算子')

PicEdge6=edge(Pic2,'roberts');
subplot(224)
imshow(PicEdge6);
title('roberts算子')



虽然我们从提取的结果来看,可能他们的差别不是很明显,但是这几个算子的基本原理还是有区别的,另外由于我采用的原始图片可能图中的猫和背景颜色有的部分很相似,所以会导致有些猫的边缘没有被提取出来,以后还需改善。

2011-11-05 09:54:50 Heaven13483 阅读数 6110
边缘提取以及边缘增强是不少图像处理软件都具有的基本功能,它的增强效果很明显,在用于识别的应用中,图像边缘也是非常重要的特征之一。图像边缘保留了原始图像中相当重要的部分信息,而又使得总的数据量减小了很多,这正符合特征提取的要求。在以后要谈到的霍夫变换(检测图像中的几何形状)中,边缘提取就是前提步骤。

这里我们只考虑灰度图像,用于图像识别的边缘提取比起仅仅用于视觉效果增强的边缘提取要复杂一些。要给图像的边缘下一个定义还挺困难的,从人的直观感受来说,边缘对应于物体的边界。图像上灰度变化剧烈的区域比较符合这个要求,我们一般会以这个特征来提取图像的边缘。但在遇到包含纹理的图像上,这有点问题,比如说,图像中的人穿了黑白格子的衣服,我们往往不希望提取出来的边缘包括衣服上的方格。但这个比较困难,涉及到纹理图像的处理等方法。

好了,既然边缘提取是要保留图像的灰度变化剧烈的区域,从数学上,最直观的方法就是微分(对于数字图像来说就是差分),在信号处理的角度来看,也可以说是用高通滤波器,即保留高频信号。这是最关键的一步,在此之前有时需要对输入图像进行消除噪声的处理。用于图像识别的边缘提取往往需要输出的边缘是二值图像,即只有黑白两个灰度的图像,其中


一个灰度代表边缘,另一个代表背景。此外,还需要把边缘细化成只有一个像素的宽度。总的说来边缘提取的步骤如下:

1,去噪声
2,微分运算
3,2值化处理
4,细化

第二步是关键,有不少书把第二步就直接称为边缘提取。实现它的算法也有很多,一般的图像处理教科书上都会介绍好几种,如拉普拉兹算子,索贝尔算子,罗伯特算子等等。这些都是模板运算,首先定义一个模板,模板的大小以3*3的较常见,也有2*2,5*5或更大尺寸的。运算时,把模板中心对应到图像的每一个像素位置,然后按照模板对应的公式对中心像素和它周围的像素进行数学运算,算出的结果作为输出图像对应像素点的值。

需要说明的是,模板运算是图像的一种处理手段--邻域处理,有许多图像增强效果都可以采用模板运算实现,如平滑效果,中值滤波(一种消除噪声的方法),油画效果,图像的凹凸效果等等。这些算法都比较简单,为人们常用。

关于前面提到的几种边缘提取算子(拉普拉兹算子,索贝尔算子,罗伯特算子),教科书上都有较为详细的介绍,我这里不多说了,(手头上没有教科书,也懒得翻译英文资料),如果你们有时间,可以把这些方法的具体情况仔细介绍一下。这里对拉普拉兹算子和索贝尔算子补充两句。拉普拉兹算子是2阶微分算子,也就是说,相当于求取2次微分,它的精度还算比较高,但对噪声过于敏感(有噪声的情况下效果很差)是它的重大缺点,所以这种算子并不是特别常用。索贝尔算子是最常用的算子之一(它是一种一阶算子),方法简单效果也不错,但提取出的边缘比较粗,要进行细化处理。另外,索贝尔算子
也可提取出图像边缘的方向信息来,有文章论证过,在不考虑噪声的情况下,它取得的边缘信息误差不超过7度。

顺便说一句,往往我们在进行边缘提取时只注意到位置信息,而忽略了边缘的方向。事实上,图像的边缘总有一定的走向,我们可以用边缘曲线的法线方向(和切线垂直的直线)来代表边缘点的方向。在图像识别的应用中,这个方向是非常重要的信息。


上面的几种算子是属于比较简单的方法,边缘提取的精度都不算特别高,下面介绍几种高级算法。首先是马尔(Marr)算子,马尔是计算机视觉这门学问的奠基人,很了不起,但这些理论很难懂。他提出的边缘提取方法可以看成两个步骤,一个是平滑作用来消除噪声,另一个是微分提取边缘,也可以说是由两个滤波器组成,低通滤波去除噪声,高通滤波提取边缘。人们也称这种方法为LOG滤波器,这也是根据它数学表达式和滤波器形状起的名字。也可以采用模板运算来实现这种算法,但模板的大小一般要在7*7以上,所以运算复杂程度比索贝尔算子等要大不少,运算时间当然也长许多。

另外一种非常重要的算法是坎尼 (Canny)算子,这是坎尼在1986年写的一篇论文里仔细论述的。他给出了判断边缘提取方法性能的指标。而坎尼算子也是图像处理领域里的标准方法,也可以说是默认的方法。比较奇怪的是,国内的图像处理教科书中,介绍坎尼算子的很少。本人见过的书中,郑南宁的' 计算机视觉与模式识别’(1998年),算是介绍的比较详细的。坎尼算子在使用时要提供给一些参数,用于控制算法的性能,实际上,对于不同的图像或不同的边缘提取目的,应该提供不同的参数,以达到最佳效果。它也有模板运算方法,模板的大小也比较大,和提供的参数有关,标准的大小差不多是17*17
,可以根据算子的可分离性用快速算法(否则就会慢的一塌糊涂),坎尼算子的2值化也很有特色,具有一定的智能性。

还有一种算法:Shen-Castan算子,大概可称为沈峻算子,总之是中国人的成果,效果和坎尼算子不相上下,这种算法在对边缘提取好坏的判别标准上有些不同。(这种方法我没用过,好象编起程序来,要比坎尼算子还复杂)

在实际的图像处理与识别应用中,有时需要根据被处理图像的种类以及实际目的,量身定做算法,边缘提取也是一样,但是基本原理都是一样的。
没有更多推荐了,返回首页