如何用Python从条纹图案获取表面拓扑？相位解包后步骤求助

我想用条纹图案恢复物体表面形状，现有干涉仪投射圆形条纹，相机拍摄反射图像。研究过相位解包的论文，但没搞懂解包后的相位和最终拓扑的关联。

我已经实现了测试图案生成函数：

def create_pattern_rings(dim, count_rings):
    pattern = np.zeros((dim[0], dim[1]))
    pix_per_ring = int((dim[1] / 2) / count_rings)
    for ring in range(0, count_rings):
        radii_out = (dim[1] / 2) - ring*pix_per_ring*2
        radii_in = radii_out - pix_per_ring
        for y in range(0, dim[0]):
            for x in range(0, dim[1]):
                if (np.sqrt((x-dim[1]/2)**2 + (y-dim[1]/2)**2)) < radii_out and (np.sqrt((x-dim[1]/2)**2 + (y-dim[1]/2)**2)) > radii_in:
                    pattern[y][x] = 1

然后做了FFT和相位解包：

def calc_fft(data, shift):
    data_fft = np.fft.fft2(data)
    if shift:
        data_fft = np.fft.fftshift(data_fft)
    return data_fft

def phase_unwrap(data_fft):
    data_unwraped = np.unwrap(np.angle(data_fft))
    return data_unwraped

现在不知道后续该怎么操作，求建议。

后续步骤建议 #

1. 明确相位与表面高度的物理逻辑 #

圆形条纹干涉的核心是光程差：物体表面某点的高度变化会导致反射光与参考光的光程差改变，直接反映为条纹的相位偏移。解包后的相位$\phi(x,y)$和表面高度$h(x,y)$的定量关系为：
$$h(x,y) = \frac{\lambda \cdot \phi(x,y)}{4\pi}$$
其中$\lambda$是投射条纹的光源波长。系数取$4\pi$是因为光往返物体表面，光程差是高度变化的2倍，对应相位差为$\frac{2\pi}{\lambda} \times 2h$，整理后得到上述公式。

2. 修正当前FFT与相位解包的错误 #

• 你直接对生成的条纹图案做FFT取相位是无效的：实际干涉测量中，相机拍摄的是参考光+物体反射光的干涉条纹，而非单纯的投射条纹。你的测试图案缺少干涉项，得到的相位没有物理意义。
• 2D相位解包不能仅用一维的np.unwrap：np.unwrap(np.angle(...))只沿默认轴解包，处理2D图像会出现大面积相位跳变错误，需要用2D解包逻辑。

3. 构建真实的干涉图模拟 #

要得到有意义的相位，首先需要模拟真实干涉场景，生成包含参考光的干涉图：

def generate_interference_pattern(dim, count_rings, h, lam=633e-9):
    # h: 物体表面高度矩阵，dim: 图像尺寸，lam: 光源波长（这里用633nm可见光示例）
    pattern = np.zeros((dim[0], dim[1]), dtype=np.float32)
    center_x, center_y = dim[1]/2, dim[0]/2
    ring_spacing = (dim[1]/2) / count_rings  # 相邻条纹的半径间隔
    
    for y in range(dim[0]):
        for x in range(dim[1]):
            r = np.sqrt((x-center_x)**2 + (y-center_y)**2)
            proj_phase = 2 * np.pi * r / ring_spacing  # 投射条纹的原始相位
            phase_shift = 4 * np.pi * h[y,x] / lam  # 高度变化带来的相位偏移
            # 干涉光强公式：I = I1 + I2 + 2*sqrt(I1*I2)*cos(相位和)
            I1, I2 = 1.0, 0.8  # 参考光与反射光的光强
            pattern[y,x] = I1 + I2 + 2*np.sqrt(I1*I2)*np.cos(proj_phase + phase_shift)
    return pattern

4. 从干涉图中提取载波相位 #

对干涉图做FFT后，需要分离出携带高度信息的载波分量：

1. 用fftshift将频谱移到中心位置
2. 定位对应条纹载波的峰值（通常在中心两侧对称位置）
3. 保留峰值所在区域，其余频谱置零
4. 逆FFT得到复信号，再提取包裹相位

示例代码：

def extract_carrier_phase(interf_img):
    fft_img = np.fft.fft2(interf_img)
    fft_shifted = np.fft.fftshift(fft_img)
    
    rows, cols = fft_shifted.shape
    cx, cy = rows//2, cols//2
    # 自动检测载波峰值（这里简化为手动指定，实际可通过阈值筛选最大值位置）
    peak_y = cy + int(cols * 0.1)  # 假设峰值在中心右侧10%宽度处
    mask = np.zeros_like(fft_shifted)
    # 给峰值周围保留一个区域（避免频谱泄露）
    mask[cx-15:cx+15, peak_y-15:peak_y+15] = 1
    fft_filtered = fft_shifted * mask
    
    # 逆FFT得到复信号，提取包裹相位
    ifft_filtered = np.fft.ifft2(np.fft.ifftshift(fft_filtered))
    wrapped_phase = np.angle(ifft_filtered)
    return wrapped_phase

5. 执行2D相位解包 #

使用2D维度的解包逻辑，比如借助scipy的unwrap函数指定双轴解包：

from scipy.signal import unwrap

def unwrap_2d(wrapped_phase):
    # 先沿x轴解包，再沿y轴解包，消除二维平面内的相位跳变
    unwrapped = unwrap(wrapped_phase, axis=1)
    unwrapped = unwrap(unwrapped, axis=0)
    return unwrapped

6. 相位转表面高度与误差修正 #

用步骤1的公式将解包后的相位转换为高度，同时需完成：

• 像素到物理尺寸的标定：根据相机与物体的距离、镜头参数，将像素坐标转换为实际物理坐标
• 系统误差消除：去除因干涉仪倾斜、镜头畸变带来的背景相位偏移，可通过测量平面物体的相位作为基准扣除

内容的提问来源于stack exchange，提问作者cqRinO