Canny邊緣檢測算法是澳洲科學家John F. Canny在1986年提出來的，不得不提一下的是當年John Canny本人才28歲！到今天已經30年過去了，Canny算法仍然是圖像邊緣檢測算法中最經典、有效的算法之一。
一起睹一下大家Canny的風采：
John Canny研究了最優邊緣檢測方法所需的特性，給出了評價邊緣檢測性能優劣的3個名額：
- 1 好的信噪比，即将非邊緣點判定為邊緣點的機率要低，将邊緣點判為非邊緣點的機率要低；
- 2 高的定位性能，即檢測出的邊緣點要盡可能在實際邊緣的中心；
- 3 對單一邊緣僅有唯一響應，即單個邊緣産生多個響應的機率要低，并且虛假響應邊緣應該得到最大抑制；
Canny算法就是基于滿足這3個名額的最優解實作的，在對圖像中物體邊緣敏感性的同時，也可以抑制或消除噪聲的影響。
Canny算子邊緣檢測的具體步驟如下：
- 一、用高斯濾波器平滑圖像
- 二、用Sobel等梯度算子計算梯度幅值和方向
- 三、對梯度幅值進行非極大值抑制
- 四、用雙門檻值算法檢測和連接配接邊緣
一、用高斯濾波器平滑圖像
高斯濾波是一種線性平滑濾波，适用于消除高斯噪聲，特别是對抑制或消除服從正态分布的噪聲非常有效。濾波可以消除或降低圖像中噪聲的影響，使用高斯濾波器主要是基于在濾波降噪的同時也可以最大限度保留邊緣資訊的考慮。
高斯濾波實作步驟：
1.1 彩色RGB圖像轉換為灰階圖像
邊緣檢測是基于對圖像灰階差異運算實作的，是以如果輸入的是RGB彩色圖像，需要先進行灰階圖的轉換。 RGB轉換成灰階圖像的一個常用公式是：
Gray = R*0.299 + G*0.587 + B*0.114
C++代碼實作起來也比較簡單，注意一般情況下圖像進行中彩色圖像各分量的排列順序是B、G、R。
//******************灰階轉換函數*************************
//第一個參數image輸入的彩色RGB圖像；
//第二個參數imageGray是轉換後輸出的灰階圖像；
//*************************************************************
void ConvertRGB2GRAY(const Mat &image,Mat &imageGray)
{
if(!image.data||image.channels()!=3)
{
return ;
}
imageGray=Mat::zeros(image.size(),CV_8UC1);
uchar *pointImage=image.data;
uchar *pointImageGray=imageGray.data;
int stepImage=image.step;
int stepImageGray=imageGray.step;
for(int i=0;i<imageGray.rows;i++)
{
for(int j=0;j<imageGray.cols;j++)
{
pointImageGray[i*stepImageGray+j]=0.114*pointImage[i*stepImage+3*j]+0.587*pointImage[i*stepImage+3*j+1]+0.299*pointImage[i*stepImage+3*j+2];
}
}
}
RGB原圖像： 轉換後的灰階圖：
1.2 生成高斯濾波卷積核
高斯濾波的過程是将灰階圖像跟高斯卷積核卷積，是以第一步是先要求解出給定尺寸和Sigma的高斯卷積核參數，以下代碼實作對卷積核參數求解：
//******************高斯卷積核生成函數*************************
//第一個參數gaus是一個指向含有N個double類型數組的指針；
//第二個參數size是高斯卷積核的尺寸大小；
//第三個參數sigma是卷積核的标準差
//*************************************************************
void GetGaussianKernel(double **gaus, const int size,const double sigma)
{
const double PI=4.0*atan(1.0); //圓周率π指派
int center=size/2;
double sum=0;
for(int i=0;i<size;i++)
{
for(int j=0;j<size;j++)
{
gaus[i][j]=(1/(2*PI*sigma*sigma))*exp(-((i-center)*(i-center)+(j-center)*(j-center))/(2*sigma*sigma));
sum+=gaus[i][j];
}
}
for(int i=0;i<size;i++)
{
for(int j=0;j<size;j++)
{
gaus[i][j]/=sum;
cout<<gaus[i][j]<<" ";
}
cout<<endl<<endl;
}
return ;
}
Sigma為1時，求得的3*3大小的高斯卷積核參數為：
Sigma為1,5*5大小的高斯卷積核參數為：
以下運算中Canny算子使用的是尺寸5*5，Sigma為1的高斯核。 更詳細的高斯濾波及高斯卷積核生成介紹可以移步這裡檢視：http://blog.csdn.net/dcrmg/article/details/52304446
1.3 高斯濾波
用在1.2中生成的高斯卷積核跟灰階圖像卷積，得到灰階圖像的高斯濾波後的圖像，抑制噪聲。 代碼實作：
//******************高斯濾波*************************
//第一個參數imageSource是待濾波原始圖像；
//第二個參數imageGaussian是濾波後輸出圖像；
//第三個參數gaus是一個指向含有N個double類型數組的指針；
//第四個參數size是濾波核的尺寸
//*************************************************************
void GaussianFilter(const Mat imageSource,Mat &imageGaussian,double **gaus,int size)
{
imageGaussian=Mat::zeros(imageSource.size(),CV_8UC1);
if(!imageSource.data||imageSource.channels()!=1)
{
return ;
}
double gausArray[100];
for(int i=0;i<size*size;i++)
{
gausArray[i]=0; //賦初值，空間配置設定
}
int array=0;
for(int i=0;i<size;i++)
{
for(int j=0;j<size;j++)
{
gausArray[array]=gaus[i][j];//二維數組到一維 友善計算
array++;
}
}
//濾波
for(int i=0;i<imageSource.rows;i++)
{
for(int j=0;j<imageSource.cols;j++)
{
int k=0;
for(int l=-size/2;l<=size/2;l++)
{
for(int g=-size/2;g<=size/2;g++)
{
//以下處理針對濾波後圖像邊界處理，為超出邊界的值指派為邊界值
int row=i+l;
int col=j+g;
row=row<0?0:row;
row=row>=imageSource.rows?imageSource.rows-1:row;
col=col<0?0:col;
col=col>=imageSource.cols?imageSource.cols-1:col;
//卷積和
imageGaussian.at<uchar>(i,j)+=gausArray[k]*imageSource.at<uchar>(row,col);
k++;
}
}
}
}
}
高斯濾波後的圖像：
跟原圖相比，圖像有一定程度的模糊。
- 二、用Sobel等梯度算子計算梯度幅值和方向
圖像灰階值的梯度可以使用最簡單的一階有限差分來進行近似，使用以下圖像在x和y方向上偏導數的兩個矩陣：
計算公式為：
其中f為圖像灰階值，P代表X方向梯度幅值，Q代表Y方向 梯度幅值，M是該點幅值，Θ是梯度方向，也就是角度。
2.1 Sobel卷積核計算X、Y方向梯度和梯度角
這裡使用較為常用的Sobel算子計算X和Y方向上的梯度以及梯度的方向角，Sobel的X和Y方向的卷積因子為：
更多關于Sobel算子的介紹可以移步這裡檢視：http://blog.csdn.net/dcrmg/article/details/52280768
使用Sobel卷積因子計算X、Y方向梯度和梯度方向角代碼實作：
//******************Sobel卷積因子計算X、Y方向梯度和梯度方向角********************
//第一個參數imageSourc原始灰階圖像；
//第二個參數imageSobelX是X方向梯度圖像；
//第三個參數imageSobelY是Y方向梯度圖像；
//第四個參數pointDrection是梯度方向角數組指針
//*************************************************************
void SobelGradDirction(const Mat imageSource,Mat &imageSobelX,Mat &imageSobelY,double *&pointDrection)
{
pointDrection=new double[(imageSource.rows-1)*(imageSource.cols-1)];
for(int i=0;i<(imageSource.rows-1)*(imageSource.cols-1);i++)
{
pointDrection[i]=0;
}
imageSobelX=Mat::zeros(imageSource.size(),CV_32SC1);
imageSobelY=Mat::zeros(imageSource.size(),CV_32SC1);
uchar *P=imageSource.data;
uchar *PX=imageSobelX.data;
uchar *PY=imageSobelY.data;
int step=imageSource.step;
int stepXY=imageSobelX.step;
int k=0;
int m=0;
int n=0;
for(int i=1;i<(imageSource.rows-1);i++)
{
for(int j=1;j<(imageSource.cols-1);j++)
{
//通過指針周遊圖像上每一個像素
double gradY=P[(i-1)*step+j+1]+P[i*step+j+1]*2+P[(i+1)*step+j+1]-P[(i-1)*step+j-1]-P[i*step+j-1]*2-P[(i+1)*step+j-1];
PY[i*stepXY+j*(stepXY/step)]=abs(gradY);
double gradX=P[(i+1)*step+j-1]+P[(i+1)*step+j]*2+P[(i+1)*step+j+1]-P[(i-1)*step+j-1]-P[(i-1)*step+j]*2-P[(i-1)*step+j+1];
PX[i*stepXY+j*(stepXY/step)]=abs(gradX);
if(gradX==0)
{
gradX=0.00000000000000001; //防止除數為0異常
}
pointDrection[k]=atan(gradY/gradX)*57.3;//弧度轉換為度
pointDrection[k]+=90;
k++;
}
}
convertScaleAbs(imageSobelX,imageSobelX);
convertScaleAbs(imageSobelY,imageSobelY);
}
數組指針pointDirection裡存放了每個點上的梯度方向角，以度為機關。由于atan求得的角度範圍是-π/2~π/2，為了便于計算，這裡對每個梯度角加了一個π/2，使範圍變成0~π，便于計算。
X方向梯度圖： Y方向梯度圖：
2.2 求梯度圖的幅值
求得X、Y方向的梯度和梯度角之後再來計算X和Y方向融合的梯度幅值，計算公式為：
代碼實作，較為簡單：
//******************計算Sobel的X和Y方向梯度幅值*************************
//第一個參數imageGradX是X方向梯度圖像；
//第二個參數imageGradY是Y方向梯度圖像；
//第三個參數SobelAmpXY是輸出的X、Y方向梯度圖像幅值
//*************************************************************
void SobelAmplitude(const Mat imageGradX,const Mat imageGradY,Mat &SobelAmpXY)
{
SobelAmpXY=Mat::zeros(imageGradX.size(),CV_32FC1);
for(int i=0;i<SobelAmpXY.rows;i++)
{
for(int j=0;j<SobelAmpXY.cols;j++)
{
SobelAmpXY.at<float>(i,j)=sqrt(imageGradX.at<uchar>(i,j)*imageGradX.at<uchar>(i,j)+imageGradY.at<uchar>(i,j)*imageGradY.at<uchar>(i,j));
}
}
convertScaleAbs(SobelAmpXY,SobelAmpXY);
}
求得的X和Y方向幅度和疊加了兩個方向上的幅值：
- 三、對梯度幅值進行非極大值抑制
求幅值圖像進行非極大值抑制，可以進一步消除非邊緣的噪點，更重要的是，可以細化邊緣。 抑制邏輯是：沿着該點梯度方向，比較前後兩個點的幅值大小，若該點大于前後兩點，則保留，若該點小于前後兩點，則置為0； 示意圖如下：
圖中四條虛線代表圖像中每一點可能的梯度方向，沿着梯度方向與邊界的上下兩個交點，就是需要拿來與中心點點（X0，Y0）做比較的點。交點值的計算采用插值法計算，以黃色的虛線所代表的梯度角Θ為例，交點處幅值為： P(X0,Y0)+（P(X0-1,Y0+1)-P(X0,Y0)）*tan(Θ)
四種情況下需要分别計算，代碼實作如下：
//******************局部極大值抑制*************************
//第一個參數imageInput輸入的Sobel梯度圖像；
//第二個參數imageOutPut是輸出的局部極大值抑制圖像；
//第三個參數pointDrection是圖像上每個點的梯度方向數組指針
//*************************************************************
void LocalMaxValue(const Mat imageInput,Mat &imageOutput,double *pointDrection)
{
//imageInput.copyTo(imageOutput);
imageOutput=imageInput.clone();
int k=0;
for(int i=1;i<imageInput.rows-1;i++)
{
for(int j=1;j<imageInput.cols-1;j++)
{
int value00=imageInput.at<uchar>((i-1),j-1);
int value01=imageInput.at<uchar>((i-1),j);
int value02=imageInput.at<uchar>((i-1),j+1);
int value10=imageInput.at<uchar>((i),j-1);
int value11=imageInput.at<uchar>((i),j);
int value12=imageInput.at<uchar>((i),j+1);
int value20=imageInput.at<uchar>((i+1),j-1);
int value21=imageInput.at<uchar>((i+1),j);
int value22=imageInput.at<uchar>((i+1),j+1);
if(pointDrection[k]>0&&pointDrection[k]<=45)
{
if(value11<=(value12+(value02-value12)*tan(pointDrection[i*imageOutput.rows+j]))||(value11<=(value10+(value20-value10)*tan(pointDrection[i*imageOutput.rows+j]))))
{
imageOutput.at<uchar>(i,j)=0;
}
}
if(pointDrection[k]>45&&pointDrection[k]<=90)
{
if(value11<=(value01+(value02-value01)/tan(pointDrection[i*imageOutput.cols+j]))||value11<=(value21+(value20-value21)/tan(pointDrection[i*imageOutput.cols+j])))
{
imageOutput.at<uchar>(i,j)=0;
}
}
if(pointDrection[k]>90&&pointDrection[k]<=135)
{
if(value11<=(value01+(value00-value01)/tan(180-pointDrection[i*imageOutput.cols+j]))||value11<=(value21+(value22-value21)/tan(180-pointDrection[i*imageOutput.cols+j])))
{
imageOutput.at<uchar>(i,j)=0;
}
}
if(pointDrection[k]>135&&pointDrection[k]<=180)
{
if(value11<=(value10+(value00-value10)*tan(180-pointDrection[i*imageOutput.cols+j]))||value11<=(value12+(value22-value11)*tan(180-pointDrection[i*imageOutput.cols+j])))
{
imageOutput.at<uchar>(i,j)=0;
}
}
k++;
}
}
}
進過非極大值抑制後的圖像如下：
跟Sobel梯度幅值圖像相比，去除了一些非局部極大值點，輪廓也進一步細化。
- 四、用雙門檻值算法檢測和連接配接邊緣
雙門檻值的機理是：
指定一個低門檻值A，一個高門檻值B，一般取B為圖像整體灰階級分布的70%，且B為1.5到2倍大小的A；
灰階值大于B的，置為255，灰階值小于A的，置為0；
灰階值介于A和B之間的，考察改像素點臨近的8像素是否有灰階值為255的，若沒有255的，表示這是一個孤立的局部極大值點，予以排除，置為0；若有255的，表示這是一個跟其他邊緣有“接壤”的可造之材，置為255，之後重複執行該步驟，直到考察完之後一個像素點。
4.1 雙門檻值處理
這個步驟裡處理大于高門檻值和小于低門檻值的像素點，分别置為255和0；
實作如下：
//******************雙門檻值處理*************************
//第一個參數imageInput輸入和輸出的的Sobel梯度幅值圖像；
//第二個參數lowThreshold是低門檻值
//第三個參數highThreshold是高門檻值
//******************************************************
void DoubleThreshold(Mat &imageIput,double lowThreshold,double highThreshold)
{
for(int i=0;i<imageIput.rows;i++)
{
for(int j=0;j<imageIput.cols;j++)
{
if(imageIput.at<uchar>(i,j)>highThreshold)
{
imageIput.at<uchar>(i,j)=255;
}
if(imageIput.at<uchar>(i,j)<lowThreshold)
{
imageIput.at<uchar>(i,j)=0;
}
}
}
}
這裡取低門檻值為60，高門檻值100，處理效果：
經過雙門檻值處理後，灰階值較低的點被消除掉，較高的點置為了255。上圖中仍有灰階值小于255的點，它們是介于高、低門檻值間的像素點。
4.2 雙門檻值中間像素濾除或連接配接處理
以下是C++編碼實作，其中在連接配接的操作上涉及到一個遞歸調用：
//******************雙門檻值中間像素連接配接處理*********************
//第一個參數imageInput輸入和輸出的的Sobel梯度幅值圖像；
//第二個參數lowThreshold是低門檻值
//第三個參數highThreshold是高門檻值
//*************************************************************
void DoubleThresholdLink(Mat &imageInput,double lowThreshold,double highThreshold)
{
for(int i=1;i<imageInput.rows-1;i++)
{
for(int j=1;j<imageInput.cols-1;j++)
{
if(imageInput.at<uchar>(i,j)>lowThreshold&&imageInput.at<uchar>(i,j)<255)
{
if(imageInput.at<uchar>(i-1,j-1)==255||imageInput.at<uchar>(i-1,j)==255||imageInput.at<uchar>(i-1,j+1)==255||
imageInput.at<uchar>(i,j-1)==255||imageInput.at<uchar>(i,j)==255||imageInput.at<uchar>(i,j+1)==255||
imageInput.at<uchar>(i+1,j-1)==255||imageInput.at<uchar>(i+1,j)==255||imageInput.at<uchar>(i+1,j+1)==255)
{
imageInput.at<uchar>(i,j)=255;
DoubleThresholdLink(imageInput,lowThreshold,highThreshold); //遞歸調用
}
else
{
imageInput.at<uchar>(i,j)=0;
}
}
}
}
}
濾除或連接配接後的最終效果：
完整工程代碼
到這裡，Canny算子檢測邊緣的四個步驟就全部完成了，以下是整個C++工程完整代碼，有興趣可以浏覽一下：
#include "core/core.hpp"
#include "highgui/highgui.hpp"
#include "imgproc/imgproc.hpp"
#include "iostream"
#include "math.h"
using namespace std;
using namespace cv;
//******************灰階轉換函數*************************
//第一個參數image輸入的彩色RGB圖像；
//第二個參數imageGray是轉換後輸出的灰階圖像；
//*************************************************************
void ConvertRGB2GRAY(const Mat &image,Mat &imageGray);
//******************高斯卷積核生成函數*************************
//第一個參數gaus是一個指向含有N個double類型數組的指針；
//第二個參數size是高斯卷積核的尺寸大小；
//第三個參數sigma是卷積核的标準差
//*************************************************************
void GetGaussianKernel(double **gaus, const int size,const double sigma);
//******************高斯濾波*************************
//第一個參數imageSource是待濾波原始圖像；
//第二個參數imageGaussian是濾波後輸出圖像；
//第三個參數gaus是一個指向含有N個double類型數組的指針；
//第四個參數size是濾波核的尺寸
//*************************************************************
void GaussianFilter(const Mat imageSource,Mat &imageGaussian,double **gaus,int size);
//******************Sobel算子計算梯度和方向********************
//第一個參數imageSourc原始灰階圖像；
//第二個參數imageSobelX是X方向梯度圖像；
//第三個參數imageSobelY是Y方向梯度圖像；
//第四個參數pointDrection是梯度方向數組指針
//*************************************************************
void SobelGradDirction(const Mat imageSource,Mat &imageSobelX,Mat &imageSobelY,double *&pointDrection);
//******************計算Sobel的X和Y方向梯度幅值*************************
//第一個參數imageGradX是X方向梯度圖像；
//第二個參數imageGradY是Y方向梯度圖像；
//第三個參數SobelAmpXY是輸出的X、Y方向梯度圖像幅值
//*************************************************************
void SobelAmplitude(const Mat imageGradX,const Mat imageGradY,Mat &SobelAmpXY);
//******************局部極大值抑制*************************
//第一個參數imageInput輸入的Sobel梯度圖像；
//第二個參數imageOutPut是輸出的局部極大值抑制圖像；
//第三個參數pointDrection是圖像上每個點的梯度方向數組指針
//*************************************************************
void LocalMaxValue(const Mat imageInput,Mat &imageOutput,double *pointDrection);
//******************雙門檻值處理*************************
//第一個參數imageInput輸入和輸出的的Sobel梯度幅值圖像；
//第二個參數lowThreshold是低門檻值
//第三個參數highThreshold是高門檻值
//******************************************************
void DoubleThreshold(Mat &imageIput,double lowThreshold,double highThreshold);
//******************雙門檻值中間像素連接配接處理*********************
//第一個參數imageInput輸入和輸出的的Sobel梯度幅值圖像；
//第二個參數lowThreshold是低門檻值
//第三個參數highThreshold是高門檻值
//*************************************************************
void DoubleThresholdLink(Mat &imageInput,double lowThreshold,double highThreshold);
Mat imageSource;
Mat imageGray;
Mat imageGaussian;
int main(int argc,char *argv[])
{
imageSource=imread(argv[1]); //讀入RGB圖像
imshow("RGB Image",imageSource);
ConvertRGB2GRAY(imageSource,imageGray); //RGB轉換為灰階圖
imshow("Gray Image",imageGray);
int size=5; //定義卷積核大小
double **gaus=new double *[size]; //卷積核數組
for(int i=0;i<size;i++)
{
gaus[i]=new double[size]; //動态生成矩陣
}
GetGaussianKernel(gaus,5,1); //生成5*5 大小高斯卷積核，Sigma=1；
imageGaussian=Mat::zeros(imageGray.size(),CV_8UC1);
GaussianFilter(imageGray,imageGaussian,gaus,5); //高斯濾波
imshow("Gaussian Image",imageGaussian);
Mat imageSobelY;
Mat imageSobelX;
double *pointDirection=new double[(imageSobelX.cols-1)*(imageSobelX.rows-1)]; //定義梯度方向角數組
SobelGradDirction(imageGaussian,imageSobelX,imageSobelY,pointDirection); //計算X、Y方向梯度和方向角
imshow("Sobel Y",imageSobelY);
imshow("Sobel X",imageSobelX);
Mat SobelGradAmpl;
SobelAmplitude(imageSobelX,imageSobelY,SobelGradAmpl); //計算X、Y方向梯度融合幅值
imshow("Soble XYRange",SobelGradAmpl);
Mat imageLocalMax;
LocalMaxValue(SobelGradAmpl,imageLocalMax,pointDirection); //局部非極大值抑制
imshow("Non-Maximum Image",imageLocalMax);
Mat cannyImage;
cannyImage=Mat::zeros(imageLocalMax.size(),CV_8UC1);
DoubleThreshold(imageLocalMax,90,160); //雙門檻值處理
imshow("Double Threshold Image",imageLocalMax);
DoubleThresholdLink(imageLocalMax,90,160); //雙門檻值中間門檻值濾除及連接配接
imshow("Canny Image",imageLocalMax);
waitKey();
system("pause");
return 0;
}
//******************高斯卷積核生成函數*************************
//第一個參數gaus是一個指向含有N個double類型數組的指針；
//第二個參數size是高斯卷積核的尺寸大小；
//第三個參數sigma是卷積核的标準差
//*************************************************************
void GetGaussianKernel(double **gaus, const int size,const double sigma)
{
const double PI=4.0*atan(1.0); //圓周率π指派
int center=size/2;
double sum=0;
for(int i=0;i<size;i++)
{
for(int j=0;j<size;j++)
{
gaus[i][j]=(1/(2*PI*sigma*sigma))*exp(-((i-center)*(i-center)+(j-center)*(j-center))/(2*sigma*sigma));
sum+=gaus[i][j];
}
}
for(int i=0;i<size;i++)
{
for(int j=0;j<size;j++)
{
gaus[i][j]/=sum;
cout<<gaus[i][j]<<" ";
}
cout<<endl<<endl;
}
return ;
}
//******************灰階轉換函數*************************
//第一個參數image輸入的彩色RGB圖像；
//第二個參數imageGray是轉換後輸出的灰階圖像；
//*************************************************************
void ConvertRGB2GRAY(const Mat &image,Mat &imageGray)
{
if(!image.data||image.channels()!=3)
{
return ;
}
imageGray=Mat::zeros(image.size(),CV_8UC1);
uchar *pointImage=image.data;
uchar *pointImageGray=imageGray.data;
int stepImage=image.step;
int stepImageGray=imageGray.step;
for(int i=0;i<imageGray.rows;i++)
{
for(int j=0;j<imageGray.cols;j++)
{
pointImageGray[i*stepImageGray+j]=0.114*pointImage[i*stepImage+3*j]+0.587*pointImage[i*stepImage+3*j+1]+0.299*pointImage[i*stepImage+3*j+2];
}
}
}
//******************高斯濾波*************************
//第一個參數imageSource是待濾波原始圖像；
//第二個參數imageGaussian是濾波後輸出圖像；
//第三個參數gaus是一個指向含有N個double類型數組的指針；
//第四個參數size是濾波核的尺寸
//*************************************************************
void GaussianFilter(const Mat imageSource,Mat &imageGaussian,double **gaus,int size)
{
imageGaussian=Mat::zeros(imageSource.size(),CV_8UC1);
if(!imageSource.data||imageSource.channels()!=1)
{
return ;
}
double gausArray[100];
for(int i=0;i<size*size;i++)
{
gausArray[i]=0; //賦初值，空間配置設定
}
int array=0;
for(int i=0;i<size;i++)
{
for(int j=0;j<size;j++)
{
gausArray[array]=gaus[i][j];//二維數組到一維 友善計算
array++;
}
}
//濾波
for(int i=0;i<imageSource.rows;i++)
{
for(int j=0;j<imageSource.cols;j++)
{
int k=0;
for(int l=-size/2;l<=size/2;l++)
{
for(int g=-size/2;g<=size/2;g++)
{
//以下處理針對濾波後圖像邊界處理，為超出邊界的值指派為邊界值
int row=i+l;
int col=j+g;
row=row<0?0:row;
row=row>=imageSource.rows?imageSource.rows-1:row;
col=col<0?0:col;
col=col>=imageSource.cols?imageSource.cols-1:col;
//卷積和
imageGaussian.at<uchar>(i,j)+=gausArray[k]*imageSource.at<uchar>(row,col);
k++;
}
}
}
}
}
//******************Sobel算子計算X、Y方向梯度和梯度方向角********************
//第一個參數imageSourc原始灰階圖像；
//第二個參數imageSobelX是X方向梯度圖像；
//第三個參數imageSobelY是Y方向梯度圖像；
//第四個參數pointDrection是梯度方向角數組指針
//*************************************************************
void SobelGradDirction(const Mat imageSource,Mat &imageSobelX,Mat &imageSobelY,double *&pointDrection)
{
pointDrection=new double[(imageSource.rows-1)*(imageSource.cols-1)];
for(int i=0;i<(imageSource.rows-1)*(imageSource.cols-1);i++)
{
pointDrection[i]=0;
}
imageSobelX=Mat::zeros(imageSource.size(),CV_32SC1);
imageSobelY=Mat::zeros(imageSource.size(),CV_32SC1);
uchar *P=imageSource.data;
uchar *PX=imageSobelX.data;
uchar *PY=imageSobelY.data;
int step=imageSource.step;
int stepXY=imageSobelX.step;
int k=0;
int m=0;
int n=0;
for(int i=1;i<(imageSource.rows-1);i++)
{
for(int j=1;j<(imageSource.cols-1);j++)
{
//通過指針周遊圖像上每一個像素
double gradY=P[(i-1)*step+j+1]+P[i*step+j+1]*2+P[(i+1)*step+j+1]-P[(i-1)*step+j-1]-P[i*step+j-1]*2-P[(i+1)*step+j-1];
PY[i*stepXY+j*(stepXY/step)]=abs(gradY);
double gradX=P[(i+1)*step+j-1]+P[(i+1)*step+j]*2+P[(i+1)*step+j+1]-P[(i-1)*step+j-1]-P[(i-1)*step+j]*2-P[(i-1)*step+j+1];
PX[i*stepXY+j*(stepXY/step)]=abs(gradX);
if(gradX==0)
{
gradX=0.00000000000000001; //防止除數為0異常
}
pointDrection[k]=atan(gradY/gradX)*57.3;//弧度轉換為度
pointDrection[k]+=90;
k++;
}
}
convertScaleAbs(imageSobelX,imageSobelX);
convertScaleAbs(imageSobelY,imageSobelY);
}
//******************計算Sobel的X和Y方向梯度幅值*************************
//第一個參數imageGradX是X方向梯度圖像；
//第二個參數imageGradY是Y方向梯度圖像；
//第三個參數SobelAmpXY是輸出的X、Y方向梯度圖像幅值
//*************************************************************
void SobelAmplitude(const Mat imageGradX,const Mat imageGradY,Mat &SobelAmpXY)
{
SobelAmpXY=Mat::zeros(imageGradX.size(),CV_32FC1);
for(int i=0;i<SobelAmpXY.rows;i++)
{
for(int j=0;j<SobelAmpXY.cols;j++)
{
SobelAmpXY.at<float>(i,j)=sqrt(imageGradX.at<uchar>(i,j)*imageGradX.at<uchar>(i,j)+imageGradY.at<uchar>(i,j)*imageGradY.at<uchar>(i,j));
}
}
convertScaleAbs(SobelAmpXY,SobelAmpXY);
}
//******************局部極大值抑制*************************
//第一個參數imageInput輸入的Sobel梯度圖像；
//第二個參數imageOutPut是輸出的局部極大值抑制圖像；
//第三個參數pointDrection是圖像上每個點的梯度方向數組指針
//*************************************************************
void LocalMaxValue(const Mat imageInput,Mat &imageOutput,double *pointDrection)
{
//imageInput.copyTo(imageOutput);
imageOutput=imageInput.clone();
int k=0;
for(int i=1;i<imageInput.rows-1;i++)
{
for(int j=1;j<imageInput.cols-1;j++)
{
int value00=imageInput.at<uchar>((i-1),j-1);
int value01=imageInput.at<uchar>((i-1),j);
int value02=imageInput.at<uchar>((i-1),j+1);
int value10=imageInput.at<uchar>((i),j-1);
int value11=imageInput.at<uchar>((i),j);
int value12=imageInput.at<uchar>((i),j+1);
int value20=imageInput.at<uchar>((i+1),j-1);
int value21=imageInput.at<uchar>((i+1),j);
int value22=imageInput.at<uchar>((i+1),j+1);
if(pointDrection[k]>0&&pointDrection[k]<=45)
{
if(value11<=(value12+(value02-value12)*tan(pointDrection[i*imageOutput.rows+j]))||(value11<=(value10+(value20-value10)*tan(pointDrection[i*imageOutput.rows+j]))))
{
imageOutput.at<uchar>(i,j)=0;
}
}
if(pointDrection[k]>45&&pointDrection[k]<=90)
{
if(value11<=(value01+(value02-value01)/tan(pointDrection[i*imageOutput.cols+j]))||value11<=(value21+(value20-value21)/tan(pointDrection[i*imageOutput.cols+j])))
{
imageOutput.at<uchar>(i,j)=0;
}
}
if(pointDrection[k]>90&&pointDrection[k]<=135)
{
if(value11<=(value01+(value00-value01)/tan(180-pointDrection[i*imageOutput.cols+j]))||value11<=(value21+(value22-value21)/tan(180-pointDrection[i*imageOutput.cols+j])))
{
imageOutput.at<uchar>(i,j)=0;
}
}
if(pointDrection[k]>135&&pointDrection[k]<=180)
{
if(value11<=(value10+(value00-value10)*tan(180-pointDrection[i*imageOutput.cols+j]))||value11<=(value12+(value22-value11)*tan(180-pointDrection[i*imageOutput.cols+j])))
{
imageOutput.at<uchar>(i,j)=0;
}
}
k++;
}
}
}
//******************雙門檻值處理*************************
//第一個參數imageInput輸入和輸出的的Sobel梯度幅值圖像；
//第二個參數lowThreshold是低門檻值
//第三個參數highThreshold是高門檻值
//******************************************************
void DoubleThreshold(Mat &imageIput,double lowThreshold,double highThreshold)
{
for(int i=0;i<imageIput.rows;i++)
{
for(int j=0;j<imageIput.cols;j++)
{
if(imageIput.at<uchar>(i,j)>highThreshold)
{
imageIput.at<uchar>(i,j)=255;
}
if(imageIput.at<uchar>(i,j)<lowThreshold)
{
imageIput.at<uchar>(i,j)=0;
}
}
}
}
//******************雙門檻值中間像素連接配接處理*********************
//第一個參數imageInput輸入和輸出的的Sobel梯度幅值圖像；
//第二個參數lowThreshold是低門檻值
//第三個參數highThreshold是高門檻值
//*************************************************************
void DoubleThresholdLink(Mat &imageInput,double lowThreshold,double highThreshold)
{
for(int i=1;i<imageInput.rows-1;i++)
{
for(int j=1;j<imageInput.cols-1;j++)
{
if(imageInput.at<uchar>(i,j)>lowThreshold&&imageInput.at<uchar>(i,j)<255)
{
if(imageInput.at<uchar>(i-1,j-1)==255||imageInput.at<uchar>(i-1,j)==255||imageInput.at<uchar>(i-1,j+1)==255||
imageInput.at<uchar>(i,j-1)==255||imageInput.at<uchar>(i,j)==255||imageInput.at<uchar>(i,j+1)==255||
imageInput.at<uchar>(i+1,j-1)==255||imageInput.at<uchar>(i+1,j)==255||imageInput.at<uchar>(i+1,j+1)==255)
{
imageInput.at<uchar>(i,j)=255;
DoubleThresholdLink(imageInput,lowThreshold,highThreshold); //遞歸調用
}
else
{
imageInput.at<uchar>(i,j)=0;
}
}
}
}
}