%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Unwrapping phase based on Ghiglia and Romero (1994) based on weighted and unweighted least-square method
% URL: https://doi.org/10.1364/JOSAA.11.000107
% Inputs:
% * psi: wrapped phase from -pi to pi
% * weight: weight of the phase (optional, default: all ones)
% Output:
% * phi: unwrapped phase from the weighted (or unweighted) least-square phase unwrapping
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function phi = phase_unwrap(psi, weight)
if (nargin < 2) % unweighted phase unwrap
% get the wrapped differences of the wrapped values
dx = [zeros([size(psi,1),1]), wrapToPi(diff(psi, 1, 2)), zeros([size(psi,1),1])];
dy = [zeros([1,size(psi,2)]); wrapToPi(diff(psi, 1, 1)); zeros([1,size(psi,2)])];
rho = diff(dx, 1, 2) + diff(dy, 1, 1);
% get the result by solving the poisson equation
phi = solvePoisson(rho);
else % weighted phase unwrap
% check if the weight has the same size as psi
if (~all(size(weight) == size(psi)))
error('Argument error: Size of the weight must be the same as size of the wrapped phase');
end
% vector b in the paper (eq 15) is dx and dy
dx = [wrapToPi(diff(psi, 1, 2)), zeros([size(psi,1),1])];
dy = [wrapToPi(diff(psi, 1, 1)); zeros([1,size(psi,2)])];
% multiply the vector b by weight square (W^T * W)
WW = weight .* weight;
WWdx = WW .* dx;
WWdy = WW .* dy;
% applying A^T to WWdx and WWdy is like obtaining rho in the unweighted case
WWdx2 = [zeros([size(psi,1),1]), WWdx];
WWdy2 = [zeros([1,size(psi,2)]); WWdy];
rk = diff(WWdx2, 1, 2) + diff(WWdy2, 1, 1);
normR0 = norm(rk(:));
% start the iteration
eps = 1e-8;
k = 0;
phi = zeros(size(psi));
while (~all(rk == 0))
zk = solvePoisson(rk);
k = k + 1;
if (k == 1) pk = zk;
else
betak = sum(sum(rk .* zk)) / sum(sum(rkprev .* zkprev));
pk = zk + betak * pk;
end
% save the current value as the previous values
rkprev = rk;
zkprev = zk;
% perform one scalar and two vectors update
Qpk = applyQ(pk, WW);
alphak = sum(sum(rk .* zk)) / sum(sum(pk .* Qpk));
phi = phi + alphak * pk;
rk = rk - alphak * Qpk;
% check the stopping conditions
if ((k >= numel(psi)) || (norm(rk(:)) < eps * normR0)) break; end;
end
end
end
function phi = solvePoisson(rho)
% solve the poisson equation using dct
dctRho = dct2(rho);
[N, M] = size(rho);
[I, J] = meshgrid([0:M-1], [0:N-1]);
dctPhi = dctRho ./ 2 ./ (cos(pi*I/M) + cos(pi*J/N) - 2);
dctPhi(1,1) = 0; % handling the inf/nan value
% now invert to get the result
phi = idct2(dctPhi);
end
% apply the transformation (A^T)(W^T)(W)(A) to 2D matrix
function Qp = applyQ(p, WW)
% apply (A)
dx = [diff(p, 1, 2), zeros([size(p,1),1])];
dy = [diff(p, 1, 1); zeros([1,size(p,2)])];
% apply (W^T)(W)
WWdx = WW .* dx;
WWdy = WW .* dy;
% apply (A^T)
WWdx2 = [zeros([size(p,1),1]), WWdx];
WWdy2 = [zeros([1,size(p,2)]); WWdy];
Qp = diff(WWdx2,1,2) + diff(WWdy2,1,1);
end
没有合适的资源?快使用搜索试试~ 我知道了~
基于加权最小二乘算法实现解包裹附matlab代码
共6个文件
png:4个
m:2个
1.该资源内容由用户上传,如若侵权请联系客服进行举报
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
版权申诉
0 下载量 41 浏览量
2022-10-26
19:02:40
上传
评论 2
收藏 935KB ZIP 举报
温馨提示
1.版本:matlab2019a,不会运行可私信 2.领域:【解包裹】 3.内容:基于加权最小二乘算法实现解包裹附matlab代码 4.适合人群:本科,硕士等教研学习使用
资源推荐
资源详情
资源评论
收起资源包目录
【解包裹】基于加权最小二乘算法实现解包裹附matlab代码.zip (6个子文件)
11.png 227KB
1.png 242KB
2.png 225KB
test_phase_unwrap.m 2KB
3.png 243KB
phase_unwrap.m 4KB
共 6 条
- 1
资源评论
天天Matlab科研工作室
- 粉丝: 3w+
- 资源: 7258
下载权益
C知道特权
VIP文章
课程特权
开通VIP
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功