道格拉斯-普克算法（Douglas–Peucker）考虑投影在线段之外的实现道格拉斯-普克算法(Douglas–Peuc

[本文主要参考]

如果原理已经很清楚了，可以直接跳过看代码

算法描述

关键点：如何判断点到线段的投影在线段之外

要判断一个点 $P$ 在二维平面上对某条线段 $AB$ 的垂直投影是否在线段 $AB$ 上，可以通过以下步骤进行计算：

定义线段：假设线段 $AB$ 的两个端点分别是 $A(x_1, y_1)$ 和 $B(x_2, y_2)$ ，我们要找的是点 $P(x_p, y_p)$ 对这条线段的垂直投影。
参数方程：线段可以表示为参数方程形式： $x = x_1 + t(x_2 - x_1)$ $y = y_1 + t(y_2 - y_1)$ 其中 $t$ 是参数，当 $t=0$ 时，得到点 $A$ ，当 $t=1$ 时，得到点 $B$ 。因此， $t$ 的有效范围是 [0, 1]。
点到直线的距离公式：可以使用点到直线的距离公式来找到点 $P$ 到线段 $AB$ 所在直线的最短距离，这个距离就是点 $P$ 到直线 $AB$ 的垂直距离。如果这个距离上的点（即垂直投影）位于线段 $AB$ 上，则 $t$ 应该在 [0, 1] 范围内。
求解 t：通过向量运算可以得到 $t$ 的值： $t = \frac{(x_p - x_1)(x_2 - x_1) + (y_p - y_1)(y_2 - y_1)}{(x_2 - x_1)^2 + (y_2 - y_1)^2}$
判断投影位置：如果 $t < 0$ 或者 $t > 1$ ，则说明点 $P$ 的垂直投影不在线段 $AB$ 上。如果 $0 \leq t \leq 1$ ，则垂直投影在线段 $AB$ 上。
计算投影点坐标：投影点 $Q$ 的坐标可以通过 $t$ 计算得出： $Q(x_q, y_q) = (x_1 + t(x_2 - x_1), y_1 + t(y_2 - y_1))$

关于t的公式

对于参数t，其实是用向量投影来推倒的（高中知识）

要求解 $t$ ，我们可以通过点 $P$ 到线段 $AB$ 的最近点 $Q$ 的位置来确定。这个最近点实际上是 $P$ 到直线 $AB$ 的垂直投影点。为了找到这个点，我们需要确保从 $P$ 到 $Q$ 的向量与从 $A$ 到 $B$ 的向量垂直，并且 $Q$ 在线段 $AB$ 上或其延长线上。

给定点 $P(x_p, y_p)$ ，线段 $AB$ 的端点 $A(x_1, y_1)$ 和 $B(x_2, y_2)$ ，我们可以将线段 $AB$ 表示为一个向量 $\vec{AB}$ ，以及一个从 $A$ 到 $P$ 的向量 $\vec{AP}$ 。

$\vec{AB} = (x_2 - x_1, y_2 - y_1)$ $\vec{AP} = (x_p - x_1, y_p - y_1)$

投影 $P$ 到 $AB$ 的公式基于向量的点积（标量积）。点积的几何意义是在一条直线上投影另一个向量的长度乘以这个向量的长度。因此，我们可以利用点积来找出 $AP$ 在 $AB$ 方向上投影的大小：

$\text{proj}_{\vec{AB}} \vec{AP} = \frac{\vec{AP} \cdot \vec{AB}}{\|\vec{AB}\|^2} \vec{AB}$

其中 $\vec{AP} \cdot \vec{AB}$ 是两个向量的点积， $\|\vec{AB}\|$ 是向量 $\vec{AB}$ 的模长。

$\vec{AP} \cdot \vec{AB} = (x_p - x_1)(x_2 - x_1) + (y_p - y_1)(y_2 - y_1)$ $\|\vec{AB}\|^2 = (x_2 - x_1)^2 + (y_2 - y_1)^2$

将这些值代入点积公式中，我们可以求得 $t$ ：

$t = \frac{\vec{AP} \cdot \vec{AB}}{\|\vec{AB}\|^2} = \frac{(x_p - x_1)(x_2 - x_1) + (y_p - y_1)(y_2 - y_1)}{(x_2 - x_1)^2 + (y_2 - y_1)^2}$

这就是求解 $t$ 的公式。注意，如果 $t$ 的值小于 0 或大于 1，那么点 $P$ 的垂直投影不在线段 $AB$ 上。如果 $t$ 在区间 [0, 1] 内，则点 $P$ 的垂直投影在线段 $AB$ 上。

代码实现（C语言）

//  Douglas-Peuker  实际上还是递归 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct{
	double x;
	double y;
}Point;

// 由于代码会顺序简化线段，可以设置一个全局变量将结果依次放入该数组中 
Point p[10];
int n = 0; 

// 拿到较小的那个
double min(double a, double b){  
    return (a > b) ? b : a;  
}  

// 获取两点之间距离 
double getDistance(Point *p1, Point *p2){
	return sqrt(pow(p1->x - p2->x, 2) + pow(p1->y - p2->y, 2));
}

// 海伦公式由三边求三角形面积 
double helen(double a, double b, double c){
	double p = (a + b + c) / 2;
	return sqrt(p*(p-a)*(p-b)*(p-c));
}

// 获得点到线段的距离（根据点的投影在不在线段上分为两种情况） 
double getPointToLine(Point *p, Point *p1, Point *p2){
	// 首先判断投影情况
	double t = ((p->x - p1->x)*(p2->x - p1->x) + (p->y - p1->y)*(p2->y - p1->y)) / ( pow(p2->x - p1->x,2) + pow(p2->y - p1->y,2) );
	
	double a = getDistance(p,p1);
	double b = getDistance(p,p2);
	double c = getDistance(p1,p2);
	 
	if(t < 0 || t > 1){
		// t的取值已经隐含了p离p1近还是离p2近，这里待优化 
		return min(a, b);
	}else{
		double s = helen(a,b,c);
		return 2 * s / c;
	}
}

// 限差为D 
void DouglasPeuker(Point *oldPoints, int oldCount, double D){
	Point head = oldPoints[0];
	Point tail = oldPoints[oldCount-1];
	
	double dmax = 0.0;
	int index = 0;
	
	for(int i=1; i<oldCount-1; i++){
		double d = getPointToLine(&oldPoints[i], &head, &tail);
		if(d > dmax){
			dmax = d;
			index = i;
		}
	}
	
	if(dmax >= D){
		Point* part1 = (Point*)malloc(sizeof(Point)*(index+1));
		Point* part2 = (Point*)malloc(sizeof(Point)*(oldCount-index));
		
		int j = 0; // 记录part2应当将元素插入的位置 
		for(int i=0; i<oldCount; i++){
			if(i<index){
				part1[i] = oldPoints[i];
			}else if(i>index){
				part2[j++] = oldPoints[i];
			}else{
				part1[i] = oldPoints[i];
				part2[j++] = oldPoints[i];
			}
		}
		DouglasPeuker(part1, index+1, D);
		DouglasPeuker(part2, oldCount-index, D);
	
	}else{
		p[n++] = tail;
	}
}

int main(){
	Point points[10] = {{1.0, 1.0}, {2.0, 2.0}, {3.0, 4.0}, {4.0, 1.0}, {5.0, 0.0}, {6.0, 3.0}, {10.0, 11.0}, {8.0, 2.0}, {9.0, 1.0}, {10.0, 6.0}};

	// 第一个节点肯定不会被简化 
	p[n++] = points[0];
	DouglasPeuker(points, 10, 1.0);

	for(int i=0; i<n; i++){
		printf("(%.1lf,%.1lf) ", p[i].x, p[i].y);
	}

	return 0;
}