diff --git a/gps_selector.py b/gps_selector.py index 66d3e5c..dbf2dc7 100644 --- a/gps_selector.py +++ b/gps_selector.py @@ -231,6 +231,6 @@ if __name__ == "__main__": # 使用示例 selector = GPSSelector( image_dir=r"G:\error_data\20240930091614\project\images", - output_dir=r"C:\datasets\ODM_output\error1_L" + output_dir=r"G:\test_data\error1_ao" ) selector.run() diff --git a/grid.py b/grid.py new file mode 100644 index 0000000..2494c32 --- /dev/null +++ b/grid.py @@ -0,0 +1,297 @@ +import csv +import numpy as np +import matplotlib.pyplot as plt +import math + +from shapely.geometry import box, MultiPoint +from shapely.ops import unary_union +from scipy.spatial import cKDTree +from utils.gps_extractor import GPSExtractor + + +# ---------------------- overlap 截断为不超过 10% ---------------------- +def clamp_overlap(overlap): + if overlap < 0: + return 0.0 + elif overlap > 0.1: + return 0.1 + else: + return overlap + + +# ====================== 1) 生成可用矩形并记录其覆盖点集 ====================== + +def generate_rectangles_with_point_indices(points, w, h, overlap=0.1, min_points=800): + """ + 在 bounding box 内，以 (w, h) + overlap 布置网格，生成所有矩形。 + 过滤：只保留"矩形内点数 >= min_points"的矩形。 + + 返回: + rect_info_list: list of (rect_polygon, covered_indices) + - rect_polygon: Shapely Polygon + - covered_indices: 一个 set，表示该矩形覆盖的所有点索引 + """ + overlap = clamp_overlap(overlap) + if len(points) == 0: + return [] + + minx, miny = np.min(points, axis=0) + maxx, maxy = np.max(points, axis=0) + + # 特殊情况：只有一个点或非常小范围 -> 很难满足 800 点 + if abs(maxx - minx) < 1e-15 and abs(maxy - miny) < 1e-15: + return [] + + # 步长 + step_x = w * (1 - overlap) + step_y = h * (1 - overlap) + + x_coords = np.arange(minx, maxx + step_x, step_x) + y_coords = np.arange(miny, maxy + step_y, step_y) + + # 建立 KDTree，加速查找 + tree = cKDTree(points) + rect_info_list = [] + + for x in x_coords: + for y in y_coords: + rect_poly = box(x, y, x + w, y + h) + rx_min, ry_min, rx_max, ry_max = rect_poly.bounds + + cx = (rx_min + rx_max) / 2 + cy = (ry_min + ry_max) / 2 + r = math.sqrt((rx_max - rx_min) ** 2 + (ry_max - ry_min) ** 2) / 2 + + candidate_ids = tree.query_ball_point([cx, cy], r) + if not candidate_ids: + continue + + covered_set = set() + for idx_pt in candidate_ids: + px, py = points[idx_pt] + if rx_min <= px <= rx_max and ry_min <= py <= ry_max: + covered_set.add(idx_pt) + + # 如果覆盖的点数不足 min_points，就不保留 + if len(covered_set) < min_points: + continue + + rect_info_list.append((rect_poly, covered_set)) + + return rect_info_list + + +# ====================== 2) 贪心算法选取子集覆盖所有点 ====================== +def cover_all_points_greedy(points, rect_info_list): + """ + 给定所有点 points 以及 "可用矩形+覆盖点集合" rect_info_list， + 要求: + - 选出若干矩形，使得所有点都被覆盖 (每个点至少属于1个选中矩形) + - 最终并集面积最小 (做近似贪心) + + 返回: + chosen_rects: 最终选出的矩形列表 (每个是 shapely Polygon) + """ + + n = len(points) + all_indices = set(range(n)) # 所有点的索引 + uncovered = set(all_indices) # 尚未被覆盖的点索引 + + chosen_rects = [] + union_polygon = None # 当前已选矩形的并集 + + # 如果没有任何矩形可用，就直接失败 + if not rect_info_list: + return [] + + # 为了在贪心过程中快速评估"新矩形带来的额外并集面积"， + # 我们每次选择矩形后更新 union_polygon，然后比较 union_polygon.union(new_rect).area - union_polygon.area + # 但 union_polygon 初始为 None + + while uncovered: + best_gain = 0 + best_new_area = float('inf') + best_rect = None + best_covered_new = set() + + for rect_poly, covered_set in rect_info_list: + # 计算能覆盖多少"尚未覆盖"的点 + newly_covered = uncovered.intersection(covered_set) + if not newly_covered: + continue + + # 计算额外增加的并集面积 + if union_polygon is None: + # 第一次选，union_polygon 为空 => new_area = rect_poly.area + new_area = rect_poly.area + area_increase = new_area + else: + # 计算 union_polygon ∪ rect_poly 的面积 + test_union = union_polygon.union(rect_poly) + new_area = test_union.area + area_increase = new_area - union_polygon.area + + # 贪心策略：最大化 (覆盖点数量) / (面积增量) + # 或者 equivalently, (覆盖点数量) 多、(面积增量) 小 都是好 + ratio = len(newly_covered) / max(area_increase, 1e-12) + + # 我们要找到 ratio 最大的那个 + if ratio > best_gain: + best_gain = ratio + best_new_area = area_increase + best_rect = rect_poly + best_covered_new = newly_covered + + if best_rect is None: + # 没有可选的矩形能覆盖任何剩余点 => 失败 (无法覆盖所有点) + return [] + + # 选中 best_rect + chosen_rects.append(best_rect) + uncovered -= best_covered_new + + # 更新并集 + if union_polygon is None: + union_polygon = best_rect + else: + union_polygon = union_polygon.union(best_rect) + + return chosen_rects + + +# ====================== 3) 主流程: 离散搜索 (w,h) + 贪心覆盖 ====================== +def find_optimal_rectangles_cover_all_points( + points, + base_w, + base_h, + overlap=0.1, + steps=5, + min_points=800 +): + """ + 在 [0.5*base_w,1.5*base_w] x [0.5*base_h,1.5*base_h] 的离散区间枚举 (w,h)， + - 生成可用矩形(≥800 点)的列表 + - 用贪心算法选出子集来覆盖所有点 + - 计算选中矩形的并集面积 + 选出面积最小的方案并返回 + """ + overlap = clamp_overlap(overlap) + n = len(points) + if n == 0: + return [], (base_w, base_h), 0.0 # 没有点就不用覆盖了 + + w_candidates = np.linspace(0.3 * base_w, 2 * base_w, steps) + h_candidates = np.linspace(0.3 * base_h, 2 * base_h, steps) + + best_rects = [] + best_area = float('inf') + best_w, best_h = base_w, base_h + + for w in w_candidates: + for h in h_candidates: + rect_info_list = generate_rectangles_with_point_indices(points, w, h, overlap, min_points) + if not rect_info_list: + # 说明没有任何矩形能达到 "≥800点" + continue + + # 用贪心覆盖所有点 + chosen_rects = cover_all_points_greedy(points, rect_info_list) + if not chosen_rects: + # 无法覆盖所有点 + continue + + # 计算并集面积 + union_poly = unary_union(chosen_rects) + area_covered = union_poly.area + + if area_covered < best_area: + best_area = area_covered + best_rects = chosen_rects + best_w, best_h = w, h + + return best_rects, (best_w, best_h), best_area + + +# ====================== 4) 读取 CSV + 可视化 ====================== +def plot_image_points_cover_all_min_area( + image_dir, # 新参数：图片文件夹路径 + base_rect_width=0.001, + base_rect_height=0.001, + overlap=0.1, + steps=5, + min_points=800 +): + """ + 从图片文件夹读取GPS坐标: + 1) 使用 GPSExtractor 从图片中提取GPS坐标 + 2) 在 [0.5*base_w,1.5*base_w] x [0.5*base_h,1.5*base_h] 离散搜索 (w,h) + 3) 对每个 (w,h), 先生成所有"含≥800点"的矩形 => 再用贪心覆盖所有点 => 计算并集面积 + 4) 最小并集面积的方案即近似最优解 + 5) 最终用该方案可视化 + """ + overlap = clamp_overlap(overlap) + + # 使用 GPSExtractor 读取图片GPS坐标 + extractor = GPSExtractor(image_dir) + gps_df = extractor.extract_all_gps() + + if gps_df.empty: + print("未能从图片中提取到GPS坐标。") + return + + points = np.column_stack((gps_df['lon'], gps_df['lat'])) # (N, 2), [x=lon, y=lat] + n = len(points) + if n == 0: + print("No points extracted from images.") + return + + # 贪心 + 离散搜索 + chosen_rects, (best_w, best_h), best_area = find_optimal_rectangles_cover_all_points( + points, + base_w=base_rect_width, + base_h=base_rect_height, + overlap=overlap, + steps=steps, + min_points=min_points + ) + + if not chosen_rects: + print(f"无法找到满足 '每个矩形≥{min_points}点' 且覆盖所有点 的方案，试着调大尺寸/步数/overlap。") + return + + # 可视化 + plt.figure(figsize=(10, 8)) + # 画点 + plt.scatter(points[:, 0], points[:, 1], c='red', s=10, label='Points') + + # 画矩形 + for i, rect in enumerate(chosen_rects): + if rect.is_empty: + continue + x, y = rect.exterior.xy + plt.fill(x, y, edgecolor='green', fill=False, alpha=0.3, + label='Chosen Rectangles' if i == 0 else "") + + plt.title( + f"Cover All Points, Each Rect≥{min_points} pts, Minimal Union Area\n" + f"base=({base_rect_width:.6f} x {base_rect_height:.6f}), overlap≤{overlap}, steps={steps}\n" + f"best (w,h)=({best_w:.6f},{best_h:.6f}), union area={best_area:.6f}, #rect={len(chosen_rects)}" + ) + plt.xlabel("Longitude") + plt.ylabel("Latitude") + plt.legend() + plt.grid(True) + plt.show() + + +# ------------------ 测试入口 ------------------ +if __name__ == "__main__": + image_dir = r"C:\datasets\134\code\images" # 替换为你的图片文件夹路径 + plot_image_points_cover_all_min_area( + image_dir, + base_rect_width=0.01, + base_rect_height=0.01, + overlap=0.05, # 会被截断到 0.1 + steps=40, + min_points=100 + )