imphase

As all this may not be easy for people unused to these
concepts of phase shifting and phase unwrapping,
here is a very detailled (and long) example:
all the steps (1 to 8) above are illustrated.
I found more useful for this example to show how to simulate 
several objects

//begin
stacksize(4e7);//much memory needed to treat pictures

nb=6;//number of black fringes on the reference plane
step=0.02;//definition of the image
theta=%pi/6;//fringes are projected on the object
//with this angle (compared to the vertical)
//(nb:this script doesn't consider shadows which can occur in real cases) 

xmax=5.3;
x=0:step:xmax;
p=(xmax/nb)/cos(theta);//distance between two black fringes
//projected on a plane perpendicular to the camera
//(=interfringe seen from the camera)
y=ones(x);
plan=(x'*y)';//matrix defining the reference plane

//Now we choose an altitude for each point of the plane:
//Below is the simulation of 3 common objects:
//a gaussian, a wedge and a pyramid.
//just uncomment the one you want (here: pyramid)

//Simulation of a gaussian object:
//size_gauss=1.5;//parameter to have a gaussian object large or narrow
//amplitude=5;//maximal height of the object
//relief_line=exp((-(x-2.5).^2)/size_gauss);//these functions can be adapted to 
//relief_column=exp((-(x-2.5).^2)/size_gauss);// simulate any object you want
//z=amplitude*relief_column'*relief_line;//compute altitude

//Simulation of a prism (or a wedge): because of the dicontinuities 
//the unwrapping process will make an error in shape reconstitution
//relief_line=ones(x); 
//start_point=56;//the prism begins on column 56
//relief_line(start_point:$)=2*(x(start_point:$)-x(start_point))...
//                                   +relief_line(start_point-1);
//relief_column=ones(x);
//z=relief_column'*relief_line;//compute altitude
//z(1:54,:)=ones(z(1:54,:));//the prism is only between
//z(100:$,:)=ones(z(100:$,:));//lines 55 to 99

//Simulation of a pyramid
relief_line=ones(x); 
start_point=floor(size(x,'c')/2);//the pyramid is centered in the picture
relief_line(1:start_point)=2*(x(1:start_point)); 
relief_line(start_point:$)=-2*(x(start_point:$)-x(start_point))...
                                    +relief_line(start_point-1); 
relief_column=ones(x);
relief_column(1:start_point)=2*(x(1:start_point)); 
relief_column(start_point:$)=-2*(x(start_point:$)-x(start_point))...
                                    +relief_column(start_point-1); 
//compute altitude:
z=relief_column'*ones(relief_line)+ones(relief_column)'*relief_line;

//phase calculation
phref=2*%pi*(plan/p);//phase of reference (phase of the points of
//the reference plane)
phobj=2*%pi*(plan+z*tan(theta))/p;//phase of the points of the object

//simulate 4 images shifted with (%pi/2) for the reference plane
//and calculate the phased image
phi=phref;
image=zeros(phi);

for a=1:4
dephi=(a-1)*%pi/2;//phase shift=(%pi/2)
lum=127.5*(1+cos(phi+dephi));//calculate luminance of the image
image(:,:,a)=lum;
end

//phase calculation for the reference plane
ph1=imphase('bucket4a',0,image(:,:,1),image(:,:,2),image(:,:,3),image(:,:,4));

//simulate 4 images shifted with (%pi/2) with fringes projected on the object
phi=phobj;
image=zeros(phi);

for a=1:4
dephi=(a-1)*%pi/2;//phase shift=(%pi/2)
lum=127.5*(1+cos(phi+dephi));//calculate luminance of the image
image(:,:,a)=lum;
end

//phase calculation for the object
ph2=imphase('bucket4a',0,image(:,:,1),image(:,:,2),image(:,:,3),image(:,:,4));

//display the original object on which we projected fringes
xset("window",0);clf;
plot3d1(1:8:size(z,'r'),1:8:size(z,'c'),z(1:8:$,1:8:$));
xtitle("original object: relief");

//display one of the images with fringes
xset("window",1);clf;imshow(image(:,:,1)/max(image(:,:,1));
xtitle(['one of the original images:';'friges are projected on the object']);

//display the phased image of the object
xset("window",2);clf;
imshow((ph2+%pi)*(1/(2*%pi)));xtitle('image of phase for the object');

//Now we verify if what we've done is correct:  
//we exploit the phased images to find again the profile of the object:
//1st step is phase unwrapping:
disp('phase unwrapping: be a bit patient');
//unwrap the phase corresponding to the referenc plane:
uwphref=unwrapl((ph1+%pi)*(127.5/%pi));
//unwrap the phase corresponding to the object:
uwphobj=unwrapl((ph2+%pi)*(127.5/%pi));

//2nd step:the difference of unwrapped phases is 
//proportionnal to the altitude z
uwpht=uwphobj-uwphref;
//3rd step: display
xset("window",3);clf;plot3d1(1:size(uwpht,'r'),1:size(uwpht,'c'),uwpht);
xtitle(['phase reconstituted from the wrapped phases';
        'this phase is proportionnal to altitude']);
    
//end