代码存档_古地图电子化_平面二维点集之间的局部普氏重叠度计算及可视化

点击查看代码

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import LineCollection
import matplotlib.font_manager as fm

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'SimSun', 'KaiTi', 'FangSong']
plt.rcParams['axes.unicode_minus'] = False

def calculate_procrustes_distance(points1, points2):
    """计算两个点集之间的普氏距离"""
    # 将点集转换为numpy数组
    points1 = np.array(points1)
    points2 = np.array(points2)
    
    # 计算平移
    centroid1 = np.mean(points1, axis=0)
    centroid2 = np.mean(points2, axis=0)
    points1_centered = points1 - centroid1
    points2_centered = points2 - centroid2
    
    # 计算缩放
    scale = np.sqrt(np.sum(points1_centered**2) / np.sum(points2_centered**2))
    points2_scaled = points2_centered * scale
    
    # 计算旋转
    H = points1_centered.T @ points2_scaled
    U, S, Vt = np.linalg.svd(H)
    R = Vt.T @ U.T
    
    # 确保R是正确的维度
    if points1.shape[1] == 2:
        # 2D情况
        R = np.array([[1, 0], [0, 1]]) if np.linalg.det(R) < 0 else R
    
    # 计算普氏距离
    points2_rotated = points2_scaled @ R.T
    d = np.sqrt(np.sum((points1_centered - points2_rotated)**2))
    return d

def partial_procrustes_overlap(points1, points2):
    """计算两个点集之间的部分普氏重叠"""
    # 将点集转换为numpy数组
    points1 = np.array(points1)
    points2 = np.array(points2)
    
    # 计算普氏距离
    d = calculate_procrustes_distance(points1, points2)
    
    # 计算部分重叠
    # 这里假设部分重叠是基于某个阈值来定义的
    threshold = 0.1
    overlap = 1 - (d / np.sqrt(np.sum(points1**2) + np.sum(points2**2)))
    
    # 如果重叠小于阈值，则认为没有部分重叠
    if overlap < threshold:
        return 0
    else:
        return overlap

def visualize_points(points1, points2, title="点集可视化"):
    """可视化两个点集"""
    points1 = np.array(points1)
    points2 = np.array(points2)
    
    plt.figure(figsize=(12, 5))
    
    # 子图1：原始点集
    plt.subplot(1, 2, 1)
    plt.scatter(points1[:, 0], points1[:, 1], color='red', label='点集1', s=100)
    plt.scatter(points2[:, 0], points2[:, 1], color='blue', label='点集2', s=100)
    
    # 连接点形成多边形
    if len(points1) > 2:
        plt.plot(np.append(points1[:, 0], points1[0, 0]), 
                np.append(points1[:, 1], points1[0, 1]), 'r--', alpha=0.5)
    if len(points2) > 2:
        plt.plot(np.append(points2[:, 0], points2[0, 0]), 
                np.append(points2[:, 1], points2[0, 1]), 'b--', alpha=0.5)
    
    plt.title("原始点集")
    plt.xlabel("X")
    plt.ylabel("Y")
    # 计算并显示欧式距离
    if len(points1) == len(points2):
        euclidean_distance = np.sqrt(np.sum((points1 - points2) ** 2))
        plt.text(0.02, 0.98, f'欧式距离: {euclidean_distance:.3f}', 
                transform=plt.gca().transAxes, fontsize=10, 
                verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    plt.legend(['点集1', '点集2'])
    plt.grid(True)
    
    # 子图2：重叠后的点集
    plt.subplot(1, 2, 2)
    
    # 计算重叠变换
    centroid1 = np.mean(points1, axis=0)
    centroid2 = np.mean(points2, axis=0)
    points1_centered = points1 - centroid1
    points2_centered = points2 - centroid2
    
    scale = np.sqrt(np.sum(points1_centered**2) / np.sum(points2_centered**2))
    points2_scaled = points2_centered * scale
    
    H = points1_centered.T @ points2_scaled
    U, S, Vt = np.linalg.svd(H)
    R = Vt.T @ U.T
    
    # 确保R是正确的维度
    if points1.shape[1] == 2:
        R = np.array([[1, 0], [0, 1]]) if np.linalg.det(R) < 0 else R
    
    points2_transformed = (points2_scaled @ R.T) + centroid1
    
    plt.scatter(points1[:, 0], points1[:, 1], color='red', label='点集1 (参考)', s=100)
    plt.scatter(points2_transformed[:, 0], points2_transformed[:, 1], 
               color='green', label='点集2 (重叠后)', s=100, alpha=0.7)
    
    # 连接点形成多边形
    if len(points1) > 2:
        plt.plot(np.append(points1[:, 0], points1[0, 0]), 
                np.append(points1[:, 1], points1[0, 1]), 'r--', alpha=0.5)
    if len(points2_transformed) > 2:
        plt.plot(np.append(points2_transformed[:, 0], points2_transformed[0, 0]), 
                np.append(points2_transformed[:, 1], points2_transformed[0, 1]), 'g--', alpha=0.5)
    
    plt.title("重叠后点集")
    plt.xlabel("X")
    plt.ylabel("Y")
    # 计算并显示普氏距离
    if len(points1) == len(points2):
        procrustes_distance = calculate_procrustes_distance(points1, points2)
        plt.text(0.02, 0.98, f'普氏距离: {procrustes_distance:.3f}', 
                transform=plt.gca().transAxes, fontsize=10, 
                verticalalignment='top', bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.5))
    
    plt.legend(['点集1 (参考)', '点集2 (重叠后)'])
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()

def visualize_error_vector_field(points1, points2, scale_factor=1.0):
    """可视化全局误差向量场
    
    参数:
    points1: 参考点集 (目标点集)
    points2: 待对齐点集 (源点集)
    scale_factor: 箭头缩放因子，默认1.0
    """
    points1 = np.array(points1)
    points2 = np.array(points2)
    
    # 计算普氏变换参数
    centroid1 = np.mean(points1, axis=0)
    centroid2 = np.mean(points2, axis=0)
    points1_centered = points1 - centroid1
    points2_centered = points2 - centroid2
    
    scale = np.sqrt(np.sum(points1_centered**2) / np.sum(points2_centered**2))
    points2_scaled = points2_centered * scale
    
    H = points1_centered.T @ points2_scaled
    U, S, Vt = np.linalg.svd(H)
    R = Vt.T @ U.T
    
    if points1.shape[1] == 2:
        R = np.array([[1, 0], [0, 1]]) if np.linalg.det(R) < 0 else R
    
    # 计算变换后的点集
    points2_transformed = (points2_scaled @ R.T) + centroid1
    
    # 计算残差向量 e_i = (R @ X_i + t) - Y_i
    # 这里Y_i是参考点，X_i是源点
    residual_vectors = points2_transformed - points1
    
    # 计算残差大小
    residual_magnitudes = np.sqrt(np.sum(residual_vectors**2, axis=1))
    max_magnitude = np.max(residual_magnitudes)
    
    # 创建图形
    plt.figure(figsize=(16, 10))
    
    # 绘制变换后的点集
    plt.scatter(points1[:, 0], points1[:, 1], color='red', label='参考点集 (Y)', s=120, zorder=5)
    plt.scatter(points2_transformed[:, 0], points2_transformed[:, 1], 
               color='green', label='对齐后点集 (R·X+t)', s=120, alpha=0.7, zorder=4)
    
    # 绘制残差向量场 - 使用更鲜艳的颜色方案
    colors = plt.cm.plasma(residual_magnitudes / max_magnitude if max_magnitude > 0 else np.zeros_like(residual_magnitudes))
    
    # 设置最小透明度确保可见性
    alphas = np.clip(0.7 + 0.3 * (residual_magnitudes / max_magnitude if max_magnitude > 0 else np.zeros_like(residual_magnitudes)), 0.7, 1.0)
    
    for i, (start, vec, mag) in enumerate(zip(points1, residual_vectors, residual_magnitudes)):
        if mag > 0:  # 只绘制非零向量
            plt.arrow(start[0], start[1], 
                     vec[0] * scale_factor, vec[1] * scale_factor,
                     head_width=0.06 * max_magnitude * scale_factor, 
                     head_length=0.12 * max_magnitude * scale_factor,
                     fc=colors[i], ec='black', alpha=alphas[i], width=0.015, linewidth=0.5)
    
    # 连接点形成多边形
    if len(points1) > 2:
        plt.plot(np.append(points1[:, 0], points1[0, 0]), 
                np.append(points1[:, 1], points1[0, 1]), 'r--', alpha=0.3, linewidth=2)
        plt.plot(np.append(points2_transformed[:, 0], points2_transformed[0, 0]), 
                np.append(points2_transformed[:, 1], points2_transformed[0, 1]), 'g--', alpha=0.3, linewidth=2)
    
    # 添加颜色条
    if max_magnitude > 0:
        sm = plt.cm.ScalarMappable(cmap=plt.cm.plasma, 
                                 norm=plt.Normalize(vmin=0, vmax=max_magnitude))
        sm.set_array([])
        cbar = plt.colorbar(sm, ax=plt.gca(), label='残差大小', shrink=0.8)
        cbar.set_label('残差大小', fontsize=12)
    
    # 添加统计信息
    plt.title('全局误差向量场可视化\n残差向量: e_i = (R·X_i + t) - Y_i', fontsize=14, pad=20)
    plt.xlabel('X 坐标', fontsize=12)
    plt.ylabel('Y 坐标', fontsize=12)
    
    # 添加文本信息
    avg_residual = np.mean(residual_magnitudes)
    max_residual = np.max(residual_magnitudes)
    info_text = f'平均残差: {avg_residual:.3f}\n最大残差: {max_residual:.3f}\n箭头缩放: {scale_factor}x'
    plt.text(0.02, 0.98, info_text, transform=plt.gca().transAxes, fontsize=10,
             verticalalignment='top', bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))
    
    plt.legend(loc='upper right')
    plt.grid(True, alpha=0.3)
    plt.axis('equal')
    
    plt.tight_layout()
    plt.show()
    
    return residual_vectors, residual_magnitudes

def visualize_comparison(points1, points2):
    """可视化比较原始点集和重叠后的点集"""
    # 计算距离和重叠度
    distance = calculate_procrustes_distance(points1, points2)
    overlap = partial_procrustes_overlap(points1, points2)
    
    # 转换为numpy数组
    points1_array = np.array(points1)
    points2_array = np.array(points2)
    
    # 计算原始欧式距离
    original_euclidean = np.sqrt(np.sum((points1_array - points2_array) ** 2))
    
    # 计算变换后的点集用于可视化
    centroid1 = np.mean(points1_array, axis=0)
    centroid2 = np.mean(points2_array, axis=0)
    points1_centered = points1_array - centroid1
    points2_centered = points2_array - centroid2
    
    scale = np.sqrt(np.sum(points1_centered**2) / np.sum(points2_centered**2))
    points2_scaled = points2_centered * scale
    
    H = points1_centered.T @ points2_scaled
    U, S, Vt = np.linalg.svd(H)
    R = Vt.T @ U.T
    
    if points1_array.shape[1] == 2:
        R = np.array([[1, 0], [0, 1]]) if np.linalg.det(R) < 0 else R
    
    points2_transformed = (points2_scaled @ R.T) + centroid1
    
    print(f"普氏距离: {distance:.4f}")
    print(f"部分普氏重叠度: {overlap:.9f}")
    print(f"原始欧式距离: {original_euclidean:.4f}")
    
    # 创建比较图形
    plt.figure(figsize=(15, 6))
    
    # 原始点集
    plt.subplot(1, 2, 1)
    plt.scatter(points1_array[:, 0], points1_array[:, 1], color='red', label='点集1', s=100)
    plt.scatter(points2_array[:, 0], points2_array[:, 1], color='blue', label='点集2', s=100)
    
    # 连接点形成多边形
    if len(points1_array) > 2:
        plt.plot(np.append(points1_array[:, 0], points1_array[0, 0]), 
                np.append(points1_array[:, 1], points1_array[0, 1]), 'r--', alpha=0.5)
    if len(points2_array) > 2:
        plt.plot(np.append(points2_array[:, 0], points2_array[0, 0]), 
                np.append(points2_array[:, 1], points2_array[0, 1]), 'b--', alpha=0.5)
    
    plt.title(f"原始点集\n欧式距离: {original_euclidean:.3f}")
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.legend(['点集1', '点集2'])
    plt.grid(True)
    
    # 重叠后点集
    plt.subplot(1, 2, 2)
    plt.scatter(points1_array[:, 0], points1_array[:, 1], color='red', label='点集1 (参考)', s=100)
    plt.scatter(points2_transformed[:, 0], points2_transformed[:, 1], 
               color='green', label='点集2 (重叠后)', s=100, alpha=0.7)
    
    # 连接点形成多边形
    if len(points1_array) > 2:
        plt.plot(np.append(points1_array[:, 0], points1_array[0, 0]), 
                np.append(points1_array[:, 1], points1_array[0, 1]), 'r--', alpha=0.5)
    if len(points2_transformed) > 2:
        plt.plot(np.append(points2_transformed[:, 0], points2_transformed[0, 0]), 
                np.append(points2_transformed[:, 1], points2_transformed[0, 1]), 'g--', alpha=0.5)
    
    plt.title(f"重叠后点集\n普氏距离: {distance:.3f}")
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.legend(['点集1 (参考)', '点集2 (重叠后)'])
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    # 示例使用
    points1 = [[2.1, 3.8], [5.9, 7.2], [8.3, 4.1], [1.5, 9.7], [6.4, 2.5], [10.2, 6.8], [3.7, 1.2]]
    points2 = [[2.3, 4.0], [6.1, 7.1], [8.5, 4.3], [1.7, 9.5], [6.2, 2.7], [10.0, 7.0], [3.9, 1.4]]
    
    print("=== 原始比较可视化 ===")
    visualize_comparison(points1, points2)
    
    print("=== 全局误差向量场可视化 ===")
    residual_vectors, residual_magnitudes = visualize_error_vector_field(points1, points2, scale_factor=5)
    
    print("残差向量统计:")
    print(f"平均残差大小: {np.mean(residual_magnitudes):.4f}")
    print(f"最大残差大小: {np.max(residual_magnitudes):.4f}")
    print(f"残差标准差: {np.std(residual_magnitudes):.4f}")
    
    # 也可以尝试不同的缩放因子
    #print("\n=== 误差向量场（箭头缩放2倍）===")
    #visualize_error_vector_field(points1, points2, scale_factor=2.0)