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

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

• 数学原理

离散函数导数

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

Laplace算子的差分形式

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

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

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

写成filter mask的形式如下，

注：

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

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

图像的Laplace操作

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

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

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

其中的参数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型的图像。

仿真结果：

原图：

Laplace操作结果：

• 基于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算子的操作结果

以灰色为主色调的显示结果就是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);

仿真结果：

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