%This program implements fringe pattern demodulation using the continuous 
%2D Paul 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

%you need to download the paul2d.m function from our website to run this program

clear; close all; clc;
S=256;
x=-S+1:1:S;
[X, Y] = meshgrid(x,x);
shape = peaks(2*S);
figure(1)
mesh(shape)

fo = 1./16;
fringes = cos(2.*pi.*fo.*X + 3.*shape);
figure(2)
colormap(gray(256))
imagesc(fringes)

wav = cwt2d(fft2(fringes), 'paul', 2:2:20, 0:pi/6:pi/2, 4); 
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))

wrapped_phase = angle(b);
unwrapped_phase = unwrap(unwrap(wrapped_phase)')'+2.*pi.*fo.*X;
figure(4)
colormap(gray(256))
imagesc(-unwrapped_phase)

figure(5)
mesh(-unwrapped_phase)