哈特曼波前探测器（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(微透镜分割后的波前') 这只是一个简单的开始。如果子镜头太少，就会有更多的波浪被丢弃，如下图所示。此外，没有考虑镜头本身的图像差异，本文只是一个简单的介绍。
用菲涅尔衍射获得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.*fftshift(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.