CN110352387B - 用于通过多深度稀疏相位恢复重建全息无透镜图像的系统和方法 - Google Patents

Abstract

一种用于无透镜成像的系统，包括与无透镜图像传感器通信的处理器。该处理器被编程用于操作该图像传感器以获得全息图。该处理器进一步被编程用于使用稀疏性假设，由该全息图生成焦深z处的重建图像X和相位W。

Description

用于通过多深度稀疏相位恢复重建全息无透镜图像的系统和 方法

相关申请的交叉引用

本申请要求2016年11月4日提交的标题为“SYSTEM AND METHOD FORRECONSTRUCTION OF HOLOGRAPHIC LENS-FREE IMAGES BY MULTI-DEPTH SPARSE PHASERECOVERY(用于通过多深度稀疏相位恢复重建全息无透镜图像的系统和方法)”的美国专利申请号62/417,708的优先权，该申请的整个披露内容通过援引并入本文。

背景技术

本披露涉及无透镜成像，更具体地涉及由全息图重建图像。

无透镜成像是基于数字全息术原理的显微技术。在数字全息术中，使用相干光源(例如，激光)来照射对象样本。在光穿过样本时，样本对光进行衍射，并且所得的衍射图样由比如电荷耦合装置(CCD)等图像传感器记录。在记录了所得的衍射图样之后，通过使用衍射过程的数学模型对衍射图样进行后处理来重建样本的图像。

与常规显微术相比，全息无透镜显微术具有几个优点。首先，由于系统中没有透镜，与传统显微镜相比，成像系统的总成本和物理尺寸可以大大减小。另外，无透镜显微术还允许与具有相同放大率的常规显微镜相比具有更宽的成像视场。最后，因为通过对所记录的衍射图样进行后处理来生成样本的图像，所以不需要操作者手动地聚焦系统，因为可以通过后处理自动调整焦深。

图1描绘了在光穿过对象时由对象引起的相干光的衍射。原则上，如果可以记录用相干光源(例如，激光器)照射的对象的精确衍射图样，则衍射图样将包含足够的信息以便在任何期望的焦深处重建样本的完美图像。然而不幸的是，由于光的波性质，对象的完整衍射图样是由空间中给定位置处的电磁波的幅度和相位参数化的复值实体(其中假设波的频率等于照射的相干光)。由于比如CCD等成像传感器典型地只能记录衍射图样的幅度而不能记录相位的事实，因此这对于在实践中应用无透镜成像提出了重大挑战。因此，不能记录完美重建样本的图像所必需的大部分信息，这在试图用常规方法重建样本的图像时表现为显著的伪影。另外，无透镜成像的传统方法典型地仅尝试在单个焦深处重建样本的图像，但是如果样本包含在多个焦深处的对象，则来自离焦对象的衍射将破坏在给定焦深处的重建图像。

发明内容

本发明披露的技术通过利用用于由潜在的多个焦深上记录的衍射图样有效地重建样本的图像、同时生成对缺失相位信息的估计的系统和方法，解决了现有方法的上述两个缺点：差的图像重建和对单个焦深的限制。使用所披露的技术导致与传统图像重建技术相比显著改善的图像质量，允许同时重建样本的三维体积，并且提供找到包含对象的焦深的稳健手段，从而消除了对手动调节焦深的需要。

本发明披露的技术与先前的方法相比是有利的，因为它提供了强大的经验收敛性并且有效地解决了大规模图像的图像重建问题。另外，本发明的方法允许由单个记录的全息图有效地重建3维体积，并且消除了对具有变化的光学波长的光源的需要。

附图说明

为了更全面地理解本披露的性质和目的，应当参考以下结合附图的详细说明，在附图中：

图1是由对象引起的光衍射的示例；

图2是使用不同重建算法的全血样本的一组图像重建，示出了大视场(左)和图像的放大部分(右)，其中顶画面示出使用基本重建算法的图像并且底画面示出了使用本披露的实施例重建的图像；

图3是示出了根据本披露的实施例的方法的图表；

图4是示出了在焦深上重建的3维体积的量值的图表，其中y轴示出了针对在[650，1150]微米范围内的100个均匀间隔的焦深，重建图像的l₁范数‖X_i‖₁，并且竖直线描绘通过手动聚焦图像获得的焦深；

图5是根据本披露的另一个实施例的系统的图；

图6描绘了由根据本披露的另一个实施例的系统获取的全息图的局部重建；

图7描绘了由根据本披露的另一个实施例的系统获取的全息图的远程重建；并且

图8是示出了根据本披露的另一个实施例的方法的图表。

具体实施方式

全息无透镜成像是基于在光已经穿过样本之后记录光的衍射图样并试图反转光衍射过程以重建样本图像的原理。不幸的是，所记录的衍射图样不包含足以完全重建样本图像的信息，因此许多现有的重建技术存在大量的伪影和降低的图像质量。本披露提供了用于与现有技术相比以显著改善的图像质量重建全息无透镜图像的技术，并且允许同时在多个焦深处重建图像的可能性。

在描述对本发明披露的方法和系统进行描述所必需的数学背景之前，描述数学符号。矩阵用大写字母(X)表示，向量用小写粗体字母(x)表示，并且标量用小写字母(x)表示。复数集合表示为

并且实数集合表示为

在给定向量x的情况下，x[i]表示x的第i个条目。同样，在给定矩阵X的情况下，X[i，j]表示X的第i行和第i列的条目。复数x的角度表示为

幅度表示为|x|，并且共轭表示为

针对矩阵(或向量)

|X|表示包含X的条目的绝对值的m×n矩阵，

表示包含X的条目的角度的m×n矩阵，并且

表示包含X的条目的共轭的m×n矩阵。

在给定矩阵(或向量)

的情况下，X的l₁范数表示为||X||₁并且定义为

X的弗罗比尼斯范数表示为||X||_F并且定义为

针对两个矩阵

和

X⊙Y表示其条目的按条目逐个作出的乘积，并且<X，Y>表示两个矩阵的点积。矩阵X的二维离散傅里叶变换表示为

并且离散傅里叶逆变换表示为

两个信号的二维卷积将表示为X*Y。所有零或一的矩阵(或向量)将分别表示为0或1，其中可以从上下文中清楚地看到大小。

全息成像

全息成像过程是基于光学衍射的过程。衍射理论的全面回顾超出了本文档的范围，但是对于全息术中使用的典型距离非常精确的常用近似是将衍射过程建模为二维卷积。具体地，针对平面X₀处的光学波前，将波前向前传播到距离z那么远的平面X_z所产生的波前由以下方程给出

X_z＝T(z)*X₀ (1)

其中T(z)表示对在距离z上的光衍射进行建模的传递函数。针对T(z)可以做出各种选择，这取决于在衍射过程的模型中选择的特定近似(例如，菲涅耳、惠更斯、远场)。在此使用的近似是宽角谱(WAS)近似，其将频域中的传递函数定义为

其中k＝2π/γ是波长为γ的光的波数，并且(k_x，k_y)分别表示水平方向和竖直方向上的空间频率。应当注意，WAS传递函数具有几个容易验证的特性，这些特性将在本发明披露的技术中被利用。

1.T(z₁)*T(z₂)＝T(z₁+z₂)。

2.T(0)*X＝X。

3.

4.线性算子

是酉算子。

还应当记得，通过使用傅里叶变换的特性可以有效地计算X与T(z)的卷积，

使用WAS近似作为光衍射的模型，如果给出从图像传感器记录的衍射图样

则可以尝试在给定的焦深

处找到样本的对应图像

满足条件H＝T(z)*X。因此，重建样本图像的估计的简单方式是求解以下形式的最小二乘问题

应当注意，因为T(z)满足以上列举的特性，并且弗罗比尼斯范数对酉运算不变的事实使得以上问题的最优解可以容易地以封闭形式计算为X^*＝T(-z)*H。然而，应当记得，图像传感器只能记录衍射图样的幅度而不能记录相位(即，在实践中，H在实际上应当是复值时却是实值)。这种限制导致重建图像中的严重伪影(通常称为双图像伪影)，这些伪影典型地表征为从对象发出的波状伪影。图2(顶画面)示出使用这种方法重建的人血的示例图像；应当注意，重建图像中存在大量的波状伪影。

解决这个问题的一种可能性是修改(4)中的问题，以解释应当只关心匹配估计全息图的量值的事实，这可以通过求解以下形式的最小二乘问题来在数学上表达

应当注意，以上问题可以用以下形式等效地表达

其中通过注意到(6)中针对X的任何值、W的最优解由

给出，并将W^*的这个值代入(6)得到(5)，可以看出(5)和(6)之间的等价关系。

虽然将(4)修改为形式(5)和(6)已经解释了记录的全息图H中缺少相位信息的事实，但是应当注意，系统现在是不明确的。这可以通过以下方式看出：考虑(6)并且注意到由于T(z)是酉算子，因此针对W的任何选择，可以生成重建图像X，满足条件H⊙W＝T(z)*X。换句话说，估计的相位项W在这个模型中完全是任意的，因此为了找到有意义的解，必须对问题做出额外的假设。

稀疏相位恢复

回想上一节，由于(6)中无透镜重建问题的性质，系统是不确定的，并且需要额外的假设才能找到有意义的解。具体地，重建图像X应当是“稀疏的”，并且每当样本中的对象仅占据具有视场的像素的一部分时，稀疏性的假设是合理的，其中许多像素等于背景强度。通过向(6)添加正则化项以鼓励解，在数学上并入稀疏性假设，其中重构图像x是稀疏的：

应当注意，测量信号的稀疏性的方法有很多，但是在本讨论中使用l₁范数，因为它具有鼓励稀疏解的期望特性，同时仍然是凸函数并且有利于有效优化。另外，稀疏度的典型测量要求大多数条目统一地都为0，而此处如果像素不包含对象，则像素的值将等于照明光的背景强度。因此，通过向模型添加附加项以使用

标量项捕获(假设为恒定的)照明光来解释非零背景。

求解稀疏相位恢复模型

虽然(7)中给出的模型基于无透镜成像重建问题的性质具有许多理论上的调整，但不幸的是，由于约束条件|W|＝1，优化问题是非凸的。然而，尽管存在这种挑战，在此描述了基于交替最小化的方法，该方法允许对所有变量进行封闭形式更新，是高效的，并且保证至少收敛到目标函数的纳什均衡点。首先，因为与T(z)的卷积针对任何z都是酉算子，所以(7)可以使用以下关系等效地重新表述

其中，等式来自弗罗比尼斯范数的酉不变性、上述T(z)算子的特性以及(2)中T(z)的定义。

根据以上关系，在其他变量保持恒定的情况下，现在可以在(7)中导出对每个变量的闭合形式更新。具体地，具有以下等式，

其中SFT_λ{·}表示复杂的软阈值算子，该算子由下式给出：

使用这些更新的方程，然后可以使用方法1中描述的对变量的一系列交替更新由记录的衍射图样有效地重建图像。虽然由于关于W的约束条件，由(7)描述的优化问题总体上是非凸的，但是在大约仅10-15次迭代内已经观察到方法1的非常强的收敛性。此外，由于方法1中的主要计算负担在于计算傅立叶变换和按元素逐个作出的乘积，因此通过在图形处理单元(GPU)上执行计算，显著加速了计算。图2(底画面)示出使用所提出的方法重建的人血图像。应当注意，基本重建方法中的波状伪像已经完全消除，并且图像中红细胞的对比度已经显著增加，其中使用基本重建方法检测不到的许多细胞现在变得清晰可见。

参考图3，本披露可以体现为用于全息重建图像的方法100。方法100包括获得103全息图。使用稀疏性假设在指定的焦深处生成106重建图像和相位。例如，可以通过求解109满足条件|W|＝1的

来生成106重建图像X和相位W。如上所述，求解109方程(7)的方法是通过使用交替最小化来更新相位W、背景水平μ和重建图像X。

在特定实施例中，通过设定112重建图像X、估计相位W和与背景的拟合度μ的初始值来生成106重建图像X和相位W。计算115传递函数T(z)的傅立叶变换，设定

通过计算118为

来更新背景水平。通过计算

来更新相位。通过

来计算124更新的Q，并且通过计算127X＝SFT_λ(Q-μexp(-izk))来更新重建图像。针对μ、W、Q和X中的每一个，重复130更新步骤以确定重建图像。可以将这些步骤重复130达预定的次数。可以通过执行以上方法100直到重建图像具有令人期望的质量来手动地确定次数。在另一个实施例中，重复130这些步骤，直到μ、W和/或X的更新值的变化小于迭代中的阈值。例如，可以重复130这些步骤，直到在这些步骤的迭代中，重建图像的更新值的变化小于10％、5％、3％、2％、或1％或任何其他百分比值。阈值可以是百分比变化或实际值。可以基于例如与执行附加迭代的成本(例如，计算时间等)平衡的重建图像的质量来选择阈值。

在另一个方面，本披露可以体现为用于无透镜成像的系统10(参见图5)。系统10具有与无透镜图像传感器12通信的处理器14。在一些实施例中，系统10可以包括无透镜图像传感器12，用于获得全息图像。图像传感器12可以是例如有源像素传感器、电荷耦合装置(CCD)、CMOS有源像素传感器等。系统10可以进一步包括光源16，比如相干光源。图像传感器12配置为与光源16协作以获得全息图像。

处理器14被编程用于执行本披露的任何方法。例如，处理器14可以被编程用于操作图像传感器12以获得全息图，并且使用稀疏性假设在一定焦深处生成重建图像和相位。处理器14可以通过例如使用交替最小化更新相位、背景水平和重建图像以求解方程(7)来生成重建图像。

参考图6-7，系统10可以配置用于“局部”重建，例如，其中图像传感器12和处理器14构成系统10。系统10可以进一步包括光源16，用于照射样本。例如，光源16可以是相干光源，比如提供相干光的激光二极管。系统10可以进一步包括样本成像室18，该样本成像室配置为在获取全息图期间容纳样本。在其他实施例(例如，如图7所描绘的实施例)中，系统20配置用于“远程”重建，其中处理器24与图像传感器分离并且通过例如有线或无线网络连接、闪存驱动器等从图像传感器接收信息。

处理器可以与存储器通信和/或包括存储器。存储器可以是例如随机存取存储器(RAM)(例如，动态RAM、静态RAM)、闪存存储器、可移除存储器等。在一些实例中，与执行在此描述的操作相关联的指令(例如，操作图像传感器、生成重建图像)可以存储在存储器和/或存储介质(在一些实施例中，其包括存储指令的数据库)内并且这些指令在处理器处执行。

在一些实例中，处理器包括一个或多个模块和/或部件。由处理器执行的每个模块/部件可以是基于硬件的模块/部件(例如，现场可编程门阵列(FPGA)、专用集成电路(ASIC)、数字信号处理器(DSP))、基于软件的模块(例如，存储在存储器和/或数据库中和/或在处理器处执行的计算机代码的模块)和/或基于硬件的模块和基于软件的模块的组合的任何组合。由处理器执行的每个模块/部件能够执行如在此描述的一个或多个特定功能/操作。在一些实例中，在处理器中包括和执行的模块/部件可以是例如过程、应用程序、虚拟机和/或某种其他硬件或软件模块/部件。处理器可以是配置为运行和/或执行这些模块/部件的任何合适的处理器。处理器可以是配置为运行和/或执行一组指令或代码的任何合适的处理装置。例如，处理器可以是通用处理器、中央处理单元(CPU)、加速处理单元(APU)、现场可编程门阵列(FPGA)、专用集成电路(ASIC)、数字信号处理器(DSP)等。

在此描述的一些实例涉及具有非暂时性计算机可读介质(也可以称为非暂时性处理器可读介质)的计算机存储产品，其上具有用于执行各种计算机实现的操作的指令或计算机代码。计算机可读介质(或处理器可读介质)在其本身不包括暂时性传播信号(例如，在比如空间或电缆等传输介质上承载信息的传播电磁波)的意义上是非暂时性的。介质和计算机代码(也可以称为代码)可以是设计和构造用于一个或多个特定目的的那些介质和计算机代码。非暂时性计算机可读介质的示例包括但不限于：比如硬盘、软盘和磁带等磁存储介质；比如光盘/数字视频光盘(CD/DVD)、光盘只读存储器(CD-ROM)和全息装置等光盘存储介质；比如光盘等磁光存储介质；载波信号处理模块；以及专门配置为存储和执行程序代码的硬件装置，比如专用集成电路(ASIC)、可编程逻辑装置(PLD)、只读存储器(ROM)和随机存取存储器(RAM)装置。在此描述的其他实例涉及计算机程序产品，该计算机程序产品可以包括例如在此讨论的指令和/或计算机代码。

计算机代码的示例包括但不限于微代码或微指令、机器指令(比如由编译器产生)、用于产生web服务的代码，以及包含由计算机使用注释器执行的更高级指令的文件。例如，可以使用Java、C++、.NET或其他编程语言(例如，面向对象的编程语言)和开发工具来实现实例。计算机代码的附加示例包括但不限于控制信号、加密代码和压缩代码。

多深度重建

到目前为止，讨论主要涉及在单个焦深处重建图像。然而，全息成像相对于常规显微术的主要优点之一是有可能重建样本的整个3维体积而不仅仅是一个焦深处的单个图像。一种可能性是使用上一节中描述的算法独立地重建多个图像，同时改变焦深，但不幸的是，如果样本包含多个焦深处的对象，则来自离焦对象的衍射图样将破坏在任何给定焦深处的重建。另外，即使在仅需要单个焦深处的一个图像的情况下，仍然需要确定正确的焦深，这手动完成的话可能是繁琐的。为了解决这些问题，本节描述了前一节中描述的模型的扩展，以重建样本的3维体积。

返回到(7)中的单焦深模型，所需的唯一修改是代替重建单个图像

重建一系列图像X_i，i＝1，...，D，其中每个X_i图像对应于指定深度z[i]处的图像。更正式地，在给定期望重建深度的向量

的情况下，则可以使用以下模型找到图像

这个模型基本上与(7)中使用的用于在单个焦深处重建图像的模型相同，但是进一步扩展到完整的3维体积。不幸的是，通过将模型扩展到完整的3维体积，不再可能像(8)中那样重新表述目标函数以导出X的闭合形式更新。相反，可以使用混合算法，其中混合算法利用交替最小化来更新W和μ变量并且利用近端梯度下降步骤来更新X。另外，大多数中间变量留在傅立叶域中以促进计算。该方法的步骤在方法2中进行描述。图4示出针对在[650，1150]微米范围内的100个均匀间隔的焦深，作为指定焦深的函数的重建图像X_i的量值。应当注意，具有最大量值的图像深度对应于通过手动聚焦重建深度而发现的焦深，从而展示了利用所提出的方法在3维体积上重建图像如何稳健地且自动地恢复样本内对象的焦深。

参考图8，在另一个方面，本披露可以体现为用于在多个焦深处全息重建图像的方法200。方法200包括获得全息图以及使用稀疏性假设由全息图生成206一系列重建图像，其中每个图像对应于指定的焦深。例如，可以通过求解209满足条件|W|＝1的

来生成206重建图像X_i和相位W。如上所述，求解209方程(13)的方法是通过使用交替最小化来更新相位W和背景水平μ并且使用近端梯度下降步骤来更新X。

在特定实施例中，通过设定212重建图像矩阵X_i(其中i＝1，...，D)、估计相位W和背景水平μ的初始值来生成206重建图像X和相位W。计算215传递函数T(z[i])的傅里叶变换(其中z[i]是焦深的向量)，设定

并且R＝G-S。通过计算218为μ＝(G[1，1]-S[1，1])/(mn).来更新背景水平。通过计算

来更新相位。通过

来计算224更新的G，并且通过R＝G-S计算227更新的R。根据

和

针对图像矩阵中的每个图像来计算230更新的图像。通过

求解计算233更新S。针对μ、W、G、R、X_i和S中的每一个，重复236更新步骤以确定重建图像。可以将这些步骤重复236达预定的次数。可以通过执行以上方法200直到重建图像具有令人期望的质量来手动地确定次数。在另一个实施例中，重复236这些步骤，直到μ，W和/或X_i的更新值的变化小于迭代中的阈值。例如，可以重复236这些步骤，直到在这些步骤的迭代中，重建图像的更新值的变化小于10％、5％、3％、2％、或1％或任何其他百分比值。阈值可以是百分比变化或实际值。可以基于例如与执行附加迭代的成本(例如，计算时间等)平衡的重建图像的质量来选择阈值。

尽管已经关于一个或多个特定实施例描述了本披露，但是应当理解，在不脱离本披露的精神和范围的情况下，可以做出本披露的其他实施例。以下是旨在说明本披露的实施例的非限制性示例权利要求。

Claims (14)

1.一种用于样本的无透镜成像的系统，该系统包括：

图像传感器；

相干光源，该相干光源被配置为照射样本；

处理器，该处理器与图像传感器通信，该处理器被编程用于：

操作该图像传感器以获得由相干光源照射的样本的全息图H；

使用重建图像X中的对象的稀疏性假设，由该全息图生成距图像传感器的距离z处的样本的重建图像X和相位W；以及

其中，通过求解满足条件

来生成距图像传感器的距离z处的该重建图像X和相位W，其中T(z)是对在z处穿过样本的光衍射进行建模的传递函数，μ是非零背景的标量解释，1是一的矩阵，λ是稀疏性参数，算子⊙表示按条目逐个作出的乘积，算子*表示二维卷积，和F范数是弗罗比尼斯范数。

2.如权利要求1所述的系统，其中，该图像传感器是有源像素传感器。

3.如权利要求1所述的系统，其中，该图像传感器是CCD。

4.如权利要求1所述的系统，其中，该图像传感器是CMOS有源像素传感器。

5.一种用于样本的图像的全息重建的方法，包括：

使用无透镜成像来获得样本的全息图H；

使用重建图像X中的对象的稀疏性假设，由该全息图生成距图像传感器的距离z处的样本的重建图像X和相位W；以及

其中，通过求解满足条件|W|＝1的

来生成距图像传感器的距离z处的该重建图像X和相位W，其中T(z)是对在z处穿过样本的光衍射进行建模的传递函数，μ是非零背景的标量解释，1是一的矩阵，λ是稀疏性参数，算子⊙表示按条目逐个作出的乘积，算子*表示二维卷积，和F范数是弗罗比尼斯范数。

6.如权利要求5所述的方法，其中，通过使用交替最小化更新W、μ和X来求解

7.如权利要求5所述的方法，其中，通过以下方式来求解

(a)设定该重建图像X、估计相位W和与背景的拟合度μ的初始值；

(b)计算传递函数T(z)的傅立叶变换，并且设定

(c)计算

(f)计算X＝SFT_λ(Q-μexp(-izk))；并且

(g)重复步骤(c)-(f)以确定重建图像X。

8.如权利要求7所述的方法，其中，将步骤(c)-(f)重复预定次数。

9.如权利要求7所述的方法，其中，这些初始值是X＝0、W＝1和μ＝0。

10.如权利要求5所述的方法，其中，T(z)被定义为使得：

其中k＝2π/γ是用于获取该全息图的具有波长γ的光的波数，并且(k_x，k_y)分别表示水平方向和竖直方向上的空间频率。

11.如权利要求5所述的方法，该方法包括：

使用稀疏性假设，由该全息图生成一系列重建图像X_i和相位W，其中，每个图像对应于距图像传感器指定距离z[i]处的图像，其中，所述重建图像X形成一系列重建图像X_i中的重建图像。

12.如权利要求9所述的方法，其中，通过求解满足条件|W|＝1的

来生成这些重建图像X_i和相位W，

其中T(z[i])是对在z[i]处穿过样本的光衍射进行建模的传递函数，μ是是非零背景的标量解释，1是一的矩阵，D是一系列重建图像X_i中的重建图像的数量，λ是稀疏性参数，算子⊙表示按条目逐个作出的乘积，算子*表示二维卷积，和F范数是弗罗比尼斯范数。

13.如权利要求12所述的方法，其中，通过使用交替最小化更新W和μ并且使用近端梯度下降步骤更新X来求解满足条件|W|＝1的

14.如权利要求13所述的方法，其中，通过以下方式求解满足条件|W|＝1的

(a)将重建图像矩阵X_i、估计相位W和与背景参数的拟合度μ设定为相应的初始值，其中i＝1，...，D；

(b)计算传递函数T(z[i])的傅立叶变换，并且设定

并且R＝G-S；

(c)计算

(d)计算

(e)计算

(f)计算R＝G-S；

(g)计算R[1，1]＝R[1，1]-mnμ；

(h)针对每个图像，计算

知

(i)计算

并且

(j)重复步骤(b)-(i)以获得重建图像X_i。

Images