-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathReto_fase (5).m
108 lines (93 loc) · 2.96 KB
/
Reto_fase (5).m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
clear;
clc;
I = double((imread('mask_1.bmp')));
I_1 = double(rgb2gray(imread('DSC_0005.JPG')));
I_2 = double(rgb2gray(imread('DSC_0006.JPG')));
I_3 = double(rgb2gray(imread('DSC_0007.JPG')));
I_4 = double(rgb2gray(imread('DSC_0008.JPG')));
I_5 = double(rgb2gray(imread('DSC_0009.JPG')));
I_6 = double(rgb2gray(imread('DSC_0010.JPG')));
I_7 = double(rgb2gray(imread('DSC_0011.JPG')));
I_8 = double(rgb2gray(imread('DSC_0012.JPG')));
I_1=I.*I_1;
I_2=I.*I_2;
I_3=I.*I_3;
I_4=I.*I_4;
I_5=I.*I_5;
I_6=I.*I_6;
I_7=I.*I_7;
I_8=I.*I_8;
Phi1= atan2(-(I_2-I_4),(I_1-I_3));
Phi2= atan2(-(I_6-I_8),(I_5-I_7));
imshow(Phi1)
figure(1);
F1= unwrap_TIE(Phi1);
surf(F1)
shading interp, colormap gray
figure(2);
F2 = unwrap_TIE(Phi2);
surf(F2)
shading interp, colormap gray
figure(3);
x=linspace(0,13,3072);
y=linspace(0,15,4608);
[X,Y] = meshgrid(y,x);
Fase_Desnv_1= F2-F1;
Fase_Desnv=Fase_Desnv_1;
surf(X,Y,Fase_Desnv)
view(2), shading interp, colormap gray
xlabel('X (cm)')
ylabel('Y (cm)')
figure(4)
x=linspace(0,13,3072);
y=linspace(0,15,4608);
[X,Y] = meshgrid(y,x);
z=-((Fase_Desnv./2*pi).*(0.5)./tand(31));
z1=z./9;
Z=z1-min(min(z1));
surf(X,Y,Z);
shading interp, colormap gray
xlabel('X (cm)')
ylabel('Y (cm)')
zlabel('Z (cm)')
function [phase_unwrap,N]=Unwrap_TIE_DCT_Iter(phase_wrap)
phi1 = unwrap_TIE(phase_wrap);
phi1=phi1+mean2(phase_wrap)-mean2(phi1); %adjust piston
K1=round((phi1-phase_wrap)/2/pi); %calculate integer K
phase_unwrap=phase_wrap+2*K1*pi;
residue=wrapToPi(phase_unwrap-phi1);
phi1=phi1+unwrap_TIE(residue);
phi1=phi1+mean2(phase_wrap)-mean2(phi1); %adjust piston
K2=round((phi1-phase_wrap)/2/pi); %calculate integer K
phase_unwrap=phase_wrap+2*K2*pi;
residue=wrapToPi(phase_unwrap-phi1);
N=0;
while sum(sum(abs(K2-K1)))>0
K1=K2;
phic=unwrap_TIE(residue);
phi1=phi1+phic;
phi1=phi1+mean2(phase_wrap)-mean2(phi1); %adjust piston
K2=round((phi1-phase_wrap)/2/pi); %calculate integer K
phase_unwrap=phase_wrap+2*K2*pi;
residue=wrapToPi(phase_unwrap-phi1);
N=N+1;
end
end
function [phase_unwrap]=unwrap_TIE(phase_wrap)
psi=exp(1i*phase_wrap);
edx = [zeros([size(psi,1),1]), wrapToPi(diff(psi, 1, 2)), zeros([size(psi,1),1])];
edy = [zeros([1,size(psi,2)]); wrapToPi(diff(psi, 1, 1)); zeros([1,size(psi,2)])];
lap = diff(edx, 1, 2) + diff(edy, 1, 1); %calculate Laplacian using the finite difference
rho=imag(conj(psi).*lap); % calculate right hand side of Eq.(4) in the manuscript
phase_unwrap = solvePoisson(rho);
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