Marr-Hildreth邊緣檢測器C++實現
Marr-Hiddreth是基於以下兩個事實的:
1 灰度變換與影象尺寸無關,因此邊緣檢測可以使用不同的尺寸運算元
2 灰度的突然變換會在一階導數中導致波峰或波谷,或在二階導數中等效的引起零交叉
用於邊緣檢測的運算元應該有兩個明顯的特點:
1 它應該是一個能計算影象中每一點處的一階導數或二階導數的數字近似的微分運算元
2 它應該能被調整以便在任何期望的尺寸上起作用,因此大運算元也可用於檢測模糊邊緣,小運算元可用於檢測銳度集中的精細細節
滿足這些條件的最令人滿意的運算元是高斯拉普拉斯運算元(LoG)
(1)
(2)
圖1
圖1(a)到(c)顯示了一個LoG的負函式的三維圖,影象和剖面,LoG的零交叉出現在x2+y2=
處,圖(d)顯示的是5*5的模板,它近似於(a)的形狀,該近似不唯一
任意尺寸的模板可以通過式(2)取樣並標定係數以使係數之和為0來生成,生成LoG濾波器的的一種更有效的方式是以n*n對式(1)進行取樣,然後將結果陣列與一個拉普拉斯模板進行卷積.
Marr-Hildreth演算法是由LoG濾波器與一副影象的卷積組成
(3)
1 對式(1)進行取樣得到n*n的高斯低通濾波器對輸入影象濾波
2 計算第一步得到的影象的拉普拉斯
3 找到步驟二所得的影象的零交叉
1.3 對高斯函式進行取樣的細節
二維高斯函式的公式:
(4)
首先我們確定
和高斯核矩陣的大小,一般我們取高斯核矩陣的大小為6
+1, 比如,
=4,size=25;
之後我們對size*size的矩陣進行迴圈,對應的左邊i,j作為x,y進行計算,不過要注意的是矩陣的中心位置是對應的x,y對應的是(0,0),所以我們需要適當的變換
再計算出矩陣所有位置對應的f(x,y)之後,因為高斯值均為小數,為了規範,我們將其轉換為整數,將(0,0)為值的值設為1,其他位置看與(0,0)位置的比值進行設定,轉化為整數後,再對核矩陣進行均值化
//構造高斯核矩陣
void GaussianKernel(cv::Mat &ss,int size,double delta){
int centerX=size/2;
int centerY=size/2;
int x=0;
int y=0;
cv::Mat gaussiankernel;
gaussiankernel=cv::Mat_<double>(size,size);
double cc=0.5/(3.1415926*delta*delta);
for(int i=0;i<size;i++){
for(int j=0;j<size;j++){
x=std::abs(i-centerX);
y=std::abs(j-centerY);
double dd=cc*exp(-((pow(x,2)+pow(y,2))/(2*delta*delta)));
gaussiankernel.at<double>(i,j)=dd;
}
}
cv::Mat ss1;
ss1=cv::Mat_<int>(size,size);
int sum=0;
//將高斯變換的小數型別轉為整型
double ff=gaussiankernel.at<double>(0,0);
for(int i=0;i<size;i++){
for(int j=0;j<size;j++){
int cc=cvRound(gaussiankernel.at<double>(i,j)/ff);
ss1.at<int>(i,j)=cc;
sum+=cc;
}
}
//均值化
for(int i=0;i<size;i++){
for(int j=0;j<size;j++){
int cc=ss1.at<int>(i,j);
ss.at<float>(i,j)=(float)(cc)/sum;
}
}
}1.4 尋找零交叉的細節
在濾波後的影象g(x,y)的任意畫素p處,尋找零交叉的一種方法是:使用以p為中心的一個3*3鄰域,p點處的零交叉意味著至少有兩個相對的鄰域畫素符號不同,有四種測試的情況:左/右,上/下和兩個對角,設定一個閾值,不僅相對鄰域的符號不同,而且他們的差值的絕對值還必須超過這個閾值.
void linyu(cv::Mat &src,cv::Mat &result,double threshold){
for (int y = 1; y < src.rows - 1; ++y)
{
for (int x = 1; x < src.cols; ++x)
{
// 鄰域判定
if ((src.at<float>(y - 1, x) *
src.at<float>(y + 1, x) < 0)&&(std::abs(src.at<float>(y - 1, x)-
src.at<float>(y + 1, x))>threshold))
{
result.at<uchar>(y, x) = 255;
}
if ((src.at<float>(y, x - 1) *
src.at<float>(y, x + 1) < 0)&&(std::abs(src.at<float>(y , x-1)-
src.at<float>(y , x+1))>threshold))
{
result.at<uchar>(y, x) = 255;
}
if ((src.at<float>(y + 1, x - 1) *
src.at<float>(y - 1, x + 1) < 0)&&(std::abs(src.at<float>(y + 1, x-1)-
src.at<float>(y - 1, x+1))>threshold))
{
result.at<uchar>(y, x) = 255;
}
if ((src.at<float>(y - 1, x - 1) *
src.at<float>(y + 1, x + 1) < 0)&&(std::abs(src.at<float>(y - 1, x-1)-
src.at<float>(y + 1, x+1))>threshold))
{
result.at<uchar>(y, x) = 255;
}
}
}
}1.5 opencv實現步驟
1 將圖片進行灰度化
2 將圖片轉換到畫素範圍[0,1]
3 生成高斯核矩陣並對影象進行平滑
4 對平滑後的影象進行拉普拉斯變換
5 進行零交叉查詢
6 得到結果
1.6 完整程式碼
#include <iostream>
#include<opencv2/opencv.hpp>
using namespace cv;
//將影象的動態範圍轉到[0,1]
void ImageAdjust(cv::Mat &src,cv::Mat &dst){
for(int i=0;i<src.rows;i++){
for(int j=0;j<src.cols;j++){
int g=src.at<uchar>(i,j);
dst.at<float>(i,j)=g/255.;
}
}
}
//構造高斯核矩陣
void GaussianKernel(cv::Mat &ss,int size,double delta){
int centerX=size/2;
int centerY=size/2;
int x=0;
int y=0;
cv::Mat gaussiankernel;
gaussiankernel=cv::Mat_<double>(size,size);
double cc=0.5/(3.1415926*delta*delta);
for(int i=0;i<size;i++){
for(int j=0;j<size;j++){
x=std::abs(i-centerX);
y=std::abs(j-centerY);
double dd=cc*exp(-((pow(x,2)+pow(y,2))/(2*delta*delta)));
gaussiankernel.at<double>(i,j)=dd;
}
}
cv::Mat ss1;
ss1=cv::Mat_<int>(size,size);
int sum=0;
//將高斯變換的小數型別轉為整型
double ff=gaussiankernel.at<double>(0,0);
for(int i=0;i<size;i++){
for(int j=0;j<size;j++){
int cc=cvRound(gaussiankernel.at<double>(i,j)/ff);
ss1.at<int>(i,j)=cc;
sum+=cc;
}
}
for(int i=0;i<size;i++){
for(int j=0;j<size;j++){
int cc=ss1.at<int>(i,j);
ss.at<float>(i,j)=(float)(cc)/sum;
}
}
}
//拉普拉斯變換
void Laplacian(cv::Mat &src, cv::Mat &dst)
{
cv::Mat element;
element=(cv::Mat_<int>(3,3)<<-1,-1,-1,-1,8,-1,-1,-1,-1);
for (int y = 1; y < src.rows-1; ++y) {
for (int x =1; x < src.cols-1; ++x) {
double value=0;
int i=0;
for (int m =y-1; m <= y+1; ++m) {
for (int n = x-1; n <= x+1; ++n) {
int m1=i/3;
int n1=i%3;
int t=element.at<int>(m1,n1);
i++;
float b=src.at<float>(m,n);
value +=t*b;
}
}
dst.at<float>(y,x)=value;
}
}
}
void linyu(cv::Mat &src,cv::Mat &result,double threshold){
for (int y = 1; y < src.rows - 1; ++y)
{
for (int x = 1; x < src.cols; ++x)
{
// 鄰域判定
if ((src.at<float>(y - 1, x) *
src.at<float>(y + 1, x) < 0)&&(std::abs(src.at<float>(y - 1, x)-
src.at<float>(y + 1, x))>threshold))
{
result.at<uchar>(y, x) = 255;
}
if ((src.at<float>(y, x - 1) *
src.at<float>(y, x + 1) < 0)&&(std::abs(src.at<float>(y , x-1)-
src.at<float>(y , x+1))>threshold))
{
result.at<uchar>(y, x) = 255;
}
if ((src.at<float>(y + 1, x - 1) *
src.at<float>(y - 1, x + 1) < 0)&&(std::abs(src.at<float>(y + 1, x-1)-
src.at<float>(y - 1, x+1))>threshold))
{
result.at<uchar>(y, x) = 255;
}
if ((src.at<float>(y - 1, x - 1) *
src.at<float>(y + 1, x + 1) < 0)&&(std::abs(src.at<float>(y - 1, x-1)-
src.at<float>(y + 1, x+1))>threshold))
{
result.at<uchar>(y, x) = 255;
}
}
}
}
int main() {
cv::Mat src=cv::imread("../3.png",1);
//影象轉為灰度圖
cv::Mat gray(src.size(),CV_32FC1);
cv::cvtColor(src,gray,cv::COLOR_BGR2GRAY);
cv::imshow("gray",gray);
cv::imwrite("../gray.png",gray);
//將影象的灰度範圍轉到[0.1]
cv::Mat gray01;
gray01=cv::Mat::zeros(gray.size(),CV_32FC1);
ImageAdjust(gray,gray01);
cv::imshow("gray01",gray01);
//生成高斯核矩陣並對影象進行平滑
int size=13;//高斯核矩陣的大小
double delta=2;
cv::Mat gaussiankernel;
cv::Mat gaussianblur;
gaussiankernel=cv::Mat_<float>(size,size);//(cv::Size(size,size),CV_64FC1);
GaussianKernel(gaussiankernel,size,delta);
cv::filter2D(gray01,gaussianblur,-1,gaussiankernel);
//對進行平滑之後的影象進行拉普拉斯變換
cv::Mat Laplace;
Laplace=cv::Mat::zeros(gaussianblur.size(),CV_32FC1);
Laplacian(gaussianblur,Laplace);
cv::imshow("Laplace",Laplace);
//找到最大值最小值,為設定閾值做參考
double max=0;
double min=0;
cv::minMaxLoc(Laplace,&min,&max);
//鄰域判定,零交叉查詢
cv::Mat result;
result=cv::Mat::zeros(Laplace.size(),CV_8U);
linyu(Laplace,result,0.02);
cv::imshow("result",result);
cv::imwrite("../result.png",result);
cv::waitKey(0);
return 0;
}1.7 實現效果
原圖
灰度影象
拉普拉斯變換影象
result