The code does not return the unwrapped phase. Why?

4 次查看（过去 30 天）
vishav Rattu
vishav Rattu 2016-11-27
评论： vishav Rattu 2016-11-28
if true
% code
phi=zeros(512,512);
I=zeros(512,512);
R=ones(512); %reference beam with all ones
for i=1:512
for k=1:512
phi(i,k)=0.0005*(((i-256)^2)+((k-256)^2));
end
end
m=1.28; % radians per pixel
Y=zeros(512,512);
for i=1:512
for k=1:512
Y(i,k)=m*k;
end
end
I=1+0.5*cos(phi +Y); %ZERO PHASE
J=1+0.5*cos(phi+ Y +(pi/2)); %PI/2 PHASE
K=1+0.5*cos(phi+ Y +(pi)); %PI PHASE
L=1+0.5*cos(phi+ Y +(3*pi/2)); %3PI/2 PHASE
P=(L-J)./(I-K); %in order to calculate the phase
ophi=atan2((P-tan(Y)),(1+P.*tan(Y))); %ophi is the calculated phase which lies between -pi/2 and pi/2
end
%in order to unwrap the phase,I use the Itoh algorithm
%by first sequentially unwrapping the all rows,one at a time.
ophi_unwrapped=zeros(512,512);
i=1;
for i=1:512
ophi_unwrapped(i,:) = unwrap(ophi(i,:));
end
%Then sequentially unwrap all the columns one at a time
i=1;
for i=1:512
ophi_unwrapped(:,i) = unwrap(ophi(:,i));
end
surf(ophi_unwrapped);
  2 个评论
vishav Rattu
vishav Rattu 2016-11-27
Correction: ophi lies between -pi to pi. Also, I should get ophi to be equal to phi as defined above.
vishav Rattu
vishav Rattu 2016-11-28
Made some changes to the program. The 'ophi' value obtained from the method below actually lies between (-pi/2 to pi/2). Now in order to unwrap this ophi, I use the Itoh algorithm again. Please keep in mind that after unwrapping we should get phi which is a paraboloid. Still the unwrapped phase is not as desired The unwrapping algorithm is shown below:
if true
% code
end
I=1+0.5*cos(phi); %ZERO PHASE
J=1+0.5*cos(phi+(pi/2)); %PI/2 PHASE
K=1+0.5*cos(phi+(pi)); %PI PHASE
L=1+0.5*cos(phi+(3*pi/2)); %3PI/2 PHASE
P=(L-J)./(I-K); %in order to calculate the phase
ophi=atan(P);
ophi_unwrapped=zeros(512,512);
k=0; % initialize k to 0
i=1; % initialize the counter to 1
alpha=0.5; % set alpha to the desired Tolerance. The maximum increase in phi is 0.26 radians per pixel
for i = 1:(size(ophi,1)-1)
ophi_unwrapped(i,:)=ophi(i,:)+(pi*k); % add 2*pi*k to ui
if((abs(ophi(i+1,:)-ophi(i,:)))>(abs(alpha))) %if diff is greater than alpha, increment or decrement k
if ophi(i+1,:)<ophi(i,:) % if the phase jump is negative, increment k
k=k+1;
else % if the phase jump is positive, decrement k
k=k-1;
end
end
end
ophi_unwrapped((i+1),:)=ophi(i+1,:)+(pi*k); % add pi*k to the last element of the input
k=0; % initialize k to 0
i=1; % initialize the counter to 1
alpha=0.5; % set alpha to the desired Tolerance.
for i = 1:(size(ophi,2)-1)
ophi_unwrapped(:,i)=ophi_unwrapped(:,i)+(pi*k); % add pi*k to ui
if((abs(ophi_unwrapped(:,i+1)-ophi_unwrapped(:,i)))>(abs(alpha))) %if diff is greater than alpha, increment or decrement k
if ophi_unwrapped(:,i+1)<ophi_unwrapped(:,i) % if the phase jump is negative, increment k
k=k+1;
else % if the phase jump is positive, decrement k
k=k-1;
end
end
end
ophi_unwrapped((i+1),:)=ophi_unwrapped(i+1,:)+(pi*k); % add pi*k to the last element of the input
figure
surf(ophi_unwrapped,'FaceColor','interp', 'EdgeColor','none','FaceLighting','phong')
view(-30,70), camlight left, axis tight
title('Wrapped phase image plotted as a surface')
xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians')

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Daniel kiracofe
Daniel kiracofe 2016-11-27
Not 100% sure I have the whole answer, but here's part of it. Look at these two lines side by side:
%by first sequentially unwrapping the all rows,one at a time.
ophi_unwrapped(i,:) = unwrap(ophi(i,:));
%Then sequentially unwrap all the columns one at a time
ophi_unwrapped(:,i) = unwrap(ophi(:,i));
The second line overwrites the results of the first line. Just blows them away. Its as if the first line had never been executed. Perhaps what you wanted was more like this:
%by first sequentially unwrapping the all rows,one at a time.
ophi_unwrapped_rows(i,:) = unwrap(ophi(i,:));
%Then sequentially unwrap all the columns one at a time
ophi_unwrapped(:,i) = unwrap( ophi_unwrapped_rows(:,i));
That gives you an intermediate variable to store the results of the row unwrapping, which you can then subsequently process for the columns.
  2 个评论
vishav Rattu
vishav Rattu 2016-11-27
I totally forgot that. Thanks mate.
But as you said, it only solves a part of the problem. The ophi_unwrapped phase should be paraboloid and not a tent like structure.
Daniel kiracofe
Daniel kiracofe 2016-11-27
well, there I don't know if I can help you. Presumably you have an error in P or Y or some other variable. I don't know what all of the intermediate things are supposed to be. I would suggest to just plot surf(P) and surf(Y) and surf(phi), etc, etc, and check all of your variables.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Transforms 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by