The code does not return the unwrapped phase. Why?

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);

[vishav Rattu]

Correction: ophi lies between -pi to pi. Also, I should get ophi to be equal to phi as defined above.

[vishav Rattu]

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')

Answer 1

[Daniel kiracofe]

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.

[vishav Rattu]

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]

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.