算法导论第十八章 计算几何：算法中的空间艺术
第十八章 计算几何：算法中的空间艺术
“几何学是描绘宇宙秩序的永恒诗篇。” —— 约翰内斯·开普勒
计算几何将数学的优雅与算法的实用性完美结合，在计算机图形学、机器人导航和地理信息系统中扮演着关键角色。本章将带您探索几何问题的算法解决方案，从基础的点线关系到复杂的空间剖分，揭示算法如何理解和操纵我们的几何世界。
18.1 几何基础：点、线和多边形
18.1.1 几何对象的表示
在计算几何中，我们使用简洁的数学结构表示几何对象：
// 二维点
typedef struct {double x;double y;
} Point;// 线段
typedef struct {Point start;Point end;
} Segment;// 多边形
typedef struct {int num_vertices;Point *vertices; // 顶点数组
} Polygon;
18.1.2 基本几何运算
点积与叉积
// 点积：a·b = |a||b|cosθ
double dot_product(Point a, Point b) {return a.x * b.x + a.y * b.y;
}// 叉积：a×b = |a||b|sinθ
double cross_product(Point a, Point b) {return a.x * b.y - a.y * b.x;
}
方向判断
// 判断点c在ab线段的哪一侧
int direction(Point a, Point b, Point c) {double cross = cross_product((Point){b.x - a.x, b.y - a.y}, (Point){c.x - a.x, c.y - a.y});if (cross > 0) return 1; // 逆时针方向（左侧）if (cross < 0) return -1; // 顺时针方向（右侧）return 0; // 共线
}
18.2 凸包问题：寻找最小包围壳
18.2.1 凸包应用场景
| 领域 | 应用 | 重要性 |
|---|---|---|
| 计算机视觉 | 物体轮廓识别 | ⭐⭐⭐⭐ |
| 路径规划 | 机器人导航区域 | ⭐⭐⭐⭐ |
| GIS系统 | 地图区域简化 | ⭐⭐⭐⭐ |
| 碰撞检测 | 快速排斥测试 | ⭐⭐⭐⭐ |
18.2.2 Graham扫描算法
算法步骤：
- 找到最低点（y最小，相同则x最小）
- 按极角排序其他点
- 扫描点集，维护凸包栈
// 比较函数：按极角排序
int compare_angles(const void *a, const void *b, void *pivot_ptr) {Point *p0 = pivot_ptr;Point *p1 = (Point *)a;Point *p2 = (Point *)b;double cross = cross_product((Point){p1->x - p0->x, p1->y - p0->y},(Point){p2->x - p0->x, p2->y - p0->y});if (cross > 0) return -1;if (cross < 0) return 1;// 共线时距离较近者优先double dist1 = (p1->x - p0->x)*(p1->x - p0->x) + (p1->y - p0->y)*(p1->y - p0->y);double dist2 = (p2->x - p0->x)*(p2->x - p0->x) + (p2->y - p0->y)*(p2->y - p0->y);return (dist1 < dist2) ? -1 : 1;
}// Graham扫描算法
Polygon graham_scan(Point *points, int n) {// 找到最低点int min_index = 0;for (int i = 1; i < n; i++) {if (points[i].y < points[min_index].y || (points[i].y == points[min_index].y && points[i].x < points[min_index].x)) {min_index = i;}}// 交换最低点到数组首位Point temp = points[0];points[0] = points[min_index];points[min_index] = temp;// 按极角排序qsort_r(points + 1, n - 1, sizeof(Point), compare_angles, &points[0]);// 构建凸包Point *hull = (Point *)malloc(n * sizeof(Point));int hull_size = 0;hull[hull_size++] = points[0];hull[hull_size++] = points[1];hull[hull_size++] = points[2];for (int i = 3; i < n; i++) {while (direction(hull[hull_size-2], hull[hull_size-1], points[i]) <= 0) {hull_size--; // 移除非凸点}hull[hull_size++] = points[i];}Polygon convex_hull;convex_hull.num_vertices = hull_size;convex_hull.vertices = hull;return convex_hull;
}
18.2.3 Jarvis步进算法
算法思想：礼品包装法
- 找到最左点作为起点
- 重复寻找下一个凸包点（极角最小）
- 直到回到起点
Polygon jarvis_march(Point *points, int n) {if (n < 3) return (Polygon){0, NULL};// 找到最左点int leftmost = 0;for (int i = 1; i < n; i++) {if (points[i].x < points[leftmost].x) {leftmost = i;}}int *hull_indices = (int *)malloc(n * sizeof(int));int hull_size = 0;int current = leftmost;do {hull_indices[hull_size++] = current;int next = (current + 1) % n;for (int i = 0; i < n; i++) {if (direction(points[current], points[next], points[i]) < 0) {next = i;}}current = next;} while (current != leftmost);// 创建凸包多边形Point *hull = (Point *)malloc(hull_size * sizeof(Point));for (int i = 0; i < hull_size; i++) {hull[i] = points[hull_indices[i]];}free(hull_indices);Polygon convex_hull;convex_hull.num_vertices = hull_size;convex_hull.vertices = hull;return convex_hull;
}
18.2.4 凸包算法比较
| 算法 | 时间复杂度 | 空间复杂度 | 适用场景 |
|---|---|---|---|
| Graham扫描 | O(n log n) | O(n) | 点集较大 |
| Jarvis步进 | O(nh) | O(n) | 凸包点少(h小) |
| QuickHull | O(n log n) | O(n) | 平均性能好 |
| Chan算法 | O(n log h) | O(n) | 理论最优 |
18.3 最近点对问题：分治的艺术
18.3.1 问题描述
给定平面上的n个点，找到距离最近的两个点。
18.3.2 分治算法实现
// 计算两点距离
double distance(Point a, Point b) {double dx = a.x - b.x;double dy = a.y - b.y;return sqrt(dx*dx + dy*dy);
}// 暴力求解（小规模）
double brute_force(Point *points, int n, Point *pair) {double min_dist = INFINITY;for (int i = 0; i < n; i++) {for (int j = i+1; j < n; j++) {double dist = distance(points[i], points[j]);if (dist < min_dist) {min_dist = dist;pair[0] = points[i];pair[1] = points[j];}}}return min_dist;
}// 分治算法
double closest_pair(Point *points_x, Point *points_y, int n, Point *pair) {if (n <= 3) {return brute_force(points_x, n, pair);}int mid = n / 2;Point mid_point = points_x[mid];// 按y排序的点分成左右两部分Point *left_y = (Point *)malloc(mid * sizeof(Point));Point *right_y = (Point *)malloc((n - mid) * sizeof(Point));int left_idx = 0, right_idx = 0;for (int i = 0; i < n; i++) {if (points_y[i].x <= mid_point.x && left_idx < mid) {left_y[left_idx++] = points_y[i];} else {right_y[right_idx++] = points_y[i];}}// 递归求解左右子集Point left_pair[2], right_pair[2];double left_dist = closest_pair(points_x, left_y, mid, left_pair);double right_dist = closest_pair(points_x + mid, right_y, n - mid, right_pair);double min_dist = left_dist;if (right_dist < min_dist) {min_dist = right_dist;pair[0] = right_pair[0];pair[1] = right_pair[1];} else {pair[0] = left_pair[0];pair[1] = left_pair[1];}// 检查跨分割线点对Point *strip = (Point *)malloc(n * sizeof(Point));int strip_size = 0;for (int i = 0; i < n; i++) {if (fabs(points_y[i].x - mid_point.x) < min_dist) {strip[strip_size++] = points_y[i];}}// 检查带状区域内点对for (int i = 0; i < strip_size; i++) {for (int j = i+1; j < strip_size && (strip[j].y - strip[i].y) < min_dist; j++) {double dist = distance(strip[i], strip[j]);if (dist < min_dist) {min_dist = dist;pair[0] = strip[i];pair[1] = strip[j];}}}free(left_y);free(right_y);free(strip);return min_dist;
}
18.3.3 算法可视化
18.4 线段相交检测
18.4.1 应用场景
- 电路设计：检查导线交叉
- 游戏开发：碰撞检测
- CAD系统：几何约束验证
- 地图导航：路径交叉检测
18.4.2 线段相交算法
// 检查点c是否在线段ab上
bool on_segment(Point a, Point b, Point c) {return (c.x <= fmax(a.x, b.x) && (c.x >= fmin(a.x, b.x)) &&(c.y <= fmax(a.y, b.y)) &&(c.y >= fmin(a.y, b.y));
}// 检查两线段是否相交
bool segments_intersect(Point p1, Point p2, Point p3, Point p4) {int d1 = direction(p3, p4, p1);int d2 = direction(p3, p4, p2);int d3 = direction(p1, p2, p3);int d4 = direction(p1, p2, p4);// 一般相交情况if (((d1 > 0 && d2 < 0) || (d1 < 0 && d2 > 0)) &&((d3 > 0 && d4 < 0) || (d3 < 0 && d4 > 0))) {return true;}// 特殊情况：共线点if (d1 == 0 && on_segment(p3, p4, p1)) return true;if (d2 == 0 && on_segment(p3, p4, p2)) return true;if (d3 == 0 && on_segment(p1, p2, p3)) return true;if (d4 == 0 && on_segment(p1, p2, p4)) return true;return false;
}
18.4.3 扫描线算法
高效检测多线段相交：
typedef struct {double x;int type; // 0:起点, 1:终点, 2:交点Segment *segment;
} Event;void sweep_line_intersection(Segment *segments, int n) {Event *events = (Event *)malloc(2 * n * sizeof(Event));int event_count = 0;// 创建事件：起点和终点for (int i = 0; i < n; i++) {// 确保起点在终点左侧if (segments[i].start.x > segments[i].end.x) {Point temp = segments[i].start;segments[i].start = segments[i].end;segments[i].end = temp;}events[event_count++] = (Event){segments[i].start.x, 0, &segments[i]};events[event_count++] = (Event){segments[i].end.x, 1, &segments[i]};}// 按x坐标排序事件qsort(events, event_count, sizeof(Event), event_compare);// 扫描线状态：当前活跃线段Segment **active_segments = (Segment **)malloc(n * sizeof(Segment *));int active_count = 0;// 处理事件for (int i = 0; i < event_count; i++) {Event e = events[i];if (e.type == 0) { // 线段起点// 将新线段插入活跃集active_segments[active_count++] = e.segment;// 检查新线段与活跃线段的交点for (int j = 0; j < active_count - 1; j++) {if (segments_intersect(e.segment->start, e.segment->end,active_segments[j]->start, active_segments[j]->end)) {printf("线段 %d 和 %d 相交\n", e.segment->id, active_segments[j]->id);}}} else if (e.type == 1) { // 线段终点// 从活跃集中移除for (int j = 0; j < active_count; j++) {if (active_segments[j] == e.segment) {// 将最后一个元素移到当前位置active_segments[j] = active_segments[active_count - 1];active_count--;break;}}}}free(active_segments);free(events);
}
18.5 Voronoi图与Delaunay三角剖分
18.5.1 Voronoi图：空间分割的艺术
每个点拥有一个区域，区域内的点距离该点比距离其他点更近。
graph TDA[点集P] --> B[平面分割]B --> C[每个点对应一个区域]C --> D[区域边界是垂直平分线]D --> E[应用：最近邻查询]
18.5.2 Delaunay三角剖分
Voronoi图的对偶图，满足：
- 外接圆不包含其他点（空圆特性）
- 最大化最小角（避免瘦长三角形）
typedef struct {Point vertices[3]; // 三角形顶点Triangle *neighbors[3]; // 相邻三角形
} Triangle;// 增量算法框架
void delaunay_triangulation(Point *points, int n) {// 创建超级三角形包含所有点Triangle super_triangle = create_super_triangle(points, n);// 初始化三角剖分TriangleList *triangles = create_triangle_list();add_triangle(triangles, super_triangle);// 逐点插入for (int i = 0; i < n; i++) {Point p = points[i];TriangleList *bad_triangles = find_bad_triangles(triangles, p);Polygon polygon = create_polygon_hole(bad_triangles);// 移除坏三角形remove_triangles(triangles, bad_triangles);// 重新三角剖分triangulate_polygon(p, polygon, triangles);}// 移除超级三角形相关三角形remove_super_triangle(triangles, super_triangle);
}
18.5.3 应用对比
| 结构 | 应用领域 | 优势 |
|---|---|---|
| Voronoi图 | 路径规划 | 自然分割空间 |
| Delaunay三角剖分 | 有限元分析 | 高质量网格生成 |
| 两者结合 | 地形建模 | 精确表示地形特征 |
| 无线网络规划 | 优化基站覆盖 |
18.6 完整实现：凸包与点对可视化
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>// [点、线段、多边形定义...]
// [几何运算函数...]
// [Graham扫描实现...]
// [最近点对实现...]// 生成随机点集
Point* generate_random_points(int n, double min_x, double max_x, double min_y, double max_y) {Point *points = (Point *)malloc(n * sizeof(Point));srand(time(NULL));for (int i = 0; i < n; i++) {points[i].x = min_x + (double)rand() / RAND_MAX * (max_x - min_x);points[i].y = min_y + (double)rand() / RAND_MAX * (max_y - min_y);}return points;
}// 可视化点集和凸包
void visualize_points(Point *points, int n, Polygon convex_hull, Point *closest_pair) {printf("点集可视化:\n");printf("点数量: %d\n", n);// 打印点坐标for (int i = 0; i < n; i++) {printf("点 %d: (%.2f, %.2f)\n", i, points[i].x, points[i].y);}// 打印凸包printf("\n凸包顶点 (%d 个):\n", convex_hull.num_vertices);for (int i = 0; i < convex_hull.num_vertices; i++) {printf("顶点 %d: (%.2f, %.2f)\n", i, convex_hull.vertices[i].x, convex_hull.vertices[i].y);}// 打印最近点对if (closest_pair) {printf("\n最近点对:\n");printf("点1: (%.2f, %.2f)\n", closest_pair[0].x, closest_pair[0].y);printf("点2: (%.2f, %.2f)\n", closest_pair[1].x, closest_pair[1].y);printf("距离: %.4f\n", distance(closest_pair[0], closest_pair[1]));}
}// 简单ASCII可视化
void ascii_visualization(Point *points, int n, Polygon convex_hull, Point *closest_pair) {const int width = 50;const int height = 20;char grid[height][width];// 初始化网格for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++) {grid[y][x] = ' ';}}// 绘制点for (int i = 0; i < n; i++) {int x = (int)((points[i].x - 0) / 100 * width);int y = (int)((points[i].y - 0) / 100 * height);if (x >= 0 && x < width && y >= 0 && y < height) {grid[y][x] = '.';}}// 绘制凸包for (int i = 0; i < convex_hull.num_vertices; i++) {int x1 = (int)((convex_hull.vertices[i].x - 0) / 100 * width);int y1 = (int)((convex_hull.vertices[i].y - 0) / 100 * height);int x2 = (int)((convex_hull.vertices[(i+1)%convex_hull.num_vertices].x - 0) / 100 * width);int y2 = (int)((convex_hull.vertices[(i+1)%convex_hull.num_vertices].y - 0) / 100 * height);// 简单线段绘制int dx = abs(x2 - x1);int dy = abs(y2 - y1);int steps = dx > dy ? dx : dy;for (int s = 0; s <= steps; s++) {double t = (double)s / steps;int x = (int)(x1 + t * (x2 - x1));int y = (int)(y1 + t * (y2 - y1));if (x >= 0 && x < width && y >= 0 && y < height) {grid[y][x] = '#';}}}// 绘制最近点对if (closest_pair) {int x1 = (int)((closest_pair[0].x - 0) / 100 * width);int y1 = (int)((closest_pair[0].y - 0) / 100 * height);int x2 = (int)((closest_pair[1].x - 0) / 100 * width);int y2 = (int)((closest_pair[1].y - 0) / 100 * height);if (x1 >= 0 && x1 < width && y1 >= 0 && y1 < height) {grid[y1][x1] = 'A';}if (x2 >= 0 && x2 < width && y2 >= 0 && y2 < height) {grid[y2][x2] = 'B';}}// 打印网格printf("\nASCII 可视化:\n");for (int y = height - 1; y >= 0; y--) {printf("%2d |", y);for (int x = 0; x < width; x++) {putchar(grid[y][x]);}printf("\n");}printf(" ");for (int x = 0; x < width; x += 5) {printf("+----");}printf("\n ");for (int x = 0; x < width; x += 5) {printf("%-5d", x);}printf("\n");
}int main() {// 生成随机点集int n = 20;Point *points = generate_random_points(n, 0, 100, 0, 100);// 计算凸包Polygon convex_hull = graham_scan(points, n);// 计算最近点对Point closest_pair[2];Point *points_y = (Point *)malloc(n * sizeof(Point));memcpy(points_y, points, n * sizeof(Point));qsort(points, n, sizeof(Point), compare_x);closest_pair(points, points_y, n, closest_pair);// 可视化结果visualize_points(points, n, convex_hull, closest_pair);ascii_visualization(points, n, convex_hull, closest_pair);// 清理内存free(convex_hull.vertices);free(points);free(points_y);return 0;
}
18.7 计算几何在游戏开发中的应用
18.7.1 碰撞检测系统
// 多边形碰撞检测
bool polygon_collision(Polygon poly1, Polygon poly2) {// 分离轴定理(SAT)for (int i = 0; i < poly1.num_vertices; i++) {Point edge = (Point){poly1.vertices[(i+1)%poly1.num_vertices].x - poly1.vertices[i].x,poly1.vertices[(i+1)%poly1.num_vertices].y - poly1.vertices[i].y};// 计算法向量Point normal = {-edge.y, edge.x};// 投影多边形1double min1 = INFINITY, max1 = -INFINITY;for (int j = 0; j < poly1.num_vertices; j++) {double projection = dot_product(normal, poly1.vertices[j]);min1 = fmin(min1, projection);max1 = fmax(max1, projection);}// 投影多边形2double min2 = INFINITY, max2 = -INFINITY;for (int j = 0; j < poly2.num_vertices; j++) {double projection = dot_product(normal, poly2.vertices[j]);min2 = fmin(min2, projection);max2 = fmax(max2, projection);}// 检查分离轴if (max1 < min2 || max2 < min1) {return false;}}return true;
}
18.7.2 视线计算
// 计算两点间视线是否被阻挡
bool line_of_sight(Point start, Point end, Segment *obstacles, int num_obstacles) {for (int i = 0; i < num_obstacles; i++) {if (segments_intersect(start, end, obstacles[i].start, obstacles[i].end)) {return false;}}return true;
}
18.8 总结与下一章预告
计算几何展示了算法与数学的完美融合：
- 凸包算法：空间边界的高效计算
- 最近点对：分治策略的经典应用
- 线段相交：几何关系的精确判定
- 空间剖分：Voronoi图与Delaunay三角剖分
下一章预告：并行算法
第十九章我们将探索并行算法的世界：
- 并行计算模型：PRAM、BSP、MapReduce
- 并行算法设计模式：分治、流水线、MapReduce
- 经典并行算法：并行排序、并行搜索
- 并行图算法：BFS、最短路径
- GPU并行计算：CUDA/OpenCL基础
随着多核处理器和分布式系统的发展，并行算法已成为解决大规模问题的关键。下一章将揭示如何设计高效并行算法，充分利用现代计算资源。
计算几何不仅是计算机科学的基石，更是连接数学与应用的桥梁。从游戏中的碰撞检测到地图导航的最短路径，从机器人运动规划到有限元分析，几何算法以优雅的数学和高效的实现塑造着我们的数字世界。