%% 下面计算 相位导数方差

% 假设 phase_map 是你的相位图二维矩阵
% K 是窗口的大小
k = 3; % 请使用实际的窗口大小替换

% 计算 x 和 y 方向的偏导
[dx, dy] = gradient(wrappedPhase);
Ksq = k^2; % 计算 K^2


half_k = floor(k / 2);

% 初始化结果矩阵
result = zeros(size(wrappedPhase));

% 镜像填充
padded_dx = padarray(dx, [half_k half_k], 'symmetric');
padded_dy = padarray(dy, [half_k half_k], 'symmetric');

% 遍历图像中的每个像素
for m = 1 : size(wrappedPhase, 1)
    for n = 1 : size(wrappedPhase, 2)
        
        % 选取窗口内的偏导数
        sub_dx = padded_dx(m : m+k-1, n : n+k-1);
        sub_dy = padded_dy(m : m+k-1, n : n+k-1);
        
        % 计算平均值
        mean_dx = mean(sub_dx(:));
        mean_dy = mean(sub_dy(:));
        
        % 计算各个像素偏导数与均值的差
        delta_x = sub_dx - mean_dx;
        delta_y = sub_dy - mean_dy;
        
        % 计算方差并存储到结果矩阵的对应位置
        result(m, n) = sqrt(sum(delta_x(:).^2) + sum(delta_y(:).^2)) / Ksq;
    end
end