%Purpose: reconstruction the complex amplitude of the object wave from an off-axis hologram
%Inputs: %'I_R', 'I_O', the recorded diffraction pattern of structured illumination;
              %'Ratio',resolution/(2pixels);
              %'w0', the pre-set carrier-frequency of structured illumination; 
%Outputs: %'Re_O',Reconstructed object wave;
                %'Carrier', w0x+1i*w0y; carrier-frequency of the illumination; 
                
close all; clear all; clc; %Initializing;

Path_exper_folder='G:\Experiment\2021_12_20_DHM'; %Destination Folder of hologram;
Path_function_folder='D:\matlab\Work\DHM\DHM_reconstruction'; %Functions folder;

Wave=561*1e-9;  %wavelength (m);
k=2*pi/Wave;       %wave vector;
Pixel=1.85*1e-6;   %pixel size(m);

%Read DH hologram: 
[File_name,Path_name] = uigetfile(strcat(Path_exper_folder,'\','*.TIFF'), 'Hologram'); 
FileToRead=strcat(Path_name,File_name), %file information of DH hologram; 
I_Holo=imread(FileToRead); I_Holo=double(I_Holo); %load DH hologram; 

%Coordinates of object space: 
[y_length,x_length]=size(I_Holo); %The size of hologram; 
x=linspace(-x_length*Pixel/2,x_length*Pixel/2,x_length);  %x coordinate(micro); 
y=linspace(-y_length*Pixel/2,y_length*Pixel/2,y_length); %y coordinate(micro); 
[xx,yy]=meshgrid(x,y); %Two-dimensional coordinates; 
figure(1); imagesc(x*1e6,y*1e6,I_Holo); colormap('bone'); xlabel('x (\mum)'); ylabel('y (\mum)'); colorbar; 

%Coordinates of frequency space: 
u=linspace(-1/(2*Pixel),1/(2*Pixel),x_length);  %Spectrum coordinate; 
v=linspace(-1/(2*Pixel),1/(2*Pixel),y_length);  %Range: 0~2pi;  
[uu,vv]=meshgrid(u,v); %Generate the two-dimensional coordinates; 
 
%Digital reference wave: 
cd (Path_function_folder); 
[w0x,w0y]=Carrier_frequency_detection(I_Holo,Pixel); %Extract the carrier-frequencies; 
Sim_R=exp(-1i*2*pi*(w0x*xx+w0y*yy)); %Digital reference wave; 
Freq_O=fftshift(fft2(fftshift(I_Holo.*Sim_R))); %Selected spectrum; 
 
Ratio=0.2; 
Filter_mask=double( sqrt(uu.^2+vv.^2)<max(u)*Ratio ); %Filtering mask to select the +1st spectrum; 
Re_O=fftshift(ifft2(fftshift(Freq_O.*Filter_mask))); %Reconstructed complex amplitude; 
 
figure(2); subplot(1,2,1); imagesc(x*1e6,y*1e6, abs(Re_O));   xlabel('x (\mum)'); ylabel('y (\mum)'); colormap('bone');  colorbar; title('Reconstructed amplitude (A.U.)');
               subplot(1,2,2); imagesc(x*1e6,y*1e6, angle(Re_O)); xlabel('x (\mum)'); ylabel('y (\mum)'); colorbar; 
               title('Wrapped phase (rad)'); 
 
Pha=Phase_unwrapping_volkovt(angle(Re_O) ); %Phase unwrapping; 
 
figure(3); subplot(1,2,1); imagesc(x*1e6,y*1e6, abs(Re_O)); colormap('bone'); xlabel('x (\mum)'); ylabel('y (\mum)'); colorbar;  title('Reconstructed amplitude (A.U.)');
               subplot(1,2,2); imagesc(x*1e6,y*1e6,Pha); xlabel('x (\mum)'); ylabel('y (\mum)'); colorbar; 
               title('Unwrapped phase (rad)');