# dingwei

**Repository Path**: hejinguo/dingwei

## Basic Information

- **Project Name**: dingwei
- **Description**: 定位算法Java实现，概化算法，外阔算法，拟合算法
- **Primary Language**: Java
- **License**: GPL-2.0
- **Default Branch**: master
- **Homepage**: https://gitee.com/cai_jin_hang/dingwei
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 3
- **Created**: 2023-10-16
- **Last Updated**: 2023-10-16

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# <center>  运动轨迹的概化算法

## <center> 1. 问题背景
在软件开发中，需要较为准确地给出某个建筑物或特定区域的轮廓，在调查了现有的库后，发现不能很好解决需求，所以考虑开发出一套算法。

## <center> 2. 需求说明
人携带手机绕建筑物一周，同时手机开GPS定位，结合相应的小程序，小程序每间隔t毫秒记录一次人此时的经纬度，就可以获得人的运动轨迹，人的轨迹近似成建筑物的轮廓。此时的数据是一个很长的经纬度数组，比如：[(39.916527, 116.397128), (39.916521, 116.397124), ...]。但是我们需要将轨迹拟合成多边形，只需要多边形的顶点经纬度即可。并求出精度范围，将拟合的多边形向外等距扩张得到新多边形以避免误差带来的问题。在之后的使用中，再根据扩张后的多边形顶点经纬度和任意给定的一点P，判断点P是否在建筑物内部，若在建筑物外部，再求出点P和建筑物的最短距离。

## <center> 3. 需求分析
根据以上描述，需要给出以下算法：
1. 经纬度与XY坐标相互转换，要求两次转换之后数值不变或满足精度要求。
2. 将有序离散点列分段直线拟合成多边形。
3. 将多边形向外等距扩张得到新多边形，求出新多边形的顶点坐标。
4. 判断点与多边形的位置关系。
5. 求出点与多边形的最短距离。

## <center>  4. 算法设计
### 4.1 经纬度与坐标转换
#### 4.1.1 经纬度转XY
经纬度本质上是球面上一点P与球心O的连线与赤道的夹角和与回归线的夹角。所以问题抽象成：给定球面上一些位置非常相近的点的球极坐标：
```python
pointListLL = [(lo_1, la_1), (lo_2, la_2), ... , (lo_n, la_n)]
```
将 pointListLL 近似看成在同一个平面$\pi$上，在平面$\pi$上找一个坐标系$XOY_{\pi}$，并在$XOY_{\pi}$内表示pointListLL。显然，平面$\pi$可以是过某一个点的球的切面，然后将其他点投影到$\pi$上即可。具体做法如下：
将第一个点在球面上平移到回归线和赤道的交点处(0, 0)，其他所有点也随之保距平移，得到平移后的点列：
```python
pointListLLTrans[i] = (pointListLL[i][0] - pointListLL[0][0], pointListLL[i][1] - pointListLL[0][1])
```
根据球极坐标的表达式：
$$
f:\left(\theta, \varphi\right)^T \rightarrow \left(x, y, z\right)^T
$$
$$
\begin{equation}
f:\left\{
\begin{aligned}
x&=R\sin\theta\cos\varphi \\
y&=R\sin\theta\sin\varphi \\
z&=R\cos\theta \nonumber
\end{aligned}
\right.
\end{equation}
$$
显然映射$f$是双射，且$f$连续，所以$f$是同胚。又有$\sin\theta\sim\theta, (\theta\rightarrow 0)$，$\cos\theta\approx 1, (\theta\rightarrow 0)$，所以$(x,y,z)$可以近似表示成：
$$
\begin{equation}
\left\{
\begin{aligned}
x&=R\theta \\
y&=R\theta\varphi \\
z&=R \nonumber
\end{aligned}
\right.
\end{equation}
$$
球面$E$上的一点$P$的邻域为$B(P;\varepsilon_1)$，过点$P$的球$E$的切面为$\pi$；根据微分几何，存在$\varepsilon_2$，使得$B(P;\varepsilon_1)\cap E$与$B(P;\varepsilon_2)\cap \pi$同胚。此时得到平面上的点列的坐标：
```python
pointListLLTransXY[i] = K * (pointListLLTrans[i][0], pointListLLTrans[i][0]*pointListLLTrans[i][1])
```
这里K选取一个充分大的数。至此，经纬度转XY坐标完成。

#### 4.1.2 XY转经纬度：
```python
pointListLLTrans[i] = 1/K * (pointListLLTransXY[i][0], pointListLLTransXY[i][1]/pointListLLTransXY[i][0])
pointListLL[i] = (pointListLLTrans[i] + pointListLL[0][0], pointListLLTrans[i] + pointListLL[0][1])
```

### 4.2 多边形的分段直线拟合
分解成三个算法：
* 等距取点
* 选取拐点
* 相邻拐点之间的点进行直线拟合

#### 4.2.1 等距取点
为了避免连续多次出现同一个点的情况，需要每隔 DISTANCE_TWO_POINT 取一个点。

#### 4.2.2 选取拐点
采用 Douglas-Peucker 算法：
(1) 在曲线首尾两点A，B之间连接一条直线AB，该直线为曲线的弦；
(2) 得到曲线上离该直线段距离最大的点C，计算其与AB的距离d；
(3) 比较该距离与预先给定的阈值threshold的大小，如果小于threshold，则该直线段作为曲线的近似，该段曲线处理完毕。
(4) 如果距离大于阈值，则用C将曲线分为两段AC和BC，并分别对两段取信进行1~3的处理。
(5) 当所有曲线都处理完毕时，依次连接各个分割点形成的折线，即可以作为曲线的近似。
![](./readmepic/a.png)

#### 4.2.3 直线拟合
基于垂线段总和最小对最小二乘法的改进，我们假定所求回归直线为：
$$
y=a_{3} x+b_{3}
$$
将我们所有的样本点向回归直线做垂线段，并将此平方总和作为度量指标来替代经典二乘法的距离。据此想法，我们定义S函数如下：
$$
S\left(a_{3}, b_{3}\right)=\sum_{i=1}^{m} \frac{\left(a_{3} x_{i}-y_{i}+b_{3}\right)^{2}}{a_{3}^{2}+1}
$$
而我们所要求的即为满足如下的函数$S$的最小值时的$a_3$和$b_3$根据$S$函数的定义我们知道$S$是一个非负函数，关于$a_3, b_3$有偏导数令其偏导为零有：
$$
\frac{\partial S\left(a_{3}, b_{3}\right)}{\partial a_{3}}=\sum_{i=1}^{m} \frac{2 x_{i}\left(a_{3} x_{i}-y_{i}+b_{3}\right)\left(a_{3}^{2}+1\right)-2 a_{3}\left(a_{3} x_{i}-y_{i}+b_{3}\right)^{2}}{a_{3}^{2}+1}=0
$$
$$
\frac{\partial S\left(a_{3}, b_{3}\right)}{\partial b_{3}}=\sum_{i=1}^{m} \frac{2\left(a_{3} x_{i}-y_{i}+b_{3}\right)}{a_{3}^{2}+1}=0
$$
即：
$$
\sum_{i=1}^{m} a_{3}\left(x_{i}^{2}-y_{i}^{2}\right)+\sum_{i=1}^{m}\left(a_{3}^{2}-1\right) x_{i} y_{i}+\sum_{i=1}^{m} b_{3}\left(1-a_{3}^{2}\right) x_{i}+\sum_{i=1}^{m} 2 a_{3} b_{3} y_{i}=a_{3} b_{3}^{2}
$$
$$
\sum_{i=1}^{m}\left(a_{3} x_{i}-y_{i}+b_{3}\right)=0
$$
$$
a_{3} \bar{x}-\bar{y}+b_{3}=0
$$
$$
a_{3}^{2}(\overline{x y}-\bar{x} \bar{y})+a_{3}\left(\overline{x^{2}-y^{2}}-\bar{x} \bar{x}+\bar{y} \bar{y}\right)+(\bar{x} \bar{y}-\overline{x y})=0
$$
$$
a_{3}=\frac{-\left(\overline{x^{2}-y^{2}}-\bar{x} \bar{x}+\bar{y} \bar{y}\right) \pm \sqrt{\left(\overline{x^{2}-y^{2}}-\bar{x} \bar{x}+\bar{y} \bar{y}\right)^{2}+4(\overline{x y}-\bar{x} \bar{y})^{2}}}{2(\overline{x y}-\bar{x} \bar{y})}((\overline{x y}-\bar{x} \bar{y}) \neq 0)
$$
$$
b_{3}=\bar{y}-a_{3} \bar{x}
$$

### 4.3 多边形外扩
![](./readmepic/b.png)
$P_1, P_2, P_3 $是多边形按顺序的三个点，现要将边$P_1P_2$向外扩张$d$，边$P_2P_3$也向外扩张$d$，得到点$P_4$，对应于点$P_2$。作线段$P_1P_2$的垂线，在垂线上选取一点$C$，使得$\left|P_2C\right|=d$；作线段$P_3P_2$的垂线，在垂线上选取一点$D$，使得$\left|P_2D\right|=d$；过$C$点作$P_1P_2$的平行线，过$D$点作$P_2P_3$的平行线，两条平行线相交于$P_4$。延长$P_3P_2$交$P_4C$于$A$，延长$P_1P_2$交$P_4D$于$B$。令向量$e_1$是P2A上的单位向量，$e_2$是$P_2B$上的单位向量。设$P_2A=ke_1$。
$$\begin{aligned}
&\because & P_2C \perp AP_4 \\
&& P_2D \perp  BP_4 \\
&& \left|P_2C\right|=\left|P_2D\right|=d\\
&\therefore & P_2P_4是\angle AP_4B的角平分线 \\
&\therefore & 平行四边形AP_4BP_2是菱形\\
&\therefore &P_2B=ke_2
\end{aligned}$$
$$\begin{aligned}
&\because S_{AP_4BP_2}&=\left|P_2A \times P_2B \right|=k^2\left| e_1 \times e_2\right|\\
&&=2S_{AP_4P_2}=kd \\
&\therefore k=\frac{d}{\left|e_1\times e_2\right|}\\
&e_1=\frac{P_3P_2}{\left|P_3P_2\right|} \\
&P_2P_4=k(e_1+e_2)
\end{aligned}$$
令$P_i$坐标是$(x_i, y_i)$, $i=1, 2, 3, 4$，则：
$$\begin{aligned}
x_4&=x_2+\frac{d}{\left| (x_2-x_3)(y_2-y_1)-(x_2-x_1)(y_2-y_3)\right|}\left((x_2-x_1)\sqrt{(x_2-x_3)^2+(y_2-y_3)^2}+(x_2-x_3)\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}\right) \\
y_4&=y_2+\frac{d}{\left| (x_2-x_3)(y_2-y_1)-(x_2-x_1)(y_2-y_3)\right|}\left((y_2-y_1)\sqrt{(x_2-x_3)^2+(y_2-y_3)^2}+(y_2-y_3)\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}\right)
\end{aligned}$$
如果要得到内缩的点的坐标，只需要将$k$加负号即可。

### 4.4 判断点与多边形的位置关系
#### 4.4.1 射线法
射线法是解决这一问题的最优方法，射线法的思想是：以目标点为端点引一条射线，计算这条射线和多边形所有边的交点数目。如果交点个数为奇数，则点在多边形部，反之则在多边形外部。图例说明，如下图所示：
![](./readmepic/c.png)
所谓射线法，关键在于单向发射，为简化问题，以水平线为例，程序实现中也是这么处理的。O点向右发出射线，与多边形的交点是B、C、D，向左发出射线，交点是A，均为奇数个。P点在多边形外，无论想哪方向发出摄像，都有2个交点。
![](./readmepic/d.png)
对于带内岛的形状，射线法同样适用，如上图所示。在实际应用中，射线法会有很多特殊情况需要讨论，全部都讨论会比较复杂，但结论是一样的。这里不做过多讨论了，不过可以给大家结论：射线法适用于所有类型的多边形进行点与多边形关系的判断，且实现相对简单，速度较快，是工程应用的不二之选。要注意的是，计算中所有数值都要选择浮点数类型，以保证计算精度。

#### 4.4.2 关于外扩算法的补充
因为不清楚多边形的点是凸点还是凹点，所以将所有点全部外扩是不行的，如果是凸点，就用外扩算法，如果是凹点，就用内缩算法。$P_{i-1}, P_i, P_{i+1}$是多边形的相邻的三个点，连接$P_{i-1},P_{i+1}$，得到新的多边形，再判断$P_i$与新多边形的位置关系，如果在内部，则为凹点，如果在外部，则为凸点。

### 4.5 求出点与多边形的最短距离
此时假设点P在多边形的外部，则P到多边形的最短距离就是P到多边形的各个边的距离的最小值，注意，多边形的边是线段，所以最短距离是到线段的最短距离，而不是到边所在的直线的最短距离。下面加以证明：
存在性：若 $F \subset \mathbf{R}^{n}$ 是非空闭集, 且 $x_{0} \in \mathbf{R}^{n}$, 则存在 $y_{0} \in F$, 有
$$
\left|x_{0}-y_{0}\right|=d\left(x_{0}, F\right) .
$$
证明： 作闭球 $\bar{B}=\bar{B}\left(x_{0}, \delta\right)$, 使得 $\bar{B} \cap F$ 不是空集. 显然
$$
d\left(x_{0}, F\right)=d\left(x_{0}, \bar{B} \cap F\right) .
$$
$\bar{B} \cap F$ 是有界闭集, 而 $\left|x_{0}-y\right|$ 看作定义在 $\bar{B} \cap F$ 上的 $y$ 的函数是连续的, 故它在 $\bar{B} \cap F$ 上达到最小值, 即存在 $y_{0} \in \bar{B} \cap F$, 使得
$$
\left|x_{0}-y_{0}\right|=\inf \left\{\left|x_{0}-y\right|: y \in \bar{B} \cap F\right\},
$$
从而有 $\left|x_{0}-y_{0}\right|=d\left(x_{0}, F\right)$.

再证明$y_0\in \partial F$：
由$F$是闭集可得$F=\partial F \cup \mathring{F}$，若$y_0\in \mathring{F}$，连接$x_0y_0$交$\partial F$于$y_0^{'}$，则$y_0^{'}\in F$， 且$|x_0y_0|>|x_0y_0^{'}|$，所以$y_0\in \partial F$

然而若$y_0$在边$P_iP_{i+1}$上但是不在顶点上时，显然有有$x_0y_0 \perp P_iP_{i+1}$，至此，结论成立。

## <center>  5. 调整参数

## <center>  6. 实际效果

## <center>  7. 算法评价
### 7.1 复杂度
|  算法   | 时间复杂度  | 空间复杂度 |
|  ----  | ----  | ---- |
| 算法4.2.1  | $O(n)$ | $O(n)$ |
| 算法4.2.2  | $O(n^2)$ | $O(n^2)$ |
| 算法4.3  | $O(n^2)$ | $O(n^2)$ |
| 算法4.4.1  | $O(n)$ | $O(n)$ |
| 算法4.5 | $O(n)$ | $O(n)$ |
|整体算法| $O(n^2)$ | $O(n^2)$ |

### 7.2 算法不足
算法对于有自交的运动轨迹和有重叠的运动轨迹不能很好的适应，所以要求携带手机的人在采点时不要原地绕圈，不要突然向相反方向运动，且绕建筑物一圈后恰好回到出发点。
在判断点与多边形位置关系的算法中，射线不能通过多边形的顶点，否则算法得到错误的结果。

## <center> 8. 参考文献
[1]钟业勋，胡宝清著. 数理地图学 地图学及其数学原理 第2版. 北京：测绘出版社, 2017.10.
[2]陈明.Douglas-Peucker算法在电子海图岛屿轮廓提取的应用[J].电子技术,2020,49(05):36-37.
[3]朱春龙,杨诚芳.建立线性相关方程的最小距离平方和法[J].人民长江,2001(09)
[4]向俊,王静,夏幼明.判断点与多边形拓扑关系的改进算法[J].计算机工程与设计,2014,35(05):1732-1737.DOI:10.16208/j.issn1000-7024.2014.05.049.
:47-49.DOI:10.16232/j.cnki.1001-4179.2001.09.020.
[5]周民强编著. 实变函数 第2版. 北京：北京大学出版社, 1995.06.
[6] A·C·米先柯（Мищенко А. С.），A·T·福明柯（Фоменко А. Т.）著；张爱和译. 微分几何与拓扑学简明教程. 北京：高等教育出版社, 2006.01.