加载中…
个人资料
彩泓
彩泓
  • 博客等级:
  • 博客积分:0
  • 博客访问:20,618
  • 关注人气:19
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
相关博文
推荐博文
谁看过这篇博文
加载中…
正文 字体大小:

图像处理中各种微分算子

(2009-07-23 13:01:16)
标签:

it

分类: 图像处理

   微分算子在图像处理中扮演重要的角色,其算法实现简单,而且边缘检测的效果又较好,因此这些基本的微分算子是学习图像处理过程中的必备方法,下面着重讨论几种常见的微分算子。

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);
 
  }
 }
}

0

阅读 评论 收藏 转载 喜欢 打印举报/Report
前一篇:高斯滤波器
  • 评论加载中,请稍候...
发评论

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

    < 前一篇高斯滤波器
      

    新浪BLOG意见反馈留言板 电话:4000520066 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有