如何计算图像的曲率

最新推荐文章于 2023-07-29 19:33:44 发布

Ein027

最新推荐文章于 2023-07-29 19:33:44 发布

阅读量1.3w

点赞数 3

分类专栏： Image Processing 文章标签：图像处理 Image 数学

Image Processing 专栏收录该内容

6 篇文章 2 订阅

订阅专栏

作者：龚元浩
链接：https://zhuanlan.zhihu.com/p/23132541
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

一般的数字图像都是一个二维离散函数 I(x_i,y_j) ，图像处理的任务通常是想得到一个新的图像 U(x_i,y_j) 。这个满足特定的性质，比方说去模糊、去雾、超分辨率、去噪、分割等等。以去模糊为例，得到的图像要比原始图像清晰。这个过程我们通常要建立一个数学模型 ${\cal E}(U)$ ，使得该模型最小值对应的即为所求。显然，模型的建立对于最终的结果至关重要。另一方面，模型也要易于求解才行。所以说，模型和求解方法是一对孪生姐妹。

最为通用的求解方法是梯度下降法，即: $\frac{\partial U^t}{\partial t}=-\frac{\partial {\cal E}}{\partial U}$

但梯度下降法也是最为低效率的方法之一。原因是多方面的。首先，它没有考虑 ${\cal E}(U)$ 的可导性。对于不可导的模型可能需要做光滑逼近。其次，在时间步长的选取上它需要满足数值稳定性要求，即便是选取高阶精度的时间离散化公式（梯形公式、龙哥库塔方法等等）。再次，它没有考虑模型 ${\cal E}(U)$ 的具体形式，所以虽然非常通用，但对于具体问题又非常没有效率。

另外一种求解方法是从如下必要条件构造迭代公式： $\frac{\partial {\cal E}}{\partial U}=0$

通常，这种方法比梯度下降法更有效率，也需要了解模型 ${\cal E}(U)$ 的具体形式。但构造这种方法一般比较困难，得到的结果很难看出其对应的物理意义。

前文基于曲率的图像处理介绍的曲率滤波方法针对曲率正则项设计了一种全新的优化方法，能够高效地减少图像的曲率而不需要计算曲率。详细的介绍见曲率滤波博士论文第六章。针对图像的曲率优化问题，这一算法是目前最高效的优化方法（没有之一）。

那么，计算图像的曲率是不是就不必要了呢？其实也不是。最明显的例子还是曲率滤波。虽然我们能优化曲率而不需要计算曲率，我们还是需要观察模型 ${\cal E}(U)$ 以确保收敛。这就需要计算曲率正则项的能量。如何计算曲率呢，尤其是针对离散的数字图像？

让我们以平均曲率为例。我们考虑曲面 $\Psi(\vec{x})=(\vec{x},U(\vec{x}))$ 。我们很容易计算它的第一和第二基本形式，从而得到它的平均曲率公式：

$H=\frac{(1+U_x^2)U_{yy}-2U_xU_yU_{xy} +(1+U_y^2)U_{xx}}{2(1+U_x^2+U_y^2)^{\frac{3}{2}}}$

大部分跟平均曲率相关的论文都通过离散化这个公式来计算平均曲率。另外一种计算方法是通过二次曲面拟合来避免离散化，即

$U(x,y)\approx f(x,y)\equiv C_5x^2+C_4y^2+C_3xy+C_2x+C_1y+ C_0$

然后通过最小二乘法先确定系数 C_i ，之后将 C_i 代入上式，得到：

$H\approx \frac{(1+C_2^2)C_4-C_2C_1C_3 +(1+C_1^2)C_5}{(1+C_1^2+C_2^2)^{\frac{3}{2}}}$

这种方法的计算量比较大，而且没有用到图像regular sampling的性质，所以在图像处理中用得不多，而是在三维点云和三角网格等不规则采样的处理中比较常见。

读到这里，聪明的读者会问：有没有不需要离散化而且利用规则采样的计算公式呢？笔者在构造平均曲率滤波的过程中，不经意地发现了一种全新的计算公式（曲率滤波博士论文的公式6.12）。希望有经验的读者斧正：

$H\approx \left( \begin{array}{ccc} \frac{-1}{16} & \frac{5}{16} &\frac{-1}{16}\\ \frac{5}{16} & -1 &\frac{5}{16}\\ \frac{-1}{16} & \frac{5}{16} &\frac{-1}{16}\\ \end{array}\right)\ast U$

这个公式是一个线性卷积操作，跟上面两个非线性的计算公式完全不同。这个公式是如何构造出来的呢？

在解释这个公式之前，让我们回忆一下微分几何中大名鼎鼎的欧拉理论（Euler Theorem, 1760，维基百科，好吧其实我当初看的是陈省身的《微分几何初步》）： $\kappa_{\theta}=\kappa_1\cos^2\theta+\kappa_2\sin^2\theta$ ，

其中 $\kappa_{\theta}$ 是方向曲率， $\kappa_1$ 和 $\kappa_2$ 为主曲率， $\theta$ 是到主平面的夹角。我们可以对 $\theta$ 求积分，得到

$\begin{split}\int\limits_{-\pi}^{\pi}\kappa_{\theta}\mathrm{d}\theta=&\kappa_1\int\limits_{-\pi}^{\pi}\cos^2\theta\mathrm{d}\theta+\kappa_2\int\limits_{-\pi}^{\pi}\sin^2\theta\mathrm{d}\theta\\=&\kappa_1\pi+\kappa_2\pi=\pi(\kappa_1+\kappa_2)=2\pi H\end{split}\,.$

从而，我们有：

$H=\frac{1}{2\pi}\int\limits_{-\pi}^{\pi}\kappa_{\theta}\mathrm{d}\theta$

注意，这个公式是数学上严格成立的，所以不是约等于符号，是严格等于。这个公式说，平均曲率是所有方向曲率的平均。这也是它为什么被称为“平均”曲率的原因（不光指它是两个主曲率的平均）。笔者在曲率滤波博士论文中计算各个方向上的方向曲率 $\kappa_{\theta}$ ，然后相加，再除以 $2\pi$ ，最终得到上述线性卷积公式。显然，方向曲率是被离散逼近的，所以上述线性计算公式是约等于符号。这里我们利用了局部窗口中方向曲率是有限的、离散的这一特性。我们的计算公式是线性的，充分利用了方向的离散型，从而避免了使用连续公式然后离散化带来的误差（计算是直接从离散图像来的）。当然，曲率滤波博士论文中方向曲率的计算是一个线性近似公式，然后假设近似误差在多次平均的情况下相互抵消（如果你想到了大数定理，说明你的数学比我好）。感兴趣的读者可以看看曲率滤波博士论文第六章。

那么线性公式和非线性公式哪个的计算精度更高呢？笔者在曲率滤波博士论文中做了简单的测试，用两个公式分别计算一个球图像的曲率，然后跟真实值比较。得到的结论是线性公式更为准确。当然，这个实验也只是针对球的图像，并不是所有图像都会有相同的结论。不过我们在500幅自然图像（BSDS500）上的对比结果表明，两者估计的平均曲率差别不大。但是从计算量和公式的优美程度上来讲，笔者还是倾向于线性卷积公式。

上图左边为原始图像，中间为非线性经典公式计算的结果，右边为线性卷积公式计算的结果。

线性卷积公式跟经典的非线性公式相比，另外一个优势是它对完全没有任何要求，而非线性公式则要求图像至少二次可导地光滑。这也导致非线性公式在大梯度附近不能正确地估计平均曲率。显然，没有光滑性要求的线性公式更适合估计平均曲率，尤其是在大梯度区域。

一个是非线性的经典公式，一个是线性的卷积公式，两者都可以用来估计平均曲率。你会选择哪一个呢？

最后，再次对比两个公式：

$H=\frac{(1+U_x^2)U_{yy}-2U_xU_yU_{xy} +(1+U_y^2)U_{xx}}{2(1+U_x^2+U_y^2)^{\frac{3}{2}}}$ ，

$H\approx \left( \begin{array}{ccc} \frac{-1}{16} & \frac{5}{16} &\frac{-1}{16}\\ \frac{5}{16} & -1 &\frac{5}{16}\\ \frac{-1}{16} & \frac{5}{16} &\frac{-1}{16}\\ \end{array}\right)\ast U$ （曲率滤波博士论文中的公式6.12）

丙申年秋

参考文献：