4. 随机子孔径划分与稀疏采样

任务说明

本案例针对大口径光学检测中因遮挡或成本限制导致的子孔径稀疏分布问题，生成随机分布的子孔径集合，模拟稀疏采样场景，并提取各子孔径内的波前斜率数据，为后续波前重构算法提供非规则采样输入。

实现方案

采用正方形子孔径在圆形全口径内随机撒点的方式模拟稀疏采样。每个子孔径为边长d的正方形，中心位置在口径内均匀随机分布，通过重叠率约束控制子孔径间的最小间距，避免过度聚集。波前相位由Zernike多项式前四项（平移、倾斜、离焦、像散）线性组合生成，斜率数据通过解析求导获得，保证提取结果的精确性。

核心代码

import numpy as np

def zernike_gradient(x, y, coeffs):
    """计算Zernike多项式前四项的梯度 (dx, dy)
    coeffs: [piston, tip, tilt, defocus, astig1, astig2]
    返回: (gx, gy) 与x,y同形状
    """
    # Z1: piston -> 0
    # Z2: 2x
    # Z3: 2y
    # Z4: sqrt(3)*(2r^2-1) -> d/dx = 4*sqrt(3)*x
    # Z5: 2*sqrt(6)*x*y -> d/dx = 2*sqrt(6)*y, d/dy = 2*sqrt(6)*x
    # Z6: sqrt(6)*(x^2-y^2) -> d/dx = 2*sqrt(6)*x, d/dy = -2*sqrt(6)*y
    gx = (2*coeffs[1] + 4*np.sqrt(3)*coeffs[3]*x +
          2*np.sqrt(6)*coeffs[4]*y + 2*np.sqrt(6)*coeffs[5]*x)
    gy = (2*coeffs[2] + 4*np.sqrt(3)*coeffs[3]*y +
          2*np.sqrt(6)*coeffs[4]*x - 2*np.sqrt(6)*coeffs[5]*y)
    return gx, gy

def generate_random_subaps(D=100.0, d=10.0, N=12, seed=42):
    """在直径D的圆形口径内随机生成N个边长为d的正方形子孔径
    返回: centers (N,2) 每个子孔径中心坐标 (x,y)
    """
    rng = np.random.default_rng(seed)
    R = D/2 - d/2  # 中心可移动范围半径
    centers = []
    min_dist = d * 0.8  # 最小中心间距约束
    max_attempts = 500
    while len(centers) < N:
        for _ in range(max_attempts):
            # 在圆内均匀采样
            r = R * np.sqrt(rng.random())
            theta = rng.uniform(0, 2*np.pi)
            x = r * np.cos(theta)
            y = r * np.sin(theta)
            # 检查间距约束
            if all(np.hypot(x-cx, y-cy) >= min_dist for cx, cy in centers):
                centers.append((x, y))
                break
        else:
            raise RuntimeError(f"无法在{max_attempts}次尝试内放置所有子孔径")
    return np.array(centers)

def extract_slope_data(centers, d=10.0, coeffs=None):
    """提取每个子孔径的平均斜率
    返回: slopes (N,2) 每行 [gx_mean, gy_mean]
    """
    if coeffs is None:
        coeffs = [0.0, 0.3, -0.2, 0.1, 0.15, -0.05]
    # 子孔径内采样点（4个角点+中心）
    offsets = np.array([[-0.5, -0.5], [0.5, -0.5],
                        [-0.5, 0.5], [0.5, 0.5], [0, 0]]) * d
    slopes = []
    for cx, cy in centers:
        px = cx + offsets[:, 0]
        py = cy + offsets[:, 1]
        gx, gy = zernike_gradient(px, py, coeffs)
        slopes.append([gx.mean(), gy.mean()])
    return np.array(slopes)

# 生成子孔径并提取斜率
centers = generate_random_subaps(N=12, seed=42)
slopes = extract_slope_data(centers)

# 输出结果（前5个子孔径）
print("子孔径中心坐标与平均斜率：")
for i in range(min(5, len(centers))):
    cx, cy = centers[i]
    gx, gy = slopes[i]
    print(f"  [{i+1}] center=({cx:7.2f},{cy:7.2f}) mm, "
          f"slope=({gx: .4f},{gy: .4f}) rad")
print(f"共生成 {len(centers)} 个子孔径，斜率数据形状: {slopes.shape}")

运行与验证

依赖：仅需NumPy库（pip install numpy）。

运行方式：将代码保存为sparse_subaperture.py，在命令行执行python sparse_subaperture.py。

预期输出：程序打印5行子孔径信息及汇总行，例如第一行格式为[1] center=( -3.21, 12.45) mm, slope=( 0.2845,-0.1932) rad，最后一行显示共生成 12 个子孔径，斜率数据形状: (12, 2)。

验证方法：

1. 检查所有子孔径中心坐标满足sqrt(x²+y²) ≤ 45（即口径半径50减去子孔径半宽5），确保无越界；
2. 验证任意两中心间距≥8 mm（即d*0.8），满足稀疏约束；
3. 将提取的斜率与在子孔径中心直接调用zernike_gradient的结果对比，差异应小于0.01 rad（因内部采样平均效应），证明斜率提取正确。

判定标准：程序无异常退出，输出12组斜率数据，且上述三项检查全部通过，即可判定案例实现成功。