The code does not return the unwrapped phase. Why?

3 Ansichten (letzte 30 Tage)
vishav Rattu
vishav Rattu am 27 Nov. 2016
Kommentiert: vishav Rattu am 28 Nov. 2016
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 Kommentare
vishav Rattu
vishav Rattu am 27 Nov. 2016
Correction: ophi lies between -pi to pi. Also, I should get ophi to be equal to phi as defined above.
vishav Rattu
vishav Rattu am 28 Nov. 2016
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')

Melden Sie sich an, um zu kommentieren.

Antworten (1)

Daniel kiracofe
Daniel kiracofe am 27 Nov. 2016
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 Kommentare
vishav Rattu
vishav Rattu am 27 Nov. 2016
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 am 27 Nov. 2016
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.

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Graphics Performance finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by