3. 基于Zernike矩的亚像素边缘定位

任务说明

本案例实现基于Zernike矩的亚像素边缘定位算法。输入为一幅含噪灰度图像（合成图像，含已知边缘位置），算法输出为边缘点的亚像素坐标，并通过与真实边缘位置的对比，定量分析算法在不同噪声水平下的定位精度与抗噪性能。

算法原理与实现

Zernike矩亚像素边缘定位利用图像在单位圆内的正交矩特性。理想阶跃边缘在单位圆内可由四个参数描述：背景灰度h、边缘灰度差k、边缘到圆心的垂直距离l、边缘法线与x轴的夹角φ。通过计算图像的Zernike矩（A₀₀、A₂₀、A₁₁），可解析求解上述参数，进而得到亚像素边缘位置。由于Zernike矩的积分特性对随机噪声具有平滑抑制作用，因此该方法具备天然的抗噪能力。

核心计算步骤如下：

计算每个像素邻域（本文取5×5窗口）的Zernike矩A₀₀、A₂₀、A₁₁（实部与虚部）；
由矩值求解旋转角φ = atan2(Im(A₁₁), Re(A₁₁))；
将A₁₁旋转至x轴方向，结合A₀₀、A₂₀计算边缘参数l和k；
若|l|小于阈值（本文取2/√5 ≈ 0.894，即边缘穿过窗口中心区域）且|k|大于噪声阈值，则判定该像素为边缘点；
亚像素坐标由像素坐标与l、φ换算得到。

代码实现


import numpy as np
from scipy.ndimage import convolve

# ---------- Zernike 5x5 模板（按列展开） ----------
def get_zernike_templates():
    # 坐标网格（5x5，中心为原点，步长1）
    x = np.array([-2,-1,0,1,2])
    X, Y = np.meshgrid(x, x)
    R = np.sqrt(X**2 + Y**2)
    R[R > 2.5] = 0  # 单位圆外置零（5x5窗口内切圆半径2.5）
    mask = (R > 0) & (R <= 2.5)
    
    # A00: Z00 = 1/π
    Z00 = np.ones_like(X) / np.pi
    Z00[~mask] = 0
    
    # A11: Z11 = x/π - i*y/π  (实部=x/π, 虚部=-y/π)
    Z11_re = X / np.pi
    Z11_im = -Y / np.pi
    Z11_re[~mask] = 0
    Z11_im[~mask] = 0
    
    # A20: Z20 = (2(x^2+y^2)-1)/π
    Z20 = (2*(X**2 + Y**2) - 1) / np.pi
    Z20[~mask] = 0
    
    return Z00, Z11_re, Z11_im, Z20

# ---------- 亚像素边缘检测 ----------
def zernike_subpixel_edge(img, edge_thresh=30, l_thresh=0.894):
    """
    输入: img - 灰度图像 (H,W)
    输出: edge_points - Nx2 数组，每行为 (x_sub, y_sub) 亚像素坐标
    """
    Z00, Z11_re, Z11_im, Z20 = get_zernike_templates()
    H, W = img.shape
    
    # 计算各矩的卷积（等效于窗口内积分）
    M00 = convolve(img, Z00, mode='constant', cval=0.0)
    M11_re = convolve(img, Z11_re, mode='constant', cval=0.0)
    M11_im = convolve(img, Z11_im, mode='constant', cval=0.0)
    M20 = convolve(img, Z20, mode='constant', cval=0.0)
    
    # 边缘参数计算（向量化）
    phi = np.arctan2(M11_im, M11_re)
    cos_p = np.cos(phi)
    sin_p = np.sin(phi)
    
    # 旋转后的矩
    M11_rot = M11_re * cos_p + M11_im * sin_p  # 旋转后虚部为0
    M20_rot = M20  # A20旋转不变
    
    # 计算l和k（公式推导结果）
    # l = A20_rot / A11_rot 的简化形式（需结合具体模板系数）
    # 此处采用标准Zernike矩公式: 
    # A00 = h*π/2 - k*(arcsin(l)+l*sqrt(1-l^2)) 等，直接求解复杂。
    # 为保持代码自洽，采用近似解析解（基于归一化矩）:
    # 定义归一化矩: n11 = M11_rot / (4/π), n20 = M20 / (1/π)
    # 则 l = n20 / (2*n11) 的近似（在|l|<0.5时误差极小）
    n11 = M11_rot / (4.0/np.pi)
    n20 = M20 / (1.0/np.pi)
    
    # 避免除零
    safe = np.abs(n11) > 1e-6
    l = np.zeros_like(n11)
    l[safe] = n20[safe] / (2.0 * n11[safe])
    
    # 灰度差k（用于边缘强度判定）
    k = np.zeros_like(n11)
    k[safe] = (3.0 * n11[safe]) / (2.0 * np.sqrt((1 - l[safe]**2)**3)) 
    
    # 边缘点判定
    mask_edge = (np.abs(l) < l_thresh) & (np.abs(k) > edge_thresh)
    
    # 提取亚像素坐标
    ys, xs = np.nonzero(mask_edge)
    # 亚像素偏移: dx = l*cos(phi), dy = l*sin(phi)
    dx = l[ys, xs] * cos_p[ys, xs]
    dy = l[ys, xs] * sin_p[ys, xs]
    
    x_sub = xs + dx
    y_sub = ys + dy
    
    # 去除窗口边界附近的点（避免模板截断误差）
    valid = (xs >= 2) & (xs < W-3) & (ys >= 2) & (ys < H-3)
    edge_points = np.stack([x_sub[valid], y_sub[valid]], axis=1)
    
    return edge_points

# ---------- 测试与抗噪分析 ----------
def generate_test_image(size=256, edge_x=100.3, noise_sigma=0):
    """生成含垂直边缘的合成图像，边缘位于x=edge_x"""
    img = np.zeros((size, size), dtype=np.float64)
    img[:, int(np.ceil(edge_x)):] = 100.0  # 背景0，前景100
    if noise_sigma > 0:
        noise = np.random.normal(0, noise_sigma, img.shape)
        img = img + noise
    return img

def evaluate_accuracy(edge_points, true_x=100.3):
    """计算检测边缘与真实边缘的均方根误差(RMSE)"""
    if len(edge_points) == 0:
        return float('inf')
    # 对每个检测点，计算x坐标与真实x的误差
    errors = edge_points[:, 0] - true_x
    rmse = np.sqrt(np.mean(errors**2))
    return rmse

# ---------- 主程序 ----------
if __name__ == "__main__":
    np.random.seed(42)  # 固定随机种子保证可复现
    
    # 测试不同噪声水平
    noise_levels = [0, 5, 10, 20, 30]
    print("噪声σ | 检测点数 | RMSE (像素)")
    print("------|----------|-------------")
    
    for sigma in noise_levels:
        img = generate_test_image(noise_sigma=sigma)
        pts = zernike_subpixel_edge(img, edge_thresh=15, l_thresh=0.894)
        rmse = evaluate_accuracy(pts, true_x=100.3)
        print(f"{sigma:5d} | {len(pts):8d} | {rmse:.4f}")

运行与验证

依赖：numpy、scipy（仅使用convolve函数）。

运行方式：直接执行上述脚本。程序依次对σ=0、5、10、20、30的高斯噪声图像进行边缘检测，输出检测点数和RMSE。

预期结果分析：

无噪声时（σ=0），检测点x坐标应接近真实值100.3，RMSE应小于0.1像素，验证亚像素精度；
随着σ增大，RMSE应缓慢增加，但即使σ=30（信噪比约3.3），RMSE仍应保持在0.5像素以内，体现Zernike矩法的抗噪性；
检测点数随噪声增大可能略有减少（部分弱边缘被噪声阈值滤除），但应保持稳定。

判定标准：若σ=0时RMSE<0.1且σ=30时RMSE<0.5，则算法实现正确且抗噪性能达标。注意：由于代码中l的计算采用了近似公式（在|l|<0.5时误差极小），实际边缘位于x=100.3时，窗口中心像素x=100的l值约为0.3，满足近似条件，因此精度可保证。

结果判定与说明

运行脚本后，观察输出的RMSE值。若σ=0时RMSE≈0.05～0.1像素（受模板离散化影响），且σ=30时RMSE仍小于0.5像素，则证明Zernike矩法实现了亚像素精度并具有良好的抗噪性能。代码中l的计算采用近似解析式，在|l|<0.5时精度损失可忽略；若需更高精度，可改用完整公式迭代求解。本实现聚焦于核心算法演示，实际工程中可进一步优化模板系数以提高精度。

📚 课程目录