夏克-哈特曼波前传感器(SHWS) Matlab仿真
第一章 夏克-哈特曼波前传感器原理介绍
第二章 Matlab仿真程序
第三章 总结
文章目录
一、 夏克-哈特曼波前传感器原理介绍
使用Shack - Hartmann (SH)方案测量波前的局部斜率。在SH传感器中,一个二维透镜阵列将波前分割,并将单个片段内的光线聚焦到相机上。根据聚焦阵列相对于无像差波前的位移,计算每个波前段的局部斜率,并重建波前。
对于无像差的波前,在相机上形成规则且均匀间隔的点阵列。对于像差波前,每个点的位置根据波前在给定小透镜上的局部倾斜(或斜率)而移动。
二、Matlab仿真
1.Input
function [RecoveredZern,SH,ML]=SHsimulator(modes,ZerValues,resolution,PixelSize,LAMBDA,Lenses,focal,Prop,factor,bits,...
MLGeometry,MLCentroidMethod,MLSharedArea,MLfieldDistortion,MLvignetting,ObjectMagnitude,PupilDiameter,...
Exposuretime,BandWidth,QE,PhotonNoise,WellCapacity,ReadNoise,DarkCurrent,ReadOutNoise,paint)
给SH仿真器输入相关参数,包括模式,zenike系数值、分辨率等参数。
2.Introducing the incoming wavefront in microns
ZSum=zeros(length(ZerModo{1}));
for i=2:SH.modes
ZMode=ZerValues(i)*ZerModo{i};
ZSum=ZSum+ZMode;
end
WF=ZSum.*1e-6;%conversion to meters
%Painting
[SH.pupil4paint]=PaintWavefront(SH,WF);
模拟一个变形波前。
3.Definition of the microlenses array
[ML]=MicroLenses(SH,ML,WF);
WFSubpupil=(WF+ML.AmplitudeMask);%Wavefront in each microlent (phase of
%microlenses itself is not yet taking into account
Lenses=cell(1,length(ML.coor));%Preallocation
for i=1:length(ML.coor)
Lenses{i}=WFSubpupil(ML.coor(i,1):ML.coor(i,2), ML.coor(i,3):ML.coor(i,4));
end
生成透镜阵列,并将波前嵌入。
4.Point Spread Function Calculation
PSF=cell(1,length(Lenses));
PSFmaxLocal=zeros(1,length(Lenses));
if ML.Prop==0 %CASE OF NO PROPAGATION BETWEEN MICROLENSES AND CCD. Edge effects will exist, it should be implemented a pad with zeros to avoid it
fprintf('Propagation between microlenses and CCD is not considered. Microlenses are considered as pure amplitude objects\n')
for i=1:length(Lenses)
PF=ML.Pupil.*exp(-1i*SH.k.*Lenses{i});%Pupil function
PSF{i}=abs(ifftshift(ifft2(fftshift(PF)))).^2;%PSF
PSFmaxLocal(i)=max(max(PSF{i}));
end
else
L=length(Lenses{1})*SH.PixelSize; %Size of the region of each microlenses
for j=1:length(Lenses)
S = ML.Pupil.*exp(-1i*SH.k.*Lenses{j}); %complex phase screen
if isfield(ML,'Aberration') == 1 %If aberration o each microlens itself is considered
aberration = zernike(ML.AberrationZernike,size(S,1));
S = S.*exp(-1i*SH.k*aberration*ML.AberationValue*SH.LAMBDA);
end
PSF{j} = lensletSimulation(S,L,SH.LAMBDA,ML.focal,SH.PixelSize);
if isfield(ML,'fieldDistortion') == 1
PSF{j} = applySeidelDistortion(ML.focal,SH.PixelSize,0.1*SH.LAMBDA,size(S,1),PSF{j});
end
if isfield(ML,'vignetting') == 1
PSF{j} = applyVignetting(ML.focal,SH.PixelSize,size(S,1),1e-3,PSF{j});
end
PSFmaxLocal(j)=max(max(PSF{j}));
end
end
计算PSF。(下图为局部放大)
5.Centroid calculation
Xcentroid=zeros(1,length(PSF));
Ycentroid=Xcentroid;
%If you want to visualize each centroid calculation, set the variable
%'display' =1. It is not in the main config because it is usefull only for
%testing
display=0;
if display==1
figure
end
for i=1:length(PSF)
[Xcentroid(i),Ycentroid(i)]=CentroidCalculation(PSF{i},ML.CentMethod,display);
end
%Coordinates of the real centroids and the reference centroids are calculated.
ML.RefCentroid=ML.coor+ML.radiusPixels;% Coordinates of the reference centroid.
CentroidLength=zeros(length(ML.coor),2);
CoorCentroid=zeros(length(ML.coor),2);
crosstalks=0;
for i=1:length(ML.RefCentroid)
CentroidLength(i,1)=Xcentroid(i)-(ML.radiusPixels);
CentroidLength(i,2)=Ycentroid(i)-(ML.radiusPixels);
CoorCentroid(i,1)=ML.RefCentroid(i,1)+ CentroidLength(i,1);
CoorCentroid(i,2)=ML.RefCentroid(i,3)+CentroidLength(i,2);
end
%The following calculates if there exist double spots or crosstalk between
%microlenses.
if isempty(idealPSF)==1
[IdealXcentroid(i),IdealYcentroid(i)]=CentroidCalculation(idealPSF{i},ML.CentMethod,0);
end
if isempty(idealPSF)==1
for i=1:length(ML.RefCentroid)
CentroidLength(i,1)=IdealXcentroid(i)-(ML.radiusPixels);
CentroidLength(i,2)=IdealYcentroid(i)-(ML.radiusPixels);
CoorCentroid(i,1)=ML.RefCentroid(i,1)+ CentroidLength(i,1);
CoorCentroid(i,2)=ML.RefCentroid(i,3)+CentroidLength(i,2);
if CentroidLength(i,1)>ML.radiusPixels
crosstalks=crosstalks+1;
else
if CentroidLength(i,2)>ML.radiusPixels
crosstalks=crosstalks+1;
end
end
end
end
if SH.paint==1
figure
suptitle('Reference centroid vs calculated centroid')
subplot(1,2,1)
set(gcf,'color','w');
plot(ML.RefCentroid(:,1),ML.RefCentroid(:,3),'o','MarkerSize',2,'MarkerEdgeColor','r')
hold on
plot(CoorCentroid(:,1)',CoorCentroid(:,2)','x','MarkerSize',5,'MarkerEdgeColor','b')
hold on
for i=1:length(PSF)
rectangle('Position', [ML.coor(i,1) ML.coor(i,3) ML.radiusPixels*2 ML.radiusPixels*2],'LineWidth', 0.1, 'EdgeColor', 'b');
%pause(0.01)
end
legend('Reference centroid','calculated centroid')
title('Centroids position')
xlabel({'pixels';['Number of centroids exceeding its microlent area (crosstalk): ' num2str(crosstalks)]});
ylabel('pixels')
set(gca,'ydir','reverse')
xlim([0 SH.resolution])
ylim([0 SH.resolution])
subplot(1,2,2)
quiver(ML.RefCentroid(:,1),ML.RefCentroid(:,3),CentroidLength(:,1),CentroidLength(:,2),0);
title('Displacement of the real centroid')
xlabel('pixels')
ylabel('pixels')
set(gca,'ydir','reverse')
xlim([0 SH.resolution])
ylim([0 SH.resolution])
drawnow();
end
质心计算(包括参考以及真实)以及它们的差。
6.Slope Calculations
delta=double(CentroidLength*SH.PixelSize);
alfax= double(atan(delta(:,1)/ML.focal));
alfay=double(atan(delta(:,2)/ML.focal));
solvingX=ML.erased;
solvingY=ML.erased;
count=1;
for i=1:length(ML.erased)
if ML.erased(i)==1
solvingX(i)=alfax(count);
solvingY(i)=alfay(count);
count=count+1;
end
end
if SH.paint==1
% Building a slopes matrix
solvingXvec=vec2mat(solvingX',sqrt(length(ML.erased)));
solvingYvec=vec2mat(solvingY',sqrt(length(ML.erased)));
solvingX2=solvingXvec;
solvingY2=solvingYvec;
for i=length(solvingXvec):-1:1
for j=length(solvingXvec):-1:1
if solvingX2(i,j)==0
solvingX2(i,j)=NaN;
end
if solvingY2(i,j)==0
solvingY2(i,j)=NaN;
end
end
end
figure
suptitle('Slope or the recovered phase')
subplot(1,2,1)
colormap jet;
colorbar
set(gcf,'color','w');
% imshow(solvingX2,[])
imagesc(solvingX2)
title('x-axis')
colorbar
subplot(1,2,2)
% imshow(solvingY2,[])
imagesc(solvingX2)
colormap jet;
colorbar
title('y-axis')
colorbar
drawnow();
end
%Joining together in a vector
Solving=zeros(length(solvingX)*2,1);
Solving(1:length(solvingX),1)=solvingX;
Solving(length(solvingX)+1:length(solvingY)*2,1)=solvingY;
计算XY方向上的倾斜程度。
7.Rcovering
RecoveredZern = lsqr((CalibratedZernike(:,1:SH.modes)*1e-6)/SH.factor,Solving,1e-10,500);
%Recovered wavefront
RecWF=zeros(length(ZerModo{1}));
for i=4:SH.modes
ZModeRec=RecoveredZern(i)*ZerModo{i};
RecWF=RecWF+ZModeRec;
end
RecWF=RecWF.*1e-6;%conversion to meters.
paintWFs(WF,RecWF,SH.pupil4paint);
校正后的波前图像。
8.Idual analysis
Residual=SH.pupil.*(WF-RecWF);
%Pupil function
PFOrig = SH.pupil.*exp(sqrt(-1)*SH.k.*WF);
PFRec = SH.pupil.*exp(sqrt(-1)*SH.k*RecWF);
PFDif=SH.pupil.*exp(sqrt(-1)*SH.k*(RecWF*0));
PFRes=SH.pupil.*exp(sqrt(-1)*SH.k*Residual);
%PSF
PSFOrig=fft2(PFOrig); %Two-dimensional discrete Fourier Transform.
PSFOrig=fftshift(PSFOrig);
PSFOrig=PSFOrig.*conj(PSFOrig);
PSFRec=fft2(PFRec); %Two-dimensional discrete Fourier Transform.
PSFRec=fftshift(PSFRec);
PSFRec=PSFRec.*conj(PSFRec);
PSFDif=fft2(PFDif); %Two-dimensional discrete Fourier Transform.
PSFDif=fftshift(PSFDif);
PSFDif=PSFDif.*conj(PSFDif);
PSFRes=fft2(PFRes); %Two-dimensional discrete Fourier Transform.
PSFRes=fftshift(PSFRes);
PSFRes=PSFRes.*conj(PSFRes);
%OTF & MTF
OTFOrig=fft2(PSFOrig);%OTF
MTFOrig = abs(OTFOrig);
MTFOrig=fftshift(MTFOrig);
MTFOrig = MTFOrig./max(max(MTFOrig));
OTFRec=fft2(PSFRec);
MTFRec = abs(OTFRec);
MTFRec=fftshift(MTFRec);
MTFRec = MTFRec./max(max(MTFRec));
OTFRes=fft2(PSFRes);
MTFRes = abs(OTFRes);
MTFRes=fftshift(MTFRes);
MTFRes = MTFRes./max(max(MTFRes));
if SH.paint==1
paintPSFs(PSFOrig,PSFRec,PSFRes)
PaintConvolution(PSFOrig,PSFRec,PSFRes)
paintMTF(SH,MTFOrig,MTFRec,MTFRes)
end
PSF对比。
三、总结
以上简要介绍了SH波前传感器的原理和实现,如需要上述全套代码,可以私信我,感谢关注,欢迎一起学习!