% Name: calculate_contrast.m
% Created 12/2015
%
% Description: Calculate pattern modulation contrast from a TIFF stack of
%  3 phase images.
%
%  contrast - Ratio of the total signal in the first to the total signal 
%             in the zeroth-order Fourier components [O'Holleran2012]
%
% Laurence Young <ljy24@cam.ac.uk>
clear all

filename = 'img.tif';
nPhases = 3;    % Number of phase images in the stack

imgStack = zeros(width, width, 3);
k0 = zeros(width,width,1);
kp1 = zeros(width,width,1);

for m = 1:nPhases
    imgStack(:,:,m) = double(imread(filename,m));
end

x0=imgStack(:,:,1);
x1=imgStack(:,:,2);
x2=imgStack(:,:,3);

first_order = sqrt((x0-x1).^2 + (x1-x2).^2 + (x2-x0).^2 );
contrast = (sqrt(2)/(width*width)) ./sum(sum(x0+x1+x2)) .* ...
               sum(sum(first_order));
