本发明属于傅里叶叠层成像技术领域，具体涉及一种基于相位迭代最小化的傅里叶叠层 成像图像重建方法。
背景技术：
自第一台显微镜诞生到现在，光学显微成像技术发展迅速，同时人们对图像质量的要求 也越来越高。人们既希望获得高分辨的图像，又希望获得大视场，然而在普通光学成像设备 中，视场和放大倍数是无法兼具的。大的视场必然导致放大倍数不足，虽然能够观察物体较 大的形貌，但无法看清物体更细微的部分；反之，高放大倍数必然导致视场过小，不利于观 察物体的形貌。从理论上讲，增大数值孔径可以提高显微镜的分辨率，但实际的效果并不是 特别明显，主要有如下方法：第一种方法是增大工作介质的折射率，一般通过注入高折射率 的液体来实现，但所存在的问题在于，液体中会混入杂质、液体溢出还有热胀冷缩等，而且 所提高的数值孔径也相当有限。另一种方法是使用屈光度更大的透镜来组成物镜，这种方法 因物镜焦距变短会使得操作距离有限，不是很方便。其他还有通过增大物镜孔径的方式来提 高数值孔径，需要对物镜进行设计还有矫正畸变，过程十分复杂。
傅里叶叠层成像技术(Fourier ptychographic microscopy，简称FPM)的提出开拓了成像 设备的新格局，其意义在于：高分辨率、大视场及定量相位成像。最初的实验设备是一个方 形的LED阵列作为光源，依次点亮每个位置的LED单元进行成像获得一系列低分辨率的子 图像，再将全部子图像通过图像重建方法生成高分辨率的图像。由于每次只点亮一个LED单 元，在采集子图像的过程中耗时较长，而且采集的子图像数据量也过于庞大，重建过程相对 缓慢。另外，傅里叶叠层成像中存在许多未知系统误差，例如LED单元照明亮度不一致、物 镜的像差等，这些误差会影响傅里叶叠层成像重建图像的质量。因此，为了进一步提高其成 像性能，国内外研究学者提出了一系列改进方法，可分为两类：第一类改进方法主要为了减 少系统成像时间，例如，Lei等(Multiplexed coded illumination for Fourier Ptychography with an LED array microscope)提出使用复合编码的方法，同时点亮多个LED单元成像。第二类 改进方法研究新的傅里叶叠层成像图像重建方法，从而提高重构精度和重构分辨率。例如， Bian等提出了一种基于Wirtinger Flow的傅里叶叠层成像图像重建方法，从而使得在FPM中 通过低信噪比子图像重建高分辨率图像成为可能。
本发明将相位恢复作为傅里叶叠层成像图像重建优化的目标函数，从而提出一种基于相 位迭代最小化的傅里叶叠层成像图像重建方法。
技术实现要素：
本发明提出一种基于相位迭代最小化的傅里叶叠层成像图像重建方法，属于傅里叶叠层 成像方法，解决在傅里叶叠层成像中如何通过所获取的一系列低分辨率子图像有效重建高分 辨率图像的问题，可以应用于各种傅里叶叠层成像系统中。
本发明的一种基于相位迭代最小化的傅里叶叠层成像图像重建方法，具体包括如下步骤：
步骤S1，生成一系列低分辨率子图像对应的幅度子图像；
步骤S2，获取重建算法参数，包括每幅低分辨率幅度子图像生成的频谱在高分辨率重建 频谱中的位置和瞳孔函数，建立低分辨率幅度子图像频谱和高分辨率重建图像频谱之间的联 系；
步骤S3，将相位恢复作为优化目标之一，生成图像重建目标优化的目标函数，并通过迭 代最小化进行求解，最终获得高分辨率重建图像频谱，具体实现如下，
步骤S31，将每幅低分辨率幅度子图像对应的二维矩阵进行矢量化操作，生成一维列向量；
步骤S32，假设经过步骤S31矢量化后的低分辨率幅度子图像幅度一维列向量表示为 bi,i∈[1,L]，其中L表示低分辨率幅度子图像的总数目，相应的相位一维列向量表示为pi， 低分辨率幅度子图像对应高分辨率图像频谱中的单元表示为zi，将图像重建目标优化的目标 函数f表示为：
其中，A(zi)＝Vec(F-1[Mat(zi)×P])，P为瞳孔函数；Mat(·)代表矢量化的逆过 程，即将一维列向量变为二维矩阵，F-1(·)代表傅里叶逆变换，Vec(·)代表二维矩阵的矢 量化运算，⊙代表点乘运算；
步骤S33，采用迭代最小化算法，分别对zi和pi进行迭代求解，具体如下，
k代表迭代次数，Δz代表梯度下降步长，fz代表f关于zi的一阶偏导数，./代表点除 运算；
步骤S34，Δz采用下式进行计算，
min(x,y)运算代表计算x和y中的最小值，um为阈值常数，τ0为比例常数， τ＝k×a，a为常数；
步骤S35，基于式(3)，可以采用下式计算，
A*(·)＝Vec(F[Mat(·)]×P*)，A*是A的共轭转置矩阵，P*为P的共轭转置矩阵， F[·]代表傅里叶变换；
步骤S36，由zi重建z，每个zi对应于z中的一部分，假设zi的频谱在高分辨率重建频 谱中的位置左上角坐标为(q1,q2)，则根据步骤S2中低分辨率幅度子图像频谱和高分辨率重建 图像频谱之间的联系，得到如下关系式，
Mat(z)(q1:q1+n1,q2:q2+n2)⊙L＝Mat(zi)×P (6)
其中，Mat(z)(1:n1,1:n2)中1:n1表示从第一行至第n1行中的所有像素，1:n2表示从第一列 至第n2列中的所有像素，L代表一个二值模板信息，其分辨率等于P，L由P生成，若坐标(i,j) 处的值P(i,j)≠0，则L(i,j)＝1，否则，L(i,j)＝0；
步骤S37，重复步骤S33至步骤S36，迭代次数为T；
步骤S4，对高分辨率重建图像频谱进行傅里叶反变换，将结果的模作为高分辨率重建后 的图像。
而且，步骤S32中，根据目标函数，初始的pi设置为一个全0的一维列向量，初始zi的值 通过如下方式生成：计算||bi||2,i∈[1,L]中的最大值q，||·||2代表二范数；计算q对应的频谱值 q'，进一步将q'放大至与高分辨率图像频谱相同大小，将放大后的结果中与zi空间位置对应 的部分作为zi的初始值。
而且，步骤S36中，相邻zi的频谱在高分辨率重建频谱中存在一定的重叠，对重叠部分 的处理方式如下，
设P1和P2分别代表两个相邻Mat(z1)和Mat(z2)对应的区域，假设两个区域对应的重叠 部分中为P3，那么Mat(z)在坐标为(i,j)处的频谱值p(i,j)计算过程如下：
(1)如果(i,j)属于P1但不属于P3，则p(i,j)等于P1对应坐标为(i,j)处的频谱值；
(2)如果(i,j)属于P2但不属于P3，则p(i,j)等于P2对应坐标为(i,j)处的频谱值；
(3)如果(i,j)属于P3，则p(i,j)等于P1对应坐标为(i,j)处的频谱值和P2对应坐标为(i,j)处 的频谱值相加取平均。
而且，所述步骤S37中迭代次数为100。
而且，所述步骤1中通过宽视场、高分辨率傅里叶折叠成像显微镜中的硬件系统生成一系 列低分辨率子图像对应的幅度子图像。
与现有技术相比，本发明的优点和有益效果：本发明方法通过建立低分辨率幅度子图像 和高分辨率重建图像之间的联系，然后将相位恢复作为优化目标之一，生成图像重建目标优 化的目标函数，并通过迭代最小化进行求解，从而将重建问题转化为多个步骤的迭代运算， 求出高分辨率重建图像频谱，与现有方法相比，本发明通过相同数目的低分辨率子图像进行 高分辨率图像重建，具有更好的重建强度图像质量。
附图说明
附图1为本发明实施例流程图；
附图2为本发明实施例中低分辨率幅度子图像频谱与高分辨率重建图像频谱之间的关系 示意图；
附图3为本发明实施例中步骤S36加权平均方法示意图；
附图4为本发明实施例中低分辨率幅度子图像频谱恢复高分辨率重建图像频谱示意图；
附图5为本发明实施例中原始高分辨率强度图像；
附图6为本发明实施例中原始高分辨率相位图像；
附图7为本发明实施例1中仿真的瞳孔函数；
附图8本发明实施例1中重建结果，(a)重建的高分辨率强度图像，(b)重建的高分辨 率相位图像；
附图9本发明实施例中基于WFP算法生成的高分辨率强度图像和高分辨率相位图像，(a) 基于WFP生成的高分辨率强度图像，(b)基于WFP生成高分辨率相位图像；
附图10为本发明实施例2中仿真的瞳孔函数；
附图11为本发明实施例2中重建结果，(a)重建的高分辨率强度图像，(b)重建的高分 辨率相位图像。
具体实施方式
下面结合附图和实施例对本发明的技术方案作进一步说明。
如图1，本发明实施例的技术方案可通过以下步骤实现：
S1.通过宽视场、高分辨率傅里叶折叠成像显微镜(Wide-field,high-resolution Fourier ptychographic microscopy)中的硬件系统生成一系列低分辨率子图像对应的幅度子图像，其中， 低分辨率子图像可以是实际拍摄得到的图像，也可以通过模拟生成；
每幅低分辨率子图像对应的幅度子图像具有相同的生成方法。假设第i幅低分辨率子图像 空间坐标(x,y)对应的像素强度为I(x,y),则对应的幅度子图像空间坐标(x,y)对应的像素幅 度I'(x,y),计算方法为：
S2.获取重建算法参数，包括每幅低分辨率幅度子图像生成的频谱在高分辨率重建频谱中 的位置，瞳孔函数，建立低分辨率幅度子图像和高分辨率重建图像之间的联系；
假设低分辨幅度子图像的分辨率为(n1,n2)，重建高分辨率图像的分辨率为(N1,N2)，其 中，每幅低分辨率幅度子图像生成的频谱在高分辨率重建频谱中的位置是根据生成每幅子图 像对应的LED在LED阵列中的空间几何位置计算得到，具体计算方法可参见文献(Wide-field, high-resolution Fourier ptychographic microscopy，宽视场、高分辨率傅里叶折叠成像显微镜)， 本发明不予撰述；瞳孔函数以二维矩阵P表示，在理想状况下可认为是一个二值化模板。实 际应用中可以通过视场角、放大倍数、像元尺寸等信息，采用(Wide-field,high-resolution Fourier ptychographic microscopy，宽视场、高分辨率傅里叶折叠成像显微镜)方法生成所需的瞳孔函 数。以第i幅低分辨率幅度子图像为例，介绍高分辨率重建图像频谱G与低分辨率幅度子图像 频谱F之间的关系。假设第i幅低分辨率幅度子图像的频谱在高分辨率重建频谱中的位置左上 角坐标为(q1,q2)，则具有如下关系式：
G(q1:q1+n1,q2:q2+n2)⊙L＝F(1:n1,1:n2)×P (1)
其中，F(1:n1,1:n2)中1:n1表示从第一行至第n1行中的所有像素，1:n2表示从第一列至第 n2列中的所有像素，L代表一个二值模板信息，其分辨率等于P，L由P生成，若坐标(i,j)处 的值P(i,j)≠0，则L(i,j)＝1，否则，L(i,j)＝0，其示意图如附图2所示。
S3.将相位恢复作为优化目标之一，生成图像重建目标优化的目标函数，并通过迭代最 小化进行求解，从而将重建问题转化为多个步骤的迭代运算，求出高分辨率重建图像频谱， 具体实现如下：
S31:将每幅低分辨率幅度子图像(由步骤S1产生)进行矢量化操作，生成一维列向量， 方法如下：
假设每幅低分辨率幅度子图像(二维矩阵M)可以表示为
mij代表二维矩阵M中坐标(i,j)处的值，即幅度值，M和N分别指的是二维矩阵M的行 数和列数。
则经过矢量化操作后，对应的一维列向量M'可以表示为：
M'＝[m11,m21,…,mM1,m12,…,mM2,…,m1N,…,mMN]T
S32:假设经过步骤S31矢量化后的低分辨率幅度子图像幅度一维列向量表示为 bi,i∈[1,L]，其中L表示低分辨率幅度子图像的总数目，相应的相位一维列向量表示为pi(初始pi可设置为一个全0的一维列向量)，其中bi和pi与M'的关系如下：bi和pi分别是 第i幅低分辨率幅度子图像对应的幅度一维列向量和相位一维列向量。M'与M对应，用于 描述如何通过二维矩阵生成所需的一维列向量。对应高分辨率图像频谱中的单元表示为zi (zi和pi的生成过程见公式(3))，根据步骤S2的关系式，可以将图像重建目标优化的目标函 数f表示为：
其中，A(zi)＝Vec(F-1[Mat(zi)×P])，P为瞳孔函数。Mat(·)代表矢量化的逆 过程，即将一维列向量变为二维矩阵。F-1(·)代表傅里叶逆变换。Vec(·)代表二维矩阵的 矢量化运算。⊙代表点乘运算。
初始zi的值可通过如下方式产生：1)计算||bi||2,i∈[1,L]中的最大值q，||·||2代表二范数。 计算q对应的频谱值q'，进一步将q'放大至与高分辨率图像频谱相同大小。那么，将该放大 后的结果中与zi空间位置对应的部分作为zi的初始值。
S33:为了求解式(2)，采用迭代最小化算法，分别对zi和pi进行迭代求解，具体如下：
k代表迭代次数，Δz代表梯度下降步长，fz代表f关于zi的一阶偏导数，./代表点除 运算。
S34:Δz可以采用下式进行计算：
min(x,y)运算代表计算x和y中的最小值，um为阈值常数，τ0为比例常数， τ＝k×a，a为常数，可以采用文献(Phase retrieval via wirtinger flow:Theory and algorithms) 中的方法，对相关参数进行取值，本实施例中um＝0.4，τ0＝330，a＝1。
S35:基于式(3)，可以采用下式计算：
A*(·)＝Vec(F[Mat(·)]×P*)，A*是A的共轭转置矩阵，P*为P的共轭转置矩阵， F[·]代表傅里叶变换。
S36:由zi重建z，其示意图如附图4所示。每个zi对应于z中的一部分，假设zi的频谱 在高分辨率重建频谱中的位置左上角坐标为(q1,q2)，则根据式(1)具有如下关系式：
Mat(z)(q1:q1+n1,q2:q2+n2)⊙L＝Mat(zi)×P (6)
考虑到相邻zi的频谱在高分辨率重建频谱中会有一定的重叠，因此需要对重叠部分进行 加权平均，其示意图如附图3所示。这里以P1和P2两个区域为例，P1和P2分别代表两个相邻 Mat(z1)和Mat(z2)对应的区域，假设两个区域对应的重叠部分中为P3，那么Mat(z)在坐标为 (i,j)处的频谱值p(i,j)计算过程如下：
(1)如果(i,j)属于P1但不属于P3，则p(i,j)等于P1对应坐标为(i,j)处的频谱值；
(2)如果(i,j)属于P2但不属于P3，则p(i,j)等于P2对应坐标为(i,j)处的频谱值；
(3)如果(i,j)属于P3，则p(i,j)等于P1对应坐标为(i,j)处的频谱值和P2对应坐标为(i,j) 处的频谱值相加取平均，本实施例中两者赋予相同的权重。；
S37:重复步骤S33至步骤S36，迭代次数为T(T为给定常数，例如100)。
S4.对高分辨率重建图像频谱进行傅里叶反变换，将结果的模作为高分辨率重建后的图 像。
实施例1：
1、仿真数据采用附图5作为原始高分辨率强度图像X(输入是强度图像，按照步骤S1生 成低分辨幅度子图像)，附图6作为原始高分辨率相位图像Y，分辨率均为512×512。基于附 图5和附图6构建仿真输入高分辨率图像Z＝X(cosY+i×sinY)；
2、基于傅里叶叠层成像原理，通过Z生成一系列低分辨率子图像，分辨率均为64×64， 合计225子图像，相邻低分辨率子图像频谱对应的高分辨率图像频谱位置有一定区域重叠，迭 代次数为100次；
3、仿真的瞳孔函数如附图7所示，其大小为64×64。以中心(假设其空间坐标为(mi,mj)) 为圆点，半径r＝26为半径画圆。对于瞳孔函数坐标为(x,y)处的值P(x,y)计算过程如下：
如果则P(x,y)＝1；否则，P(x,y)＝0。；
4、重建得到的高分辨率强度图像X'和高分辨率相位图像Y'(分辨率均为512×512)如 附图8所示；
5、同时通过WFP算法(Wirtinger Flow)对上述仿真数据进行实验验证，得到的高分辨率 强度图像和高分辨率相位图像如附图9所示。对比附图8(a)、附图9(a)及附图5，可以看出附图 8(a)与附图5更加的相似，附图9(a)相对较为模糊，纹理及边缘均不太清晰。因此本发明方法 具有更好的重建图像质量。
实施例2：
1、采用与实施例1相同的仿真输入高分辨率图像Z＝X(cosY+i×sinY)；
2、基于傅里叶叠层成像原理，通过Z生成一系列低分辨率子图像，分辨率均为128×128， 合计36子图像，相邻低分辨率子图像频谱对应的高分辨率图像频谱位置有一定区域重叠，迭 代次数为100次；
3、仿真的瞳孔函数如附图10所示，其大小为128×128。以中心(假设其空间坐标为 (mi,mj))为圆点，半径r＝51为半径画圆。对于瞳孔函数坐标为(x,y)处的值P(x,y)计算过 程如下：
如果则P(x,y)＝1；否则，P(x,y)＝0。
4、重建得到的高分辨率强度图像X'和高分辨率相位图像Y'(分辨率均为512×512) 如附图11所示。对比附图11和附图9，可以发现附图11中重建强度图像的边缘及细节更为清晰， 具有更好的重建强度图像质量。这说明低分辨率子图像分辨率的提升有益于提升重建强度图 像质量。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技 术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不 会偏离本发明的精神或者超越所附权利要求书所定义的范围。