夏克-哈特曼波前传感器(SHWS) Matlab仿真

夏克-哈特曼波前传感器(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波前传感器的原理和实现，如需要上述全套代码，可以私信我，感谢关注，欢迎一起学习！