道格拉斯-普克算法 - 把一堆復雜的線條變得簡單，同時盡量保持原來的樣子

flyfish

道格拉斯-普克算法（Douglas-Peucker Algorithm解決的問題其實很日常：把一堆復雜的線條（比如地圖上的道路、河流，或者GPS記錄的軌跡）變得簡單，同時盡量保持原來的樣子。

舉個例子，假設用GPS記錄了一條徒步路線，每走1米就記一個點，最后生成了1000個點的折線。但其實很多相鄰的點幾乎在一條直線上，完全沒必要都保留——存起來占空間，畫出來也累贅。這時候這個算法就派上用場了：它能自動刪掉那些“多余”的點，比如直線段中間的點，只留下關鍵的拐點，讓線條變簡單，但看起來還是你走的那條路。

大概是上世紀70年代初。1972年有個叫烏爾斯·拉默的人先提出了類似思路，1973年道格拉斯和普克兩個人又完善了這個方法，所以后來就用他們的名字命名了。

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import io# 設置中文字體，確保中文正常顯示
plt.rcParams["font.family"] = ["SimHei", "sans-serif"]
plt.rcParams['axes.unicode_minus'] = False  # 解決負號顯示問題def point_to_segment_dist(point, start, end):"""計算點到線段的垂直距離"""if np.allclose(start, end):return np.linalg.norm(point - start)# 計算線段的單位向量line_vec = end - startline_len = np.linalg.norm(line_vec)unit_line_vec = line_vec / line_len# 計算從起點到點的向量point_vec = point - start# 計算投影長度projection_length = np.dot(point_vec, unit_line_vec)# 如果投影長度超出線段范圍，則計算到端點的距離if projection_length < 0:return np.linalg.norm(point - start)elif projection_length > line_len:return np.linalg.norm(point - end)# 計算投影點projection = start + projection_length * unit_line_vec# 計算點到投影點的距離return np.linalg.norm(point - projection)def douglas_peucker(points, epsilon, frames=None):"""道格拉斯-普克算法實現，并記錄每一步的處理過程用于可視化參數:points: 待簡化的點集epsilon: 距離閾值frames: 存儲每一步處理結果的列表"""if len(points) < 3:if frames is not None:frames.append({'points': points.copy(),'start_idx': 0,'end_idx': len(points) - 1,'max_dist_idx': None,'is_terminal': True})return pointsstart_point = points[0]end_point = points[-1]# 計算所有中間點到線段的距離distances = []for i in range(1, len(points) - 1):dist = point_to_segment_dist(points[i], start_point, end_point)distances.append((dist, i))if not distances:if frames is not None:frames.append({'points': points.copy(),'start_idx': 0,'end_idx': len(points) - 1,'max_dist_idx': None,'is_terminal': True})return np.array([start_point, end_point])# 找到最大距離的點max_dist, max_dist_idx = max(distances, key=lambda x: x[0])# 記錄當前步驟if frames is not None:frames.append({'points': points.copy(),'start_idx': 0,'end_idx': len(points) - 1,'max_dist_idx': max_dist_idx if max_dist > epsilon else None,'is_terminal': max_dist <= epsilon})# 如果最大距離大于閾值，則遞歸處理if max_dist > epsilon:left_points = points[:max_dist_idx + 1]right_points = points[max_dist_idx:]left_simplified = douglas_peucker(left_points, epsilon, frames)right_simplified = douglas_peucker(right_points, epsilon, frames)# 合并結果（去掉重復的分割點）return np.vstack([left_simplified[:-1], right_simplified])else:# 所有點都足夠接近線段，直接返回起點和終點return np.array([start_point, end_point])def create_frames(points, epsilon):"""創建算法執行過程的幀列表"""frames = []douglas_peucker(points, epsilon, frames)return framesdef draw_frame(frame, frames, points, epsilon, step_num, total_steps, figsize=(10, 6)):"""繪制單個幀"""fig, (ax1, ax2) = plt.subplots(1, 2, figsize=figsize)# 左側子圖：當前處理的線段ax1.set_title("當前處理的線段")ax1.set_xlabel("X")ax1.set_ylabel("Y")ax1.set_xlim(points[:, 0].min() - 1, points[:, 0].max() + 1)ax1.set_ylim(points[:, 1].min() - 1, points[:, 1].max() + 1)# 繪制原始曲線（半透明）ax1.plot(points[:, 0], points[:, 1], 'b-', alpha=0.3, label='原始曲線')# 繪制當前處理的線段current_points = frame['points']start_idx = frame['start_idx']end_idx = frame['end_idx']start_point = current_points[start_idx]end_point = current_points[end_idx]ax1.plot([start_point[0], end_point[0]], [start_point[1], end_point[1]], 'r-', linewidth=2, label='當前線段')# 繪制最大距離點if frame['max_dist_idx'] is not None:max_dist_idx = frame['max_dist_idx']max_point = current_points[max_dist_idx]# 繪制最大距離點ax1.plot(max_point[0], max_point[1], 'go', markersize=8, label='最遠點')# 計算并繪制垂直線projection = compute_projection(max_point, start_point, end_point)ax1.plot([max_point[0], projection[0]], [max_point[1], projection[1]], 'g--', linewidth=1, label='垂直距離')# 顯示距離值dist = point_to_segment_dist(max_point, start_point, end_point)mid_x = (max_point[0] + projection[0]) / 2mid_y = (max_point[1] + projection[1]) / 2ax1.text(mid_x, mid_y, f'd={dist:.2f}', ha='center', va='bottom', bbox=dict(facecolor='white', alpha=0.8))ax1.legend()# 右側子圖：累積簡化結果ax2.set_title("累積簡化結果")ax2.set_xlabel("X")ax2.set_ylabel("Y")ax2.set_xlim(points[:, 0].min() - 1, points[:, 0].max() + 1)ax2.set_ylim(points[:, 1].min() - 1, points[:, 1].max() + 1)# 繪制原始曲線（半透明）ax2.plot(points[:, 0], points[:, 1], 'b-', alpha=0.3, label='原始曲線')# 收集所有已處理的線段simplified_points = []# 從第一步到當前步，收集所有終端節點（即已處理完的線段）for f in frames[:step_num + 1]:if f['is_terminal']:p = f['points']simplified_points.append(p[0])simplified_points.append(p[-1])# 去重并排序（按X坐標）if simplified_points:simplified_points = np.array(simplified_points)_, idx = np.unique(simplified_points[:, 0], return_index=True)simplified_points = simplified_points[np.sort(idx)]# 繪制簡化后的曲線ax2.plot(simplified_points[:, 0], simplified_points[:, 1], 'r-', linewidth=2, label='簡化曲線')ax2.plot(simplified_points[:, 0], simplified_points[:, 1], 'ro', markersize=5)ax2.legend()# 添加標題和步驟信息if frame['is_terminal']:step_text = f"步驟 {step_num + 1}/{total_steps}: 所有點距離 ≤ ε，保留首尾點"else:step_text = f"步驟 {step_num + 1}/{total_steps}: 保留最遠點，分割曲線"fig.suptitle(f"道格拉斯-普克算法 (閾值 ε = {epsilon})", fontsize=14)plt.figtext(0.5, 0.01, step_text, ha="center", fontsize=12)# 保存當前幀為圖像buf = io.BytesIO()plt.savefig(buf, format='png', bbox_inches='tight')buf.seek(0)image = Image.open(buf)plt.close(fig)return imagedef compute_projection(point, start, end):"""計算點在直線上的投影"""if np.allclose(start, end):return startline_vec = end - startpoint_vec = point - startline_len_sq = np.sum(line_vec ** 2)# 計算投影系數t = np.dot(point_vec, line_vec) / line_len_sq# 限制投影在端點之間t = max(0, min(1, t))return start + t * line_vecdef create_gif(points, epsilon, output_path='douglas_peucker.gif', duration=1000):"""創建道格拉斯-普克算法執行過程的GIF動畫"""# 生成所有幀frames = create_frames(points, epsilon)# 繪制每一幀并保存為GIFimages = []for i, frame in enumerate(frames):image = draw_frame(frame, frames, points, epsilon, i, len(frames))images.append(image)# 保存為GIFimages[0].save(output_path,save_all=True,append_images=images[1:],duration=duration,loop=0  # 0表示無限循環)print(f"GIF動畫已保存至: {output_path}")return output_path# 示例：創建一個帶噪聲的正弦曲線并生成GIF
if __name__ == "__main__":# 生成示例數據np.random.seed(42)  # 設置隨機種子，確保結果可重現x = np.linspace(0, 10, 100)y = np.sin(x) + np.random.normal(0, 0.3, size=len(x))  # 添加隨機噪聲points = np.column_stack([x, y])# 設置閾值epsilon = 0.5# 創建GIF動畫create_gif(points, epsilon, output_path='douglas_peucker.gif', duration=1000)