目的：实现区域法的Matlab代码
1.理论介绍
波前重建是一个积分过程，它将一系列独立的斜率或梯度测量值转换成一个平滑变化的三维表面，定义要校正的波前误差。波前梯度（斜率）测量不可避免地受到随机噪声的干扰，这是因为光的量子性质和探测过程中电子的添加。由于波前上任意两点之间存在多条梯度（斜率）路径，因此没有唯一的波前精确满足所有测量的梯度，因此采用了统计解。通常的标准是最小化重建波前和单个梯度测量之间的均方误差。
波前传感器进行区域测量的所有自适应光学系统都需要进行波前重建，即区域法波前重构即从波前探测器测得的波前斜率中恢复出波前的相位值。重建器的几何结构（即波前评估节点与测量区域之间的关系）取决于所用传感器的类型。剪切干涉仪和夏克哈特曼传感器等梯度传感器使用不同的探测器配置。根据子孔径和测量斜率之间的关系分为：休晋模型（Hudgin）[1]、绍契威尔模型（Southwell）[2]和弗雷德模型（Fried）[3]，如下图所示。这里只介绍用于夏克-哈特曼波前探测器的Fried和Southwell模型的重构算法。
1.1 Fried模型，对于N×N 微透镜阵列构成的网格：
这种配置主要用于Shack-Hartmann传感器，其中在相同的子孔径中测量和梯度。Fried曾对这种结构进行过分析。尽管它使用单个探测器阵列，但这种配置由两个独立的网格组成，当测量的梯度分解为45°分量，连接对角节点时，可以看出这一点。两个独立网格的存在需要某种方法来建立它们之间的关系。两个网格之间存在位移，产生棋盘形波前图案，波前传感器对此不敏感。这不是一件小事，产生的错误可能很难消除。
以3*3的子透镜为例，解释上式，对其展开
改写为矩阵形式：
记作：S=AΦ，A根据N而确定。S是波前探测器测得的斜率，Φ是重构矩阵的相位值。
其与变形镜驱动器的匹配方式如下图所示：
Fried配置最常用于Shack–Hartmann波前传感器，在相同的子孔径中测量和梯度。其中驱动器与四个相邻透镜的角交点对齐。上右图显示了97个驱动器DM的配置。驱动器周围的透镜对Fried配置中的驱动器位移非常敏感，使校准更容易。然而，梯度测量可以分解为45°分量，连接对角节点。这将导致两个独立的交错网络。
1.2 Hudgin 模型
在 Hudgin 模型中，根据两个相邻网格点之间的中点测量波前斜率，x、y 方向上每对相邻网格点之间的相位可以估计为：
1.3 而Southwell模型
对于Shack Hartmann传感器，默认几何结构称为Southwell配置，其特征是波前样本与波前斜率测量值一致。这种配置最直观，因为它被设计用来估计每个子孔径的波前误差。然而，波前斜率与子孔径中心处的期望波前值相关时，需要采用间接计算路线。必须首先将Southwell配置中的波前斜率数据转换为Hudgin配置，如图2.11b中四个子孔径的补丁所示。在当地，计算只是相邻斜坡的简单平均值（连接两个节点的斜率（梯度）是在以这些节点为中心的区域中测量的两个梯度的平均值）。
连接两个节点的斜率（梯度）是在以这些节点为中心的区域中测量的两个梯度的平均值：
其中，Sx和Sy是代表局部波前斜率的标量，m和n是子孔径的序号。在Hudgin配置下（上式和Hudgin的公式联合），得到斜率可以直接与采样波前关系：
如果我们把它的极限设为d，当d接近零时，它就变成了导数的标准定义。为了估计整个波前，必须对整个瞳孔进行计算。上式记作：
连续区间 V.S. 离散空间比较：
即，对于N×N 微透镜阵列构成的网格，ds分隔的两个相邻相位值之差表示：
以3*3的子透镜为例，对其展开解释：
改成矩阵形式：
这就是为啥重构矩阵C中元素只有0、0.5，而矩阵E
记作: CS=EΦ, C和E根据N的个数来确定，矩阵C中元素只有0、0.5；而矩阵E的元素只有0、1、-1。S是波前探测器测得的斜率，Φ是重构矩阵的相位值。
其与变形镜驱动器的匹配方式如下图所示：
在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.