距离及其近似计算

距离及其近似计算

---

1.导读
本文简要推导在一些需要计算距离，计算资源又不提供开平方运算的情况下，如何逼近常用的欧式距离？
在泛函分析中，距离通常由范数诱导，因此在信号处理中常用的距离为基于$L^p$范数的距离，这里给出$L^p$范数的定义。

$$d(p_1,p_2)=\sqrt[p]{(x_1-x_2)^p+(y_1-y_2)^p)}$$
通常的欧式距离，即为取$p=2$的一个特例
$$d(p_1,p_2)=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2)}$$
在ISP的pipeline或者其他通过ASIC逻辑电路来计算距离的计算系统中，开方运算功耗较大，需要较大的面积，因此常常采用基于$L^1$范数的距离度量，也即棋盘距离：
$$d(p_1,p_2)=|x_1-x_2|+|y_1-y_2|$$
对于lens shading correction或者沿径向去噪等对距离要求较高应用，需要对欧式距离有一个很高精度的逼近，因此需要考虑对欧式距离进一步逼近。

2.具体内容
这里提供两种对欧式距离逼近的策略。
2.1 采用两个正方形及其旋转组合成的八边形逼近
八边形示意图如下：

原理与解释：

$$\begin{align*} d(p_1, p_2) &= \lVert x_1-x_2 \rVert_{oct} \\ &= \left\{ \begin{array}{ll} &|x_1-x_2|+|y_1-y_2| &, |x_1-x_2|>|y_1-y_2| tan22.5 \quad\&\&\quad |y_1-y_2|>|x_1-x_2| tan22.5 \\ &\sqrt{2}max(|x_1-x_2|,|y_1-y_2|) &, otherwise \end{array} \right. \end{align*}$$

推导过程：

如上图所示：八边形分为蓝色和红色部分，以第一象限为例：蓝色部分点在直线$|x|+|y|=\sqrt{2}l$，蓝色部分到八边形中心的距离为：

$$\begin{align*} d_{blue} & =|x-0|+|y-0| \\ & = |x|+|y| \\ & = \sqrt{2}l \end{align*}$$
红色部分：

$$\left\{ \begin{array}{ll} |x|&= l &, when|x|>|y| \\ |y|&= l &, when|y|>|x| \end{array} \right.$$

要使红色部分到中心的距离和蓝色部分到中心的距离相等

$$\begin{align*} d_{red} & = d_{blue} \\ & = \sqrt{2}l \\ & = \sqrt{2}max(|x|, |y|) \end{align*}$$
从而八边形距离公式如本部分开始所述。

示例代码（MATLAB）

function d = L2_distance_oct(x1, y1, x2, y2)

dx_abs = abs(x1 - x2);
dy_abs = abs(y1 - y2);

if ( dx_abs > dy_abs * tand(22.5) ) && ( dy_abs > dx_abs * tand(22.5) )
  d = dx_abs + dy_abs;
else
  d = sqrt(2) * max(dx_abs, dy_abs);
end


end

2.2 斜面距离逼近

斜面距离示意图如下：

原理与解释：
斜面距离是对邻域中欧式距离的整数近似。从像素$p$到其4-邻域的像素只需要水平移动或者垂直移动（称为a-move）。因为所有移动根据对称或旋转都是相等的，对离散距离的唯一可能的定义就是$d_4$距离，其中$a=1$。从像素$p$到其8-邻域的像素不仅有水平移动或垂直移动，还有对角移移动（称为b-move）。同时考虑这两种移动，可把斜面距离记为$d_{a,b}$。最自然的$b$值为$\sqrt{2}a$，但是为了计算的简单和减少存储量，将$a$值和$b$值都取整数。最常用的一组值是$a=3$和$b=4$。这组值可如下获得。考虑像素$p_1$和$p_2$之间在水平方向相差的像素个数为$n_x$，在垂直方向相差像素个数为$n_y$（不失一般性，设$n_x > n_y$），则像素$p_1$和$p_2$之间的斜面距离为

$$D(p_1,p_2)=(n_x-n_y)a+n_yb$$
它与欧式距离之差为：
$$\Delta D(p_1,p_2)=\sqrt{n_x^2+n_y^2}-[(n_x-n_y)a+n_yb]$$
如果取$a=1$，$b=\sqrt{2}$,将$\Delta D(p_1,p_2)$对$n_y$求导数得
$$\Delta D^{'}(p_1,p_2)=\frac{n_y}{\sqrt{n_x^2+n_y^2}}-(\sqrt{2}-1)$$
令导数为零，计算$\Delta D(p_1,p_2)$的极值可得
$$n_y=\sqrt{(\sqrt{2}-1)/2}n_x$$
即两距离之差$\Delta D(p_1,p_2)$在满足上式时区最大值$(\sqrt{(\sqrt{2}-1)/2}-1)n_x \approx -0.09n_x$(此时直线与横轴间夹角为$24.5^{\circ}$)。进一步可证明当$b=1/2^{1/2}+(2^{1/2}-1)=1.351$时$\Delta D(p_1,p_2)$的最大值可达到最小。因为$4/3 \approx 1.33$，随意斜面距离取$a=3$和$b=4$
此时斜面距离的计算公式为：
$$\begin{align*} D(p_1,p_2)= \left\{ \begin{array}{ll} (n_x-n_y)×3+n_y×4 &, \quad when\quad n_x>n_y \\ (n_y-n_x)×3+n_x×4 &, \quad otherwise \end{array} \right. \end{align*}$$
示例代码(MATLAB)

function d = L2_distance_slope(x1, y1, x2, y2)

dx_abs = abs(x1 - x2);
dy_abs = abs(y1 - y2);

if dx_abs > dy_abs
  d = (dx_abs - dy_abs) * 3 + 4 * dy_abs;
else
  d = (dy_abs - dx_abs) * 3 + 4 * dx_abs; 
end

end

综上，本文提供了两种对$L^1$的距离的菱形区域的改进方案，为需要对图像沿半径进行区域划分提供了不同的选择，根据计算资源选择对应的方案。

参考文献：
1.图像工程（中册），第三版，p10