目的：实现区域法的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中解决propertiesxslope;    % x方向上的斜率n ;        % 子透镜的个数；不是有效子透镜yslope;    % y方向上的斜率S;         % S矩阵 S E C 是Southwell重构波前中的过渡矩阵E;         % E矩阵C;         % C矩阵SouthWF;   % Southwell重构出来的波前mask;      % 根据x或者y方向上的斜率定义出maskFriedWF;   % 采用Fried重构出的波前A;         % Fried重构波前中，从斜率到波前相位的过渡矩阵Fried_mask ;% Fried的maskendmethods  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);endendmethods(Static) %% function E = getE(n)E     = zeros(2*n*(n-1),n*n); for i = 1:nfor 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;endendend%%function C = getC(n)C = zeros(2*n*(n-1),2*n*n);for i = 1:nfor 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;endendend%%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:Nif mod (i+num , m) == 0num = num+1;end% Ax MATRIX WILL HAVE THE FORMAx(i,i+num)     = -1; %  -1 -1  0  1  1  0  0  0 0Ax(i,i+num+1)   = -1; %   0 -1 -1  0  1  1  0  0 0Ax(i,i+num+m)   =  1; %   0  0  0 -1 -1  0  1  1 0Ax(i,i+num+m+1) =  1; %   0  0  0  0 -1 -1  0  1 1% Ay MATRIX WILL HAVE THE FORMAy(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 0Ay(i,i+num+m)   = -1; %   0  0  0 -1  1  0 -1  1 0Ay(i,i+num+m+1) =  1; %   0  0  0  0 -1  1  0 -1 1endA         = 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]);endfunction Ainv = SVDgetInv(A)% 本代码用于求矩阵的伪逆， 也可以直接用pinv求得[u,d, v] = svd(A);   % 奇异值分解Dinv     = zeros(M); for i=1:Mif d(i,i)<10^-6Dinv(i, i) = 0;elseDinv(i, i) = 1/d(i,i);endend 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中有效子透镜的masknElements            = 2*nLenslet+1; % Linear number of lenslet+actuatorvalidLensletActuator = zeros(nElements);index                = 2:2:nElements; % Lenslet indexvalidLensletActuator(index,index) = validLenslet;for xLenslet = indexfor yLenslet = indexif validLensletActuator(xLenslet,yLenslet)==1xActuatorIndice = [xLenslet-1,xLenslet-1,...xLenslet+1,xLenslet+1];yActuatorIndice = [yLenslet-1,yLenslet+1,...yLenslet+1,yLenslet-1];validLensletActuator(xActuatorIndice,yActuatorIndice) = 1;endendendindex = 1:2:nElements; % Actuator indexval   = logical(validLensletActuator(index,index));% figure(10); imagesc(val)endend
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.

区域法(zonal method)重构波前相关推荐

哈特曼波前传感器区域法重构算法实例
一.引言目前哈特曼波前传感器波前重建的算法主要有两种,即区域波前重建法和模式波前重建法. 而区域法重构按求解方法又可分为快速傅里叶变换法,迭代求解法,矩阵向量法等. 本文将分析一下迭代法求解法重构波 ...
matlab区域法直线度,用最小区域法求直线度误差的探讨
庸 [1l } 程技术学院学报 1991年第 l期 7OURNAL OF TANGSf iAN INSTITUTE OF TECHNOLOGY 1．1991 用最小区域法求直线度误差 ...
最小误差阈值分割 matlab,原创：最小包容区域法处理圆度误差的程序算法
希望可以加分,这能鼓励我们对生活工作中遇到的问题应用matlab来解决. 在工作中遇到这样的一个问题:当用三坐标测量圆特征,得到平面上离散分布的若干点,这些点近似分布在圆周上.如何评价该圆特征的圆度? ...
结构光学习 | 波前重构技术
引子光波波前误差是影响发射激光束的质量或光学成像质量最主要的因素,但在自适应光学系统中,一般不能直接获得光波波前误差的数据以进行校正,而只能测得离散的波前斜率或离焦面上的光强分布,这就需要从上述离散 ...
基于Matlab传统有限差分的最小二乘积分波前重构算法
三维面形重建是利用测得的表面梯度(斜率)或法向矢量数据重建物体表面形貌的过程,是光学测量.机器视觉中的一个经典问题,是哈特曼波前检测.干涉测量.光栅投影测量等技术的关键步骤之一. 本文开发了一种基于传 ...
R语言：接受拒绝法（Acceptance-Rejection Method）生成随机数
接受拒绝法Acceptance-Rejection Method 原理解释具体步骤连续情况举例离散情况举例计算机能够通过调用自带函数,实现对均匀分布等常见分布的采样.但在实际情况中分布通常由其 ...
牛顿法求解1-100的平方根python_使用牛顿-拉弗森法定义平方根函数（Newton-Raphson method Square Root Python）...
牛顿法(Newton's method)又称为牛顿-拉弗森法(Newton-Raphson method),是一种近似求解实数方程式的方法.(注:Joseph Raphson在1690年出版的< ...
基于Matlab有限差分的高阶迭代最小二乘积分的波前重构算法
该算法根据x与y方向波前斜率,可对其波前进行重构. 一.算法验证首先,分别得到x方向与y方向的波前斜率,分别如下图所示: 图1 X方向波前斜率图2 Y方向波前斜率基于有限差分的高阶迭代最小二乘积 ...
代码重构(二)：类重构规则
在上篇博客<代码重构(一):函数重构规则(Swift版)>中,详细的介绍了函数的重构规则,其中主要包括:Extract Method, Inline Method, Inline Temp ...

区域法(zonal method)重构波前