图像处理中的一阶偏导数和二阶偏导数

http://blog.csdn.net/xiaofengsheng/article/details/6023368

http://www.cnblogs.com/german-iris/p/4840647.html

 1. 一阶差分:

 

2. 二阶偏导数的推导和近似:

 

3. 上式以点(i+1,j)为中心,用i代换i+1可得以(i,j)为中心的二阶偏导数则有:

 

4. 同理:

 

5. 进而可推导:

 

6. 这样我们就可以很好的运用其他的一阶偏导的定义,如SIFT特征OpenCV实现版本中的一阶以及二阶偏导:

[cpp]  view plain  copy
  1. /* 
  2. Computes the partial derivatives in x, y, and scale of a pixel in the DoG 
  3. scale space pyramid. 
  4.  
  5. @param dog_pyr DoG scale space pyramid 
  6. @param octv pixel's octave in dog_pyr 
  7. @param intvl pixel's interval in octv 
  8. @param r pixel's image row 
  9. @param c pixel's image col 
  10.  
  11. @return Returns the vector of partial derivatives for pixel I 
  12.     { dI/dx, dI/dy, dI/ds }^T as a CvMat* 
  13. */  
  14. static CvMat* deriv_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c )  
  15. {  
  16.     CvMat* dI;  
  17.     double dx, dy, ds;  
  18.   
  19.     dx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) -  
  20.         pixval32f( dog_pyr[octv][intvl], r, c-1 ) ) / 2.0;  
  21.     dy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) -  
  22.         pixval32f( dog_pyr[octv][intvl], r-1, c ) ) / 2.0;  
  23.     ds = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) -  
  24.         pixval32f( dog_pyr[octv][intvl-1], r, c ) ) / 2.0;  
  25.   
  26.     dI = cvCreateMat( 3, 1, CV_64FC1 );  
  27.     cvmSet( dI, 0, 0, dx );  
  28.     cvmSet( dI, 1, 0, dy );  
  29.     cvmSet( dI, 2, 0, ds );  
  30.   
  31.     return dI;  
  32. }  
  33.   
  34.   
  35.   
  36. /* 
  37. Computes the 3D Hessian matrix for a pixel in the DoG scale space pyramid. 
  38.  
  39. @param dog_pyr DoG scale space pyramid 
  40. @param octv pixel's octave in dog_pyr 
  41. @param intvl pixel's interval in octv 
  42. @param r pixel's image row 
  43. @param c pixel's image col 
  44.  
  45. @return Returns the Hessian matrix (below) for pixel I as a CvMat* 
  46.  
  47.     / Ixx  Ixy  Ixs / <BR> 
  48.     | Ixy  Iyy  Iys | <BR> 
  49.     / Ixs  Iys  Iss / 
  50. */  
  51. static CvMat* hessian_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c )  
  52. {  
  53.     CvMat* H;  
  54.     double v, dxx, dyy, dss, dxy, dxs, dys;  
  55.   
  56.     v = pixval32f( dog_pyr[octv][intvl], r, c );  
  57.     dxx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) +   
  58.             pixval32f( dog_pyr[octv][intvl], r, c-1 ) - 2 * v );  
  59.     dyy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) +  
  60.             pixval32f( dog_pyr[octv][intvl], r-1, c ) - 2 * v );  
  61.     dss = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) +  
  62.             pixval32f( dog_pyr[octv][intvl-1], r, c ) - 2 * v );  
  63.     dxy = ( pixval32f( dog_pyr[octv][intvl], r+1, c+1 ) -  
  64.             pixval32f( dog_pyr[octv][intvl], r+1, c-1 ) -  
  65.             pixval32f( dog_pyr[octv][intvl], r-1, c+1 ) +  
  66.             pixval32f( dog_pyr[octv][intvl], r-1, c-1 ) ) / 4.0;  
  67.     dxs = ( pixval32f( dog_pyr[octv][intvl+1], r, c+1 ) -  
  68.             pixval32f( dog_pyr[octv][intvl+1], r, c-1 ) -  
  69.             pixval32f( dog_pyr[octv][intvl-1], r, c+1 ) +  
  70.             pixval32f( dog_pyr[octv][intvl-1], r, c-1 ) ) / 4.0;  
  71.     dys = ( pixval32f( dog_pyr[octv][intvl+1], r+1, c ) -  
  72.             pixval32f( dog_pyr[octv][intvl+1], r-1, c ) -  
  73.             pixval32f( dog_pyr[octv][intvl-1], r+1, c ) +  
  74.             pixval32f( dog_pyr[octv][intvl-1], r-1, c ) ) / 4.0;  
  75.   
  76.     H = cvCreateMat( 3, 3, CV_64FC1 );  
  77.     cvmSet( H, 0, 0, dxx );  
  78.     cvmSet( H, 0, 1, dxy );  
  79.     cvmSet( H, 0, 2, dxs );  
  80.     cvmSet( H, 1, 0, dxy );  
  81.     cvmSet( H, 1, 1, dyy );  
  82.     cvmSet( H, 1, 2, dys );  
  83.     cvmSet( H, 2, 0, dxs );  
  84.     cvmSet( H, 2, 1, dys );  
  85.     cvmSet( H, 2, 2, dss );  
  86.   
  87.     return H;  
  88. }  

 

参考:

(1)http://hi.baidu.com/shareshow/blog/item/34abdf544725cf54d109069b.html

(2)SIFT的OpenCV实现




OpenCV-跟我一起学数字图像处理之拉普拉斯算子

Laplace算子和Sobel算子一样,属于空间锐化滤波操作。起本质与前面的Spatial Filter操作大同小异,下面就通过Laplace算子来介绍一下空间锐化滤波,并对OpenCV中提供的Laplacian函数进行一些说明。

  • 数学原理

离散函数导数

离散函数的导数退化成了差分,一维一阶差分公式和二阶差分公式分别为,

CodeCogsEqn

CodeCogsEqn(2)

Laplace算子的差分形式

分别对Laplace算子x,y两个方向的二阶导数进行差分就得到了离散函数的Laplace算子。

在一个二维函数f(x,y)中,x,y两个方向的二阶差分分别为,

CodeCogsEqn(3)

CodeCogsEqn(4)

所以Laplace算子的差分形式为,

CodeCogsEqn(5)

写成filter mask的形式如下,

0 1 0
1 -4 1
0 1 0
注意该mask的特点,mask在上下左右四个90度的方向上结果相同,也就是说在90度方向上无方向性。为了让该mask在45度的方向上也具有该性质,对该filter mask进行扩展定义为,
1 1 1
1 -8 1
1 1 1
 

注:

有时我们也会见到不同于上述结果的Laplace算子的filter mask,

0 -1 0
-1 4 -1
0 -1 0
 
-1 -1 -1
-1 8 -1
-1 -1 -1

其原因是在定义二阶导数的时候采用了相反的定义,这个无关紧要,但是要注意,当用Laplace算子滤波后的图像与原图叠加时,混合操作是加还是减因上述的定义而异。

图像的Laplace操作

如同本文开始时说的那样,将Laplace算子写成filter mask后,其操作大同小异于其他的空间滤波操作。将filter mask在原图上逐行移动,然后mask中数值与其重合的像素相乘后求和,赋给与mask中心重合的像素,对图像的第一,和最后的行和列无法做上述操作的像素赋值零,就得到了拉普拉斯操作结果。

拉普拉斯操作结果与原图的混合

因为Laplace算子是二阶导数操作,其在强调图像素中灰度不连续的部分的同时也不在强调灰度值连续的部分。这样会产生一个具有很明显的灰度边界,但是没有足够特征的黑色背景。背景特征可以通过原图像与Laplace算子操作后的图像混合恢复。用公式,

CodeCogsEqn(6)

其中的参数c的取值和上面的两种mask定义有关,当mask中心的数值取正时c=-1,相反c=1;

  • 基于OpenCV的Laplace算子的计算

OpenCV中Laplacian函数可以实现对图像的Laplace操作,具体用法如下,

Laplacian( src_gray, dst, ddepth, kernel_size, scale, delta, BORDER_DEFAULT );

参数意义为,

  1. src_gray,输入图像
  2. dst,Laplace操作结果
  3. ddepth,输出图像深度,因为输入图像一般为CV_8U,为了避免数据溢出,输出图像深度应该设置为CV_16S
  4. kernel_size,filter mask的规模,我们的mask时3x3的,所以这里应该设置为3
  5. scale,delta,BORDER_DEFAULT,默认设置就好

基于OpenCV的Laplace算子仿真代码段如下,

复制代码
//load the Original Image and get some informations
Mat src = imread("012.jpg",0);
namedWindow("OriginalImage");
imshow("OriginalImage",src);
CV_Assert(src.depth() == CV_8U);

//OpenCV solution - Laplacian
Mat dst,abs_dst_laplace;
Laplacian(src,dst,CV_16S,3);
convertScaleAbs(dst,abs_dst_laplace);

//show the result
namedWindow("result_laplacian");
imshow("result_laplacian",abs_dst_laplace);
复制代码

其中convertScaleAbs函数功能是将CV_16S型的输出图像转变成CV_8U型的图像。

仿真结果:

原图:

original

Laplace操作结果:

abs_dst_laplae

  • 基于mask operation原理仿真

Laplace算子滤波仿真

根据数学原理中介绍的算法,编写相应代码,进行相关仿真。其中对Laplace操作结果进行了图像拉伸显示,因为Laplace操作结果的像素值范围可能落在了[0,255]之外,而计算机在显示的时候将赋值全部置为0,大于255的像素全部显示成255。

代码段如下,

复制代码
//get some informations of original image
int nr = src.rows;
int nc = src.cols;
int n = nr*nc;
int arr[9] = {0};

//scan the whole pixels of original image 
//and do Laplacian Operation
int* table_lap = new int[n];
int* table_orig = new int[n];
int l;
for (int i=0;i<n;i++)
{
    table_lap[i] = 0;
    table_orig[i] = 0;
}
for (int i=1;i<nr-1;i++)
{
    const uchar* previous = src.ptr<uchar>(i-1);
    const uchar* current = src.ptr<uchar>(i);
    const uchar* next = src.ptr<uchar>(i+1);
    for (int j=1;j<nc-1;j++)
    {
        for (int k=0;k<3;k++)
        {
            arr[k] = previous[j+k-1];
            arr[k+3] = current[j+k-1];
            arr[k+6] = next[j+k-1];
        }
        l = nc*i+j;        //calculate the location in the table of current pixel
        Lmaskoperation(table_lap,arr,l);
        table_orig[l] = arr[4];
    }
}

//pixels scale
uchar* La_scaled = new uchar[n];
table_scale(table_lap,La_scaled,n);

//padding values
Mat LaResult_own;
LaResult_own.create(src.size(),src.type());
uchar* p = NULL;
for (int i=0;i<nr;i++)
{
    p = LaResult_own.ptr<uchar>(i);
    for (int j=0;j<nc;j++)
    {
        l = nc*i+j;
        p[j] = La_scaled[l];
    }
}

//show results
namedWindow("LaResult_own");
imshow("LaResult_own",LaResult_own);
复制代码

其中Lmaskoperation是我写的mask为Laplace mask的mask operation操作函数,函数段如下,

复制代码
//**********************//
//Laplacian mask operation
//**********************//
void Lmaskoperation(int* table,int* arr,int l)
{
    int tmp[9] = {-1,-1,-1,-1,8,-1,-1,-1,-1};
    for (int i=0;i<9;i++)
    {
        table[l] = table[l] + tmp[i]*arr[i];
    }
}
复制代码

tabel_scale函数就是我写的图像拉伸函数,将Laplace操作结果拉伸到[0,255],具体函数段如下,

复制代码
//*****************************//
//scale the pixels to [0 255]
//*****************************//
void table_scale(int* table,uchar* result,int n)
{
    int min = table[0];
    int max = table[0];
    for (int i=0;i<n;i++)
    {
        if(min>table[i])
        {
            min = table[i];
        }
        if(max<table[i])
        {
            max = table[i];
        }
    }
    for (int i=0;i<n;i++)
    {
        result[i] = (uchar)(255*(table[i]-min)/(max-min));
    }
}
复制代码

仿真结果,拉伸后Laplace算子的操作结果

LaResult_own

以灰色为主色调的显示结果就是Laplace算子操作拉伸后显示的一大特点。

Laplace滤波图像与原图像的混合

我使用的mask中心值为正,所以混合操作需要原图减去Laplace滤波图像,代码段如下,

复制代码
//blending with the original image using Eq g(x,y)=f(x,y)+c*Lap(x,y)
int* table_blend = new int[n];
for(int i=0;i<n;i++)
{
    table_blend[i] = table_orig[i] - table_lap[i];
    if(table_blend[i]<0)
    {
        table_blend[i] = 0;
    }
    else if (table_blend[i]>255)
    {
        table_blend[i] = 255;
    }
}

//padding values to blending result
Mat Blresult;
Blresult.create(src.size(),src.type());
for (int i=0;i<nr;i++)
{
    p = Blresult.ptr<uchar>(i);
    for(int j=0;j<nc;j++)
    {
        l = nc*i+j;
        p[j] = table_blend[l];
    }
}

//show blending result
namedWindow("blending result_laplacian");
imshow("blending result_laplacian",Blresult);
复制代码

仿真结果:

blending result_laplacian

最后给出冈萨雷斯在介绍Laplacian时所给素材的仿真结果

blending result_laplacian



评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值