圖像處理中各種邊緣檢測的微分算子簡單比較(Sobel，Robert， Prewitt，Laplacian，Canny)

http://blog.csai.cn/user1/44205/archives/2008/23882.html

不同圖像灰度不同，邊界處一般會有明顯的邊緣，利用此特征可以分割圖像。需要說明的是：邊緣和物體間的邊界並不等同，邊緣指的是圖像中像素的值有突變的地方，而物體間的邊界指的是現實場景中的存在於物體之間的邊界。有可能有邊緣的地方並非邊界，也有可能邊界的地方並無邊緣，因為現實世界中的物體是三維的，而圖像只具有二維信息，從三維到二維的投影成像不可避免的會丟失一部分信息；另外，成像過程中的光照和噪聲也是不可避免的重要因素。正是因為這些原因，基於邊緣的圖像分割仍然是當前圖像研究中的世界級難題，目前研究者正在試圖在邊緣提取中加入高層的語義信息。

在實際的圖像分割中，往往只用到一階和二階導數，雖然，原理上，可以用更高階的導數，但是，因為噪聲的影響，在純粹二階的導數操作中就會出現對噪聲的敏感現象，三階以上的導數信息往往失去了應用價值。二階導數還可以說明灰度突變的類型。在有些情況下，如灰度變化均勻的圖像，只利用一階導數可能找不到邊界，此時二階導數就能提供很有用的信息。二階導數對噪聲也比較敏感，解決的方法是先對圖像進行平滑濾波，消除部分噪聲，再進行邊緣檢測。不過，利用二階導數信息的算法是基於過零檢測的，因此得到的邊緣點數比較少，有利於後繼的處理和識別工作。

各種算子的存在就是對這種導數分割原理進行的實例化計算，是為了在計算過程中直接使用的一種計算單位；

 

Roberts算子：邊緣定位准，但是對噪聲敏感。適用於邊緣明顯且噪聲較少的圖像分割。Roberts邊緣檢測算子是一種利用局部差分算子尋找邊緣的算子,Robert算子圖像處理後結果邊緣不是很平滑。經分析，由於Robert算子通常會在圖像邊緣附近的區域內產生較寬的響應，故采用上述算子檢測的邊緣圖像常需做細化處理，邊緣定位的精度不是很高。

Prewitt算子：對噪聲有抑制作用，抑制噪聲的原理是通過像素平均，但是像素平均相當於對圖像的低通濾波，所以Prewitt算子對邊緣的定位不如Roberts算子。

Sobel算子：Sobel算子和Prewitt算子都是加權平均，但是Sobel算子認為，鄰域的像素對當前像素產生的影響不是等價的，所以距離不同的像素具有不同的權值，對算子結果產生的影響也不同。一般來說，距離越遠，產生的影響越小。

Isotropic Sobel算子：加權平均算子，權值反比於鄰點與中心點的距離，當沿不同方向檢測邊緣時梯度幅度一致，就是通常所說的各向同性。

在邊沿檢測中，常用的一種模板是Sobel 算子。Sobel 算子有兩個，一個是檢測水平邊沿的；另一個是檢測垂直平邊沿的 。Sobel算子另一種形式是各向同性Sobel(Isotropic Sobel)算子，也有兩個，一個是檢測水平邊沿的 ，另一個是檢測垂直平邊沿的 。各向同性Sobel算子和普通Sobel算子相比，它的位置加權系數更為准確，在檢測不同方向的邊沿時梯度的幅度一致。由於建築物圖像的特殊性，我們可以發現，處理該類型圖像輪廓時，並不需要對梯度方向進行運算，所以程序並沒有給出各向同性Sobel算子的處理方法。
    由於Sobel算子是濾波算子的形式，用於提取邊緣，可以利用快速卷積函數，簡單有效，因此應用廣泛。美中不足的是，Sobel算子並沒有將圖像的主體與背景嚴格地區分開來，換言之就是Sobel算子沒有基於圖像灰度進行處理，由於Sobel算子沒有嚴格地模擬人的視覺生理特征，所以提取的圖像輪廓有時並不能令人滿意。  在觀測一幅圖像的時候，我們往往首先注意的是圖像與背景不同的部分，正是這個部分將主體突出顯示，基於該理論，我們可以給出閾值化輪廓提取算法，該算法已在數學上證明當像素點滿足正態分布時所求解是最優的。

上面的算子是利用一階導數的信息，屬於梯度算子范疇。

Laplacian算子：這是二階微分算子。其具有各向同性，即與坐標軸方向無關，坐標軸旋轉後梯度結果不變。但是，其對噪聲比較敏感，所以，圖像一般先經過平滑處理，因為平滑處理也是用模板進行的，所以，通常的分割算法都是把Laplacian算子和平滑算子結合起來生成一個新的模板。

Laplacian算子一般不以其原始形式用於邊緣檢測，因為其作為一個二階導數，Laplacian算子對噪聲具有無法接受的敏感性；同時其幅值產生算邊緣，這是復雜的分割不希望有的結果；最後Laplacian算子不能檢測邊緣的方向；所以Laplacian在分割中所起的作用包括：（1）利用它的零交叉性質進行邊緣定位；（2）確定一個像素是在一條邊緣暗的一面還是亮的一面；一般使用的是高斯型拉普拉斯算子（Laplacian of a Gaussian,LoG)，由於二階導數是線性運算，利用LoG卷積一幅圖像與首先使用高斯型平滑函數卷積改圖像，然後計算所得結果的拉普拉斯是一樣的。所以在LoG公式中使用高斯函數的目的就是對圖像進行平滑處理，使用Laplacian算子的目的是提供一幅用零交叉確定邊緣位置的圖像；圖像的平滑處理減少了噪聲的影響並且它的主要作用還是抵消由Laplacian算子的二階導數引起的逐漸增加的噪聲影響。

 

 

 

 

 

 

 

 

 

微分算子在圖像處理中扮演重要的角色，其算法實現簡單，而且邊緣檢測的效果又較好，因此這些基本的微分算子是學習圖像處理過程中的必備方法，下面著重討論幾種常見的微分算子。

1.Sobel

 其主要用於邊緣檢測，在技術上它是以離散型的差分算子，用來運算圖像亮度函數的梯度的近似值，缺點是Sobel算子並沒有將圖像的主題與背景嚴格地區分開來，換言之就是Sobel算子並沒有基於圖像灰度進行處理，由於Sobel算子並沒有嚴格地模擬人的視覺生理特征，所以提取的圖像輪廓有時並不能令人滿意，算法具體實現很簡單，就是3*3的兩個不同方向上的模板運算，這裡不再寫出。

2.Robert算子

根據任一相互垂直方向上的差分都用來估計梯度，Robert算子采用對角方向相鄰像素只差

3.Prewitt算子

   該算子與Sobel算子類似，只是權值有所變化，但兩者實現起來功能還是有差距的，據經驗得知Sobel要比Prewitt更能准確檢測圖像邊緣。

4.Laplacian算子

   拉普拉斯算子是一種二階微分算子，若只考慮邊緣點的位置而不考慮周圍的灰度差時可用該算子進行檢測。對於階躍狀邊緣，其二階導數在邊緣點出現零交叉，並且邊緣點兩旁的像素的二階導數異號。

5.Canny算子

該算子功能比前面幾種都要好，但是它實現起來較為麻煩，Canny算子是一個具有濾波，增強，檢測的多階段的優化算子，在進行處理前，Canny算子先利用高斯平滑濾波器來平滑圖像以除去噪聲，Canny分割算法采用一階偏導的有限差分來計算梯度幅值和方向，在處理過程中，Canny算子還將經過一個非極大值抑制的過程，最後Canny算子還采用兩個閾值來連接邊緣。

下面算法是基於的算法不可能直接運行，只是我把Canny的具體實現步驟寫了出來，若需用還要自己寫。

該算子具體實現方法：

// anny.cpp: implementation of the Canny class.
//
//////

#include "anny.h"
#include "math.h"
//#include "algorithms.h"
//#include "algorithm.h"
#include "stdlib.h"
//#include "maths.h"
//using namespace std;
//////
// Construction/Destruction
//////

Canny::Canny(int PicHeight,int PicWidth,double ** PicData,double PicSigma,double  PicRatioLow,double PicRatioHigh)
{
 iHeight=PicHeight;
 iWidth=PicWidth;
 iData=PicData;
 sigma=PicSigma;
 dRatioLow=PicRatioLow;
 dRatioHigh=PicRatioHigh;
}

Canny::~Canny()
{

}

void Canny::CannyArith(int **iEdgePoint)
{
 int i;
 int **iGradX ;                       // 指向x方向導數的指針
    int **iGradY ;                         // 指向y方向導數的指針
    int **iExtent ;                      // 梯度的幅度
 iGradX=new int *[iHeight];
 for(i=0;i<iHeight;i++) 
  iGradX[i]=new int[iWidth];
 iGradY=new int *[iHeight];
 for(i=0;i<iHeight;i++) 
  iGradY[i]=new int[iWidth];
 iExtent=new int *[iHeight];
 for(i=0;i<iHeight;i++) 
  iExtent[i]=new int[iWidth];
 
 // 對原圖象進行濾波
         GaussionSmooth();
// 計算X,Y方向上的方向導數
    DirGrad(iGradX,iGradY); 
    // 計算梯度的幅度
    GradExtent(iGradX,iGradY,iExtent);
    // 應用non-maximum 抑制
   NonMaxSuppress(iExtent,iGradX,iGradY,iEdgePoint);
 // 應用Hysteresis，找到所有的邊界
    Hysteresis(iExtent,iEdgePoint);
 // 釋放內存
 for(i=0;i<iHeight;i++)  
        delete  []*(iGradX+i);  
    delete   iGradX;
 for(i=0;i<iHeight;i++)  
        delete  []*(iGradY+i);  
    delete   iGradY;
 for(i=0;i<iHeight;i++)  
        delete  []*(iExtent+i);  
    delete   iExtent; 
 
 
}

 
void Canny::GaussionSmooth()
{
 int i,j,k;                             //循環變量
 int iWindowSize;                       //記錄模板大小的變量
 int iHalfLen;                          //模板大小的一半
 double *pdKernel;                         //模板各點的權值
 double dDotMul;                        //模板與對應像素點的卷積和
 double dWeightSum;                     //模板的權值累加和
 double **dTemp;                         //記錄圖像數據的中間變量
 //開辟空間
 dTemp=new double *[iHeight];           
 for(i=0;i<iHeight;i++)
  dTemp[i]=new double[iWidth];
 //獲得模板長度和模板的各個權值
 MakeGauss(&pdKernel,&iWindowSize);
 //得到模板的一半長度
 iHalfLen=iWindowSize/2;
 //對圖像對水方向根據模板進行平滑
 for(i=0;i<iHeight;i++)
 {
  for(j=0;j<iWidth;j++)
  {
   dDotMul=0;
   dWeightSum=0;
   for(k=(-iHalfLen);k<=iHalfLen;k++)
   {
    if((k+j>=0)&&(k+j<iWidth))
    {
     dDotMul+=iData[i][j+k]*pdKernel[k+iHalfLen];
     dWeightSum+=pdKernel[k+iHalfLen];

    }
   }
   dTemp[i][j]=dDotMul/dWeightSum;
  }
 }
 //對圖像垂直方向上根據模板的轉置進行平滑(注意圖像數據是在水平平滑之後進行的)
 for(i=0;i<iWidth;i++)
 {
  for(j=0;j<iHeight;j++)
  {
   dDotMul=0;
   dWeightSum=0;
   for(k=(-iHalfLen);k<=iHalfLen;k++)
   {
    if((k+j>=0)&&(k+j<iHeight))
    {
     dDotMul+=dTemp[j+k][i]*pdKernel[k+iHalfLen];
     dWeightSum+=pdKernel[k+iHalfLen];
     
    }
   }
   iData[j][i]=dDotMul/dWeightSum;
  }
 }
 //空間釋放
    delete []pdKernel;
 pdKernel=NULL;
 for(i=0;i<iHeight;i++)  
        delete  []*(dTemp+i);  
    delete   dTemp;  
   

}

void Canny::MakeGauss(double **pdKernel,int *iWindowSize)
{
 int i;                             //循環變量
 int nCenter;                       //確定高斯模板的一半長度
 double dDistance;                  //一維高斯模板各點離中心點的距離
 double PI=3.1415926;               //圓周率              
 double dValue;                     //中間變量，記錄高斯模板各點的權值(未經歸一化)
 double dSum=0;                     //中間變量，記錄高斯模板各點權值的總和
 *iWindowSize=int(1+2*int(3*sigma+0.5));    //確定一維高斯模板長度，根據概率論的知識，選取[-3*sigma, 3*sigma]以內的數據。
 nCenter=(*iWindowSize)/2;          //得到一半長度
 *pdKernel=new double[*iWindowSize];//開辟記錄各點權值的空間

 //利用高斯分布函數(正太分布)確定各點的權值，主要是根據高斯分布離中心點的距離越遠，所取的值就越小，這與圖像有些
 //相似，離中心點越遠，對中心點的影響就越小。
 for(i=0;i<(*iWindowSize);i++)
 {
  dDistance=double(i-nCenter);
  //高斯分布函數求值
  dValue=exp((-1/2)*dDistance*dDistance/(sigma*sigma))/(sqrt(2*PI)*sigma);
  (*pdKernel)[i]=dValue;
  dSum+=dValue;
  
 }
 //歸一化(因為要不改變原圖像的灰度區域，就必須保證各權值之和為1
 for(i=0;i<(*iWindowSize);i++)
 {
  (*pdKernel)[i] /= dSum;
 }
}

void Canny::DirGrad(int **iGradX,int **iGradY)
{
 int i,j,temp1,temp2;
 //水平方向的方向導數(下面都是用min和max對邊界值做了相應的處理)
 for(i=0;i<iHeight;i++)
 {
  for(j=0;j<iWidth;j++)
  {

   if(iWidth-1<j+1)
    temp1=iWidth-1;
   else
    temp1=j+1;
   if(0<j-1)
    temp2=j-1;
   else
    temp2=0;
   
   iGradX[i][j]=int(iData[i][temp1]-iData[i][temp2]);
  }
 }
 //垂直方向的方向導數
 for(i=0;i<iWidth;i++)
 {
  for(j=0;j<iHeight;j++)
  {
   if(iHeight-1<j+1)
    temp1=iHeight-1;
   else
    temp1=j+1;
   if(0<j-1)
    temp2=j-1;
   else
    temp2=0;
   iGradY[j][i]=int(iData[temp1][i]-iData[temp2][i]);
  }
 }
}

 
void Canny::GradExtent(int **iGradX,int **iGradY,int **iExtent)
{
 int i,j;
 double iTemp1,iTemp2;
 for(i=0;i<iHeight;i++)
 {
  for(j=0;j<iWidth;j++)
  {
   iTemp1=iGradX[i][j]*iGradX[i][j];
   iTemp2=iGradY[i][j]*iGradY[i][j];
   iExtent[i][j]=int(sqrt(iTemp1+iTemp2)+0.5);
  }
 }
}

 
void Canny::NonMaxSuppress(int **iExtent,int **iGradX,int **iGradY,int **dUnchRst)
{
 int i,j;
 int gx,gy;                     //記錄像素點X,Y 方向的方向導數值
 int g1,g2,g3,g4;               //各個領域的梯度值
 double weight;                    //比重
 double dTemp1,dTemp2,dTemp;       //中間變量
 //處理邊緣值(邊緣點不可能是邊界點
  for(i=0;i<iHeight;i++)
 {
  dUnchRst[i][0]=0;
  dUnchRst[i][iWidth-1]=0;
 }
 for(j=0;j<iWidth;j++)
 {
  dUnchRst[0][j]=0;
  dUnchRst[iHeight-1][j]=0;
 }
 //標記有可能是邊界點的像素點
 for(i=1;i<iHeight-1;i++)
 {
  for(j=1;j<iWidth-1;j++)
  {
   //梯度值是0的像素點不可能是邊界點
   if(iExtent[i][j]==0)
    dUnchRst[i][j]=0;
   else
   {
    dTemp=iExtent[i][j];
    gx=iGradX[i][j];
    gy=iGradY[i][j];
    //下面都是判斷當前像素點的梯度值和其領域像素點的梯度值，如大於就有可能是邊界點，如小於就不可能是邊界點
    if(abs(gy)>abs(gx))
    {
                       weight=double(abs(gx)/abs(gy));
        g2=iExtent[i-1][j];
        g4=iExtent[i+1][j];
        if(gx*gy>0)
        {
         g1=iExtent[i-1][j-1];
         g3=iExtent[i+1][j+1];
        }
        else
        {
         g1=iExtent[i-1][j+1];
         g3=iExtent[i+1][j-1];
        }
    }
    else
    {
     weight=double(abs(gy)/abs(gx));
     g2=iExtent[i][j+1];
     g4=iExtent[i][j-1];
     if(gx*gy>0)
     {
      g1=iExtent[i+1][j+1];
      g3=iExtent[i-1][j-1];
     }
     else
     {
      g1=iExtent[i-1][j+1];
      g3=iExtent[i+1][j-1];
     }
    }
    dTemp1=weight*g1+(1-weight)*g2;
    dTemp2=weight*g3+(1-weight)*g4;
    //當大於的時候就有可能是邊界點
    if(dTemp>=dTemp1&&dTemp>=dTemp2)
    {
     dUnchRst[i][j] = 128 ;
    }
    else
    {
     
     dUnchRst[i][j]=0 ;
    }
   }

  } 
 }
}

 

void Canny::Hysteresis(int **iExtent,int **iEdgePoint)
{
 int i,j;
 int iThreHigh;
 int iThreLow;
 SetThreshold(iExtent,&iThreHigh,&iThreLow,iEdgePoint);
 for(i=0;i<iHeight;i++)
 {
  for(j=0;j<iWidth;j++)
  {
   if((iEdgePoint[i][j]==128)&&(iExtent[i][j]>=iThreHigh))
   {
    iEdgePoint[i][j]=255;
    TraceEdge(i,j,iThreLow,iEdgePoint,iExtent);

   }

  }
 }
 // 那些還沒有被設置為邊界點的象素已經不可能成為邊界點
 for(i=0;i<iHeight;i++)
 {
  for(j=0;j<iWidth;j++)
  {
   
   if(iEdgePoint[i][j]!=255)
   {
    // 設置為非邊界點
    iEdgePoint[i][j] = 0 ;
   }
  }
  }

}

void Canny::SetThreshold(int **iExtent,int *iThreHigh,int *iThreLow,int **iEdgePoint)
{

 int i,j,k;
 int GradHist[1024];                     //統計梯度直方圖的數據，梯度最大值不可能超過1024
 int iEdgeNum;                           //邊界點的數量
 int iGradMax=0;                         //邊界點的梯度最大值
 int iHighCount;                         //根據iRatioHigh小於高閾值像素的個數
 //初始化
 for(i=0;i<1024;i++)
  GradHist[i]=0;
 //梯度直方圖統計
 for(i=0;i<iHeight;i++)
 {
  for(j=0;j<iWidth;j++)
  {
   if(iEdgePoint[i][j]==128)
   {
    GradHist[iExtent[i][j]]++;
   }
  }
 }
 iEdgeNum=0;
 //找出最大梯度和統計邊界點的個數
 for(i=0;i<1024;i++)
 {
  if(GradHist[i]!=0)
   iGradMax=i;
  iEdgeNum+=GradHist[i];
 }
 //獲得小於高閾值的個數
 iHighCount=int(iEdgeNum*dRatioHigh+0.5);
 k=1;
 iEdgeNum=GradHist[1];
 //求出高閾值
 while((k<=(iGradMax-1))&&(iEdgeNum<iHighCount))
 {
  k++;
  iEdgeNum+=GradHist[k];
 }
 *iThreHigh=k;
 //根據高閾值和比例關系求得低閾值
 *iThreLow=int((*iThreHigh)*dRatioLow+0.5);

}

void Canny::TraceEdge(int y,int x,int iThreLow,int **iEdgePoint,int **iExtent)
{
 // 對8鄰域象素進行查詢
 int xNb[8] = {1, 1, 0,-1,-1,-1, 0, 1} ;
 int yNb[8] = {0, 1, 1, 1,0 ,-1,-1,-1} ;
 int yy ;
 int xx ;
 int k  ;
 for(k=0;k<8;k++)
 {
  yy=y+yNb[k] ;
  xx=x+xNb[k] ;
 // 如果該象素為可能的邊界點，又沒有處理過, 並且梯度大於閾值
  if(iEdgePoint[yy][xx]==128&&iExtent[yy][xx]>=iThreLow)
  {
 // 把該點設置成為邊界點
  iEdgePoint[yy][xx]=255 ;
 // 以該點為中心進行跟蹤
      //TraceEdge(yy,xx,iThreLow,iEdgePoint,iExtent);
 
  }
 }
}

 

 

 

http://blog.sina.com.cn/s/blog_5fbeac730100ew38.html