The code does not return the unwrapped phase. Why?

7 ビュー (過去 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 月 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 件のコメント
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.

サインインしてコメントする。

Community Treasure Hunt

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

Start Hunting!

Translated by