1. 基于Zernike多项式拟合的波前重构入门

任务说明

本案例实现基于Zernike多项式的最小二乘波前重构：给定一组离散采样点上的波前斜率（x、y方向差分），通过Zernike模式法拟合出波前相位分布，并输出重构波前与残差。这是自适应光学中从波前传感器斜率数据恢复相位的最基础方法。

实现方案

采用前36项Zernike多项式（标准Noll序号，对应径向阶数n≤7），在单位圆内规则网格上构造模式矩阵，利用最小二乘求解模式系数，进而重构波前。输入为模拟生成的已知相位（由前几项Zernike叠加），取其解析差分作为“测量斜率”，验证重构精度。

代码实现

import numpy as np
from math import sqrt, atan2, factorial

def zernike_poly(n, m, rho, theta):
    """返回第(n,m)项Zernike在极坐标(rho,theta)处的值（Noll归一化）"""
    if m >= 0:
        # 径向多项式 R_n^m
        s = 0
        for k in range((n - m)//2 + 1):
            coef = (-1)**k * factorial(n-k)
            coef /= (factorial(k) * factorial((n+m)//2 - k) * factorial((n-m)//2 - k))
            s += coef * rho**(n - 2*k)
        # Noll归一化系数
        norm = sqrt(2*(n+1)) if m != 0 else sqrt(n+1)
        return norm * s * np.cos(m * theta)
    else:
        m_abs = -m
        s = 0
        for k in range((n - m_abs)//2 + 1):
            coef = (-1)**k * factorial(n-k)
            coef /= (factorial(k) * factorial((n+m_abs)//2 - k) * factorial((n-m_abs)//2 - k))
            s += coef * rho**(n - 2*k)
        norm = sqrt(2*(n+1))
        return norm * s * np.sin(m_abs * theta)

def noll_to_poly(j):
    """Noll序号j(从1开始)转(n,m)"""
    n = 0
    while (n+1)*(n+2)//2 < j:
        n += 1
    # 该n阶起始序号
    start = n*(n+1)//2 + 1
    pos = j - start
    # 同阶内排列：m = n, n-2, ..., -n
    m = n - 2*pos
    return n, m

def zernike_derivatives(x, y, j_list):
    """
    计算单位圆内点(x,y)处各Zernike项对x和y的导数
    返回 Gx, Gy: 形状(len(j_list), len(x))
    """
    rho = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    Gx = np.zeros((len(j_list), len(x)))
    Gy = np.zeros((len(j_list), len(x)))
    eps = 1e-10
    for idx, j in enumerate(j_list):
        n, m = noll_to_poly(j)
        # 数值导数（解析式复杂，用中心差分保证精度）
        h = 1e-6
        for i in range(len(x)):
            # d/dx
            xp = x[i] + h
            xm = x[i] - h
            rp = sqrt(xp**2 + y[i]**2)
            rm = sqrt(xm**2 + y[i]**2)
            tp = atan2(y[i], xp)
            tm = atan2(y[i], xm)
            zp = zernike_poly(n, m, rp, tp)
            zm = zernike_poly(n, m, rm, tm)
            Gx[idx, i] = (zp - zm) / (2*h)
            # d/dy
            yp = y[i] + h
            ym = y[i] - h
            rp = sqrt(x[i]**2 + yp**2)
            rm = sqrt(x[i]**2 + ym**2)
            tp = atan2(yp, x[i])
            tm = atan2(ym, x[i])
            zp = zernike_poly(n, m, rp, tp)
            zm = zernike_poly(n, m, rm, tm)
            Gy[idx, i] = (zp - zm) / (2*h)
    return Gx, Gy

def reconstruct_from_slopes(x, y, gx_meas, gy_meas, n_terms=36):
    """
    最小二乘波前重构
    参数:
        x, y: 采样点坐标(单位圆内)
        gx_meas, gy_meas: 测量的x/y斜率
        n_terms: 使用的Zernike项数(Noll序号1~n_terms)
    返回:
        coeffs: Zernike系数
        W_rec: 重构波前(与x同形状)
    """
    j_list = list(range(1, n_terms+1))
    Gx, Gy = zernike_derivatives(x, y, j_list)
    # 构造模式矩阵 A: 每行对应一个采样点的两个斜率方程
    n_pts = len(x)
    A = np.zeros((2*n_pts, n_terms))
    b = np.zeros(2*n_pts)
    for i in range(n_pts):
        A[2*i, :] = Gx[:, i]
        A[2*i+1, :] = Gy[:, i]
        b[2*i] = gx_meas[i]
        b[2*i+1] = gy_meas[i]
    # 最小二乘
    coeffs, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
    # 重构波前
    W_rec = np.zeros(n_pts)
    for idx, j in enumerate(j_list):
        n, m = noll_to_poly(j)
        rho = np.sqrt(x**2 + y**2)
        theta = np.arctan2(y, x)
        W_rec += coeffs[idx] * zernike_poly(n, m, rho, theta)
    return coeffs, W_rec

# ============ 测试与验证 ============
if __name__ == "__main__":
    # 1. 生成单位圆内规则网格采样点
    N = 20  # 每边采样数
    x_lin = np.linspace(-0.95, 0.95, N)
    X, Y = np.meshgrid(x_lin, x_lin)
    mask = X**2 + Y**2 <= 0.95**2  # 圆内点
    x = X[mask]
    y = Y[mask]

    # 2. 构造真实波前: Z4(散焦) + Z5(像散) + Z6(像散45°)
    true_coeffs = {4: 0.8, 5: 0.5, 6: -0.3}
    W_true = np.zeros_like(x)
    for j, c in true_coeffs.items():
        n, m = noll_to_poly(j)
        rho = np.sqrt(x**2 + y**2)
        theta = np.arctan2(y, x)
        W_true += c * zernike_poly(n, m, rho, theta)

    # 3. 计算“测量”斜率(解析差分近似)
    h = 1e-6
    gx_meas = np.zeros_like(x)
    gy_meas = np.zeros_like(x)
    for i in range(len(x)):
        # 对真实波前求数值梯度
        xp = x[i] + h
        xm = x[i] - h
        rp = sqrt(xp**2 + y[i]**2)
        rm = sqrt(xm**2 + y[i]**2)
        tp = atan2(y[i], xp)
        tm = atan2(y[i], xm)
        Wp = sum(c * zernike_poly(*noll_to_poly(j), rp, tp) for j, c in true_coeffs.items())
        Wm = sum(c * zernike_poly(*noll_to_poly(j), rm, tm) for j, c in true_coeffs.items())
        gx_meas[i] = (Wp - Wm) / (2*h)
        yp = y[i] + h
        ym = y[i] - h
        rp = sqrt(x[i]**2 + yp**2)
        rm = sqrt(x[i]**2 + ym**2)
        tp = atan2(yp, x[i])
        tm = atan2(ym, x[i])
        Wp = sum(c * zernike_poly(*noll_to_poly(j), rp, tp) for j, c in true_coeffs.items())
        Wm = sum(c * zernike_poly(*noll_to_poly(j), rm, tm) for j, c in true_coeffs.items())
        gy_meas[i] = (Wp - Wm) / (2*h)

    # 4. 执行重构
    coeffs, W_rec = reconstruct_from_slopes(x, y, gx_meas, gy_meas, n_terms=36)

    # 5. 评估
    residual = W_true - W_rec
    rms_true = np.sqrt(np.mean(W_true**2))
    rms_res = np.sqrt(np.mean(residual**2))

    print(f"采样点数: {len(x)}")
    print(f"真实波前RMS: {rms_true:.4f} rad")
    print(f"重构残差RMS: {rms_res:.4f} rad")
    print(f"相对误差: {rms_res/rms_true*100:.2f}%")
    print("\n前6项重构系数(理论值: Z4=0.8, Z5=0.5, Z6=-0.3):")
    for i in range(6):
        print(f"  Z{i+1}: {coeffs[i]:.4f}")

运行与验证

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

预期输出：程序打印采样点数、真实波前RMS、重构残差RMS、相对误差百分比以及前6项Zernike系数。由于真实波前仅由Z4、Z5、Z6构成，且斜率采样无噪声，重构出的Z4、Z5、Z6系数应分别接近0.8、0.5、-0.3，其余系数接近0，残差RMS应远小于真实波前RMS（相对误差通常在1%以下，主要来源于数值差分与网格离散误差）。

判定标准：若Z4、Z5、Z6的重构系数与理论值偏差小于0.05，且相对误差小于5%，则本案例实现正确。可自行修改true_coeffs中的系数或增加Zernike项数（如加入Z7彗差）来进一步验证算法的通用性。