% This program uses G-S algorithm to realize the design of beam shaping
% device
%% Initiation
clear;
clc;
N=512;
M=512;
f=500;
lamda=560e-6;
L0=15;
r=5;
dx0=L0/N;
dy0=L0/M;
L=lamda*f*N/L0;
dx=L/N;
dy=L/M;
% - N and M are the number of horizontal and vertical sampling points,
% respectively
% - f is the focal length, and its nuit is micrometre
% - lamda is the wave length, and its unit is micrometre
% - L0 is the size of binary optical element, and its unit is micrometre
% - r is the clear radius of the binary optical element, and its unit is
% micrometre
% - dx0 and dy0 are the horizontal and vertical sample interval in binary
% optical element, respectively, and their unit are micrometre
% - L is the size of image plane, and its unit is micrometre
% - dx and dy are the horizontal and vertical sample interval in image
% plane
%% Creating the orginal BOE palne and the ideal image plane
[m,n]=meshgrid(linspace(-N/2,N/2-1,N),linspace(-M/2,M/2-1,M));
I0=zeros(N,M);
R=sqrt((dx0*m).^2+(dy0*n).^2);
index1=find(R<=r);
I0(index1)=1;
energy0=sum(sum(I0));
chemin='.\image\';
[nom,chemin]=uigetfile([chemin,'*.*'],['Input the picture'],100,100);
[I,MAP]=imread([chemin,nom]);
I=rgb2gray(I);
I=imresize(I,[1024 1024]);
I=im2double(~imbinarize(I));
% The image sampler
I_i=zeros(N,M);
for i=1:1024/N:1024
for j=1:1024/M:1024
I_i((i-1)*N/1024+1,(j-1)*M/1024+1)=I(i,j);
end
end
% energy conservation
energy=sum(sum(I_i));
I_i=I_i*energy0/energy;
%% G-S algorithm
[I0_phase,err]=GSA_Q(I0,I_i,8,100);
%% Actual result of imaging and calculating the rms and ESEU
I1=abs(fftshift(fft2(sqrt(I0).*exp(1i*I0_phase)))/sqrt(M*N)).^2;
I1_u=I1/max(max(I1));
index2=find(I_i>0.8);
ave_I1=sum(I1(index2))/(M*N);
rms=1/(length(index2)-1)*sqrt(sum(((I1(index2)-ave_I1)/ave_I1).^2));
eta=sum(I1(index2))/sum(sum(I1));
%% Ploting the result
figure(1);
mesh(I_i);
title('The intensity distribution of ideal image plane','FontName','Times New Roman','FontSize',32)
figure(2);
mesh(I1);
title('The intensity distribution of actual image plane','FontName','Times New Roman','FontSize',32)
figure(3);
bar3(I0_phase);
title('The phase distribution','FontName','Times New Roman','FontSize',32)
figure(4);
imshow(I_i);
title('The ideal image plane','FontName','Times New Roman','FontSize',32)
figure(5)
imshow(I1_u);
title('The actual image plane','FontName','Times New Roman','FontSize',32)
figure(6)
imshow(I0_phase);
title('The phase of the binanry optical element','FontName','Times New Roman','FontSize',32)
figure(7)
stairs(I0_phase(N/2,:));
title('The phase of the binanry optical element (Central section)','FontName','Times New Roman','FontSize',32)
xlabel('\itPeriods','FontName','Times New Roman','FontSize',20);
ylabel('\itPhase/rad','FontName','Times New Roman','FontSize',20);
axis([1 M -4 4]);
figure(8)
plot(err,'markersize',2);
xlabel('The number of iterations','FontName','Times New Roman','FontSize',20);
ylabel('error (F-norm)','FontName','Times New Roman','FontSize',20);
axis([1 200 0 inf]);
title('The relationship between error and the number of iterations','FontName','Times New Roman','FontSize',32);
评论0
最新资源