CN101865673A - 一种微观光场采集与三维重建方法及装置 - Google Patents

Abstract

本发明的实施例提出了一种微观光场采集与三维重建方法，包括通过控制微观样本沿光轴方向的微量移动，采集并记录微观样本的光场信息以得到图像栈；根据得到的图像栈及所述图像栈对应的Poisson模型构建统计理论模型，并根据压缩感知理论模型得到微观样本的梯度；和根据统计理论模型和压缩感知理论模型建立联合模型，并对联合模型进行迭代求解以获得微观样本的三维结构。本发明的实施例还提出了一种微观光场采集与三维重建装置，包括微观光场采集模块和三维重建模块。本发明实施例提供的方法和装置普适性更强，与现有方法相比效果更显著，并且操作更灵活。

Description

一种微观光场采集与三维重建方法及装置

技术领域

本发明涉及计算机应用技术和信号处理领域，具体涉及一种微观光场采集与三维重建方法及装置。

背景技术

以三维形式观察、分析和理解细胞等微观样本的组织结构是国内外专家学者的研究热点，微观光场采集与三维重建可以实现微观样本的三维成像及可视化。

目前，可以实现微观样本三维成像及可视化的光学成像系统主要包括体式显微成像系统、3D-X射线显微成像系统和共焦显微成像系统。其中，体式显微成像系统利用双目镜筒中左右两光束不并行产生的体式角，为左右两眼提供具有立体感的图像，该种显微成像系统产生的放大倍率不如常规显微成像系统，主要用于工业用途，且由于集中适应性冲突(convergence-accommodation conflict)而容易引起人体视觉疲劳。

3D-X射线显微成像系统需要使用X射线做光源，对样本损伤很大，且由于需要样本在工作台上旋转，采集几百到上千幅图像用于三维建模，一般需要半小时到一小时左右，再加上后期三维建模时间，整体耗时很长。

共焦显微成像系统是目前最为普遍使用的科研型显微镜，其成像系统利用物镜和探测器之间的针孔装置，确保只有聚焦平面的光子才能通过针孔，顺利到达探测器，极大减少了非聚焦平面的光线对所成图像的模糊影响。由于共焦装置具有很强的层析能力，因而可以直接获得清晰的三维微观样本图像，并且可以检测到样本内部很深的位置，尤其适合检测晶胚等厚组织样本。

但是共焦显微成像系统有很多限制：

1)价格及其昂贵，一般为几十万美元。

2)针孔装置使得探测器探测到的光子数量及其微少，因而图像对比度远不如宽场显微成像系统。

3)针孔装置的使用决定每次只能对样本的一个点进行聚焦成像，即使有更加先进的针孔装置出现，可以同时对多个点进行成像，效率也同样非常低下。

4)共焦显微成像系统以激光作为光源，会使承载在样本中的荧光探针发生光漂白和光裂解，且严重影响活体样本的存活能力。

针对现有三维显微成像系统难以兼顾图像高质量、样本无损伤、处理高效率、设备低成本和人体舒适感等多种需求的现状，设计并实现能同时满足这些需求的三维显微成像系统，具有非常重要的意义。

发明内容

本发明的目的旨在至少解决上述技术缺陷之一，特别针对利用微观光场的透射属性和宽场显微光学成像系统，多截面高密度采集微观光场，获得二维图像栈，利用三维重建方法，再现微观样本的三维结构本，提出了一种微观光场采集与三维重建方法及装置。

为实现上述目的，本发明实施例的一方面公开了一种微观光场采集与三维重建方法，包括如下步骤：

通过控制微观样本沿光轴方向的微量移动，采集并记录所述微观样本的光场信息以得到图像栈；

根据得到的所述图像栈及所述图像栈对应的Poisson模型构建统计理论模型，并根据微观样本的梯度场稀疏特性构建压缩感知理论模型；

根据所述统计理论模型和压缩感知理论模型建立联合模型，并利用所述联合模型得到微观样本重建公式，对该公式进行迭代求解以获得微观样本的三维结构。

本发明实施例的另一方面公开了一种微观光场采集与三维重建装置，包括微观光场采集模块和三维重建模块。

其中，所述微观光场采集模块，用于通过控制微观样本沿光轴方向的微量移动，采集并记录所述微观样本的光场信息以得到图像栈，所述微观光场采集模块包括工作台、宽场显微镜、载物台运动控制器、光源和制冷数码电荷耦合元件；

所述三维重建模块，用于根据得到的所述图像栈及所述图像栈对应的Poisson模型构建统计理论模型，并根据微观样本的梯度场稀疏特性构建压缩感知理论模型；根据所述统计理论模型和压缩感知理论模型建立联合模型，并利用所述联合模型得到微观样本重建公式，对该公式进行迭代求解以获得微观样本的三维结构。

本发明实施例提供的微观光场采集与三维重建方法及装置与其它直接显示样本三维结构的显微光学成像系统包括体式显微成像系统、3D-X射线成像系统和共焦成像系统相比，具有如下优点：

1)普适性更强。本发明提供的方法使用宽场显微成像系统采集微观光场，设备结构相对简单。例如无需共焦成像系统的针孔装置和激光光源，以及3D-X射线显微成像系统的X射线光源。研究者可以直接观测到微观样本的全局图像，而不是像共焦系统由局部扫描组合成全局图像，效率更高。激光光源和X射线光源一方面增加了实验条件的苛刻性，另一方面更容易引起对光照敏感的活体样本受到光损伤。

而宽场显微成像系统使用普通光源，且拍摄图像时间短，有效地克服了光漂白和光裂解，非常适合在极低光强下获得样本图像，例如承载着荧光探针的活体细胞。特别是在检测自发光样本时，利用激光光源的共焦显微成像系统根本无法成像，只能利用宽场显微成像系统获取样本图像栈，然后利用反卷积技术重建样本的三维结构。

同时，宽场显微成像系统和计算机价格低廉，因此更加适合在相关研究领域内广泛使用。因此，宽场显微成像系统无论在使用范围，还是在应用环境上，都具有更强的普适性。

2)效果更显著。本发明提供的方法对微观样本的三维重建是通过反卷积技术来实现的，反卷积可以进一步提高图像细节信息的清晰度和对比度，同时可以去除噪声，提高分辨率，尤其是轴向分辨率，而轴向分辨率的优良表现是昂贵的光学仪器如共焦或是双光子激发显微成像系统等的主要优势。经研究调查表明，宽场显微成像系统利用反卷积方法能够获得相当于甚至有时优于共焦显微成像系统成像的结果。

3)操作更灵活。其它三维光学成像系统在实验时需要连续完成，尤其对于共焦成像和X射线成像，耗时很长。而本发明提供的方法可以在宽场成像系统获取图像栈以后，随个人意志进行计算机后期处理工作，操作更灵活。同时该方法通过对微观样本进行三维建模仿真，可以满足操作人员随意调整视角进行观察的需要，观察更方便，实用性也更强。

本发明提出的上述方案，对现有系统的改动很小，不会影响系统的兼容性，而且实现简单、高效。

本发明附加的方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明上述的和/或附加的方面和优点从下面结合附图对实施例的描述中将变得明显和容易理解，其中：

图1为根据本发明实施例的微观光场采集与三维重建方法的流程框图；

图2为图1中微观光场采集与三维重建方法的流程示意图；

图3为根据本发明实施例的变步长采集微观光场示意图；

图4为根据本发明实施例的三维重建算法设计流程图；

图5为根据本发明实施例的微观样本f和点扩散函数h_k迭代过程示意图；

图6为根据本发明实施例的微观光场采集与三维重建装置的结构框图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

为了实现本发明之目的，本发明公开了一种微观光场采集与三维重建方法。该方法利用微观光场的透射属性和宽场显微光学成像系统，多截面高密度采集微观光场信息，获得二维图像栈。根据所获得的二位图像栈，在研究并估计点扩散模型基础上，设计并实现同时满足上述需求的三维重建方法，以真实再现微观样本的三维结构本。

结合图1所示，该微观光场采集与三维重建方法包括如下步骤：

S101：采集微观光场信息。

首先对微观光场进行说明，微观光场是微观折射光场的简称，指微观场景下任意入射光条件下的通用折射光场。由于显微光学成像系统是单一垂直视角，不会产生视差，无法使用宏观多视角的采集方法来采集光场。同时由于微观光场以透射光为主，穿越样本内部到达探测器，所以可以采用多截面的方式采集微观光场信息。

根据以上分析，在本步骤中，采集微观光场信息包括：通过控制微观样本沿光轴方向的微量移动，采集并记录微观样本的光场信息得到图像栈。

具体的说，结合图2所示首先利用光源向微观样本提供照明，然后控制微观样本沿光轴方向的微量移动。其中，微观样本放置于显微镜载物台上，通过载物台运动控制器控制显微镜载物台沿光轴方向的微量移动。

在理论上，图像栈中的图像越多，获取的样本信息就越完备，复原的三维结构也越接近真实状态。但是在实际应用中，过多的图像会导致计算复杂度大量增加，耗时长，尤其对于活体样本，需要在极短的时间内得到所有光场信息，因此图像数目不宜过多。然而，图像栈中二维图像数量要符合乃奎斯特采样定理，才能保证无失真的复原出样本三维结构。

为了解决上述矛盾，可以采用变步长的方式移动载物台以达到获取二维图像栈的目的。具体的说，对于复杂或是感兴趣的图像部分，选取的步长略小，对于简单或是不感兴趣的图像部分，选取的步长略大。通过上述方式，可以获取足够的光场信息用于三维重建。

为准确完备的采集微观样本微观光场，在本实施例中，采用宽场显微光学成像系统，利用变步长移动载物台的变焦方式高密度采集样本微观光场。

利用变步长方式沿光轴方向移动载物台，可以得到一系列微观样本不同焦平面的2D图像，即图像栈g。根据采集得到的图像栈g，利用制冷数码CCD(Charge-coupled Device，电荷耦合元件)记录微观光场信息。

图3示出了采用变步长方式的微观光场采集的示意图。如图3中所示，根据图像的复杂程度以及用户对该部分图形的兴趣程度，载物台沿光轴方向的微量移动步长逐渐变化，即微观样本沿光轴方向的微量移动步长逐渐变化。对于复杂或是用户感兴趣的图像部分，选取的步长略小；对于简单或是用户不感兴趣的图像部分，选取的步长略大。

S102：根据上述图像栈g，重建三维微观样本。

由于显微镜的探测器本质上是一个光子计数器，尤其是在低光照条件下，光子数一般非常少。在检测到的光子数目中，统计变量可以被Poisson噪声模型非常好的描述。另外，Poisson噪声模型可以得到更快的样本函数的收敛速度。并且Poisson噪声模型可以获得比Gaussian噪声模型更好的反卷积结果。因此，本实施例是在Poisson噪声模型的基础上进行实施。

具体的说，图4示出了三维重建微观样本的算法设计流程图，包括如下步骤：

A1：建立统计理论模型；

对于沿光轴z轴获取的一系列2D图像切片{g₁，g₂，...，g_k}，即图像栈g的图像切片，其对应的Poisson模型联合概率密度函数为：

$p\left(g_k\right|f,h_k)=\underset{k=1}{\overset{K}{\mathrm{\Π}}}\frac{{(f\⊗h_k)}^{g_k}}{g_k!}\exp (-(f\⊗h_k))---\left(1\right).$

因此，该联合密度函数的对数似然函数如下式：

$L(f,h_k)=-\underset{k}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}{h}_{k}f+\underset{k}{\mathrm{\Σ}}\left(g_k\ln \underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}{h}_{k}f\right)---\left(2\right),$

在计算联合密度函数的对数似然函数时，由于ln(g_k！)为常量所以省略。

其中，f为微观样本，K为图像栈g中的图像切片数，h_k为第k张图像栈的图像切片的点扩散函数，g_k为图像栈的2D图像切片。

A2：根据微观样本的梯度场稀疏特性构建压缩感知理论模型；

由于图像梯度

可以理解为图像g_k沿x方向与滤波器[-1 1]的卷积以及沿y方向与滤波器[-1 1]^T的卷积。设m＝[-1 1]，则

${\▿g}_{\mathrm{kx}}=g_k\⊗m=f\⊗h\⊗m=f\⊗m\⊗h={\▿f}_x\⊗h.$

其中，h为点扩散函数。由上式可以看出，微观样本f的梯度可以通过观测图像的梯度和PSF(point-spread function，点散布函数)进行复原。

由于图像梯度是稀疏的，所以根据压缩感知(compressed sensing)理论，可以通过求解如下优化问题复原微观样本梯度：

$\min {\left|\right|\▿f\left|\right|}_{l_1}s.t.{\left|\right|\▿f\⊗h_k\text{-}{\▿g}_k\left|\right|}_{l_2}={\mathrm{\σ}}^2---\left(3\right)$

其中，

为第k张图像切片的梯度场，

为微观样本的梯度，σ²为噪声方差，l₁和l₂分别表示1范数和2范数。

A3：建立联合模型。

根据统计理论模型得到的联合密度函数的对数似然函数和压缩感知模型中复原微观样本梯度，即结合式(2)和(3)，通过正则法得到最大似然估计函数：

$L_{\mathrm{cs}}(f,h_k)=L(f,h_k)-\mathrm{\α}\underset{k}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}|\▿f|-$

$\mathrm{\λ}\underset{k}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}({({\▿f}_x\⊗h_k-{\▿g}_{\mathrm{kx}})}^2+{({\▿f}_y\⊗h_k-{\▿g}_{\mathrm{ky}})}^2-{\mathrm{\σ}}^2)---\left(4\right)$

α＞0，用于权衡复原图像f的光滑性和对观测图像g的数据依赖性；λ为拉格朗日乘子。

根据上述建立的联合模型，计算f和h的迭代估计模型。具体的说，利用联合模型对f求导并令其等于0，推导出的微观样本重建公式，即f迭代估计模型。利用联合模型对h求导并令其等于0，推导出的点扩散函数重建公式，即h迭代估计模型。

然后，利用已得到的f和h的迭代估计模型，设置初始值。图5示出了微观样本f和点扩散函数h迭代示意图。结合图5所示，设置初始值包括：设置f和h的初始值。其中，点扩散函数h的初始值可以通过理论方法得到。微观样本f的初始值，通过将图像栈的数据g设置为f的初始值得到。

初始值设置成功后，对f和h的迭代估计模型并进行盲反卷积，当判断f和h的估计值符合视觉要求时，输出微观样本f的估计值和点扩散函数h的估计值。

本实施例采用盲反卷积重建样本3D结构。在三维反卷积过程中，根据成像模型

(忽略背景和噪声)可知，除了已知观测数据g，即图像栈外，还需要确定点扩散函数h，才能有效的复原微观样本f的三维结构。而盲反卷积方法进行微观样本的三维重建时，可以在估计出微观样本f的同时，估计出点扩散函数h。

具体的说，f和h的迭代估计模型并进行盲反卷积，当判断f和h符合视觉要求时，输出微观样本f的估计值和点扩散函数h的估计值，包括如下步骤：

B1：根据已设置的点扩散函数h、微观样本f的初始值和图像栈g估计新的微观样本；

B2：对步骤B1中新的微观样本的估计值设置灰度值非负限制和边界限制。根据上述约束条件的限制结果，得到新的微观样本的估计值，更新微观样本的估计值；

具体的说，灰度值非负限制：在算法迭代过程中，由于傅里叶变换或是减运算等操作，会导致像素灰度值出现负数的情况，这在实际中是不可能的。因此在迭代中，采用如强制负灰度值为零的限制方法，以保证最后结果符合实际情况。

边界限制：(i，j)表示图像x＝i，y＝j的位置坐标，1≤i≤M，1≤j≤N，M和N为图像像素的行数和列数。

对于图像边界点有：

$\left\{\begin{array}{c}u(1,j)=u(2,j)\\ u(M,j)=u(M-1,j)\\ u(i,1)=u(i,2)\\ u(i,N)=u(i,N-1)\end{array}\right.,$

满足上式要求的图像坐标，符合边界限制的条件。根据上述约束条件的限制结果，可以得到新的微观样本的估计值，

B3：判断步骤B2中新的微观样本的估计值是否符合视觉要求，如果符合视觉要求，则输出新的微观样本的估计值和点扩散函数的估计值。

如果不符合视觉要求，则根据点扩散函数h、新的微观样本的估计值和图像栈g重新估计新的点扩散函数。并对新的点扩散函数的估计值设置灰度值非负限制和边界限制，得到新的点扩散函数的估计值，更新点扩散函数的估计值，执行步骤B1。

其中，视觉要求是指就是达到用户满意的程度，即用户认为该微观样本的重建三维结构达到其要求。总之，通过添加灰度值的非负限制和边界限制，不断更新迭代，得到满意的微观样本f估计值和点扩散函数h估计值。由此，可以达到三维重建微观样本的目的。

本发明实施例还公开了一种微观光场采集与三维重建装置。图6示出了该微观光场采集与三维重建装置的结构框图。如图6中所示，该装置100包括微观光场采集模块110和三维重建模块120。

具体的说，微观光场采集模块110，用于采集微观光场信息，包括工作台、宽场显微镜、载物台运动控制器、光源和制冷数码电荷耦合元件(制冷数码CCD)。其中，宽场显微镜、载物台运动控制器和光源均固定于工作台上。首先利用光源向微观样本提供照明，然后控制微观样本沿光轴方向的微量移动。将微观样本放置于显微镜载物台上，通过载物台运动控制器控制显微镜载物台沿光轴方向的微量移动。

在理论上，图像栈中的图像越多，获取的样本信息就越完备，复原的三维结构也越接近真实状态。但是在实际应用中，过多的图像会导致计算复杂度大量增加，耗时长，尤其对于活体样本，需要在极短的时间内得到所有光场信息，因此图像数目不宜过多。然而，图像栈中二维图像数量要符合乃奎斯特采样定理，才能保证无失真的复原出样本三维结构。

为了解决上述矛盾，可以采用变步长的方式移动载物台以达到获取二维图像栈的目的。具体的说，对于复杂或是感兴趣的图像部分，选取的步长略小，对于简单或是不感兴趣的图像部分，选取的步长略大。通过上述方式，可以获取足够的光场信息用于三维重建。

为准确完备的采集微观样本微观光场，在本实施例中，微观光场采集模块110采用宽场显微光学成像系统，利用变步长移动载物台的变焦方式高密度采集样本微观光场，由此得到一系列微观样本不同焦平面的2D图像，称为图像栈g。根据采集得到的图像栈g，利用制冷数码CCD记录微观光场信息。

结合图3所示，根据图像的复杂程度以及用户对该部分图形的兴趣程度，载物台沿光轴方向的微量移动步长逐渐变化，即微观样本沿光轴方向的微量移动步长逐渐变化。对于复杂或是用户感兴趣的图像部分，选取的步长略小；对于简单或是用户不感兴趣的图像部分，选取的步长略大。

微观光场采集与三维重建装置100还包括三维重建模块120，该模块根据来自微观光场采集模块110的图像栈g，三维重建微观样本。

由于，Poisson噪声模型可以得到更快的样本函数的收敛速度。并且Poisson噪声模型可以获得比Gaussian噪声模型更好的反卷积结果。因此，本实施例是在Poisson噪声模型的基础上进行实施。

如图4所示，三维重建模块120进行三维重建微观样本的算法设计，包括：

A1：建立统计理论模型；

对于微观光场采集模块110沿光轴z轴获取的一系列2D图像切片{g₁，g₂，...，g_k}，即图像栈g的图像切片，其对应的Poisson模型联合概率密度函数为：

$p\left(g_k\right|f,h_k)=\underset{k=1}{\overset{K}{\mathrm{\Π}}}\frac{{(f\⊗h_k)}^{g_k}}{g_k!}\exp (-(f\⊗h_k))---\left(1\right).$

因此，该联合密度函数的对数似然函数如下式：

$L(f,h_k)=-\underset{k}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}{h}_{k}f+\underset{k}{\mathrm{\Σ}}\left(g_k\ln \underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}{h}_{k}f\right)---\left(2\right),$

在计算联合密度函数的对数似然函数时，由于ln(g_k！)为常量所以省略。

其中，f为微观样本，K为图像栈g中的图像切片数，h_k为第k张图像栈的图像切片的点扩散函数，g_k为图像栈的2D图像切片。

A2：微观样本的梯度场稀疏特性构建压缩感知理论模型；

由于图像梯度

可以理解为图像g_k沿x方向与滤波器[-1 1]的卷积以及沿y方向与滤波器[-1 1]^T的卷积。设m＝[-1 1]，则

${\▿g}_{\mathrm{kx}}=g_k\⊗m=f\⊗h\⊗m=f\⊗m\⊗h={\▿f}_x\⊗h.$

其中，h为点扩散函数。由上式可以看出，微观样本f的梯度可以通过观测图像的梯度和PSF(point-spread function，点散布函数)进行复原。

由于图像梯度是稀疏的，所以根据压缩感知(compressed sensing)理论，可以通过求解如下优化问题复原微观样本梯度：

$\min {\left|\right|\▿f\left|\right|}_{l_1}s.t.{\left|\right|\▿f\⊗h_k\text{-}{\▿g}_k\left|\right|}_{l_2}={\mathrm{\σ}}^2---\left(3\right)$

其中，

为第k张图像切片的梯度场，

为微观样本的梯度，σ²为噪声方差，l₁和l₂分别表示1范数和2范数。

A3：建立联合模型。

根据统计理论模型得到的联合密度函数的对数似然函数和压缩感知模型，即结合式(2)和(3)，通过正则法得到最大似然估计函数：

$L_{\mathrm{cs}}(f,h_k)=L(f,h_k)-\mathrm{\α}\underset{k}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}|\▿f|-$

$\mathrm{\λ}\underset{k}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}({({\▿f}_x\⊗h_k-{\▿g}_{\mathrm{kx}})}^2+{({\▿f}_y\⊗h_k-{\▿g}_{\mathrm{ky}})}^2-{\mathrm{\σ}}^2)---\left(4\right)$

α＞0，用于权衡复原图像f的光滑性和对观测图像g的数据依赖性；λ为拉格朗日乘子。

根据上述建立的联合模型，计算f和h的迭代估计模型。具体的说，根据上述建立的联合模型，计算f和h的迭代估计模型。具体的说，利用联合模型对f求导并令其等于0，推导出微观样本重建公式，即f迭代估计模型。利用联合模型对h求导并令其等于0，推导出点扩散函数重建公式，即h迭代估计模型。

然后，利用已得到的f和h的迭代估计模型，设置初始值。结合图5所示，设置初始值包括：设置f和h的初始值。其中，点扩散函数h的初始值可以通过理论方法得到。微观样本f的初始值，通过将图像栈的数据g设置为f的初始值得到。

初始值设置成功后，对f和h的迭代估计模型并进行盲反卷积，当判断f和h的估计值符合视觉要求时，输出微观样本f的估计值和点扩散函数h的估计值。

本实施例采用盲反卷积重建样本3D结构。在三维反卷积过程中，根据成像模型

(忽略背景和噪声)可知，除了已知观测数据g，即图像栈外，还需要确定点扩散函数h，才能有效的复原微观样本f的三维结构。而盲反卷积方法进行微观样本的三维重建时，可以在估计出微观样本f的同时，估计出点扩散函数h。

具体的说，三维重建模块120对f和h的迭代估计模型并进行盲反卷积，当判断f和h符合视觉要求时，输出微观样本f的估计值和点扩散函数h的估计值，包括如下步骤：

B1：根据已设置的点扩散函数h、微观样本f的初始值和图像栈g估计新的微观样本；

B2：对步骤B1中新的微观样本的估计值设置灰度值非负限制和边界限制。根据上述约束条件的限制结果，得到新的微观样本的估计值，更新微观样本的估计值；

具体的说，灰度值非负限制：在算法迭代过程中，由于傅里叶变换或是减运算等操作，会导致像素灰度值出现负数的情况，这在实际中是不可能的。因此在迭代中，采用如强制负灰度值为零的限制方法，以保证最后结果符合实际情况。

边界限制：(i，j)表示图像x＝i，y＝j的位置坐标，1≤i≤M，1≤j≤N，M和N为图像像素的行数和列数。

对于图像边界点有：

$\left\{\begin{array}{c}u(1,j)=u(2,j)\\ u(M,j)=u(M-1,j)\\ u(i,1)=u(i,2)\\ u(i,N)=u(i,N-1)\end{array}\right.,$

满足上式要求的图像坐标，符合边界限制的条件。根据上述约束条件的限制结果，可以得到新的微观样本的估计值，

B3：判断步骤B2中新的微观样本的估计值是否符合视觉要求，如果符合视觉要求，则输出新的微观样本的估计值和点扩散函数的估计值。

如果不符合视觉要求，则根据点扩散函数h、新的微观样本的估计值和图像栈g重新估计新的点扩散函数。并对新的点扩散函数的估计值设置灰度值非负限制和边界限制，得到新的点扩散函数的估计值，更新点扩散函数的估计值，执行步骤B1。

其中，视觉要求是指就是达到用户满意的程度，即用户认为该微观样本的重建三维结构达到其要求。总之，通过添加灰度值的非负限制和边界限制，不斯更新迭代，得到满意的微观样本f估计值和点扩散函数h估计值。由此，可以达到重建三维微观样本的目的。

本发明实施例提供的微观光场采集与三维重建方法及装置与其它直接显示样本三维结构的显微光学成像系统包括体式显微成像系统、3D-X射线成像系统和共焦成像系统相比，具有如下优点：

1)普适性更强。本发明提供的方法使用宽场显微成像系统采集微观光场，设备结构相对简单。例如无需共焦成像系统的针孔装置和激光光源，以及3D-X射线显微成像系统的X射线光源。研究者可以直接观测到微观样本的全局图像，而不是像共焦系统由局部扫描组合成全局图像，效率更高。激光光源和X射线光源一方面增加了实验条件的苛刻性，另一方面更容易引起对光照敏感的活体样本受到光损伤。

而宽场显微成像系统使用普通光源，且拍摄图像时间短，有效地克服了光漂白和光裂解，非常适合在极低光强下获得样本图像，例如承载着荧光探针的活体细胞。特别是在检测自发光样本时，利用激光光源的共焦显微成像系统根本无法成像，只能利用宽场显微成像系统获取样本图像栈，然后利用反卷积技术重建样本的三维结构。

同时，宽场显微成像系统和计算机价格低廉，因此更加适合在相关研究领域内广泛使用。因此，宽场显微成像系统无论在使用范围，还是在应用环境上，都具有更强的普适性。

2)效果更显著。本发明提供的方法对微观样本的三维重建是通过反卷积技术来实现的，反卷积可以进一步提高图像细节信息的清晰度和对比度，同时可以去除噪声，提高分辨率，尤其是轴向分辨率，而轴向分辨率的优良表现是昂贵的光学仪器如共焦或是双光子激发显微成像系统等的主要优势。经研究调查表明，宽场显微成像系统利用反卷积方法能够获得相当于甚至有时优于共焦显微成像系统成像的结果。

3)操作更灵活。其它三维光学成像系统在实验时需要连续完成，尤其对于共焦成像和X射线成像，耗时很长。而本发明提供的方法可以在宽场成像系统获取图像栈以后，随个人意志进行计算机后期处理工作，操作更灵活。同时该方法通过对微观样本进行三维建模仿真，可以满足操作人员随意调整视角进行观察的需要，观察更方便，实用性也更强。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (15)

1.一种微观光场采集与三维重建方法，其特征在于，该方法包括如下步骤：

通过控制微观样本沿光轴方向的微量移动，采集并记录所述微观样本的光场信息以得到图像栈；

根据得到的所述图像栈及所述图像栈对应的Poisson模型构建统计理论模型，并根据微观样本的梯度场稀疏特性构建压缩感知理论模型；

根据所述统计理论模型和压缩感知理论模型建立联合模型，并利用所述联合模型得到的微观样本重建公式，进行迭代求解以获得微观样本的三维结构。

2.如权利要求1所述的方法，其特征在于，所述图像栈为所述微观样本不同焦平面的2D图像。

3.如权利要求1所述的方法，其特征在于，所述统计理论模型为：

其中，K为图像栈的切片数，g_k为图像栈中的第k张图像切片，f为微观样本，h_k为第k张图像栈的图像切片的点扩散函数。

4.如权利要求3所述的方法，其特征在于，所述根据微观样本的梯度场稀疏特性构建的压缩感知理论模型：

$\min {\left|\right|\▿f\left|\right|}_{l_1}s.t.{\left|\right|\▿f\⊗h_k-{\▿g}_k\left|\right|}_{l_2}={\mathrm{\σ}}^2,$

其中，为第k张图像切片的梯度场，

为微观样本的梯度，σ²为噪声方差，l₁和l₂分别表示1范数和2范数。

5.如权利要求4所述的方法，其特征在于，所述联合模型为：

$L_{\mathrm{cs}}(f,h_k)=L(f,h_k)-\mathrm{\α}\underset{k}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}|\▿f|-\mathrm{\λ}\underset{k}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}({({\▿f}_x\⊗h_k-{\▿g}_{\mathrm{kx}})}^2+{({\▿f}_y\⊗h_k-{\▿g}_{\mathrm{ky}})}^2-{\mathrm{\σ}}^2),$

其中，α＞0，用于权衡复原图像f的光滑性和对观测图像g的数据依赖性，λ为拉格朗日乘子。

6.如权利要求5所述的方法，其特征在于，利用所述联合模型得到的微观样本重建公式，进行迭代求解以获得微观样本的三维结构，包括以下步骤：

利用所述联合模型分别对所述微观样本f和点扩散函数h进行求导，并令其等于0，得到微观样本重建公式和点扩散函数重建公式；

确定点扩散函数h的初始值，并将图像栈g设为微观样本f的初始值；

对所述f和h的迭代估计模型并进行盲反卷积，当判断所述f和h符合视觉要求时，输出微观样本f的估计值和点扩散函数h的估计值。

7.如权利要求1所述的方法，其特征在于，所述通过控制微观样本沿光轴方向的微量移动采用变步长方式。

8.一种微观光场采集与三维重建装置，其特征在于，包括微观光场采集模块和三维重建模块，

所述微观光场采集模块，用于通过控制微观样本沿光轴方向的微量移动，采集并记录所述微观样本的光场信息以得到图像栈，所述微观光场采集模块包括工作台、宽场显微镜、载物台运动控制器、光源和制冷数码电荷耦合元件；

所述三维重建模块，用于根据得到的所述图像栈及所述图像栈对应的Poisson模型构建统计理论模型，并根据微观样本的梯度场稀疏特性构建压缩感知理论模型；

根据所述统计理论模型和压缩感知理论模型建立联合模型，并利用所述联合模型得到的微观样本重建公式，进行迭代求解以获得微观样本的三维结构。

9.如权利要求8所述的装置，其特征在于，

所述工作台，用于固定所述宽场显微镜、载物台运动控制器和光源；

所述光源，用于向所述微观样本提供照明；

所述载物台运动控制器，用于控制所述微观样本沿光轴方向的微量移动，采集所述图像栈；

所述制冷数码电荷耦合元件，用于记录所述微观光场信息。

10.如权利要求8所述的装置，其特征在于，所述图像栈为所述微观样本不同焦平面的2D图像。

11.如权利要求8所述的装置，其特征在于，所述统计理论模型为：

其中，K为图像栈的切片数，g_k为图像栈中的第k张图像切片，f为微观样本，h_k为第k张图像栈的图像切片的点扩散函数。

12.如权利要求11所述的装置，其特征在于，所述根据微观样本的梯度场稀疏特性构建压缩感知理论模型：

$\min {\left|\right|\▿f\left|\right|}_{l_1}s.t.{\left|\right|\▿f\⊗h_k-{\▿g}_k\left|\right|}_{l_2}={\mathrm{\σ}}^2,$

其中，

为第k张图像切片的梯度场，

为微观样本的梯度，σ²为噪声方差，l₁和l₂分别表示1范数和2范数。

13.如权利要求12所述的装置，其特征在于，所述三维重建模块建立的联合模型为：

$L_{\mathrm{cs}}(f,h_k)=L(f,h_k)-\mathrm{\α}\underset{k}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}|\▿f|-\mathrm{\λ}\underset{k}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}\underset{x}{\mathrm{\Σ}}({({\▿f}_x\⊗h_k-{\▿g}_{\mathrm{kx}})}^2+{({\▿f}_y\⊗h_k-{\▿g}_{\mathrm{ky}})}^2-{\mathrm{\σ}}^2),$

其中，α＞0，用于权衡复原图像f的光滑性和对观测图像g的数据依赖性，λ为拉格朗日乘子。

14.如权利要求13所述的装置，其特征在于，所述三维重建模块利用所述联合模型得到的微观样本重建公式，进行迭代求解以获得微观样本的三维结构包括：

利用所述联合模型分别对所述微观样本f和点扩散函数h进行求导，并令其等于0，得到微观样本重建公式和点扩散函数重建公式；

确定点扩散函数h的初始值，并将图像栈g设为微观样本f的初始值；

对所述f和h的迭代估计模型并进行盲反卷积，当判断所述f和h符合视觉要求时，输出微观样本f的估计值和点扩散函数h的估计值。

15.如权利要求8所述的装置，其特征在于，所述通过控制微观样本沿光轴方向的微量移动采用变步长方式。

Images