scipy.spatialのVoronoiにおいて無限遠点を含む領域(ポリゴン)をどうにか描画領域内だけでも覆うようにして着色したい．

voronoi.py

# %%
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PolyCollection
from sklearn.decomposition import PCA
from scipy.spatial import Voronoi
from scipy.spatial import ConvexHull
import colorsys


def compare_reconstruction(data):
    """Compare reconstruction errors of SVD and EVD for increasing components."""
    data_mean = np.mean(data, axis=0)
    data_centered = data - data_mean

    # SVD
    U, singular_values, Vt = np.linalg.svd(data_centered, full_matrices=False)
    svd_errors = []
    for k in range(1, 8):
        approx_data_svd = U[:, :k] @ np.diag(singular_values[:k]) @ Vt[:k, :]
        error_svd = np.linalg.norm(data_centered - approx_data_svd, 'fro')
        svd_errors.append(error_svd)

    # EVD (X.T @ X)
    covariance_matrix = (1 / (data.shape[0] - 1)) * data_centered.T @ data_centered
    eigenvalues, eigenvectors = np.linalg.eigh(covariance_matrix)
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]
    evd_errors = []
    for k in range(1, 8):
        V_k = eigenvectors[:, :k]
        approx_data_evd = data_centered @ V_k @ V_k.T
        error_evd = np.linalg.norm(data_centered - approx_data_evd, 'fro')
        evd_errors.append(error_evd)

    # Plot results
    plt.plot(range(1, 8), svd_errors, label='SVD')
    plt.plot(range(1, 8), evd_errors, label='EVD (X.T @ X)')
    plt.xlabel('Number of components (k)')
    plt.ylabel('Reconstruction error (Frobenius norm)')
    plt.title('SVD vs EVD Reconstruction Error')
    plt.legend()
    plt.grid()
    plt.show()

    print("SVD errors:", svd_errors)
    print("EVD errors:", evd_errors)


# def calculate_polygon_area(polygon):
#     """Calculate the area of a polygon using the shoelace formula."""
#     num_vertices = len(polygon)
#     area = 0.5 * sum(
#         (p1[0] + p0[0]) * (p1[1] - p0[1])
#         for p0, p1 in [[polygon[i % num_vertices], polygon[(i + 1) % num_vertices]] for i in range(num_vertices)]
#     )
#     return abs(area)


def calculate_cross_product(point1, point2, point3):
    """Calculate the cross product of vectors point1->point2 and point2->point3."""
    vec1 = np.array(point2) - np.array(point1)
    vec2 = np.array(point3) - np.array(point2)
    return vec1[0] * vec2[1] - vec1[1] * vec2[0]


def find_opposite_bend_index(points):
    """Return the index of a point where the polygon bends oppositely."""
    num_points = len(points)
    prev_cross_product = None
    for i in range(num_points):
        p1 = points[i]
        p2 = points[(i + 1) % num_points]
        p3 = points[(i + 2) % num_points]
        cross_product = calculate_cross_product(p1, p2, p3)
        if prev_cross_product is not None and prev_cross_product * cross_product < 0:
            return (i + 1) % num_points
        prev_cross_product = cross_product
    return -1


def is_polygon_convex(points, verbose=False):
    """Check if the polygon defined by points is convex."""
    opposite_bend_idx = find_opposite_bend_index(points)
    if opposite_bend_idx == -1:
        return True

    num_points = len(points)
    idx = opposite_bend_idx
    p1 = points[idx]
    p2 = points[(idx + 1) % num_points]
    p3 = points[(idx + 2) % num_points]
    cross_product = calculate_cross_product(p1, p2, p3)
    if verbose:
        print(f"Non-convex at points {p1}, {p2}, {p3}: cross_product = {cross_product}")
    return False


def check_segment_intersection(point1, point2, point3, point4, verbose=False):
    """Check if segments point1-point2 and point3-point4 intersect (excluding endpoint contact)."""
    # Exclude intersection if endpoints coincide
    if np.allclose(point1, point3) or np.allclose(point1, point4) or \
       np.allclose(point2, point3) or np.allclose(point2, point4):
        if verbose:
            print(f"Endpoint coincidence detected: {point1}, {point2}, {point3}, {point4}")
        return False

    def is_counter_clockwise(a, b, c):
        cross_product = calculate_cross_product(a, b, c)
        return cross_product > 1e-8  # Adjusted threshold for numerical stability

    # Orientation test using cross products
    d1 = is_counter_clockwise(point1, point2, point3)
    d2 = is_counter_clockwise(point1, point2, point4)
    d3 = is_counter_clockwise(point3, point4, point1)
    d4 = is_counter_clockwise(point3, point4, point2)

    # General intersection condition
    intersects = d1 != d2 and d3 != d4
    if intersects and verbose:
        print(f"Intersection detected: {point1}-{point2} and {point3}-{point4}")
        print(f"ccw values: d1={d1}, d2={d2}, d3={d3}, d4={d4}")

    # Additional check for endpoints lying on segments
    if not intersects:
        def is_point_on_segment(point, seg_start, seg_end):
            px, py = point
            ax, ay = seg_start
            bx, by = seg_end
            return (min(ax, bx) <= px <= max(ax, bx) and
                    min(ay, by) <= py <= max(ay, by) and
                    abs(calculate_cross_product(seg_start, point, seg_end)) < 1e-8)

        if (is_point_on_segment(point1, point3, point4) or
            is_point_on_segment(point2, point3, point4) or
            is_point_on_segment(point3, point1, point2) or
                is_point_on_segment(point4, point1, point2)):
            if verbose:
                print(f"Endpoint on segment detected: {point1}, {point2}, {point3}, {point4}")
            return False

    return intersects


def is_polygon_simple(polygon, verbose=False):
    """Check if a polygon is simple (non-self-intersecting)."""
    num_vertices = len(polygon)
    for i in range(num_vertices):
        for j in range(i + 2, num_vertices):
            # Skip adjacent edges
            if j == (i + 1) % num_vertices or i == (j + 1) % num_vertices:
                continue
            p1 = polygon[i]
            p2 = polygon[(i + 1) % num_vertices]
            p3 = polygon[j]
            p4 = polygon[(j + 1) % num_vertices]
            if check_segment_intersection(p1, p2, p3, p4, verbose):
                if verbose:
                    print(f"Self-intersection detected between edges {p1}-{p2} and {p3}-{p4}")
                return False
    return True


# def calculate_distance_squared(point1, point2):
#     """Calculate the squared Euclidean distance between two points."""
#     return (point1[0] - point2[0])**2 + (point1[1] - point2[1])**2


def calculate_perpendicular_bisector(point1, point2):
    """Calculate the perpendicular bisector of the segment between point1 and point2."""
    x1, y1 = point1
    x2, y2 = point2
    # Midpoint
    mid_x, mid_y = (x1 + x2) / 2, (y1 + y2) / 2
    # Slope of the original line
    try:
        slope_orig = (y2 - y1) / (x2 - x1)
        slope_perp = -1 / slope_orig
    except ZeroDivisionError:
        slope_perp = None  # Vertical line case
    intercept = mid_y - (slope_perp * mid_x if slope_perp is not None else 0)
    return slope_perp, intercept


def lighten_color(rgba):
    """Convert a color to a lighter version by adjusting HSV saturation and value."""
    rgb_normalized = rgba[:3]
    h, s, v = colorsys.rgb_to_hsv(*rgb_normalized)
    s *= 0.2  # Reduce saturation
    v = 0.9   # Set high value for lightness
    return colorsys.hsv_to_rgb(h, s, v)


# def build_polygon_from_region(voronoi, region, intersection_points, order):
#     """Build a polygon by replacing infinite vertices with intersection points."""
#     polygon = []
#     for vertex_idx in region:
#         if vertex_idx != -1:
#             polygon.append(voronoi.vertices[vertex_idx])
#         else:
#             # Add intersection points based on order
#             if order[0] == 0 or len(intersection_points) > 1:
#                 polygon.append(intersection_points[order[0]])
#             if order[1] == 0 or len(intersection_points) > 1:
#                 polygon.append(intersection_points[order[1]])
#     return polygon


def calculate_circle_line_intersections(center, radius, slope, intercept):
    """Calculate intersection points of a line with a circle."""
    mx, my = center
    a = 1 + slope**2
    b = 2 * (slope * (intercept - my) - mx)
    c = mx**2 + (intercept - my)**2 - radius**2
    discriminant = b**2 - 4 * a * c
    if discriminant < 0:
        return None
    x1 = (-b + np.sqrt(discriminant)) / (2 * a)
    x2 = (-b - np.sqrt(discriminant)) / (2 * a)
    y1 = slope * x1 + intercept
    y2 = slope * x2 + intercept
    return np.array([[x1, y1], [x2, y2]])


def calculate_convex_hull(polygons):
    """Calculate the convex hull of combined polygon points."""
    combined_points = np.vstack(polygons)
    unique_points = np.unique(combined_points, axis=0)
    hull = ConvexHull(unique_points)
    return unique_points[hull.vertices]


def set_uniform_grid_ticks(ax):
    """Set uniform grid ticks with nice, rounded intervals."""
    def nice_number(value):
        """Round value to a nice number (e.g., 0.1, 0.2, 0.5, 1, 2, 5, 10)."""
        if value == 0:
            return 1.0
        power = np.floor(np.log10(value))
        normalized = value / (10 ** power)
        if normalized < 1.5:
            nice = 1
        elif normalized < 3:
            nice = 2
        elif normalized < 7:
            nice = 5
        else:
            nice = 10
        return nice * (10 ** power)

    x_ticks = ax.get_xticks()
    y_ticks = ax.get_yticks()
    x_interval = np.diff(x_ticks)[0] if len(x_ticks) > 1 else 1
    y_interval = np.diff(y_ticks)[0] if len(y_ticks) > 1 else 1
    tick_interval = min(x_interval, y_interval)
    tick_interval = nice_number(tick_interval)
    x_lim = ax.get_xlim()
    y_lim = ax.get_ylim()
    x_ticks_new = np.arange(
        np.floor(x_lim[0] / tick_interval) * tick_interval,
        np.ceil(x_lim[1] / tick_interval) * tick_interval + tick_interval,
        tick_interval
    )
    y_ticks_new = np.arange(
        np.floor(y_lim[0] / tick_interval) * tick_interval,
        np.ceil(y_lim[1] / tick_interval) * tick_interval + tick_interval,
        tick_interval
    )
    x_ticks_new = x_ticks_new[(x_ticks_new >= x_lim[0]) & (x_ticks_new <= x_lim[1])]
    y_ticks_new = y_ticks_new[(y_ticks_new >= y_lim[0]) & (y_ticks_new <= y_lim[1])]
    ax.set_xticks(x_ticks_new)
    ax.set_yticks(y_ticks_new)


def prepare_data_IMP():
    FEATURE_MAP = {
        1: [
            np.array([0.75, 0.96, 0.98, 0.75, 0.32, 0.50, 0.57]),
            np.array([0.74193548, 0.93548387, 0.9223301, 0.00970874, 0.30097087, 0.48261823, 0.5483871]),
        ],
        2: [np.array([0.88, 0.86, 0.88, 0.89, 0.59, 0.51, 0.70])],
        3: [
            np.array([0.88, 0.72, 0.92, 0.91, 0.59, 0.52, 0.40]),
            np.array([0.88, 0.72, 0.92, 0.91, 0.59, 0.52, 0.62]),
        ],
        4: [np.array([0.78, 0.96, 0.87, 0.89, 0.64, 0.55, 0.58])],
        5: [np.array([0.90, 0.73, 0.99, 0.93, 0.59, 0.50, 0.52])],
        6: [
            np.array([0.88, 0.71, 0.75, 0.80, 0.57, 0.43, 0.38]),
            np.array([0.88, 0.74, 0.71, 0.83, 0.57, 0.43, 0.38]),
        ],
        7: [np.array([0.87, 0.96, 0.88, 0.86, 0.54, 0.62, 0.75])],
        8: [np.array([0.87, 0.69, 0.86, 0.97, 0.60, 0.45, 0.40])],
        9: [
            np.array([0.86, 0.72, 0.81, 0.75, 0.57, 0.44, 0.30]),
            np.array([0.85, 0.70, 0.80, 0.77, 0.58, 0.45, 0.35]),
        ],
        0: [
            np.array([0.75, 0.84, 0.92, 0.68, 0.64, 0.48, 0.54]),
            np.array([0.75, 0.82, 0.90, 0.70, 0.62, 0.50, 0.45]),
        ],
    }

    feature_data = []
    labels = []
    for digit, templates in FEATURE_MAP.items():
        for template in templates:
            feature_data.append(template)
            labels.append(digit)

    feature_data = np.array(feature_data)
    labels = np.array(labels)

    # PCA with all components
    pca = PCA(n_components=feature_data.shape[1])  # 7 components
    pca.fit(feature_data)

    return pca, feature_data, labels


def prepare_data(use_pca=True, component_indices=[0, 1]):
    """Prepare feature data, perform PCA or SVD, and compute Voronoi diagram with specified components."""
    pca, feature_data, labels = prepare_data_IMP()

    # Select specified components
    selected_components = pca.components_[component_indices]
    projected_data_pca = feature_data @ selected_components.T

    if not use_pca:
        # SVD (manual implementation)
        data_mean = np.mean(feature_data, axis=0)
        data_centered = feature_data - data_mean
        U, singular_values, Vt = np.linalg.svd(data_centered, full_matrices=False)
        principal_components_svd = Vt[component_indices, :]  # Select specified components
        projected_data_svd = data_centered @ principal_components_svd.T
        num_samples = feature_data.shape[0]
        explained_variance_svd = (singular_values[component_indices]**2) / (num_samples - 1)
        explained_variance_ratio_svd = explained_variance_svd / np.sum(singular_values**2 / (num_samples - 1))

        # Compare PCA and SVD
        print("=== PCA vs SVD Comparison ===")
        print("PCA components_:\n", selected_components)
        print("SVD components:\n", principal_components_svd)
        print("PCA explained_variance_:\n", pca.explained_variance_[component_indices])
        print("SVD explained_variance:\n", explained_variance_svd)
        print("PCA explained_variance_ratio_:\n", pca.explained_variance_ratio_[component_indices])
        print("SVD explained_variance_ratio:\n", explained_variance_ratio_svd)
        print("PCA projected_data (first 5 points):\n", projected_data_pca[:5])
        print("SVD projected_data (first 5 points):\n", projected_data_svd[:5])
        print("Data mean (PCA mean_):\n", pca.mean_)
        print("Data mean (SVD computed):\n", data_mean)

        # Adjust signs
        for i in range(2):
            if np.allclose(selected_components[i], -principal_components_svd[i], atol=1e-6):
                principal_components_svd[i] = -principal_components_svd[i]
                projected_data_svd[:, i] = -projected_data_svd[:, i]
                print(f"Flipped sign of component {component_indices[i]+1} for consistency")

        print("Difference in components_:\n", np.abs(selected_components - principal_components_svd).max())
        print("Difference in projected_data:\n", np.abs(projected_data_pca - projected_data_svd).max())

        compare_reconstruction(feature_data)

        voronoi = Voronoi(projected_data_svd)
        return voronoi, pca, projected_data_svd, labels, component_indices
    else:
        voronoi = Voronoi(projected_data_pca)
        return voronoi, pca, projected_data_pca, labels, component_indices


def is_point_in_polygon(points, q, epsilon=1e-10):
    points = np.array(points)
    q = np.array(q)
    min_x, max_x = np.min(points[:, 0]), np.max(points[:, 0])
    min_y, max_y = np.min(points[:, 1]), np.max(points[:, 1])
    if not (min_x <= q[0] <= max_x and min_y <= q[1] <= max_y):
        return False

    def in_y_range(p1, p2, q):
        return (p1[1] > q[1]) != (p2[1] > q[1])

    def is_left_side(p1, p2, q):
        x_inters = (q[1] - p1[1]) * (p2[0] - p1[0]) / (p2[1] - p1[1]) + p1[0]
        return q[0] < x_inters - epsilon

    def is_on_vertex_or_edge(p1, p2, q):
        if np.allclose(q, p1, atol=epsilon) or np.allclose(q, p2, atol=epsilon):
            return True
        if min(p1[1], p2[1]) <= q[1] <= max(p1[1], p2[1]):
            if abs(p2[1] - p1[1]) > epsilon:
                x_inters = (q[1] - p1[1]) * (p2[0] - p1[0]) / (p2[1] - p1[1]) + p1[0]
                if abs(q[0] - x_inters) < epsilon:
                    return True
        return False

    n, inside = len(points), False
    for i in range(n):
        p1, p2 = points[i], points[(i + 1) % n]
        if is_on_vertex_or_edge(p1, p2, q):
            return True
        if abs(p1[1] - p2[1]) < epsilon:
            continue
        if in_y_range(p1, p2, q) and is_left_side(p1, p2, q):
            inside = not inside
    return inside


def closest_polygon_and_vertex(polygons, q, epsilon=1e-10):
    q = np.array(q)
    min_distance = float('inf')
    best_poly_idx = -1
    best_vertex_idx = -1
    best_is_vertex = False

    for poly_idx, points in enumerate(polygons):
        points = np.array(points)
        n = len(points)

        if is_point_in_polygon(points, q, epsilon):
            continue

        vertex_distances = np.linalg.norm(points - q, axis=1)
        min_vertex_dist = np.min(vertex_distances)
        min_vertex_idx = np.argmin(vertex_distances)

        min_edge_dist = float('inf')
        min_edge_idx = -1
        for i in range(n):
            p1 = points[i]
            p2 = points[(i + 1) % n]
            v = p2 - p1
            w = q - p1
            v_dot_v = np.dot(v, v)
            if v_dot_v < epsilon:
                continue
            t = np.dot(w, v) / v_dot_v
            if t < 0:
                dist = np.linalg.norm(q - p1)
                if dist < min_edge_dist:
                    min_edge_dist = dist
                    min_edge_idx = i
            elif t > 1:
                dist = np.linalg.norm(q - p2)
                if dist < min_edge_dist:
                    min_edge_dist = dist
                    min_edge_idx = (i + 1) % n
            else:
                projection = p1 + t * v
                dist = np.linalg.norm(q - projection)
                if dist < min_edge_dist:
                    min_edge_dist = dist
                    min_edge_idx = i

        if min_vertex_dist < min_edge_dist:
            if min_vertex_dist < min_distance:
                min_distance = min_vertex_dist
                best_poly_idx = poly_idx
                best_vertex_idx = min_vertex_idx
                best_is_vertex = True
        else:
            if min_edge_dist < min_distance:
                min_distance = min_edge_dist
                best_poly_idx = poly_idx
                best_vertex_idx = min_edge_idx
                best_is_vertex = False

    # if best_poly_idx == -1:
    #     raise ValueError("No valid polygon found (all points may be inside polygons)")

    return best_poly_idx, best_vertex_idx, best_is_vertex, min_distance


def get_unique_points(points):
    # 点をタプルに変換
    points_tuples = [tuple(p) for p in points]
    # カウント
    counts = {}
    for p in points_tuples:
        counts[p] = counts.get(p, 0) + 1
    # 元の順序で1回出現の点のみ抽出
    return [list(p) for p in points_tuples if counts[p] == 1]


def plot_voronoi(voronoi, pca, projected_data, labels, use_pca=True, component_indices=[0, 1], ax=None):
    """Plot the Voronoi diagram with finite and infinite regions for specified components."""
    is7x7 = True if ax is not None else False

    xlim = [f(projected_data[:, 0]) for f in [min, max]]
    ylim = [f(projected_data[:, 1]) for f in [min, max]]
    x_margin = 0.05 * (xlim[1] - xlim[0])
    y_margin = 0.05 * (ylim[1] - ylim[0])
    xlim = (xlim[0] - x_margin, xlim[1] + x_margin)
    ylim = (ylim[0] - y_margin, ylim[1] + y_margin)

    region_to_point = {voronoi.point_region[point_idx]: point_idx for point_idx in range(len(voronoi.points))}
    color_map = [rgba[:3] for rgba in plt.cm.tab10(np.arange(10))]

    if not is7x7:
        fig, ax = plt.subplots(figsize=(11.69, 8.27), dpi=300)
        ax.set_xlabel(f'PC{component_indices[0]+1}')
        ax.set_ylabel(f'PC{component_indices[1]+1}')
    ax.set_aspect('equal')
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.grid()

    finite_polygons = [voronoi.vertices[region] for region in voronoi.regions if len(region) > 0 and -1 not in region]
    finite_polygons_pidx = [region_to_point[ridx] for ridx, region in enumerate(
        voronoi.regions) if len(region) > 0 and -1 not in region]
    finite_colors = [
        lighten_color(color_map[labels[region_to_point[region_idx]]])
        for region_idx, region in enumerate(voronoi.regions)
        if len(region) > 0 and -1 not in region
    ]
    ax.add_collection(PolyCollection(finite_polygons, facecolors=finite_colors, edgecolors='gray'))

    infinite_ridge_pairs = [k for k, v in voronoi.ridge_dict.items() if -1 in v]
    infinite_region_polygons = {}
    center = np.mean(voronoi.points, axis=0)
    radius0 = max(np.linalg.norm(voronoi.points - center, axis=1))
    radius1 = max(np.linalg.norm(voronoi.vertices - center, axis=1))
    radius = max([radius0, radius1])

    for point_pair in infinite_ridge_pairs:
        # point間の垂直二等分線のradius上の点を求め，pointのpolygonに加える
        points = [voronoi.points[point_idx] for point_idx in point_pair]
        slope, intercept = calculate_perpendicular_bisector(points[0], points[1])
        intersection_points = calculate_circle_line_intersections(center, radius, slope, intercept)
        if intersection_points is None:
            continue
        region0 = voronoi.regions[voronoi.point_region[point_pair[0]]]
        polygon00 = [voronoi.vertices[vertex_idx] if vertex_idx != -
                     1 else intersection_points[0] for vertex_idx in region0]
        polygon01 = [voronoi.vertices[vertex_idx] if vertex_idx != -
                     1 else intersection_points[1] for vertex_idx in region0]
        region1 = voronoi.regions[voronoi.point_region[point_pair[1]]]
        polygon10 = [voronoi.vertices[vertex_idx] if vertex_idx != -
                     1 else intersection_points[0] for vertex_idx in region1]
        polygon11 = [voronoi.vertices[vertex_idx] if vertex_idx != -
                     1 else intersection_points[1] for vertex_idx in region1]
        is_convex00 = is_polygon_convex(polygon00)
        is_convex01 = is_polygon_convex(polygon01)
        is_convex10 = is_polygon_convex(polygon10)
        is_convex11 = is_polygon_convex(polygon11)
        is_convex0 = is_convex00 and is_convex10  # 0を追加で凸包
        is_convex1 = is_convex01 and is_convex11  # 1を追加で凸包
        vidx = voronoi.ridge_dict[point_pair][1]  # [0]側が-1
        innerPoint = voronoi.vertices[vidx]
        rs = np.linalg.norm(intersection_points-innerPoint, axis=1)
        if is_convex0 and not is_convex1:
            infinite_region_polygons.setdefault(point_pair[0], []).append(polygon00)
            infinite_region_polygons.setdefault(point_pair[1], []).append(polygon10)
        elif not is_convex0 and is_convex1:
            infinite_region_polygons.setdefault(point_pair[0], []).append(polygon01)
            infinite_region_polygons.setdefault(point_pair[1], []).append(polygon11)
        elif is_convex0 and is_convex1:
            if rs[0] < rs[1]:
                infinite_region_polygons.setdefault(point_pair[0], []).append(polygon00)
                infinite_region_polygons.setdefault(point_pair[1], []).append(polygon10)
            else:
                infinite_region_polygons.setdefault(point_pair[0], []).append(polygon01)
                infinite_region_polygons.setdefault(point_pair[1], []).append(polygon11)

    infinite_polygons_pidx = []
    infinite_polygons = []
    infinite_colors = []
    for point_idx, polygons in infinite_region_polygons.items():
        infinite_polygons_pidx.append(point_idx)
        infinite_polygons.append(calculate_convex_hull(polygons))
        infinite_colors.append(lighten_color(color_map[labels[point_idx]]))

    if True:
        # 四隅の点がどのpolygonにも所属していない場合に相手をみつける．
        for q in [[xlim[0], ylim[0]], [xlim[0], ylim[1]], [xlim[1], ylim[0]], [xlim[1], ylim[1]]]:
            fpinq = [is_point_in_polygon(polygon, q) for polygon in finite_polygons]
            ifpinq = [is_point_in_polygon(polygon, q) for polygon in infinite_polygons]
            if not any(fpinq) and not any(ifpinq):
                best_f_poly_idx, best_f_vertex_idx, best_f_is_vertex, min_f_distance = closest_polygon_and_vertex(
                    finite_polygons, q)
                best_if_poly_idx, best_if_vertex_idx, best_if_is_vertex, min_if_distance = closest_polygon_and_vertex(
                    infinite_polygons, q)
                if best_f_poly_idx == -1 and best_if_poly_idx == -1:
                    continue
                isBoth = best_f_poly_idx != -1 and best_if_poly_idx != -1
                if best_f_poly_idx != -1 and best_if_poly_idx == -1 or (isBoth and min_f_distance < min_if_distance):
                    if not best_f_is_vertex:
                        n = len(finite_polygons[best_f_poly_idx])
                        a = np.insert(
                            finite_polygons[best_f_poly_idx], best_f_vertex_idx, q, axis=0)
                        b = np.insert(
                            finite_polygons[best_f_poly_idx], (best_f_vertex_idx + 1) % n, q, axis=0)
                        if is_polygon_simple(a):
                            finite_polygons[best_f_poly_idx] = a
                        elif is_polygon_simple(b):
                            finite_polygons[best_f_poly_idx] = b
                elif best_f_poly_idx == -1 and best_if_poly_idx != -1 or (isBoth and min_f_distance >= min_if_distance):
                    if not best_if_is_vertex:
                        n = len(infinite_polygons[best_if_poly_idx])
                        a = np.insert(
                            infinite_polygons[best_if_poly_idx], best_if_vertex_idx, q, axis=0)
                        b = np.insert(
                            infinite_polygons[best_if_poly_idx], (best_if_vertex_idx + 1) % n, q, axis=0)
                        if is_polygon_simple(a):
                            infinite_polygons[best_if_poly_idx] = a
                        elif is_polygon_simple(b):
                            infinite_polygons[best_if_poly_idx] = b

    ax.add_collection(PolyCollection(infinite_polygons, facecolors=infinite_colors, edgecolors='gray'))

    for point_idx, point in enumerate(voronoi.points):
        label_margin = 0.1
        x, y = point
        if xlim[0] - label_margin <= x <= xlim[1] + label_margin and ylim[0] - label_margin <= y <= ylim[1] + label_margin:
            ax.text(
                x, y, str(labels[point_idx]),
                # x, y, str(point_idx),
                color=color_map[labels[point_idx]], size=15, ha='center', va='center'
            )
    if False:
        if False:
            delta = radius
            ax.set_xlim([f([xlim[0], xlim[1], center[0]-delta, center[0]+delta]) for f in [min, max]])
            ax.set_ylim([f([ylim[0], ylim[1], center[1]-delta, center[1]+delta]) for f in [min, max]])
            ax.add_patch(plt.Circle(center, radius, color='r', fill=False))
            for vidx, (x, y) in enumerate(voronoi.vertices):
                ax.text(x, y, str(vidx), color='black', size=15)
        else:
            delta = 0.1
            ax.set_xlim(xlim[0]-delta, xlim[1]+delta)
            ax.set_ylim(ylim[0]-delta, ylim[1]+delta)
            for vidx, (x, y) in enumerate(voronoi.vertices):
                x = xlim[0]-delta if x < xlim[0]-delta else (xlim[1]+delta if x > xlim[1]+delta else x)
                y = ylim[0]-delta if y < ylim[0]-delta else (ylim[1]+delta if y > ylim[1]+delta else y)
                # ax.text(x, y, str(vidx), color='black', size=15)
    set_uniform_grid_ticks(ax)

    if not is7x7:
        if use_pca:
            title = f'PCA (PC{component_indices[0]+1}+PC{component_indices[1]+1}'
            variance_ratio = pca.explained_variance_ratio_[component_indices].sum()
            title += f', Variance Ratio: {variance_ratio:.2f})'
            ax.set_title(title)
        method = 'evd' if use_pca else 'svd'
        component_str = f'pc{component_indices[0]+1}_pc{component_indices[1]+1}'
        plt.savefig(f'voronoi_{method}_{component_str}.png')
        plt.show()


def plot7x7():
    pca, feature_data, labels = prepare_data_IMP()
    fig = plt.figure(facecolor='white', figsize=(8.27*3, 8.27*3), dpi=300)
    for i in range(7):
        for j in range(7):
            component_indices = [i, j]
            ax = plt.subplot2grid((7, 7), component_indices)
            if i == j:
                v = pca.explained_variance_ratio_[i]
                ax.set_axis_off()
                ax.text(0.5, 0.5, f'{(v*100):0.1f}%', size=40, ha='center', va='center', transform=ax.transAxes)
            else:
                projected_data_pca = feature_data @ pca.components_[component_indices].T
                voronoi = Voronoi(projected_data_pca)
                plot_voronoi(
                    voronoi, pca, projected_data_pca, labels,
                    use_pca=True, component_indices=component_indices,
                    ax=ax
                )
    plt.tight_layout()
    plt.savefig(f'voronoi_7x7.svg')
    plt.show()


def main():
    """Main function to prepare data and plot Voronoi diagrams for different component pairs."""
    use_pca = True
    # Test different component combinations
    component_pairs = [
        # [0, 1],  # PC1+PC2 (default)
        # [0, 2],  # PC1+PC3
        # [0, 3],  # PC1+PC4
        # [0, 4],  # PC1+PC5
        # [0, 5],  # PC1+PC6
        # [0, 6],  # PC1+PC7
        # [1, 2],  # PC2+PC3
        # [1, 3],  # PC2+PC4
        # [1, 4],  # PC2+PC5
        # [1, 5],  # PC2+PC6
        # [1, 6],  # PC2+PC7
        # [2, 3],  # PC3+PC4
        # [2, 4],  # PC3+PC5
        # [2, 5],  # PC3+PC6
        # [2, 6],  # PC3+PC7
        [3, 4],  # PC4+PC5
        # [3, 5],  # PC4+PC6
        # [3, 6],  # PC4+PC5
        # [4, 5],  # PC4+PC5
        # [4, 6],  # PC4+PC5
        # [5, 6],  # PC4+PC5
    ]
    # component_pairs = [[i, j] for i in range(6) for j in range(i+1, 7)]
    for component_indices in component_pairs:
        voronoi, pca, projected_data, labels, component_indices = prepare_data(
            use_pca=use_pca, component_indices=component_indices
        )
        plot_voronoi(
            voronoi, pca, projected_data, labels,
            use_pca=use_pca, component_indices=component_indices
        )


if __name__ == "__main__":
    main()
    # plot7x7()