利用深度学习解包裹相位

本文采用软件为pycharm以及matlab，相关pycharm配置可以看我的另一篇文章，手把手教你配置深度学习环境。

本文通过仿真数据训练一个改进的 U‑Net 网络，将归一化的包裹相位图（范围 [-1,1]）映射为连续的相位图（弧度值），从而消除 2π 跳变，恢复真实相位。

数据生成

仿真方法：随机生成平滑的连续相位图（代表物体形貌），加上噪声后包裹得到归一化的包裹相位图。

关键函数：generate_smooth_phase（二次曲面 + 正弦扰动 + 随机高斯斑）→ add_phase_noise → wrap_phase → 归一化。

训练集：动态生成，每轮训练产生新的样本（20,000 训练 + 2,000 验证），无需存储图片。

模型架构

基础骨架：编码器-解码器结构，跳跃连接保留细节。

增强模块：

残差块（Residual Block）：缓解梯度消失，提升深层网络性能。

通道注意力（CBAM）：自适应重标定特征通道，增强对条纹边缘和相位跳变的敏感度。

输入/输出：单通道灰度图（尺寸 128×128，与训练一致）。

激活函数：ReLU，批归一化。

训练策略

损失函数：MSE + 0.1 × Gradient Loss

MSE 保证像素精度。

Gradient Loss（Sobel 边缘差）强制保留相位梯度，避免过度平滑。

优化器：Adam（lr=0.001），学习率阶梯衰减（StepLR，每 20 轮乘 0.5）。

早停：验证损失连续 7 轮不下降则停止训练，保存最佳模型 best_unwrap.pth。

加速：图像尺寸 128×128，训练样本 8000，num_workers=4，CPU 下无混合精度。
（条件允许的话用gpu加速，图像尺寸什么的话也可以往上走走）

测试流程

加载模型：读取 best_unwrap.pth，设置 eval() 模式。

数据预处理：

从 MATLAB 生成的 .mat 文件读取归一化包裹相位（变量 wrapped_norm）。

若尺寸不为 128×128，用 OpenCV 双线性插值缩放。

推理：input_tensor = torch.from_numpy(phase_norm).unsqueeze(0).unsqueeze(0).float() → 模型前向 → 输出归一化连续相位。

后处理：

用 denormalize_continuous 恢复弧度值（乘以 2π×2）。

若原始尺寸被缩放，再插值回原尺寸。

输出：同时保存为 .npy 和 .mat 文件，变量名 unwrapped_phase。

流程

写好pycharm和matlab代码，然后在pycharm中运行train代码，训练模型，再运行test代码，得到输出的.mat文件，再去matlab中运行代码即可。

代码部分

pycharm中四个代码分别为models.py，test.py，train.py，utils.py。

models.py

import torch
import torch.nn as nn

class ChannelAttention(nn.Module):
    def __init__(self, in_channels, reduction=16):
        super().__init__()
        self.avg_pool = nn.AdaptiveAvgPool2d(1)
        self.max_pool = nn.AdaptiveMaxPool2d(1)
        self.fc = nn.Sequential(
            nn.Linear(in_channels, in_channels // reduction),
            nn.ReLU(inplace=True),
            nn.Linear(in_channels // reduction, in_channels)
        )
        self.sigmoid = nn.Sigmoid()

    def forward(self, x):
        b, c, _, _ = x.size()
        avg_out = self.fc(self.avg_pool(x).view(b, c))
        max_out = self.fc(self.max_pool(x).view(b, c))
        out = self.sigmoid(avg_out + max_out).view(b, c, 1, 1)
        return x * out


class ResidualBlock(nn.Module):
    def __init__(self, in_ch, out_ch):
        super().__init__()
        self.conv1 = nn.Conv2d(in_ch, out_ch, 3, padding=1)
        self.bn1 = nn.BatchNorm2d(out_ch)
        self.relu = nn.ReLU(inplace=True)
        self.conv2 = nn.Conv2d(out_ch, out_ch, 3, padding=1)
        self.bn2 = nn.BatchNorm2d(out_ch)
        if in_ch != out_ch:
            self.shortcut = nn.Conv2d(in_ch, out_ch, 1)
        else:
            self.shortcut = nn.Identity()

    def forward(self, x):
        residual = self.shortcut(x)
        x = self.relu(self.bn1(self.conv1(x)))
        x = self.bn2(self.conv2(x))
        x = self.relu(x + residual)
        return x


class DoubleConv(nn.Module):
    def __init__(self, in_ch, out_ch):
        super().__init__()
        self.conv = nn.Sequential(
            nn.Conv2d(in_ch, out_ch, 3, padding=1),
            nn.BatchNorm2d(out_ch),
            nn.ReLU(inplace=True),
            nn.Conv2d(out_ch, out_ch, 3, padding=1),
            nn.BatchNorm2d(out_ch),
            nn.ReLU(inplace=True),
        )

    def forward(self, x):
        return self.conv(x)


class UNet(nn.Module):
    def __init__(self, in_channels=1, out_channels=1, features=[64, 128, 256, 512], use_cbam=True):
        super().__init__()
        self.use_cbam = use_cbam
        self.encoder = nn.ModuleList()
        self.pool = nn.MaxPool2d(2, 2)

        # 编码器
        for feature in features:
            block = ResidualBlock(in_channels, feature) if use_cbam else DoubleConv(in_channels, feature)
            self.encoder.append(block)
            in_channels = feature
            if use_cbam:
                self.encoder.append(ChannelAttention(feature))

        # 瓶颈
        self.bottleneck = ResidualBlock(features[-1], features[-1]*2) if use_cbam else DoubleConv(features[-1], features[-1]*2)
        if use_cbam:
            self.bottleneck_att = ChannelAttention(features[-1]*2)

        # 解码器
        self.upconvs = nn.ModuleList()
        self.decoders = nn.ModuleList()
        for feature in reversed(features):
            self.upconvs.append(nn.ConvTranspose2d(feature*2, feature, kernel_size=2, stride=2))
            decoder = ResidualBlock(feature*2, feature) if use_cbam else DoubleConv(feature*2, feature)
            self.decoders.append(decoder)
            if use_cbam:
                self.decoders.append(ChannelAttention(feature))

        self.final_conv = nn.Conv2d(features[0], out_channels, kernel_size=1)

    def forward(self, x):
        skip_connections = []
        for i, encoder in enumerate(self.encoder):
            x = encoder(x)
            if i % 2 == 0:   # 只保存残差块后的输出，跳过注意力层
                skip_connections.append(x)
            if i < len(self.encoder)-1 and (i+1) % 2 == 0:
                x = self.pool(x)

        x = self.bottleneck(x)
        if self.use_cbam:
            x = self.bottleneck_att(x)

        for idx in range(len(self.upconvs)):
            x = self.upconvs[idx](x)
            skip = skip_connections[-(idx+1)]
            if x.shape != skip.shape:
                x = nn.functional.interpolate(x, size=skip.shape[2:], mode='bilinear', align_corners=False)
            x = torch.cat((skip, x), dim=1)
            x = self.decoders[idx*2](x)          # 残差块
            if self.use_cbam:
                x = self.decoders[idx*2+1](x)    # 注意力层

        return self.final_conv(x)

test.py

import torch
import numpy as np
import cv2
import os
from pathlib import Path
from tqdm import tqdm
import scipy.io as sio
from models import UNet
from utils import denormalize_continuous

#  配置 
MODEL_PATH = "best_unwrap.pth"
OUTPUT_DIR = "./unwrapped_results"
DEVICE = torch.device("cuda" if torch.cuda.is_available() else "cpu")
TARGET_SIZE = (128, 128)        # 与训练尺寸一致
USE_CBAM = True                 # 与训练时一致

def main_single():
    # 加载模型
    model = UNet(in_channels=1, out_channels=1, use_cbam=USE_CBAM).to(DEVICE)
    model.load_state_dict(torch.load(MODEL_PATH, map_location=DEVICE))
    model.eval()

    # 加载归一化包裹相位
    mat_data = sio.loadmat(r"D:\current\6.10\wrapped_phase_norm.mat")   # 修改为实际路径
    phase_norm = mat_data['wrapped_norm'].astype(np.float32)

    h, w = phase_norm.shape
    # 调整尺寸到网络输入大小
    if (h, w) != TARGET_SIZE:
        phase_norm = cv2.resize(phase_norm, TARGET_SIZE, interpolation=cv2.INTER_LINEAR)

    # 转换为张量并确保数据类型为 float32
    input_tensor = torch.from_numpy(phase_norm).unsqueeze(0).unsqueeze(0).to(DEVICE).float()

    with torch.no_grad():
        output_norm = model(input_tensor).squeeze().cpu().numpy()
    continuous = denormalize_continuous(output_norm)  

    # 恢复原始尺寸
    if (h, w) != TARGET_SIZE:
        continuous = cv2.resize(continuous, (w, h), interpolation=cv2.INTER_LINEAR)

    os.makedirs(OUTPUT_DIR, exist_ok=True)

    # 保存为 .npy 文件
    npy_path = os.path.join(OUTPUT_DIR, "unwrapped_phase.npy")
    np.save(npy_path, continuous)
    print(f"Unwrapped phase (npy) saved to {npy_path}")

    # 同时保存为 .mat 文件，方便 MATLAB 读取
    mat_path = os.path.join(OUTPUT_DIR, "unwrapped_phase.mat")
    sio.savemat(mat_path, {"unwrapped_phase": continuous})
    print(f"Unwrapped phase (mat) saved to {mat_path}")

def main_batch():
    os.makedirs(OUTPUT_DIR, exist_ok=True)
    model = UNet(in_channels=1, out_channels=1, use_cbam=USE_CBAM).to(DEVICE)
    model.load_state_dict(torch.load(MODEL_PATH, map_location=DEVICE))
    model.eval()

    input_folder = "./wrapped_phases"   # 存放 .mat 文件的文件夹
    mat_files = list(Path(input_folder).glob("*.mat"))
    for mat_path in tqdm(mat_files):
        mat = sio.loadmat(str(mat_path))
        # 变量名统一为 'wrapped_norm'，如果不是就调整
        phase_norm = mat['wrapped_norm'].astype(np.float32)
        h, w = phase_norm.shape
        if (h, w) != TARGET_SIZE:
            phase_norm = cv2.resize(phase_norm, TARGET_SIZE, interpolation=cv2.INTER_LINEAR)
        input_tensor = torch.from_numpy(phase_norm).unsqueeze(0).unsqueeze(0).to(DEVICE).float()
        with torch.no_grad():
            output_norm = model(input_tensor).squeeze().cpu().numpy()
        continuous = denormalize_continuous(output_norm)
        if (h, w) != TARGET_SIZE:
            continuous = cv2.resize(continuous, (w, h), interpolation=cv2.INTER_LINEAR)

        # 保存为 .npy
        npy_out = Path(OUTPUT_DIR) / mat_path.with_suffix('.npy').name
        np.save(str(npy_out), continuous)

        # 同时保存为 .mat
        mat_out = Path(OUTPUT_DIR) / mat_path.with_suffix('.mat').name
        sio.savemat(str(mat_out), {"unwrapped_phase": continuous})

    print(f"All results saved to {OUTPUT_DIR}")

if __name__ == "__main__":
    main_single()   

train.py

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import Dataset, DataLoader
import numpy as np
from tqdm import tqdm
import os
from models import UNet
from utils import generate_pair

#  CPU 优化配置 
DEVICE = torch.device("cuda" if torch.cuda.is_available() else "cpu")
BATCH_SIZE = 16                    
EPOCHS = 100
LEARNING_RATE = 0.001
TRAIN_SAMPLES = 8000              
VAL_SAMPLES = 1000                
IMAGE_SIZE = (128, 128)           
MODEL_SAVE_PATH = "best_unwrap.pth"
PATIENCE = 7
USE_CBAM = True
NUM_WORKERS = 2                   

#  动态数据集 
class WrappedPhaseDataset(Dataset):
    def __init__(self, num_samples, shape=IMAGE_SIZE):
        self.num_samples = num_samples
        self.shape = shape

    def __len__(self):
        return self.num_samples

    def __getitem__(self, idx):
        phase_scale = np.random.uniform(4, 12)
        noise_std = np.random.uniform(0.05, 0.2)
        wrapped, continuous = generate_pair(shape=self.shape, phase_scale=phase_scale, noise_std=noise_std)
        return (torch.from_numpy(wrapped).unsqueeze(0).float(),
                torch.from_numpy(continuous).unsqueeze(0).float())

#  梯度损失 
def gradient_loss(pred, target):
    sobel_x = torch.tensor([[-1,0,1],[-2,0,2],[-1,0,1]], dtype=torch.float32).view(1,1,3,3).to(pred.device)
    sobel_y = torch.tensor([[-1,-2,-1],[0,0,0],[1,2,1]], dtype=torch.float32).view(1,1,3,3).to(pred.device)
    pred_dx = nn.functional.conv2d(pred, sobel_x, padding=1)
    pred_dy = nn.functional.conv2d(pred, sobel_y, padding=1)
    target_dx = nn.functional.conv2d(target, sobel_x, padding=1)
    target_dy = nn.functional.conv2d(target, sobel_y, padding=1)
    return nn.functional.l1_loss(pred_dx, target_dx) + nn.functional.l1_loss(pred_dy, target_dy)

#  早停 
class EarlyStopping:
    def __init__(self, patience=7, min_delta=0.0001):
        self.patience = patience
        self.min_delta = min_delta
        self.counter = 0
        self.best_loss = None
        self.early_stop = False

    def __call__(self, val_loss):
        if self.best_loss is None:
            self.best_loss = val_loss
        elif val_loss > self.best_loss - self.min_delta:
            self.counter += 1
            if self.counter >= self.patience:
                self.early_stop = True
        else:
            self.best_loss = val_loss
            self.counter = 0
        return self.early_stop

#  训练 
def main():
    train_dataset = WrappedPhaseDataset(TRAIN_SAMPLES, IMAGE_SIZE)
    val_dataset = WrappedPhaseDataset(VAL_SAMPLES, IMAGE_SIZE)
    train_loader = DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True, num_workers=NUM_WORKERS, pin_memory=False)
    val_loader = DataLoader(val_dataset, batch_size=BATCH_SIZE, shuffle=False, num_workers=NUM_WORKERS, pin_memory=False)

    model = UNet(in_channels=1, out_channels=1, use_cbam=USE_CBAM).to(DEVICE)
    criterion = nn.MSELoss()
    optimizer = optim.Adam(model.parameters(), lr=LEARNING_RATE)
    scheduler = optim.lr_scheduler.StepLR(optimizer, step_size=20, gamma=0.5)
    early_stopping = EarlyStopping(patience=PATIENCE)

    best_loss = float('inf')
    print(f"Training on {DEVICE}, train samples: {TRAIN_SAMPLES}, image size: {IMAGE_SIZE}, use_cbam={USE_CBAM}")

    for epoch in range(1, EPOCHS+1):
        model.train()
        train_loss = 0.0
        pbar = tqdm(train_loader, desc=f"Epoch {epoch}")
        for wrapped, cont in pbar:
            wrapped, cont = wrapped.to(DEVICE), cont.to(DEVICE)
            optimizer.zero_grad()
            output = model(wrapped)
            loss_mse = criterion(output, cont)
            loss_grad = gradient_loss(output, cont)
            loss = loss_mse + 0.1 * loss_grad
            loss.backward()
            optimizer.step()
            train_loss += loss.item()
            pbar.set_postfix({'Loss': loss.item()})
        avg_train_loss = train_loss / len(train_loader)

        # 验证
        model.eval()
        val_loss = 0.0
        with torch.no_grad():
            for wrapped, cont in val_loader:
                wrapped, cont = wrapped.to(DEVICE), cont.to(DEVICE)
                output = model(wrapped)
                loss = criterion(output, cont)
                val_loss += loss.item()
        avg_val_loss = val_loss / len(val_loader)
        print(f"Epoch {epoch} | Train Loss: {avg_train_loss:.6f} | Val Loss: {avg_val_loss:.6f}")

        if avg_val_loss < best_loss:
            best_loss = avg_val_loss
            torch.save(model.state_dict(), MODEL_SAVE_PATH)
            print(f"  -> Best model saved (loss: {best_loss:.6f})")

        if early_stopping(avg_val_loss):
            print(f"Early stopping triggered after {epoch} epochs. Best val loss: {best_loss:.6f}")
            break

        scheduler.step()

    print(f"Training finished. Best model saved to {MODEL_SAVE_PATH}")

if __name__ == "__main__":
    main()

utils.py

import numpy as np
import torch

def generate_smooth_phase(shape=(256,256), scale=8.0):
    """生成平滑连续相位图，模拟物体形貌（弧度）"""
    h, w = shape
    x = np.linspace(-1, 1, w)
    y = np.linspace(-1, 1, h)
    X, Y = np.meshgrid(x, y)
    # 二次曲面基础 + 正弦波纹 + 随机高斯斑
    phase = scale * (0.5 * (X**2 + Y**2))
    phase += scale * 0.3 * np.sin(2*np.pi*X*2) * np.cos(2*np.pi*Y*2)
    if np.random.rand() > 0.5:
        H = np.random.randn(h, w) * 0.3
        phase += scale * 0.1 * H
    return phase.astype(np.float32)


def wrap_phase(phase):
    """将相位包裹到 [-pi, pi]"""
    return ((phase + np.pi) % (2*np.pi)) - np.pi


def add_phase_noise(phase, std=0.1):
    """添加高斯噪声（弧度）"""
    noise = np.random.normal(0, std, phase.shape)
    return phase + noise


def generate_pair(shape=(256,256), phase_scale=8.0, noise_std=0.1):
    """生成一对 (包裹相位, 连续相位) 归一化后的数据"""
    continuous = generate_smooth_phase(shape, phase_scale)
    noisy = add_phase_noise(continuous, noise_std)
    wrapped = wrap_phase(noisy)
    # 连续相位归一化：假设相位范围在 [-2π*2, 2π*2] 左右
    continuous_norm = continuous / (2 * np.pi * 2)
    # 包裹相位归一化到 [-1,1]
    wrapped_norm = wrapped / np.pi
    return wrapped_norm, continuous_norm


def denormalize_continuous(norm_cont):
    """将网络输出的归一化连续相位恢复为真实弧度值"""
    return norm_cont * (2 * np.pi * 2)

MATLAB代码

matlab部分代码主要是为了看深度学习效果

%% 对比深度学习解包裹效果 + 定量分析
clc; clear; close all;

%%  文件路径
wrapped_file = 'D:\current\6.10\wrapped_phase_norm.mat';  
dl_file      = 'D:\pycharm_project\xiangweijiebaoguo\unwrapped_results\unwrapped_phase.mat';

%%  加载包裹相位
if ~exist(wrapped_file, 'file')
    error('包裹相位文件不存在: %s', wrapped_file);
end
tmp = load(wrapped_file);
var_names = fieldnames(tmp);
wrapped_var = [];
for i = 1:numel(var_names)
    if contains(var_names{i}, 'wrapped')
        wrapped_var = var_names{i};
        break;
    end
end
if isempty(wrapped_var)
    error('找不到包裹相位变量（应包含 "wrapped" 字段）');
end
wrapped_phase = tmp.(wrapped_var);
fprintf('加载包裹相位变量: %s\n', wrapped_var);
% 如果是归一化数据，转换为弧度
if max(wrapped_phase(:)) <= 1 && min(wrapped_phase(:)) >= -1
    wrapped_phase = wrapped_phase * pi;
    fprintf('检测到归一化包裹相位，已乘以 π 恢复弧度值。\n');
end

%%  加载解包裹相位 
if ~exist(dl_file, 'file')
    error('解包裹相位文件不存在: %s', dl_file);
end
tmp2 = load(dl_file);
var_names2 = fieldnames(tmp2);
unwrapped_var = [];
for i = 1:numel(var_names2)
    if contains(var_names2{i}, 'unwrapped')
        unwrapped_var = var_names2{i};
        break;
    end
end
if isempty(unwrapped_var)
    % 如果没找到包含'unwrapped'的变量，则取第一个变量
    unwrapped_var = var_names2{1};
    fprintf('未找到明确解包裹变量，默认使用第一个变量: %s\n', unwrapped_var);
end
unwrapped_phase = tmp2.(unwrapped_var);
fprintf('加载解包裹相位变量: %s\n', unwrapped_var);

% 尺寸要一致
if ~isequal(size(wrapped_phase), size(unwrapped_phase))
    error('两个相位图尺寸不一致，请检查！');
end
[rows, cols] = size(unwrapped_phase);

%%  1. 相位范围与整体统计
wrap_min = min(wrapped_phase(:));
wrap_max = max(wrapped_phase(:));
unwrap_min = min(unwrapped_phase(:));
unwrap_max = max(unwrapped_phase(:));
unwrap_range = unwrap_max - unwrap_min;
fprintf('\n========== 相位统计 ==========\n');
fprintf('包裹相位范围: [%.2f, %.2f] rad\n', wrap_min, wrap_max);
fprintf('解包裹相位范围: [%.2f, %.2f] rad, 总跨度 = %.2f rad\n', unwrap_min, unwrap_max, unwrap_range);

%%  2. 相邻像素跳变检测 
dx = abs(diff(unwrapped_phase, 1, 2));
dy = abs(diff(unwrapped_phase, 1, 1));
jump_thresh = pi * 0.8;
total_pairs = rows*(cols-1) + (rows-1)*cols;
jump_ratio = (sum(dx(:) > jump_thresh) + sum(dy(:) > jump_thresh)) / total_pairs;
fprintf('\n========== 连续性指标 ==========\n');
fprintf('相邻像素跳变 > %.1f rad 的比例: %.4f %%\n', jump_thresh, jump_ratio*100);
if jump_ratio < 0.01
    fprintf('✓ 无明显 2π 跳变残留，解包裹连续性好。\n');
elseif jump_ratio < 0.05
    fprintf('⚠ 存在少量跳变（%.2f%%），可能有局部包裹错误。\n', jump_ratio*100);
else
    fprintf('✗ 大量跳变，解包裹失败或存在严重误差。\n');
end

%%  3. 表面粗糙度（RMS）和 平滑度 
rms_unwrap = sqrt(mean(unwrapped_phase(:).^2));
diff_h = diff(unwrapped_phase, 1, 2);
diff_v = diff(unwrapped_phase, 1, 1);
smoothness = std([diff_h(:); diff_v(:)]);
fprintf('\n========== 平滑度指标 ==========\n');
fprintf('解包裹相位的 RMS 值: %.3f rad\n', rms_unwrap);
fprintf('相邻像素差值的标准差: %.4f rad (越小越平滑)\n', smoothness);
if smoothness < 0.2
    fprintf('✓ 相位非常平滑，噪声抑制良好。\n');
elseif smoothness < 0.5
    fprintf('○ 相位有一定波动，可接受。\n');
else
    fprintf('✗ 相位粗糙，可能仍有噪声或解包裹错误。\n');
end

%%  4. 梯度分析 
[Gx, Gy] = gradient(unwrapped_phase);
grad_mag = sqrt(Gx.^2 + Gy.^2);
max_grad = max(grad_mag(:));
fprintf('\n========== 梯度分析 ==========\n');
fprintf('梯度最大值: %.4f rad/像素\n', max_grad);
if max_grad < 0.3
    fprintf('✓ 全局梯度较小，相位变化平缓，形貌平滑。\n');
elseif max_grad < 1.0
    fprintf('○ 局部存在较陡峭变化（可能是真实边缘）。\n');
else
    fprintf('⚠ 出现尖锐跳变，可能存在过冲或解包裹错误。\n');
end

%%  5. 自洽性检查（重新包裹误差） 
wrapped_diff = mod(unwrapped_phase + pi, 2*pi) - pi;
error_wrap = wrapped_diff - wrapped_phase;
max_err = max(abs(error_wrap(:)));
fprintf('\n========== 自洽性检查 ==========\n');
fprintf('重新包裹后的相位与原包裹相位最大差异: %.4f rad\n', max_err);
if max_err < 0.1
    fprintf('✓ 解包裹结果与原始包裹相位自洽（无漂移）。\n');
else
    fprintf('⚠ 存在系统性偏差（可能常数偏移或尺度问题）。\n');
end

%%  最终结论 
fprintf('\n========== 综合结论 ==========\n');
if jump_ratio < 0.02 && smoothness < 0.3 && max_grad < 0.5
    fprintf('深度学习解包裹效果良好，相位连续平滑，无明显误差，可用于后续形貌还原。\n');
elseif jump_ratio < 0.1 && smoothness < 0.8
    fprintf('解包裹效果基本可用，但局部存在轻微缺陷，建议人工检查边缘区域。\n');
else
    fprintf('解包裹效果不理想，请检查训练数据或模型参数。\n');
end

%%  绘制图形
% 2D 并排显示
figure('Name', 'Phase Comparison');
subplot(1,2,1); imagesc(wrapped_phase, [-pi pi]); colormap(jet); colorbar; title('Wrapped Phase'); axis image off;
subplot(1,2,2); imagesc(unwrapped_phase, [unwrap_min unwrap_max]); colormap(jet); colorbar; title('DL Unwrapped'); axis image off;

% 重新包裹误差图
figure; imagesc(error_wrap, [-0.5 0.5]); colormap(jet); colorbar; title('Re-wrap Error (rad)'); axis image off;

% 中心行剖面
center_row = round(rows/2);
figure; plot(1:cols, wrapped_phase(center_row,:), 'r-', 1:cols, unwrapped_phase(center_row,:), 'b-');
grid on; legend('Wrapped', 'DL Unwrapped'); title(sprintf('Profile row %d', center_row));
xlabel('Pixel'); ylabel('Phase (rad)');

输出结果

这个是我的一组数据输入后的输出结果，可供参考。

加载包裹相位变量: wrapped_norm
检测到归一化包裹相位，已乘以 π 恢复弧度值。
加载解包裹相位变量: unwrapped_phase

========== 相位统计 ==========
包裹相位范围: [-3.14, 3.14] rad
解包裹相位范围: [0.74, 9.62] rad, 总跨度 = 8.88 rad

========== 连续性指标 ==========
相邻像素跳变 > 2.5 rad 的比例: 0.0000 %
✓ 无明显 2π 跳变残留，解包裹连续性好。

========== 平滑度指标 ==========
解包裹相位的 RMS 值: 5.397 rad
相邻像素差值的标准差: 0.0205 rad (越小越平滑)
✓ 相位非常平滑，噪声抑制良好。

========== 梯度分析 ==========
梯度最大值: 0.1261 rad/像素
✓ 全局梯度较小，相位变化平缓，形貌平滑。

========== 自洽性检查 ==========
重新包裹后的相位与原包裹相位最大差异: 6.2790 rad
⚠ 存在系统性偏差（可能常数偏移或尺度问题）。

========== 综合结论 ==========
深度学习解包裹效果良好，相位连续平滑，无明显误差，可用于后续形貌还原。