<< imnoise SIP - Scilab Image Processing Toolbox improfile >>

SIP - Scilab Image Processing Toolbox >> SIP - Scilab Image Processing Toolbox > imphase


Image phase calculation.

Calling Sequence

phase_value = imphase(formula_name,thresh,im1,im2,im3[,im4,im5])


a 2D array containing the calculated phase (between (-%pi) and (+%pi)).
a string, the name of the formula to use to calculate the phase. May be
3 images, phase_shift=(%pi/2)
3 images, phase_shift=(2*%pi/3)
3 images, phase_shift=(2*%pi/3)
4 images, phase_shift=(%pi/2)
5 images, phase_shift=(%pi/2)
5 images, phase_shift=(%pi/2)
if the difference for a pixel between 2 images is less than this threshold, than this pixels is considered to be belong to the background. matrix containing an image.
im1 im2...
matrices containing images.


imphase calculates the phased image resulting from the combination of several phase-shifted pictures.

This function is used in profilometry (or in interferometry) to calculate the phase on any point of an image. The elevation of a point is proportionnal to the calculated of this point.

Here is a usual method for proceeding:

1) Fringes are projected on a reference plane, perpendicular to a CCD camera.

2) A 1st image is taken with the CCD camera.

3) Fringes are deplaced from a quarter of an interfringe (in this case, phase shift = %pi/2)

4) A 2nd image is taken. And so on...

5) The phased image(phref) is calculated from these images. Several algorithms exist depending on the phase shift (often (2*%pi/3) or (%pi/2)) and the number of images available.

6) Repeat the same operations with an object instead of the reference plane. Calculate (phobj): the phased image for the object.

7) Unwrap (phref) and (phobj) in (uwphref) and (uwphobj).

8) The altitude map is proportional to (uwphobj)-(uwphref)

Using five images (when it is possible) should lead to better phase calculation.

An alternative method is to substract (phobj-phref) modulo 256 (if all images are 8bit) and only then unwrap the resulting image: in this case, there is only one unwrapping process.

Bellow are the formula used:

'bucket3a': phase shift is 90 degres ; usage: general purpose (vibration environnements) ; phase=arctan((i3-i2)/(i1-i2))

'bucket3b': phase shift is 120 degres ; usage: general purpose (vibration environnements) ; p1=2*(%pi/3);p2=2*p1 ; phase=arctan(((i3-i2)+(i1-i3)*cos(p1)+(i2-i1)*cos(p2))/ ((i1-i3)*sin(p1)+(i2-i1)*sin(p2)))

'bucket3c': phase shift is 120 degres ; usage: general purpose (vibration environnements) ; phase=arctan((sqrt(3)*(i3-i2))/(2*i1-i2-i3))

'bucket4a': phase shift is 90 degres ; usage: general purpose ; phase=arctan((i4-i2)/(i1-i3))

'bucket5a': phase shift is 90 degres ; usage: phase shift correction (Hariharan) ; phase=arctan((2*(i4-i2))/(i1-2*i3+i5))

'bucket5b': phase shift is 90 degres ; usage: enhanced phase shift error correction (Creath/Schmit) ; phase=arctan((i1-4*i2+4*i4-i5)/(i1+2*i2-6*i3+2*i4+i5))

Note about Object mask: If the value of a pixel doesn't vary more than the threshold between two images, than we consider that this pixel is in the background (and not on the object on which we project fringes), and we give it the phase = 0 radian.


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

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) 

p=(xmax/nb)/cos(theta);//distance between two black fringes
//projected on a plane perpendicular to the camera
//(=interfringe seen from the camera)
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
//start_point=56;//the prism begins on column 56
//                                   +relief_line(start_point-1);
//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
start_point=floor(size(x,'c')/2);//the pyramid is centered in the picture
//compute altitude:

//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

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

//phase calculation for the reference plane

//simulate 4 images shifted with (%pi/2) with fringes projected on the object

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

//phase calculation for the object

//display the original object on which we projected fringes
xtitle("original object: relief");

//display one of the images with fringes
xtitle(['one of the original images:';'friges are projected on the object']);

//display the phased image of the object
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:
//unwrap the phase corresponding to the object:

//2nd step:the difference of unwrapped phases is 
//proportionnal to the altitude z
//3rd step: display
xtitle(['phase reconstituted from the wrapped phases';
        'this phase is proportionnal to altitude']);


software FRINGE ANALYSIS by HOLO3 (St Louis, FRANCE)

software INTELLIWAVE by engsynthesis

"Modelisation de forme 3D par methode de moire de projection et analyse par decalage de phases", Cyril Breque and Fabrice Bremand

"Metrologie optique par decalage de phase", Yves Surrel,conservatoire national des arts et metiers

A google search with keywords: phase shifting interferometry or moire or phase unwrapping should lead to good introductory documents on the subject.



The latest version of the Scilab Image Processing toolbox can be found at


See Also

<< imnoise SIP - Scilab Image Processing Toolbox improfile >>