关于图像中的导数问题

在图像处理的过程中，对其求导数是十分常见的事情：

• 如求图像中某点的梯度（含一阶导数）
• 如求图像中某点的拉普拉斯响应（含二阶纯偏导数）
• 再如求图像中某点的 Hessian 矩阵（含二阶混合偏导数）

---

一些重要的声明

• 对于边缘点求导像素值缺省的问题，优先建议使用 复制填充 方式，其次是 零填充 方式
• 针对图像而言，它的 两个混合偏导数一定相等 （本人已推导证明）
• 在图像中， 求导数和求卷积可交换顺序 ，例如：
  • 用高斯偏导核去卷积一幅图像等价于图像先做高斯滤波再求梯度
  • 用高斯二阶偏导核（ LoG 算子）去卷积一幅图像等价于图像先做高斯滤波再求拉普拉斯响应
  • 一般而言，我们都是先高斯滤波，再进行求导数操作（而不是对高斯算子进行求导数操作），因为这样更加符合认知，不容易出错（实现难度也减小了不少）
• 一般来说，认为 x 方向就是图像 行 的方向， y 的方向就是图像 列 的方向（参考冈萨雷斯的教材）
• 一般来说，都是直接使用 相关性 代替 卷积 进行滤波，这样做是没问题的

---

求导数的具体表达式

中心差分的误差相对于前向、后向差分的误差更小，因此优先考虑使用中心差分（来求导数）

一、前向、后向差分（不涉及二阶导时可使用）

• 前向差分：

\[ \frac{\partial f}{\partial x} = f(x+1,y) - f(x,y)\;\;, \;\;\frac{\partial f}{\partial y} = f(x,y+1) - f(x,y) \]

• 后向差分：

\[ \frac{\partial f}{\partial x} = f(x,y) - f(x-1,y)\;\;, \;\;\frac{\partial f}{\partial y} = f(x,y) - f(x,y-1) \]

二、中心差分（建议使用，使用 Taylor 展开导出）

• 一阶中心差分：

\[ \frac{\partial f}{\partial x} = \frac{f(x+1,y) - f(x-1,y)}{2}\;\;, \;\;\frac{\partial f}{\partial y} = \frac{f(x,y+1) - f(x,y-1)}{2} \]

• 二阶纯偏导中心差分：

\[ \frac{\partial^{2} f}{\partial x^{2}} = f(x+1,y) + f(x-1,y) -2f(x,y)\;\;, \;\;\frac{\partial^{2} f}{\partial y^{2}} = f(x,y+1) + f(x,y-1) -2f(x,y) \]

• 二阶混合偏导中心差分（偏导必相等）：

\[ \frac{\partial^{2} f}{\partial x \partial y} = \frac{\partial^{2} f}{\partial y \partial x} = \frac{[f(x-1,y-1) + f(x+1,y+1)] - [f(x-1,y+1) + f(x+1,y-1)]}{4} \]

---

实际应用

• 梯度幅值算子：（若 f 是高斯函数，则为高斯偏导核）

使用高斯偏导核卷积图像求梯度幅值 = 图像先高斯滤波，再求梯度幅值

\[ \lvert \frac{\partial f}{\partial x} \rvert + \lvert \frac{\partial f}{\partial y} \rvert \]

• 拉普拉斯算子：（若 f 是高斯函数，则为高斯二阶偏导核）

使用高斯二阶偏导核卷积图像 = 图像先高斯滤波，再求拉普拉斯响应

\[ \frac{\partial^{2} f}{\partial x^{2}} + \frac{\partial^{2} f}{\partial y^{2}} \]

• Hessian 矩阵：

Hessian 矩阵是针对于图像中的某个具体像素而言的

它分为使用邻域信息和不使用邻域信息两种

使用邻域信息又分为盒式平均和高斯平均两种

\[ \begin{pmatrix} \frac{\partial^{2} f}{\partial x^{2}} & \frac{\partial^{2} f}{\partial x \partial y} \\ \frac{\partial^{2} f}{\partial y \partial x} & \frac{\partial^{2} f}{\partial x^{2}})\end{pmatrix} \]

图像 Hessian 矩阵的三种求法

图像中的每一个像素点都能求出它的 Hessian 矩阵

要注意一点： Harris 角点检测算法中的矩阵使用的是一阶导数，不是 Hessian 矩阵

• 不使用邻域信息的 Hessian 矩阵

close all; clear; clc;
img = imread('lena.tif');
if size(img,3) == 3
    img = rgb2gray(img);
end
img = double(img);

%% 不使用邻域信息的Hessian矩阵
fxx = imfilter(img,[1;-2;1],'replicate');
fyy = imfilter(img,[1,-2,1],'replicate');
fxy = imfilter(img,[1;0;-1]*[1,0,-1]/4,'replicate');
fyx = fxy;
H = cell(size(img));
for i = 1:numel(img)
   H{i} = [fxx(i),fxy(i);fyx(i),fyy(i)];
end

• 使用 3×3 邻域信息的 Hessian 矩阵（使用盒式核滤波）

close all; clear; clc;
img = imread('lena.tif');
if size(img,3) == 3
    img = rgb2gray(img);
end
img = double(img);

%% 使用3×3邻域信息的Hessian矩阵（使用盒式核滤波）
h = fspecial('average',3);
blur_img = imfilter(img,h,'replicate');
fxx = imfilter(blur_img,[1;-2;1],'replicate');
fyy = imfilter(blur_img,[1,-2,1],'replicate');
fxy = imfilter(blur_img,[1;0;-1]*[1,0,-1]/4,'replicate');
fyx = fxy;
H = cell(size(img));
for i = 1:numel(img)
   H{i} = [fxx(i),fxy(i);fyx(i),fyy(i)];
end

• 使用 3×3 邻域信息的 Hessian 矩阵（使用高斯核滤波）

close all; clear; clc;
img = imread('lena.tif');
if size(img,3) == 3
    img = rgb2gray(img);
end
img = double(img);

%% 使用3×3邻域信息的Hessian矩阵（使用高斯核滤波）
h = fspecial('gaussian',3,1);
blur_img = imfilter(img,h,'replicate');
fxx = imfilter(blur_img,[1;-2;1],'replicate');
fyy = imfilter(blur_img,[1,-2,1],'replicate');
fxy = imfilter(blur_img,[1;0;-1]*[1,0,-1]/4,'replicate');
fyx = fxy;
H = cell(size(img));
for i = 1:numel(img)
   H{i} = [fxx(i),fxy(i);fyx(i),fyy(i)];
end