##?GNU Octave, version 4.4.1
##Copyright (C) 2018 John W. Eaton and others.
##This is free software; see the source code for copying conditions.
##There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
##FITNESS FOR A PARTICULAR PURPOSE. For details, type 'warranty'.
##
##Octave was configured for "x86_64-w64-mingw32".
##
##Additional information about Octave is available at https://www.octave.org.
##
##Please contribute if you find this software useful.
##For more information, visit https://www.octave.org/get-involved.html
##
##Read https://www.octave.org/bugs.html to learn how to submit bug reports.
##For information about changes from previous versions, type 'news'.
% It is assumed that the images of the intensity patterns before and after the
% mirror was tilted are recorded and saved in a computer as image files jpg, png, tif, etc.
%% Reconstruction code
R=double(imread('File name'));% Read image0 (before tilt) and convert to matrix with floating points
R1=R(:,:,1);% Extract one color channel
[m,n]=size(R1);% Read matrix size
R2=zeros(m*2,n*2);%Create a matrix size which is twice that of the initial matrix
R2(m/2+1:m+m/2,n/2+1:n+n/2)=R1;%Zero padding of R1 by placing at the center of R2
Ry=double(imread('File name'));% Read image1 (after tilt) and convert to matrix with floating points
Ry1=Ry(:,:,1);% Extract one color channel
Ry2=zeros(m*2,n*2);
Ry2(m/2+1:m+m/2,n/2+1:n+n/2)=Ry1;
% Matched filter
Recon_MF=ifftshift(ifft2(fft2(Ry2).*conj((fft2(R2)))));% Reconstruct the complex image
figure % Open a new figure window
imagesc(abs(Recon_MF));% Display the absolute of the reconstructed image
% Phase-only filter
Recon_POF=ifftshift(ifft2(fft2(Ry2).*conj(exp(1i*arg(fft2(R2))))));% Reconstruct the complex image
figure % Open a new figure window
imagesc(abs(Recon_POF));% Display the absolute of the reconstructed image
% Non-linear filter
R11=fftshift(fft2(R2));% Calculate the Fourier transform
R22=fftshift(fft2(Ry2));% Calculate the Fourier transform
alpha=-0.3;% o
beta=0.8;% r
Recon_NLR=ifftshift(ifft2(((abs(R11).^alpha).*exp(1i*arg(R11))).*conj((abs(R22).^beta).*exp(1i*arg(R22)))));% Reconstruct the complex image
Recon_NLR=flipud(fliplr(Recon_NLR)); % Reorient the image
figure % Open a new figure window
imagesc(abs(Recon_NLR))% Display the absolute of the reconstructed image