Phase change/phase shift in leaf surface moisture

Question

I am writing a phase retrieval algorithm in Matlab, I have few different images that represent different milliliter of water on the leaf surface( leaf surface moisture) that needed to be inserted into the algorithm to see the phase change/phase shift that happened in each millimeter of water on the leaf surface, but I am not too sure how to see the phase change/phase shift of each image and compare it. Do I represent the phase change/phase shift in the graph or there are other methods to see the phase change/phase shift? Do I plot the phase or get a value? May I get some helps with this? Thank you.

EDIT: The topic I'm doing now is interferometry, and it produces fringes on each of the images. Each of these fringes will contain the phase values/information, what my code suppose to do is to retrieve the phase information/ values of each image and then I need to see the phase change/phase shift that happened on each image, but I'm not sure how I can do this and I fail to do so. Based on the algorithm, I need to prove that as the moisture level increase/decrease, the phase value will increase/decrease as well/ phase value will shift towards how may degree. Now I don't know am I suppose to see the phase change/phase shift of each image via a graph or other methods. I have attached an image to show the example of input images.Figure: Interferogram for different ml of water I hope that I can get some help with this. Thank you.

EDIT 2: image of a sample graph

Here is my code:

close all; clc; workspace; 
Iterations = 1100;       %Number of iterations 
p = 0.1;                 %Time to pause between display of images 
tic;                     %Timer Initialization
error = [];
%This whole part is to generate a mesh plot, i still dont know what info
%can get from here but just leave it first. 
x = linspace(-10,10,256);
y = linspace(-10,10,256);
[X,Y] = meshgrid(x,y);
x0 = 0;              % Center
y0 = 0;              % Center
sigma = 2;       % Beam Waistfringes
A = 1;               % Peak of the Beam 
%These 2 lines, res and Source, is basically input intensity. U need
%this in order to get the exponential image to process ltr.
res = ((X-x0).^2 + (Y-y0).^2)./(2*sigma^2);  
Source = A  .* exp(-res); 
%These 2 lines is for mesh plot also, ignore first
surf(Source);
shading interp 
%Reading Fringe Image (Target)
Target = imread('0ml(1).jpeg');
Target = rgb2gray(Target);
Target = im2double(Target);
%Get the size of Target
[rowsimgA, colsimgA, numberOfColorChannelsimgA] = size(Source);
[rowsimgB, colsimgB, numberOfColorChannelsimgB] = size(Target);
%imresize function is basically to resize any image to a particular size.
%In this code u can see that i resize the Source beam input to the same
%size as the target image in the format of (Target, [rows, columns])
if rowsimgB ~= rowsimgA || colsimgA ~= colsimgB
    Target = imresize(Target, [rowsimgA colsimgA]); 
end
%FFT process
A = fftshift(ifft2(fftshift(Target)));
%Initiate loop
% The loop im not too sure yet, but basically when u multiplay the source
% and the phase, u produce exponential image, and then u do fft and ifft in
% order to reconstruct the image.
for i = 1:Iterations
  B = abs(Source) .* exp(1i*angle(A)); 
  C = fftshift(fft2(fftshift(B)));
  D = abs(Target) .* exp(1i*angle(C));
  A = fftshift(ifft2(fftshift(D)));
  error = [error; (1/sqrt(rowsimgA*colsimgA)*sqrt(sum(sum(abs(C)-abs(Target)).^2)))];  
end
%Display Outputs
figure(2), imagesc(Target), colorbar, title('Original Fringe')
figure(3), i = 1:1:i; plot(i,(error')); title('Error');
%Phase Mask
figure(4), imagesc(angle(A)), title('Phase Mask');
%Last Pattern
figure(5), imagesc(abs(C)), colorbar, title('Reconstructed Image');
%Reconstructed Image (Phase Only)
figure(6), imagesc(angle(D)), colorbar, title('Reconstructed Image (Phase Only)');
%Reconstructed Image (Amplitude Only)
figure(7), imagesc(abs(D)), colorbar, title('Reconstructed Image (Amplitude Only)');
%Display Total Time Taken for Phase Retrieval
toc; 

`

Answer 1

I feel that there may be a few aspects that need clarification before addressing your question.

First, it seems that you are using the Gerchberg-Saxton algorithm which is indeed the most common phase retrieval algorithm. In its original form, this was designed to dress a specific problem that is likely very similar to your problem: you want to know the magnitude and phase of a field, but can only measure the magnitude. This is common, for example in optics, where cameras only measure intensity. The solution, as proposed by G&S, is to measure the magnitude, then apply some known transform to the field, for example a Fourier transform using a lens, and then measure the magnitude again. Turns out that if you know the magnitude of your field before and after the transform, and you know the transform itself, you can infer the phase using the G-S algorithm. I think this is what you are after.

It seems like you have based yourself on the Wikipedia page for the G-S algorithm. This is tackling a slightly different problem often found in holography. In that case, you have a Gaussian laser beam, for which you can modulate the phase in any way you please but you have no control of the amplitude. You then Fourier transform this field, e.g. using a lens, and you want to create a particular intensity pattern afterwards. Turns out that the G-S algorithm allows you to do this too, but the form of the algorithm is slightly different. This is what is presented on the Wikpedia page.

As far as you're concerned, you don't want to assume a Gaussian input field. Instead, you want to use a measured field corresponding to before the transform, and a different measured field corresponding to after the transform. It's not clear from your question whether you actually have both of these. If you don't, you may need to change your optical setup. Also, I believe this will only work if you illuminate with coherent light.

Finally, to address your actual question, you should consider what the phase means. In your case a change in phase is caused by a longer path-length through water than through air. Hence, you are measuring the `thickness' of the water. If the thickness of the water changes spatially over your image, you should plot an image of the phase, and this is indeed what you do here:

imagesc(angle(A))

Keep in mind, the phase wraps around, so if you have a thickness of water that is more than a wavelength, you will get fringes and it may be tricky to infer the actual water thickness.

EDIT (In response to OP's edit):

This is a somewhat different problem again, and I don't think the Gerchberg-Saxton algorithm will help you here. Instead, I recommend looking into off-axis holography. This involves interfering your object beam and reference beam at an angle, which you seem to have done already (?).

To calculate the phase you have to Fourier transform, spatially filter, and inverse Fourier transform back. Popoff provides an accessible explanation for this on his website: https://www.wavefrontshaping.net/post/id/12, and there are some good papers on this too, for example an overview given by Verrier and Atlan in Applied Optics 2011.