ImplicitSurfaceNetwork/application/test_circle_gwn.py

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable


def arc_gwn_batch(a, b, bulge, P):
    """
    向量化批处理接口：对 (N, 2) 的点阵 P 同时计算单段圆弧的广义绕数 W。
    """
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    P = np.asarray(P, dtype=float)

    half_ab   = 0.5 * (a + b)
    bisec_vec = 0.5 * np.array([a[1] - b[1], b[0] - a[0]])
    center    = half_ab + bisec_vec / bulge

    ca = a - center
    cb = b - center
    cp = P - center

    D = np.linalg.norm(cp, axis=1)
    r = np.linalg.norm(ca)
    theta = 2.0 * math.atan(bulge)

    W = np.empty(len(P))

    mask_center = D < 1e-12
    mask_circle = np.abs(D - r) < 1e-12
    mask_normal = ~(mask_circle)

    # W[mask_center] = theta / (2.0 * math.pi)
    W[mask_circle] = theta / (4.0 * math.pi)

    cp_n = cp[mask_normal]
    D_n  = D[mask_normal]

    sigma   = np.where(r > D_n, 1.0, -1.0)
    Rd      = r * D_n
    r_p_D   = r + D_n
    abs_r_D = np.abs(r - D_n)

    dot_a  = ca.dot(cp_n.T)
    cross1 = cp_n[:, 0] * ca[1] - ca[0] * cp_n[:, 1]
    Ga     = np.arctan2(r_p_D * cross1, abs_r_D * (Rd + dot_a))

    dot_b  = cb.dot(cp_n.T)
    cross2 = cp_n[:, 0] * cb[1] - cb[0] * cp_n[:, 1]
    Gb     = np.arctan2(r_p_D * cross2, abs_r_D * (Rd + dot_b))

    delta = Gb - Ga
    fix   = delta * theta < 0.0
    delta[fix] += math.copysign(math.pi, theta)

    two_pi_W = 0.5 * theta + sigma * delta
    W[mask_normal] = two_pi_W / (2.0 * math.pi)

    return W


def plot_multi_arc_gwn_field(arcs, margin=0.5, nx=400, ny=400):
    """
    绘制多段圆弧的广义绕数之和。

    参数
    ----
    arcs : list of tuple
        每段圆弧表示为 (a, b, bulge)，其中 a, b 为长度为 2 的序列。
    margin : float
        AABB 外扩边界。
    nx, ny : int
        采样分辨率。
    """
    arcs = [(np.asarray(a, dtype=float), np.asarray(b, dtype=float), float(bulge))
            for a, b, bulge in arcs]

    # ---- 1. 计算所有圆弧的几何信息，确定全局 AABB ----
    all_arc_x = []
    all_arc_y = []
    all_centers = []
    all_rs = []

    for a, b, bulge in arcs:
        half_ab   = 0.5 * (a + b)
        bisec_vec = 0.5 * np.array([a[1] - b[1], b[0] - a[0]])
        center    = half_ab + bisec_vec / bulge
        r         = np.linalg.norm(a - center)
        theta     = 2.0 * math.atan(bulge)

        start_angle = math.atan2(a[1] - center[1], a[0] - center[0])
        end_angle   = start_angle + theta
        t_arc       = np.linspace(start_angle, end_angle, 200)
        arc_x       = center[0] + r * np.cos(t_arc)
        arc_y       = center[1] + r * np.sin(t_arc)

        all_arc_x.append(arc_x)
        all_arc_y.append(arc_y)
        all_centers.append(center)
        all_rs.append(r)

    # 全局 AABB
    x_min = min(ax.min() for ax in all_arc_x) - margin
    x_max = max(ax.max() for ax in all_arc_x) + margin
    y_min = min(ay.min() for ay in all_arc_y) - margin
    y_max = max(ay.max() for ay in all_arc_y) + margin

    # ---- 2. 均匀采样并累加所有圆弧的 GWN ----
    xs = np.linspace(x_min, x_max, nx)
    ys = np.linspace(y_min, y_max, ny)
    X, Y = np.meshgrid(xs, ys)

    points = np.stack([X.ravel(), Y.ravel()], axis=1)
    W = np.zeros(len(points))

    for a, b, bulge in arcs:
        W += arc_gwn_batch(a, b, bulge, points)

    W = W.reshape(X.shape)

    # ---- 3. 三分区独立归一化 ----
    eps        = np.finfo(float).eps * 1e6
    bound_low  = -eps
    bound_high =  eps

    # 负区
    W_neg = W.copy()
    W_neg[W >= bound_low] = np.nan
    W_min = np.nanmin(W_neg)
    if not np.isnan(W_min):
        W_neg = (bound_low - W_neg) / (bound_low - W_min)

    # 零区
    W_zero = W.copy()
    W_zero[(W < bound_low) | (W > bound_high)] = np.nan
    W_zero = (W_zero - bound_low) / (2.0 * eps)

    # 正区
    W_pos = W.copy()
    W_pos[W <= bound_high] = np.nan
    W_max = np.nanmax(W_pos)
    if not np.isnan(W_max):
        W_pos = (W_pos - bound_high) / (W_max - bound_high)

    # ---- 4. 绘图 ----
    fig, ax = plt.subplots(figsize=(10, 8))

    im_neg  = ax.pcolormesh(X, Y, W_neg,  cmap='Blues',  vmin=0, vmax=1, shading='auto')
    im_zero = ax.pcolormesh(X, Y, W_zero, cmap='Greys',  vmin=0, vmax=1, shading='auto')
    im_pos  = ax.pcolormesh(X, Y, W_pos,  cmap='Reds',   vmin=0, vmax=1, shading='auto')

    # 叠加所有圆弧
    for idx, (a, b, bulge) in enumerate(arcs):
        ax.plot(all_arc_x[idx], all_arc_y[idx], 'k-', linewidth=2.5, zorder=4)
        ax.plot(a[0], a[1], 'go', markersize=7, zorder=5)
        ax.plot(b[0], b[1], 'mo', markersize=7, zorder=5)
        ax.plot(all_centers[idx][0], all_centers[idx][1], 'k+', markersize=8, mew=1.5, zorder=5)

    ax.set_aspect('equal')
    ax.set_xlim(x_min, x_max)
    ax.set_ylim(y_min, y_max)
    ax.set_title(f'Multi-Arc GWN Sum  ({len(arcs)} arcs)')

    # ---- 5. 三个独立 colorbar ----
    divider = make_axes_locatable(ax)

    cax_neg = divider.append_axes("left", size="5%", pad=0.4)
    cbar_neg = fig.colorbar(im_neg, cax=cax_neg, orientation='vertical')
    cbar_neg.ax.yaxis.set_ticks_position('left')
    cbar_neg.ax.yaxis.set_label_position('left')
    cbar_neg.set_label('Negative $W$ (normalized)', rotation=90, labelpad=10)
    cbar_neg.set_ticks([0, 1])
    cbar_neg.set_ticklabels([f'{bound_low:.2e}', f'{W_min:.3f}' if not np.isnan(W_min) else 'N/A'])

    cax_pos = divider.append_axes("right", size="5%", pad=0.4)
    cbar_pos = fig.colorbar(im_pos, cax=cax_pos, orientation='vertical')
    cbar_pos.set_label('Positive $W$ (normalized)', rotation=270, labelpad=15)
    cbar_pos.set_ticks([0, 1])
    cbar_pos.set_ticklabels([f'{bound_high:.2e}', f'{W_max:.3f}' if not np.isnan(W_max) else 'N/A'])

    cax_zero = divider.append_axes("bottom", size="5%", pad=0.6)
    cbar_zero = fig.colorbar(im_zero, cax=cax_zero, orientation='horizontal')
    cbar_zero.set_label(f'Near-zero $W$ (normalized)  [{bound_low:.2e}, {bound_high:.2e}]', labelpad=5)
    cbar_zero.set_ticks([0, 1])
    cbar_zero.set_ticklabels([f'{bound_low:.2e}', f'{bound_high:.2e}'])

    plt.tight_layout()
    plt.show()


# ==================== 主程序 ====================
if __name__ == '__main__':
    # 示例 1：两段圆弧拼成近似整圆
    arcs = [
        ([2.0, 0.0], [0.0, 0.0], float("inf")),   # 上半圆
        ([0.0, 0.0], [2.0, 0.0], float("inf")),   # 下半圆（注意 bulge>0 表示逆时针）
    ]
    plot_multi_arc_gwn_field(arcs, margin=0.5, nx=1024, ny=1024)