我用matlab产生了100100个随机点图案（如图，其中有若干亮点），现在我想要每一个点对应一个像素，输出为100100像素的图案。（因为我要用Gerchberg-Saxton算法恢复相位，每一个像素对应一个相位值，需要100100个相位，所以像素值也需要100100）。
程序：
clear
clc
N=101; %101为光屏的范围，即Figure中-50到50
lda=633e-9;dx=300e-9; %波长及单元结构中心间隔
x0=(-N/2+1/2)*dx:dx:(N/2-1/2)*dx;y0=x0; %超表面位置坐标,即x0=(最小值:间隔:最大值)
x1sin=linspace(-1/dx/2,1/dx/2,N)*lda;y1sin=x1sin;
% 对应论文中公式5-5和5-6,linspace(a,b,N)表示生成N个等间隔点，这些点的间距为(b-a)/(N-1),即横向空间频率；
% x1sin与y1sin用于计算衍射角，即矢量与面yoz、xoz夹角
[xx1sin,yy1sin]=meshgrid(x1sin,y1sin);
% 基于向量x1sin和向量y1xin中包含的坐标返回二维网络坐标；
% 即形成一个length(x1sin)和length(y1sin)的方形网格坐标。
trans=find((xx1sin.^2+yy1sin.^2)<1); %寻找行波
xx1_k=xx1sin(trans);yy1_k=yy1sin(trans); %求出行波对应的横向空间频率
kz=sqrt(1-xx1_k.^2-yy1_k.^2); %求出行波对应的横向空间频率（与z轴）
tic %计时开始
r=zeros(N);k=0;ks=0;
% 创建一个NxN全0矩阵，k为设置为亮的衍射级次的个数(实时随机点数),ks(可能)为统计的循环次数(没啥意义)
Nb=2022;%目标随机点数（圆圈内的），即行波中目标衍射级次的数量
while k<Nb %产生随机点
s1=unidrnd(N); %产生一组从1到N的随机数，相当于论文中mx
s2=unidrnd(N); %相当于论文中my
if (xx1sin(s1,s2).^2+yy1sin(s1,s2).^2)<1 %找出行波
r(s1,s2)=1;r(N+1-s1,N+1-s2)=1; %将衍射级次设置为亮，且满足中心对称性
k=sum(r(:)); %统计设置为亮的衍射级次的个数，相当于论文中nb
end
ks=ks+1;
end
toc %计时结束
k %在命令行窗口输出k值
ks %在命令行窗口输出ks值
rs=r(trans);
for i=1:101
for j=1:101
E2(5*i-4:5*i,5*j-4:5*j)=r(i,j);
end
end
figure(1),imagesc(linspace(-50,50,505),linspace(-50,50,505),E2);
% 绘制结果图
% imagesc(x,y,C)指定图像位置,使用x和y可指定与C(1,1)和C(m,n)对应的边角的位置。
axis square %将当前坐标系图形设置为方形
axis([-50,50,-50,50]) %设置坐标刻度范围
wordsize=20; %x,y轴文本字体大小为20
xlabel('Orders in x direction','fontsize',wordsize) %x轴说明
ylabel('Orders in y direction','fontsize',wordsize) %y轴说明
set(gca,'FontSize',20); %设置x，y轴上数字字体尺寸为20磅
colormap(gray) %设置图像的颜色为灰度渐进
saveas(gcf,'picture.png');
这是我要用的GS算法
clear all
I=imread('Picture.png');
image=double(I(:,:,1));% 将目标图像读进来
[M,N]=size(image);
phase=zeros(N);
A=fftshift(fft2(fftshift(image)));
% 先对带有线段的图像进行傅里叶变换
Source=ones(N);% 假设入射的为平面光
for n=1:200 %这个for语句里就是GS方法的步骤
B=Source.*exp(i*angle(A));
C=fftshift(fft2(fftshift(B)));
C=rot90(rot90(C));
D=image.*exp(i*angle(C) );
A=fftshift(fft2(fftshift(D)));
end
holo=angle(A);% 这个holo就是超表面所需要的相位
% 下面是再现步骤
Amp=abs(A)/max(max(abs(A)));
E=fftshift(fft2(fftshift(exp(i*holo))));
E=rot90(rot90(E));
Recon=abs(E).^2;
figure(2),imagesc(Recon);colormap(hot);axis equal;axis tight;title('模拟全息图'); %输出全息图
figure(3),imagesc(holo);colormap(hot);title('相位分布图'); %输出相位分布图
```
