波前重构 Southwell区域法 Shack-Hartmann波前探测器

内容简介：推导见正文，代码见代码，另外一个很有用的Shack-Hartmann波前探测器资源链接见其他资源

逼逼叨：做Shack-Hartmann波前探测器，需要波前重构，懒得自己写，想找个现成的，发现有限差分最小二乘居然要卖一百多巨款！！！matlab有自带的函数（ lsqr - 求解线性方程组 - 最小二乘法），矩阵自己手动写一下，发现要不了多久，放出来。

推导

目的：需要根据测量得到的波前斜率重构波前相位。
方程：写在word里懒得敲了，直接截图了

结果如下：

代码

%测量单个二次相位透镜，看不同的波前斜率的输出的位置
close all;

%% 参数设置
M_len=55;N_len=55;%单个透镜包含的像素数量
num_len_x=15;num_len_y=15;%x,y方向包含的透镜数量
M=num_len_y*M_len;
N=num_len_x*N_len;
x=(-(N-1)/2:(N-1)/2).*dx;
y=(-(M-1)/2:(M-1)/2).*dy;
[Y,X]=ndgrid(x,y);%z=0面的坐标

%% 生成相位分布，计算斜率
phase_test_set=-k0/2*(X.*X+Y.*Y)./(r_max*num_len_y*sqrt(2));
x_center=zeros(num_len_y,num_len_x);%中心位置
y_center=zeros(num_len_y,num_len_x);
%中心位置坐标
for i=1:num_len_x
    for j=1:num_len_y
        x_center(i,j)=(-(num_len_x+1)/2+j)*N_len*dx;
        y_center(i,j)=(-(num_len_y+1)/2+i)*M_len*dy;
    end
end
S_x_set=-k0*x_center./(r_max*num_len_y*sqrt(2));%中心斜率
S_y_set=-k0*y_center./(r_max*num_len_y*sqrt(2));
r2_center=x_center.*x_center+y_center.*y_center;
phi_center_set=-k0/2*(r2_center)./(r_max*num_len_y*sqrt(2));


%% 根据dW/dx,DW/dy波前重构
%认为恢复的相位的坐标在透镜的中心位置采用southwell
Px=N_len*dx;Py=M_len*dy;
%构造矩阵
%Sx和系数矩阵的上半部分
SA_x=zeros(num_len_y*(num_len_x-1),1);
B_h=zeros(num_len_y*(num_len_x-1),2*num_len_y*num_len_x);
B_l=zeros((num_len_y-1)*num_len_x,2*num_len_y*num_len_x);
for i=1:num_len_y
   for j=1:num_len_x-1
    SA_x((i-1)*(num_len_x-1)+j)=S_x_set(i,j)+S_x_set(i,j+1);
    %B的上半部分
    hang=(i-1)*(num_len_x-1)+j;%方程编号
    lie_1=(i-1)*(num_len_x)+j;%对应的phi位置
    lie_2=lie_1+1;
    B_h(hang,lie_1)=-1;
    B_h(hang,lie_2)=1;
   end
end
%Sy和系数矩阵的下半部分
SA_y=zeros(num_len_x*(num_len_y-1),1);
for i=1:num_len_y-1
   for j=1:num_len_x
    SA_y((i-1)*num_len_x+j)=S_y_set(i,j)+S_y_set(i+1,j);
    %B的下半部分
    hang=(i-1)*(num_len_x)+j;%方程编号
    lie_1=(i-1)*(num_len_x)+j;%对应的phi位置
    lie_2=i*(num_len_x)+j;
    B_l(hang,lie_1)=-1;
    B_l(hang,lie_2)=1;
   end
end
%组合SA矩阵
SA=zeros(2*num_len_y*(num_len_x-1),1);
SA(1:num_len_y*(num_len_x-1),1)=Px/2*SA_x;
SA(1+num_len_y*(num_len_x-1):2*num_len_y*(num_len_x-1),1)=Py/2*SA_y;
B=zeros((num_len_y-1)*num_len_x+num_len_y*(num_len_x-1),2*num_len_y*num_len_x);
B(1:num_len_y*(num_len_x-1),:)=B_h;
B(1+num_len_y*(num_len_x-1):2*num_len_y*(num_len_x-1),:)=B_l;

phi_matrix = lsqr(B,SA,1e-10,500);
phi_solve=zeros(num_len_y,num_len_x);
for i=1:num_len_y
   for j=1:num_len_x
     phi_solve(i,j)=phi_matrix((i-1)*(num_len_x)+j,1);
   end
end

phi_solve=phi_solve-mean(mean(phi_solve))+mean(mean(phi_center_set));%减去相差的常数
figure;
meshc(x,y,phase_test_set);%原相位
hold on;
plot3(x_center,y_center,phi_solve,"*");%恢复的相位

其他资源

matlab社区上的沙克-哈特曼模拟器
基本每一步都有，要考虑的都考虑到了，而且都写成函数了，直接调用就行，非常好用。
引用格式：SergioB (2023). Shack-Hartmann Simulator (https://github.com/SergioBonaqueGonzalez/Shack-Hartmann-Simulator), GitHub. 检索来源 2023/9/12.
资源链接：https://ww2.mathworks.cn/matlabcentral/fileexchange/63851-shack-hartmann-simulator