1.理论介绍
波前重建是一个积分过程，它将一系列独立的斜率或梯度测量值转换为光滑变化的三维表面，并定义要校正的波前误差。由于光的量子性质和探测过程中电子的添加，波前梯度(斜率)的测量不可避免地会受到随机噪声的干扰。由于波前任意两点之间有多个梯度（斜率）路径，因此没有唯一的波前精确满足所有测量梯度，因此采用了统计解。一般标准是最小化重建波前与单梯度测量之间的平方误差。 所有由波前传感器测量的自适应光学系统都需要进行波前重建，即从波前探测器测量的波前斜率中恢复波前的相位值。根据所用传感器的类型，重建器的几何结构（即波前评估节点与测量区域之间的关系）。不同的探测器配置了剪切干涉仪、夏克哈特曼传感器等梯度传感器。分为子孔径与测量斜率之间的关系：休晋模型（Hudgin）[1]、绍契威尔模型（Southwell）[2]和弗雷德模型（Fried）[3]如下图所示。这里只介绍夏克-哈特曼波前探测器Fried和Southwell模型重构算法。
1.1 Fried模型，对于N×N由微透镜阵列组成的网格：
该配置主要用于Shack-Hartmann在相同的子孔径中测量和梯度的传感器。Fried分析了这种结构。虽然采用单个探测器阵列，但当测量梯度分解为45时，该配置由两个独立的网格组成°在连接对角节点时，可以看到重量。两个独立网格的存在需要某种方法来建立它们之间的关系。两个网格之间存在位移，产生棋盘波前图案，波前传感器对此不敏感。这不是一件小事，错误可能很难消除。
改写为矩阵形式：
记作：
Fried最常用的配置Shack–Hartmann波前传感器在相同的子孔径中测量和梯度。驱动器对齐四个相邻透镜的角交点。右上图显示了97个驱动器DM的配置。驱动器周围的透镜对Fried配置中的驱动器位移非常敏感，使校准更容易。但梯度测量可分解为45°重量，连接对角节点。这将导致两个独立的交错网络。
1.2Hudgin 模型
1.3 而Southwell模型
对于Shack Hartmann默认几何结构称为传感器Southwell其特点是波前样本与波前斜率测量值一致。这种配置是最直观的，因为它被设计用来估计每个子孔径的波前误差。然而，当波前斜率与子孔径中心的预期波前值相关时，需要间接计算路线。必须首先将Southwell波前斜率数据在配置中转换为Hudgin配置，如图2.11b中四个子孔径补丁显示。在当地，只计算相邻斜坡的简单平均值（连接两个节点的斜率（梯度）是以这些节点为中心区域测量的两个梯度的平均值）。
连接两个节点的斜率(梯度)是以这些节点为中心的两个梯度的平均值:
其中，Sx和Sy标量代表局部波前斜率，m和n是子孔径的序号。Hudgin配置下(上和Hudgin公式组合)，斜率可直接与采样波前关系:
如果我们把它的极限设为d，当d接近零时，它成为导数的标准定义。为了估计整个波前，必须计算整个瞳孔。上述记录：
连续区间 V.S. 离散空间比较：
即，对于N×N由微透镜阵列组成的网格，ds分隔的两个相邻相位值之差表示：
改成矩阵形式：
这就是为啥重构矩阵C中元素只有
记作:
在Southwell配置中，波前传感器的每个透镜集中在波前校正器的一个驱动器上。这种配置使得AO系统更难校准，因为一个驱动器的移动在以执行器为中心的相应透镜中没有被感应到。该驱动器的梯度影响函数主要由其相邻透镜测量。
2.Matlba代码实现：
命名为：ZonalRecWF.m，是zonal method reconstructed wavefront的缩写哦，代码被整合成一个函数，方便自己调用而已。
classdef ZonalRecWF
% 区域法重构波前：
% 1.SouthWell模型；这里的斜率是需要转置的
% x方向：0.5*(Sx(i,j+1)+Sx(i,j)) = W(i,j+1)-W(i,j); i=1,N ;j=1,N-1; N是子透镜的个数
% y方向：0.5*(Sx(i+1,j)+Sx(i,j)) = W(i+1,j)-W(i,j); i=1,N-1;j=1,N
% 可以写为 C * s = E * w; C矩阵中的元素全是0.5，E矩阵中的元素为1和-1，两者都是稀疏矩阵
% s是x和y方向上斜率组合在一起后的结果；w是波前数据也即是要求的波前相位值；
% 2.Fried模型重构波前, W是单个透镜四个个点上的值
% W(i+0,j) W(i+0,j+1) W(i+0,j+2)
% 透镜(i+0,j+0) 透镜(i+0,j+1)
% W(i+1,j) W(i+1,j+1) W(i+1,j+2)
% 透镜(i+1,j+0) 透镜(i+1,j+1)
% W(i+2,j) W(i+2,j+1) W(i+2,j+2)
% x方向：gx(i,j) = 0.5*({W(i+1,j+1)+W(i+1,j)} - {W(i,j+1)+W(i,j)});
% y方向：gy(i,j) = 0.5*({W(i+1,j+1)+W(i,j+1)} - {W(i+1,j)+W(i,j)});
% 调用方式： www = ZonalRecWF(dZx,dZy);
% 调用之后，可能colorbar上的数值不一样，自己根据倍数自行调整；
% 暂时还不能解决Southwell重构出波前的mask问题，还需要重新从slopemmse中解决
properties
xslope; % x方向上的斜率
n ; % 子透镜的个数；不是有效子透镜
yslope; % y方向上的斜率
S; % S矩阵 S E C 是Southwell重构波前中的过渡矩阵
E; % E矩阵
C; % C矩阵
SouthWF; % Southwell重构出来的波前
mask; % 根据x或者y方向上的斜率定义出mask
FriedWF; % 采用Fried重构出的波前
A; % Fried重构波前中，从斜率到波前相位的过渡矩阵
Fried_mask ;% Fried的mask
end
methods
function obj = ZonalRecWF(xslope,yslope)
p = inputParser;
p.addOptional('xslope', @isnumeric);
p.addOptional('yslope', @isnumeric);
p.parse(xslope,yslope);
obj.xslope = xslope;
obj.yslope = yslope;
obj.mask = xslope;
obj.mask(find(obj.mask~=0)) = 1;
obj.n = length(obj.xslope);
obj.S = obj.getS(obj.xslope',obj.yslope',obj.n);
obj.E = obj.getE(obj.n);
obj.C = obj.getC(obj.n);
obj.SouthWF = pinv(obj.E) * obj.C * obj.S;
obj.SouthWF = reshape(obj.SouthWF, obj.n, obj.n).* obj.mask;
obj.SouthWF = (obj.SouthWF)';
[obj.FriedWF,obj.A] = obj.Fried_zonal(obj.xslope,obj.yslope,obj.n);
obj.Fried_mask = obj.validActuator(size(obj.xslope,1), obj.mask);
end
end
methods(Static)
%%
function E = getE(n)
E = zeros(2*n*(n-1),n*n);
for i = 1:n
for j = 1:(n-1)
E((i-1)*(n-1)+j,(i-1)*n+j) = -1;
E((i-1)*(n-1)+j,(i-1)*n+j+1) = 1;
E((n+i-1)*(n-1)+j,i+(j-1)*n) = -1;
E((n+i-1)*(n-1)+j,i+j*n) = 1;
end
end
end
%%
function C = getC(n)
C = zeros(2*n*(n-1),2*n*n);
for i = 1:n
for j = 1:(n-1)
C((i-1)*(n-1)+j,(i-1)*n+j) = 0.5;
C((i-1)*(n-1)+j,(i-1)*n+j+1) = 0.5;
C((n+i-1)*(n-1)+j,n*(n+j-1)+i) = 0.5;
C((n+i-1)*(n-1)+j,n*(n+j)+i) = 0.5;
end
end
end
%%
function S = getS(xslope,yslope,n)
S = [reshape(xslope, 1,n*n) reshape(yslope, 1, n*n)]';
end
%%
function [WaveFront,A] = Fried_zonal(xslope,yslope,n)
% 本代码是使用区域法中的Fried重构波前
% n是微透镜的个数
% 特别注意 这里的xslope,yslope不需要转置；如果，我说的是如果哦
% 如果重构的波前和原始波前不一致，则在转置；
num = 0;
m = n+1;
M = m*m;
N = n*n;
Ax = zeros(N, M);
Ay = zeros(N, M);
for i = 1:N
if mod (i+num , m) == 0
num = num+1;
end
% Ax MATRIX WILL HAVE THE FORM
Ax(i,i+num) = -1; % -1 -1 0 1 1 0 0 0 0
Ax(i,i+num+1) = -1; % 0 -1 -1 0 1 1 0 0 0
Ax(i,i+num+m) = 1; % 0 0 0 -1 -1 0 1 1 0
Ax(i,i+num+m+1) = 1; % 0 0 0 0 -1 -1 0 1 1
% Ay MATRIX WILL HAVE THE FORM
Ay(i,i+num) = -1; % -1 1 0 -1 1 0 0 0 0
Ay(i,i+num+1) = 1; % 0 -1 1 0 -1 1 0 0 0
Ay(i,i+num+m) = -1; % 0 0 0 -1 1 0 -1 1 0
Ay(i,i+num+m+1) = 1; % 0 0 0 0 -1 1 0 -1 1
end
A = cat(1, Ax, Ay); % 从斜率到波前相位的过渡矩阵
A = A*0.5;
Ainv = pinv(A); % pinv求伪逆与SVD分解求得的值一样
sx = xslope(:);
sy = yslope(:);
sxy = cat(1, sx, sy);
WaveFront = Ainv*sxy;
WaveFront = reshape(WaveFront,[m,m]);
end
function Ainv = SVDgetInv(A)
% 本代码用于求矩阵的伪逆， 也可以直接用pinv求得
[u,d, v] = svd(A); % 奇异值分解
Dinv = zeros(M);
for i=1:M
if d(i,i)<10^-6
Dinv(i, i) = 0;
else
Dinv(i, i) = 1/d(i,i);
end
end
u = u(:,1:M); % u should have been size (A) ??
ut = u';
Ainv = v*Dinv*ut; % 伪逆矩阵 即pinv(A)
end
% Fried重构的mask函数
function val = validActuator(nLenslet, validLenslet)
% nLenslet = 15;
% validLenslet = wfs.validLenslet;% SHWFS中有效子透镜的mask
nElements = 2*nLenslet+1; % Linear number of lenslet+actuator
validLensletActuator = zeros(nElements);
index = 2:2:nElements; % Lenslet index
validLensletActuator(index,index) = validLenslet;
for xLenslet = index
for yLenslet = index
if validLensletActuator(xLenslet,yLenslet)==1
xActuatorIndice = [xLenslet-1,xLenslet-1,...
xLenslet+1,xLenslet+1];
yActuatorIndice = [yLenslet-1,yLenslet+1,...
yLenslet+1,yLenslet-1];
validLensletActuator(xActuatorIndice,yActuatorIndice) = 1;
end
end
end
index = 1:2:nElements; % Actuator index
val = logical(validLensletActuator(index,index));
% figure(10); imagesc(val)
end
end
end
主函数：
clc; clear all; close all
W = ZonalRecWF(x方向的斜率，y方向的斜率); % 其都是N*N的矩阵数据
figure(1);imagesc(W.SouthWF);colormap(jet);colorbar
figure(2);imagesc(W.FriedWF .* W.Fried_mask);colormap(jet);colorbar
% 如果出错 就用：
% figure(2);imagesc(W.FriedWF );colormap(jet);colorbar
参考文献：
[1] Hudgin R H. Wave-front reconstruction for compensated imaging[J]. Journal of the Optical Society of America, 1977, 67(3): 375-378.
[2] Fried D L. Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements[J]. Journal of the Optical Society of America, 1977, 67(3): 370-375.
[3] Southwell W H. Wave-front estimation from wave-front slope measurements[J]. Journal of the Optical Society of America, 1980, 70(8): 998-1006.