图像曲率 与 黎曼几何

一、图像曲率

物体曲率图：

基于曲率的图像处理, （曲率滤波算法比传统的几何流算法快1000到10000倍） 论文查看

下图是高斯曲率滤波的算法伪代码（计算太简单了，而且优化效果比state of the art还要快三个数量级，也就是1000倍，推导过程见 论文 ）：

github代码：https://github.com/YuanhaoGong/CurvatureFilter

一般的数字图像都是一个二维离散函数 $I(x,y)$，图像处理的任务通常是想得到一个新的图像 $U(x,y)$。这个 $U(x,y)$满足特定的性质，比方说去模糊、去雾、超分辨率、去噪、分割等等。以去模糊为例，得到的图像要比原始图像清晰。

计算图像的曲率¹: 计算公式，滤波核

$$H=\frac{(1+U_x^2)U_{yy}-2U_x{U}_y{U}_{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$$

python 示例：

#encoding=utf-8
import numpy as np
import cv2
import os



def calcAndDrawHist(image, color):
    hist = cv2.calcHist([image], [0], None, [256], [0.0, 255.0])
    minVal, maxVal, minLoc, maxLoc = cv2.minMaxLoc(hist)
    histImg = np.zeros([256, 256, 3], np.uint8)
    hpt = int(0.9 * 256)

    for h in range(256):
        intensity = int(hist[h] * hpt / maxVal)
        cv2.line(histImg, (h, 256), (h, 256 - intensity), color)
    return histImg


def calculate_curvature(img):
    x, y = np.gradient(img)  # 一阶导数
    xx, xy = np.gradient(x)  # 二阶偏导数
    yx, yy = np.gradient(y)  # 二阶偏导数
    Iup = (1 + x * x) * yy - 2 * x * y * xy + (1 + y * y) * xx   # 公式的分子
    Idown = np.power((2 * (1 + x * x + y * y)), 1.5)             # 公式的分母
    final = Iup / Idown
    final = abs(final)
    final = (final - final.min()) / (final.max() - final.min())  # 将结果归一化
    final = (final * 255).astype(np.uint8)
    return final


if __name__ == '__main__':
    img_u8 = cv2.imread("len_std.jpg", 0)
    final = calculate_curvature(img_u8)
    histImg = calcAndDrawHist(final, [0, 0, 255])
    cv2.imwrite("len_std_final.png", final)
    cv2.imwrite("len_std_histImg.png", histImg)

高斯曲率滤波

from skimage import io
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import time



def update_GC(img, row, col):
    img_ij = img[row:-1:2, col:-1:2]
    img_prev = img[row - 1:-2:2, col:-1:2]
    img_next = img[row + 1::2, col:-1:2]
    img_left = img[row:-1:2, col - 1:-2:2]
    img_rigt = img[row:-1:2, col + 1::2]
    img_leUp = img[row - 1:-2:2, col - 1:-2:2]
    img_riUp = img[row - 1:-2:2, col + 1::2]
    img_leDn = img[row + 1::2, col - 1:-2:2]
    img_riDn = img[row + 1::2, col + 1::2]

    d1 = (img_prev + img_next) / 2.0 - img_ij
    d2 = (img_left + img_rigt) / 2.0 - img_ij
    d3 = (img_leUp + img_riDn) / 2.0 - img_ij
    d4 = (img_leDn + img_riUp) / 2.0 - img_ij
    d5 = (img_prev + img_leUp + img_left) / 3.0 - img_ij
    d6 = (img_prev + img_riUp + img_rigt) / 3.0 - img_ij
    d7 = (img_leDn + img_left + img_next) / 3.0 - img_ij
    d8 = (img_rigt + img_riDn + img_next) / 3.0 - img_ij

    d = d1 * (np.abs(d1) <= np.abs(d2)) + d2 * (np.abs(d2) < np.abs(d1))
    d = d * (np.abs(d) <= np.abs(d3)) + d3 * (np.abs(d3) < np.abs(d))
    d = d * (np.abs(d) <= np.abs(d4)) + d4 * (np.abs(d4) < np.abs(d))
    d = d * (np.abs(d) <= np.abs(d5)) + d5 * (np.abs(d5) < np.abs(d))
    d = d * (np.abs(d) <= np.abs(d6)) + d6 * (np.abs(d6) < np.abs(d))
    d = d * (np.abs(d) <= np.abs(d7)) + d7 * (np.abs(d7) < np.abs(d))
    d = d * (np.abs(d) <= np.abs(d8)) + d8 * (np.abs(d8) < np.abs(d))

    img_ij[...] += d


def CF(inputimg, filterType=2, total_iter=10):
    outputimg = np.copy(inputimg)
    localFunc = {
        2: update_GC,
    }
    update = localFunc.get(filterType)
    for iter_num in range(total_iter):
        update(outputimg, 1, 1)
        update(outputimg, 2, 1)
        update(outputimg, 1, 2)
        update(outputimg, 2, 2)
    return outputimg


def gauss_filter_img(img_gray):
    lena = img_gray.astype('float32', copy=True)
    filterType = 2
    lena_filtered = CF(lena, filterType)
    return lena - lena_filtered


if __name__ == '__main__':
	# read the image
	lena_original = io.imread('../images/lena.png')
	lena = lena_original.astype('float32', copy= True)
	
	# set the filter type: 0 (Total Variation), 1 (Mean Curvature), 2 (Gaussian Curvature)
	filterType = 2;
	
	# curvature filter
	lena_filtered = CF(lena, filterType)

	fig, axes = plt.subplots(nrows=1, ncols=3, figsize = (17,8))
	axes[0].imshow(lena, cmap=cm.Greys_r)
	axes[0].set_title('Original Image')
	axes[0].set_axis_off()
	axes[1].imshow(lena_filtered, cmap=cm.Greys_r)
	axes[1].set_title('Filtered Image (10 iterations)')
	axes[1].set_axis_off()
	
	axes[2].imshow(lena - lena_filtered, cmap=cm.Greys_r)
	axes[2].set_title('Difference')
	axes[2].set_axis_off()

二、图像梯度与黎曼距离

关于梯度函数或卷积核 。

协方差矩阵 $\mathbf{A}$ 和 $\mathbf{B}$ 之间的黎曼距离为 (1)，其中 $\mathbf{A}$ 和 $\mathbf{B}$ 必须为正定矩阵。
协方差矩阵 $\mathbf{A}$ 和 $\mathbf{B}$ 之间的Log det距离为 (2).
$$d_{riemann}={\left(\sum _i\log ({\lambda }_i{)}^2\right)}^{-1/2}其中:{\lambda }_i=\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{h}(\mathbf{A},\mathbf{B})\text{\hspace{1em}(1)}$$
$$d_{Log-det}=\sqrt{\left(\log (\det (\frac{\mathbf{A}+\mathbf{B}}{2}))-0.5\times \log (\det (\mathbf{A})\det (\mathbf{B}))\right)}\text{\hspace{1em}(2)}$$

黎曼几何与流形学习：https://zhuanlan.zhihu.com/p/81945115
SPD流形上常用操作：https://zhuanlan.zhihu.com/p/332451578

2.1 黎曼几何

黎曼几何 是研究内蕴几何的几何分支。通俗来讲，就是我们可能生活在弯曲的空间中，比如一只生活在二维球面的蚂蚁，作为生活在弯曲空间中的个体，我们并没有足够多的智慧去把我们的弯曲嵌入到更高维的空间中去研究，就好比蚂蚁只懂得在球面上爬，不能从“三维空间的曲面”这一观点来认识球面，因为球面就是它们的世界。因此，我们就有了内蕴几何，它告诉我们，即便是身处弯曲空间中，我们依旧能够测量长度、面积、体积等，我们依旧能够算微分、积分，甚至我们能够发现我们的空间是弯曲的！也就是说，身处球面的蚂蚁，只要有足够的智慧，它们就能发现曲面是弯曲的——跟哥伦布环球航行那样——它们朝着一个方向走，最终却回到了起点，这就可以断定它们自身所处的空间必然是弯曲的——这个发现不需要用到三维空间的知识。

黎曼度量–测地线 ²
微分几何是把曲面看成是三维空间中的二维子集，而黎曼几何则是从二维曲面本身内蕴地研究几何问题。

黎曼几何时，往“客观实体”方向思考，总是有益的。

有了度规，可以很自然地引入“测地线”这一实体。狭义来看，它就是两点间的最短线——是平直空间的直线段概念的推广（实际的测地线不一定是最短的，但我们先不纠结细节，而且这不妨碍我们理解它，因为测地线至少是局部最短的）。不难想到，只要两点确定了，那么不管使用什么坐标，两点间的最短线就已经确定了，因此这显然是一个客观实体。有一个简单的类比，就是不管怎么坐标变换，一个函数 $f(x)$ 的图像极值点总是确定的——不管你变还是不变，它就在那儿，不偏不倚。

向量与联络³

当我们在我们的位置建立起自己的坐标系后，我们就可以做很多测量，测量的结果可能是一个标量，比如温度、质量，这些量不管你用什么坐标系，它都是一样的。当然，有时候我们会测量向量，比如速度、加速度、力等，这些量都是客观实体，但因为测量结果是用坐标的分量表示的，所以如果换一个坐标，它的分量就完全不一样了。

不同位置的人可能出于各种原因，使用了不同的坐标系，因此，当我们写出一个向量 $A^{\mu }$时，严格来讲应该还要注明是在x位置测量的：$A^{\mu }(x)$，只有不引起歧义的情况下，我们才能省略它。

由于不同位置使用了不同的坐标系、坐标的单位长度不一样、坐标轴的指向不一样，等等。
因此，我们需要将一个位置的向量坐标，变换到另外一个位置的坐标中去，才能对分量进行比较。这里，我们只考虑两个距离为无穷小的位置 $x$ 和 $x+\mathrm{d}x$的变换，即怎么将位于 $x+\mathrm{d}x$处测量的向量 $A^{\mu }(x+\mathrm{d}x)$（$x+\mathrm{d}x$处有自己的坐标系），用$x$ 处的坐标系表示出来。很显然，这涉及到一个变换矩阵，因此关键是确定这个变换矩阵。

黎曼曲率⁴，

黎曼曲率提供了一种方案，让身处空间内部的人也能计算自身所处空间的弯曲程度。因此，能够身处空间之中而发现空间中的弯曲与否，是一件很了不起的事情，就好像我们已经超越了我们现有的空间，到了更高维的空间去“居高临下”那样。

如果站在更高维空间的角度看，就容易发现空间的弯曲。比如弯曲空间中有一条测地线，从更高维的空间看，它就是一条曲线，可以计算曲率等，但是在原来的空间看，它就是直的，测地线就是直线概念的一般化，因此不可能通过这种途径发现空间的弯曲性，必须有一些迂回的途径。可能一下子不容易想到，但是各种途径都殊途同归后，就感觉它是显然的了。

怎么更好地导出黎曼曲率来，使得它能够明显地反映出弯曲空间跟平直空间的本质区别呢？几种导出黎曼曲率的方式。文章中三种不同的黎曼曲率张量的导出方式，分别从三个角度表明了弯曲空间与平直空间的区别：

1. 平直空间中，协变导数的次序是可以交换的，弯曲空间则不是；
2. 平直空间中，测地线分布的均匀的、线性的，弯曲空间则不是；
3. 平直空间中，一个向量去“溜达”一圈回来之后，并没有变化，而弯曲空间中，向量去“溜达”完之后，可能就不是原来的向量了。

曲率的计数与计算(Python)⁵

1、曲率的独立分量

黎曼曲率张量是一个非常重要的张量，当且仅当它全部分量为0时，空间才是平直的。它也出现在爱因斯坦的场方程中。总而言之，只要涉及到黎曼几何，黎曼曲率张量就必然是核心内容。

已经看到，黎曼曲率张量有4个指标，这也意味着它有 $n^4$个分量，$n$ 是空间的维数。那么在2、3、4维空间中，它就有16、81、256个分量了，可见，要计算它，是一件相当痛苦的事情。幸好，这个张量有很多的对称性质，使得独立分量的数目大大减少，我们来分析这一点。

首先我们来导出黎曼曲率张量的一些对称性质，这部分内容是跟经典教科书是一致的。定义：
$$R_{\mu \alpha \beta \gamma }=g_{\mu \nu }{R}_{\alpha \beta \gamma }^{\nu }$$

定义这个量的原因，要谈及逆变张量和协变张量的区别，我们这里主要关心几何观，因此略过对张量的详细分析。这个量被称为完全协变的黎曼曲率张量，有时候也直接叫做黎曼曲率张量，只要不至于混淆，一般不做区分。通过略微冗长的代数运算（在一般的微分几何、黎曼几何或者广义相对论教材中都有），可以得到
$$\begin{array}{rl}& R_{\mu \alpha \beta \gamma }=-R_{\mu \alpha \gamma \beta }\\ & R_{\mu \alpha \beta \gamma }=-R_{\alpha \mu \beta \gamma }\\ & R_{\mu \alpha \beta \gamma }=R_{\beta \gamma \mu \alpha }\\ & R_{\mu \alpha \beta \gamma }+R_{\mu \beta \gamma \alpha }+R_{\mu \gamma \alpha \beta }=0\end{array}$$

前两个式子很快就告诉我们，如果 $\beta =\gamma$ 或 $\mu =\alpha$，那么 $R_{\mu \alpha \gamma \beta }=0$，这是反对称的结果。下面的计数表明，独立分量的数目还可以以进一步减少。我们可以将每个等式看成一个约束，有一个约束，独立分量就减少1，我们要考虑有多少个独立的约束，来得出有多少个独立分量。

首先，$R_{\mu \alpha \beta \gamma }$可以看成一个 $n\times n$ 矩阵（前两个指标），由这个矩阵的每个元素都是一个 $n\times n$ 矩阵（后两个指标），由$R_{\mu \alpha \beta \gamma }=-R_{\mu \alpha \gamma \beta }$ 可知，作为矩阵每个元素的矩阵，是一个反对称矩阵，因此只有 $n(n-1)/2$个独立分量，而由 $R_{\mu \alpha \beta \gamma }=-R_{\alpha \mu \beta \gamma }$ 知道，作为矩阵的矩阵，它也是反对称的，因此也只有 $n(n-1)/2$ 个独立分量，也就是说，根据前两个式子就可以知道总的独立分量数不超过$[n(n-1)/2]2$。

接着，第三、第四个约束可以为我们进一步减少分量数目。Rμαβγ=Rβγμα表明，如果将剩下的$[n(n-1)/2]2$个分量排成$n(n-1)/2\times n(n-1)/2$的矩阵，那么它是对称的，这时候它的独立分量数目就只有
$$\frac{1}{2}\frac{n(n-1)}{2}\left[\frac{n(n-1)}{2}+1\right]$$

最后一个约束 $R_{\mu \alpha \beta \gamma }+R_{\mu \beta \gamma \alpha }+R_{\mu \gamma \alpha \beta }=0$ 是针对四个不相同指标的，因为一旦有一对相同指标，就能够转化为前三个约束的线性组合。因此，独立分量的数目还需要减少$\left(\genfrac{}{}{0ex}{}{n}{4}\right)$：
$$\frac{1}{2}\frac{n(n-1)}{2}\left[\frac{n(n-1)}{2}+1\right]-\left(\genfrac{}{}{0ex}{}{n}{4}\right)=\frac{n^2(n^2-1)}{12}$$

于是，在2、3、4、5维空间中，黎曼曲率张量的独立分量数目分别是1、6、20、50。特别地，在二维空间中（即研究三维空间中的二维曲面之时），曲率张量只有一个独立分量 。

2、曲率张量的计算

利用Python中的符号计算库SymPy【文档】 可以帮助我们快速地计算联络系数和曲率张量。SymPy是Python的一个符号计算库，具有小巧、开源、强大等特点。SymPy内置了专门用来处理微分几何和张量的模块，这使得它可以非常方便地计算 ${\Gamma }_{\alpha \beta }^{\mu }$ 和 $R_{\alpha \beta \gamma }^{\mu }$。
安装：pip install sympy 或者 conda install -c anaconda sympy

from sympy.diffgeom import Manifold, Patch, CoordSystem
from sympy.diffgeom import TensorProduct as TP
from sympy.diffgeom import metric_to_Riemann_components as Riemann
from sympy.diffgeom import metric_to_Christoffel_2nd as Christoffel
from sympy import Symbol, Function, latex, sin

n = 2
M = Manifold('M', n)
P = Patch('P', M)

coord = CoordSystem('coord', P, ['x%s'%i for i in range(n)])
x = coord.coord_functions()
dx = coord.base_oneforms()
g = [[1, 0], [0, sin(x[0])**2]]
metric = sum([g[i][j]*TP(dx[i], dx[j]) for i in range(n) for j in range(n)])

C = Christoffel(metric)
R = Riemann(metric)

这是以二维球面坐标为例的，其中n定义了维数，g定义了度规矩阵，如果读者要修改，只需要修改这两个参数就好。里边的x变量是坐标的集合，$\mathrm{d}x$变量是微分元的集合。为了一致，这里的指标是从 $0$ 开始的，也就是遍历 $0\sim n-1$。得到结果是

$$\begin{array}{rlrlrl}[& [[0,0],& [0,-\sin (x0)\cdot \cos (x0)]],& [[0,\frac{\cos (x0)}{\sin (x0)}],& [\frac{\cos (x0)}{\sin (x0)},0]]& ]\\ [& [[[0,0],& [0,0]],& [[0,\sin (x0)\cdot \cdot 2],& [-\sin (x0)\cdot \cdot 2,0]]& ],\\ [& [[0,-1],& [1,0]],& [[0,0],& [0,0]]]& ]\end{array}$$

for i1 in range(n):
    for i2 in range(n):
        for i3 in range(n):
            if C[i1, i2, i3]:
                print 'Gamma^%s_%s%s:'%(i1, i2, i3), C[i1, i2, i3]

for i1 in range(n):
    for i2 in range(n):
        for i3 in range(n):
            for i4 in range(n):
                if R[i1, i2, i3, i4]:
                    print 'Riemann^%s_%s%s%s:'%(i1, i2, i3, i4), R[i1, i2, i3, i4]

鸣谢：

---

1. https://zhuanlan.zhihu.com/p/23132541 ↩︎

2. https://www.spaces.ac.cn/archives/3977 ↩︎

3. https://www.spaces.ac.cn/archives/3998 ↩︎

4. https://www.spaces.ac.cn/archives/4014 ↩︎

5. https://www.spaces.ac.cn/archives/4026 ↩︎