哈特曼波前探测器（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:lensfor j = 1:lenslen_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;endend
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:numZnm(:,:,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.

基于Matlab模拟哈特曼波前探测器相关推荐

模拟哈特曼波前探测器
这是模拟哈特曼波前探测器的另一份代码,这里不需要设置微透镜的参数. 波前 ==> 光点,使用的是离散DFT对应的fft(,),而不是 FFT对应的 fft2(),最大的优点是运算速度极快,其原 ...
基于哈特曼波前传感器光斑阵列图像直接斜率求解算法
一.背景介绍夏克哈特曼波前传感器可获取光斑阵列图像,通常情况下,通过计算每个光斑的质心以获取斜率信息.此外,也可以通过对光斑阵列图像进行傅里叶变换来获取斜率信息. 当然,也可直接对光斑阵列图像直接进 ...
一种针对夏克哈特曼波前传感器质心数据求解波前斜率的处理方法
一.导出质心数据针对夏克哈特曼波前传感器(型号:索雷博)导出的质心数据"Save Centroid Date",本文提供一种基于质心数据的斜率矩阵获取及波前重构方法. 图 1 哈 ...
哈特曼波前传感器区域法重构算法实例
一.引言目前哈特曼波前传感器波前重建的算法主要有两种,即区域波前重建法和模式波前重建法. 而区域法重构按求解方法又可分为快速傅里叶变换法,迭代求解法,矩阵向量法等. 本文将分析一下迭代法求解法重构波 ...
【图像重建】基于matlab布雷格曼迭代算法集合ART算法CT图像重建【含Matlab源码 1905期】
一.获取代码方式获取代码方式1: 完整代码已上传我的资源:[图像重建]基于matlab布雷格曼迭代算法集合ART算法CT图像重建[含Matlab源码 1905期] 获取代码方式2: 通过订阅紫极神光 ...
基于Matlab模拟、检测和跟踪飞机着陆进场中异常的仿真（附源码）
目录一.介绍二.生成和标记轨迹三.定义方案四.运行方案并检测异常轨道五.将跟踪异常报告与事实进行比较六.总结七.程序该示例显示了如何自动检测最终接近机场跑道的飞机的偏差和异常.在此示例 ...
【物理应用】基于matlab模拟井筒多相流【含Matlab源码 2152期】
⛄一.获取代码方式获取代码方式1: 完整代码已上传我的资源:[物理应用]基于matlab模拟井筒多相流[含Matlab源码 2152期] 点击上面蓝色字体,直接付费下载,即可. 获取代码方式2: 付 ...
【光学】基于matlab模拟拉盖尔高斯【含Matlab源码 2167期】
⛄一.获取代码方式获取代码方式1: 完整代码已上传我的资源:[光学]基于matlab模拟拉盖尔高斯[含Matlab源码 2167期] 点击上面蓝色字体,直接付费下载,即可. 获取代码方式2: 付费专 ...
基于Matlab模拟独立瑞利衰落下双分支分集接收机QPSK的误码率
基于Matlab模拟独立瑞利衰落下双分支分集接收机QPSK的误码率本文介绍了如何使用Matlab仿真独立瑞利衰落下双分支分集接收机QPSK的误码率,其中包括如何生成数据.如何进行信道仿真.如何实现Q ...

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

1.1 模拟微透镜阵列。

1.2 微透镜分割波前

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

参考文献：

基于Matlab模拟哈特曼波前探测器相关推荐

最新文章

热门文章