基于Matlab模拟哈特曼波前探测器

        哈特曼波前探测器（Shack-Hartmann Wavefront Sensor, SHWFS）[1]用于测量波前像差，本文介绍如何用Matlab模拟SHWFS工作过程，即从像差中获得光点阵图像、计算光斑质心、波前斜率、重构波前。（暂时不对其原理详细介绍）

1.1 模拟微透镜阵列。

lens      = 20;             % 子透镜个数
lens_res  = 30;             % 子透镜对应的CCD区域的分辨率
res       = lens * lens_res;% CCD的分辨率
xp        = linspace(-1,1,lens_res); 
[X,Y]     = meshgrid(xp,xp);
[rho]     = sqrt(X.^2+Y.^2);
len       = ones(lens_res,lens_res) .* (rho<=1);% 单个子透镜
mask1     = (rho<=1);
valid_num = numel(rho(rho<=1));
xp        = linspace(-1,1,res); 
[X,Y]     = meshgrid(xp,xp);
[rho]     = sqrt(X.^2+Y.^2);
mask      = ones(res,res) .* (rho<=1);% 单个子透镜
Circles   = repmat(len,lens,lens);  % 
aa        = Circles .* mask;
figure(3); subplot(121);imagesc(Circles);colorbar;
subplot(122);imagesc(aa);colorbar;% 圆形口径下的子透镜

得到的子透镜如下，但是右侧有的子透镜还是不完整，需要进一步处理。

加上这行代码后：

for i = 1:lens
    for j = 1:lens
        len_t  = aa((i-1) * lens_res + 1 : i * lens_res ,...
                    (j-1) * lens_res + 1 : j * lens_res);
        if(sum(sum(len_t)) < valid_num)
            aa((i-1) * lens_res + 1 : i * lens_res , ...
                (j-1) * lens_res + 1 : j * lens_res) = 0;
        end
    end
end
figure(3); imagesc(aa);colorbar;

有效的子透镜如下图所示。每个小黄色的圆形即代表子透镜，其元素值均为1.

 1.2 微透镜分割波前

任意选取一组Zernike系数：RealZerCoes 如下

real_coe = [ ...
    -0.0804,-1.5517,1.7531,  -0.0317,-0.1377, 0.1954,0.0167,-0.0031,0.0298,0.0410...
    -0.1640, 0.0846,-0.0254, 0.0734,-0.1656,-0.2289, 0.1348,-0.0130, 0.0624,-0.1272...
    -0.0083, 0.0949, 0.0321, 0.0319, 0.0067,-0.0526,-0.0635,-0.0137,-0.0004, 0.0381...
    -0.0893,-0.0089, 0.0125,-0.0202,-0.0010];
num = numel(real_coe);
Znm = zeros(res,res,num);
for i= 1:num
    Znm(:,:,i) = Zernike(i,res); % 参考第一篇代码
end
WF = reshape(Znm,[res*res,num]) * real_coe';
WF = reshape(WF,[res,res]);
figure(20); imagesc(WF);colorbar;colormap(jet);title('波前数据')
figure(21); subplot(121);imagesc(WF);colorbar;colormap(jet);title('波前,单位[um]')
subplot(122);imagesc(WF .* aa);colorbar;colormap(jet);title('被微透镜分割后的波前')

这里只是从简单入手，如果子透镜个数太少，就会有更多的波前没有被使用而被丢弃，如下图所示。另外，透镜本身的像差也没有考虑在内，本文只是做了简单的介绍而已。

 1.3 某个子透镜形成的PSF图像

使用菲涅尔衍射获得 PSF图像。傅里叶光学中给出了相关公式。

lambda      = 7.80e-07;
pixelSize   = 10e-6;        % 单个像素尺寸10um
k           = 2*pi/lambda;  % 波数   
focalLength = 0.01;         % 子透镜的焦距[m]
i = 10; j = 10; % 选择坐标为第(i,j)的子透镜
% 准备求第(i,j)子透镜前表面对应的PSF图像
WF_i_j      = WF .* aa * (1e-6); % 单位换算到m
WF_i_j      = WF_i_j((i-1) * lens_res + 1 : i * lens_res , ...
                    (j-1) * lens_res + 1 : j * lens_res);
WF_i_j_1 = WF_i_j;
WF_i_j_1((WF_i_j_1==0)) = nan;
figure(22); surf(WF_i_j_1);colorbar;colormap(jet);title('子透镜(i,j)前表面的波前')
SS          = mask1 .* exp(-1i*k .* WF_i_j);
% 菲涅尔衍射方法    
[M,~]   = size(SS);    % 物面复振幅分布占据像素的个数   
L1      = 3e-4;        % CCD单个像素为10um，每个子透镜对应的像素30*30，去大小为300um*300um
dx1     = L1/M;        % 物面采样间距 = CCD的单个像素尺寸[um] 10um
z       = focalLength;
L2      = lambda*z/dx1;% 像面尺寸,0.78um*10mm/10um=780mm
dx2     = lambda*z/L1; % 像面的最小尺寸,0.78um*10mm/300um=26um
x2      = -L2/2:dx2:L2/2-dx2;%  -390nm:26nm:364.4nm, 50等分
[X2,Y2] = meshgrid(x2,x2); 
c       = 1/(1i*lambda*z)*exp(1i*k/(2*z)*(X2.^2+Y2.^2)); % 系数
u2      = c.*ifftshift(fft2(fftshift(SS)))*dx1^2;
outpsf0 = (abs(u2).^2);
%{
    outpsf0 - 30*30,当时假定单个子透镜在CCD上对应的像素数为30*30，且单个像素为10um
    但是，在透镜后焦面上真实的复振幅分布的最小间距(像素)为 dx2=26um，可见最小像素变大了；
    现在要把最小像素由26um，减少到10um，采用的是 imresize()函数；经过该操作后，又出来了新的问题，
    即：像素个数由30*30 变为了78*78，只需要裁剪出中间的30*30区域即可。
%}
outpsf1  = imresize(outpsf0,(L2/length(outpsf0))/pixelSize,'bilinear'); % Scaling the pixel size of the propagated field to the correct one.
figure(5); 
subplot(121);imagesc(outpsf0);colorbar;colormap(jet); title('插值前')         
subplot(122);imagesc(outpsf1);colorbar;colormap(jet); title('插值后')         
      

某个子透镜前表面采集的畸变波前，与其对应的PSF图像如下：

而对于理想的参考波前，被采样后的波前和PSF图像如下：

因为理想的平面波成像要在正中央，此时只要对正中央采集出20*20像素（单个像素大小10微米）即可，并记下其对应的坐标

[midx, midy]= SimulatingSHWFS.CentroidCalculation( outpsf_fater,1, 0);
midx        = round(midx);
midy        = round(midy);
leftx       = midx - lens_res/2+1;
rightx      = midx + lens_res/2;
lefty       = midy - lens_res/2+1;
righty      = midy + lens_res/2;
% 确保参考波前的PSF的质心在正中心
figure(7);imagesc(outpsf_fater(leftx : rightx ,lefty : righty));
colorbar;colormap(jet); title('裁剪后的PSF')   

mid0 = round(SimulatingSHWFS.CentroidCalculation( ...
             outpsf_fater(leftx : rightx ,lefty : righty),1, 0));

 

此时得到右78*78 ==> 20*20 裁剪的左端点 和右端点，对与任意畸变波前裁剪都以这个为标准即可。

% 参考波前下 左端点为32，右端点为51；以这个为基准，对畸变波前形成的PSF裁剪
leftx    = 32;
rightx   = 51;
lefty    = 32;
righty   = 51;
figure(8);imagesc(outpsf_fater(leftx : rightx ,lefty : righty));
colorbar;colormap(jet); title('裁剪后的PSF')          
[midx0, midy0] = (SimulatingSHWFS.CentroidCalculation( ...
             outpsf_fater(leftx : rightx ,lefty : righty),1, 0))

对所有子透镜执行相同的操作，即可得到相应的光点阵图像（这里没有考虑到相机的设置问题，以及子透镜本身所带的光学像差）

下面的代码是模拟SHWFS的全部代码，上面的代码是对其做了很大的简化，看起来比较方便，至于Zernike模式法重构波前不做介绍，需要还请私聊，以及基于LCOS的重构方法都可私聊获得代码；仿真中后者重构精度更高，但是还没有在实验室验证

参考文献：

[1] Platt B C, Shack R. History and principles of Shack-Hartmann wavefront sensing[J]. Journal of Refractive Surgery, 2001, 17(5): S573-S577.