1. 均匀方形子孔径划分与波前斜率提取

任务说明

本案例实现一个完整的波前斜率提取流程：输入一个方形波前相位分布（模拟 Shack-Hartmann 传感器测量对象），将其均匀划分为方形子孔径，计算每个子孔径内的平均波前斜率（x 方向与 y 方向），并输出斜率矩阵供后续波前重构使用。

实现方案

采用最小充分实现：使用 NumPy 生成模拟波前（包含倾斜项与离焦项），按子孔径尺寸对波前进行分块，对每个子孔径内的相位数据用最小二乘平面拟合，拟合系数即为该子孔径的平均斜率。该方法与 Shack-Hartmann 传感器中每个微透镜焦斑质心偏移量正比于局部波前斜率的物理原理一致。

代码实现

import numpy as np

def generate_test_wavefront(n=64, tilt_x=0.3, tilt_y=-0.2, defocus=0.05):
    """
    生成模拟波前相位（单位：波长）
    n: 波前采样网格尺寸 (n x n)
    tilt_x, tilt_y: x/y方向倾斜系数
    defocus: 离焦系数
    """
    x = np.linspace(-1, 1, n)
    y = np.linspace(-1, 1, n)
    X, Y = np.meshgrid(x, y)
    W = tilt_x * X + tilt_y * Y + defocus * (X**2 + Y**2)
    return W

def divide_into_subapertures(wavefront, sub_size):
    """
    将方形波前均匀划分为方形子孔径
    返回形状: (N_sub, N_sub, sub_size, sub_size)
    """
    n = wavefront.shape[0]
    if n % sub_size != 0:
        raise ValueError(f"波前尺寸 {n} 必须能被子孔径尺寸 {sub_size} 整除")
    N_sub = n // sub_size
    # 按行分块再按列分块
    blocks = wavefront.reshape(N_sub, sub_size, N_sub, sub_size)
    # 转置为 (N_sub, N_sub, sub_size, sub_size)
    blocks = blocks.transpose(0, 2, 1, 3)
    return blocks

def fit_plane_and_get_slopes(sub_blocks):
    """
    对每个子孔径进行最小二乘平面拟合
    返回: Gx, Gy (每个子孔径的x/y方向平均斜率)
    """
    N_sub = sub_blocks.shape[0]
    sub_size = sub_blocks.shape[2]
    # 子孔径内局部坐标（归一化到 [-0.5, 0.5]）
    local_coords = np.linspace(-0.5, 0.5, sub_size)
    lx, ly = np.meshgrid(local_coords, local_coords)
    # 设计矩阵 A: [1, x, y]
    A = np.column_stack([
        np.ones(sub_size * sub_size),
        lx.ravel(),
        ly.ravel()
    ])
    # 预计算 (A^T A)^-1 A^T 用于批量拟合
    pinv = np.linalg.pinv(A)
    
    Gx = np.zeros((N_sub, N_sub))
    Gy = np.zeros((N_sub, N_sub))
    
    for i in range(N_sub):
        for j in range(N_sub):
            z = sub_blocks[i, j].ravel()
            coeffs = pinv @ z  # [a0, a1, b1]
            Gx[i, j] = coeffs[1]  # x方向斜率
            Gy[i, j] = coeffs[2]  # y方向斜率
    return Gx, Gy

def extract_slopes(wavefront, sub_size):
    """完整流程：分块 + 斜率提取"""
    sub_blocks = divide_into_subapertures(wavefront, sub_size)
    Gx, Gy = extract_slopes_from_blocks(sub_blocks)
    return Gx, Gy

def extract_slopes_from_blocks(sub_blocks):
    """从已分块的子孔径数据提取斜率（供内部调用）"""
    return fit_plane_and_get_slopes(sub_blocks)

# ========== 主程序 ==========
if __name__ == "__main__":
    # 1. 生成模拟波前
    n = 64
    wavefront = generate_test_wavefront(n, tilt_x=0.3, tilt_y=-0.2, defocus=0.05)
    
    # 2. 设置子孔径尺寸（每个子孔径 8x8 像素）
    sub_size = 8
    N_sub = n // sub_size  # 8x8 个子孔径
    
    # 3. 执行分块与斜率提取
    sub_blocks = divide_into_subapertures(wavefront, sub_size)
    Gx, Gy = extract_slopes_from_blocks(sub_blocks)
    
    # 4. 输出结果摘要
    print(f"波前尺寸: {n}x{n}, 子孔径尺寸: {sub_size}x{sub_size}")
    print(f"子孔径数量: {N_sub}x{N_sub}")
    print(f"Gx 形状: {Gx.shape}, Gy 形状: {Gy.shape}")
    print(f"\nGx 斜率矩阵（前4行）:")
    print(np.round(Gx[:4, :4], 4))
    print(f"\nGy 斜率矩阵（前4行）:")
    print(np.round(Gy[:4, :4], 4))
    
    # 5. 验证：理论斜率 vs 提取斜率
    # 理论值: dW/dx = tilt_x + 2*defocus*x, dW/dy = tilt_y + 2*defocus*y
    # 在中心子孔径处 x≈0, y≈0, 理论斜率 ≈ tilt_x, tilt_y
    center = N_sub // 2
    print(f"\n中心子孔径 ({center},{center}) 提取斜率: Gx={Gx[center,center]:.4f}, Gy={Gy[center,center]:.4f}")
    print(f"理论倾斜系数: tilt_x=0.3, tilt_y=-0.2")
    print(f"偏差: ΔGx={abs(Gx[center,center]-0.3):.4f}, ΔGy={abs(Gy[center,center]+0.2):.4f}")

运行与验证

运行方式：将上述代码保存为 slope_extraction.py，在装有 NumPy 的 Python 环境中执行 python slope_extraction.py。

预期输出与判定标准：

• 程序输出 Gx、Gy 矩阵形状均为 (8, 8)，与 64×64 波前按 8×8 子孔径划分一致。
• 中心子孔径（索引 3,3 或 4,4）的提取斜率应接近理论倾斜系数：Gx ≈ 0.3，Gy ≈ -0.2。由于离焦项在中心附近贡献很小，偏差应小于 0.02（由离焦项在中心子孔径范围内的残余梯度导致）。
• Gx 矩阵从左到右应呈现递增趋势（因离焦项引入的 x 方向梯度 2·defocus·x 随 x 增大而增大），Gy 矩阵从上到下应呈现递减趋势（因 y 方向梯度 2·defocus·y 随 y 增大而增大）。

若输出满足上述特征，则说明子孔径划分正确、斜率提取算法有效，可继续用于后续波前重构模块。