clear
close all
clc

% 'Dammann_numerical.mat' contains the transition points data in [Zhou1995]
% [Zhou1995] C. Zhou and L. Liu, "Numerical study of Dammann array 
% illuminators," Applied Optics, vol. 34, p. 5961, September 1995. 
load('Dammann_Numerical_Data.mat');
Dammann_numerical=Dammann_Numerical_Data.Dammann_numerical;
fprintf('Numerical data loaded.\n');

%% BEGIN of Dammann Grating Generation
fprintf('Generating Dammann Grating...\n');
X_Array_Size=6;
Y_Array_Size=6;

Filename_Prefix=sprintf('Dammann_%d-by-%d_',X_Array_Size,Y_Array_Size);

X_Transition_Points=Dammann_numerical{table2array(cell2table(Dammann_numerical(:,1)))==X_Array_Size,2};
Size_X_Transition_Points=size(X_Transition_Points,2);
Y_Transition_Points=Dammann_numerical{table2array(cell2table(Dammann_numerical(:,1)))==Y_Array_Size,2};
Size_Y_Transition_Points=size(Y_Transition_Points,2);


Dammann_Pixel_X_Size=512;
Dammann_Pixel_Y_Size=512;
Dammann_Gap=128;

Light_X_Size=ceil(1*Dammann_Pixel_X_Size);
Light_Y_Size=ceil(1*Dammann_Pixel_Y_Size);
Light_Sigma=0.5;

Dammann_1D_X=zeros(Dammann_Pixel_X_Size,Dammann_Pixel_Y_Size);
Dammann_1D_X_FINAL=zeros(Dammann_Pixel_X_Size,Dammann_Pixel_Y_Size);
Dammann_1D_Y=zeros(Dammann_Pixel_X_Size,Dammann_Pixel_Y_Size);
Dammann_1D_Y_FINAL=zeros(Dammann_Pixel_X_Size,Dammann_Pixel_Y_Size);

for i=1:1:Dammann_Pixel_X_Size
    Dammann_1D_X(:,i)=mod(i-1,Dammann_Gap);
end
for i=1:1:Dammann_Pixel_Y_Size
    Dammann_1D_Y(i,:)=mod(i-1,Dammann_Gap);
end
Dammann_1D_X=uint16(Dammann_1D_X*65536/Dammann_Gap);
Dammann_1D_Y=uint16(Dammann_1D_Y*65536/Dammann_Gap);

for ix=1:1:Size_X_Transition_Points
    bin_X(:,:,ix)=imbinarize(Dammann_1D_X,X_Transition_Points(ix));
    if mod(ix,2)
        Dammann_1D_X_FINAL=logical(Dammann_1D_X_FINAL+bin_X(:,:,ix));
    else
        Dammann_1D_X_FINAL=logical(Dammann_1D_X_FINAL-bin_X(:,:,ix));
    end
end

for iy=1:1:Size_Y_Transition_Points
    bin_Y(:,:,iy)=imbinarize(Dammann_1D_Y,Y_Transition_Points(iy));
    if mod(iy,2)
        Dammann_1D_Y_FINAL=logical(Dammann_1D_Y_FINAL+bin_Y(:,:,iy));
    else
        Dammann_1D_Y_FINAL=logical(Dammann_1D_Y_FINAL-bin_Y(:,:,iy));
    end
end

Dammann_2D=xor(Dammann_1D_X_FINAL,Dammann_1D_Y_FINAL);
Dammann_2D=double(Dammann_2D*65535);
figure
imshow(Dammann_2D);
fprintf('\tDammann Grating Generated.\n');
% END of Dammann Grating Generation


%% START Image Export
fprintf('Exporting images...\n');
%  Dammann Grating
Dammann_2D_Image=uint16(Dammann_2D);
Dammann_2D_Image_Filename=sprintf('%sPattern_Original_PX%d_PY%d_G%d',Filename_Prefix,Dammann_Pixel_X_Size,Dammann_Pixel_Y_Size,Dammann_Gap);
imwrite(Dammann_2D_Image,[pwd '\Dammann Grating\' Dammann_2D_Image_Filename '.png']);
fprintf('\tExporting finished.\n');

%% START Image Show
fprintf('Image Show...\n');
figure('Name',Dammann_2D_Image_Filename,'NumberTitle','off');
imshow(Dammann_2D_Image);
fprintf('\t-----Complete-----\n');
