%This program implements fringe pattern demodulation using the continuous 
%2D fan continuous wavelet transform algorithm 
%and employs the phase estimation method.

%Install YAWTB toolbox to run this program. See the link --
%http://rhea.tele.ucl.ac.be/yawtb/
%This program is written By Dr. Abdulbasit Abid and Dr. Munther Gdeisat 
%at The General Engineering %Research Institute at 
%Liverpool John Moores University, Liverpool (UK) on 8th Sept. 2009.
%Email m.a.gdeisat@ljmu.ac.uk

% Please see the paper entitled "Spatial carrier fringe pattern
% demodulation by use of a two-dimensional continuous wavelet transform" 
% Applied Optics, Vol. 45, Issue 34, pp. 8722-8732 doi:10.1364/AO.45.008722
% By:Munther Gdeisat, David Burton and Michael Lalor

clear; close all; clc;

%construct a simulated object
S=256;
x=-S+1:1:S;
[X, Y] = meshgrid(x,x);
aa = size(X);
noise = randn(aa(1), aa(2));
shape = peaks(2*S);
figure(1)
mesh(shape)

%construct a simulated fringe pattern
fo = 1./16;
beta = 0.5;
fringes = cos(2.*pi.*fo.*X + beta*shape)+ 1*shape + 0.*noise;
figure(2)
colormap(gray(256))
imagesc(fringes)

%This is the best method to remove the background illumination of a fringe
%pattern. Use it when it is required to do so. 
%kernel = [-1 -2 -1; 0 0 0; 1 2 1]';
%fringes = filter2(kernel, fringes);

%2D Morlet wavelet transform
%look at the fringes in a fringe pattern. Can you determine the number of
%pixels in a fringe. You can use colormap(gray(256)) and imagesc(fringes)
%commands in addition to the zoom capbility in Matlab to do this
intial_scale = 2; %minimum number of pixels in a fringe - 2
scale_step = 2;
final_scale = 20; %maximum number of pixels in a fringe + 5
scales = intial_scale:scale_step:final_scale;

initial_angle = 0;
angle_step = pi/6;
final_angle = pi/2;
angles = initial_angle:angle_step:final_angle; % angles in radians

%2D-CWT
wav = cwt2d(fft2(fringes), 'fan', scales, angles); 

%extract the maximum ridge
sizecoeff = size(wav.data);
for i=1:sizecoeff(1)
    for j=1:sizecoeff(2)
        a = reshape(wav.data(i,j,:,:), sizecoeff(3),  sizecoeff(4));
        b(i,j) = max(max(a));
    end
end
figure(3)
colormap(gray(256))
imagesc(angle(b))

%extract the wrapped phase map
wrapped_phase = angle(b);

%unwrap the phase using a simple phase unwrapper and
%remove the spatial carrier fo from the extracted phase
unwrapped_phase = unwrap(unwrap(wrapped_phase)')'- 2.*pi.*fo.*X;
figure(4)
colormap(gray(256))
imagesc(unwrapped_phase)
figure(5)
mesh(unwrapped_phase)