一种图像重构方法及系统与流程

本申请涉及图像重构技术领域，特别涉及一种图像重构方法及系统。

背景技术：

作为光学显微镜的一种，荧光显微镜是现代生命科学研究中必不可少的研究工具，然而传统光学显微镜的分辨率受到衍射极限的限制。近年来，一系列突破衍射极限的超分辨显微成像技术相继出现，例如随机光学重构显微镜、光激活定位显微镜、荧光光激活定位显微镜、受激发射损耗显微镜、结构光照明显微镜等。结构光照明显微镜(structuredilluminationmicroscope,sim)是目前主流的超分辨成像技术之一。该技术采用加载了一定频率信息的结构光对样品进行照明(以下统称为照明光)，将传统显微镜无法获取的高频信息编码到采集图像中，再通过图像处理算法对样品信息进行重构，其横向分辨率可以达到传统分辨率的两倍。

超分辨重构算法是sim的核心技术之一，其算法主要步骤包括：图像预处理、频谱分离、计算照明光空间频率、参数修正、去卷积、频谱平移以及频谱融合。影响重构图像质量的因素有很多，例如：平移相位不准确而导致的频谱分离不彻底、照明光参数的求解误差、去卷积过程引入的噪声放大、点扩散函数误差、离焦光线等。实验系统中照明光参数主要有空间频率、方向、调制度和初始相位等参数，由于实验系统的不稳定，这些参数是无法预先知道的，只能通过一系列原始测量数据对这些参数进行求解，其求解精度对最终重构图像的质量有重要影响。

目前，照明光参数的求解算法主要有以下几种方法：文献1[phase-shiftestimationinsinusoidallyilluminatedimagesforlateralsuperresolution]通过分析原始数据在频域中的峰值点的位置及其相位来确定照明光的空间频率和相位，这种方法对于高频率的照明光和信噪较大的数据并不适用；文献2[surpassingthelateralresolutionlimitbyafactoroftwousingstructuredilluminationmicroscopy]通过计算平移之后与之前的组分进行互相关，进而求解得到照明光参数，但此方法计算量大，耗时较长，不适用于sim数据的快速处理；文献3[phaseoptimisationforsturcturedilluminationmicroscopy]在文献2的基础上增加了迭代优化算法，加快了算法的计算效率，但算法较为复杂，计算精度与初始值有关，初始值不当会引起算法无法收敛至准确值。

技术实现要素：

本申请提供了一种图像重构方法及系统，旨在至少在一定程度上解决现有技术中的上述技术问题之一。

为了解决上述问题，本申请提供了如下技术方案：

一种图像重构方法，包括：

步骤a：对图像进行频谱分离，得到所述图像的频谱分量；

步骤b：对所述频谱分量进行峰值点搜索，获取照明光的方向和空间频率；

步骤c：通过线性拟合算法计算所述照明光的调制度和初始相位，通过所述照明光的调制度和初始相位对所述频谱分量进行修正，将所述频谱分量的系数进行归一化；

步骤d：对所述归一化的频谱分量进行去卷积和平移处理，并将所述去卷积和平移处理后的频谱分量进行频谱组合，得到所述图像的重构结果。

本申请实施例采取的技术方案还包括：所述步骤a还包括：对所述图像进行预处理；所述预处理包括：对所述图像进行光强修正，并对所述图像的边缘进行模糊化处理。

本申请实施例采取的技术方案还包括：在所述步骤b中，所述获取照明光的方向和空间频率具体包括：根据照明光的空间频率初始值设定滤波函数，对所述频谱分量进行滤波，对所述滤波后的频谱分量进行峰值点搜索，确定滤波后的频谱分量的最大值，并在所述最大值附近±1个像素的范围内，采用质心法确定所述峰值点位置，并精确到亚像素，进而得到所述照明光的方向和空间频率。

本申请实施例采取的技术方案还包括：在所述步骤c中，所述通过照明光的调制度和初始相位对频谱分量进行修正具体包括：

步骤c1：中心组分的匹配：中心组分的表达式为：

上述公式中，是光学传递函数，d＝1,2,3；n＝0,1,2；m＝-1,0,1，是频谱分量，假设b＝s*a，则s＝c1,0/cd,0，通过对表达式a和b进行线性拟合得到参数s，并将所有其他组分用参数s进行修正，即：

其中表示将d1用d2替换；

步骤c2：边缘组分与中心组分进行匹配：根据傅里叶频移定理：

首先将边缘频谱变换到实域，乘以相移因子再变换至频域，得到移位之后的频谱：

同理，获得平移之后的otf，即得到表达式a和b：

假设b＝s*a，通过复数线性拟合得到参数边缘组分修正为：

本申请实施例采取的技术方案还包括：在所述步骤c中，所述通过线性拟合算法计算照明光的调制度和初始相位具体为：考虑到噪声和计算误差的存在，需要对表达式a和b内的数据进行筛选。首先计算表达式a和b的有效支持域，对于中心组分修正，所述支持域取|otf|>0.1的区域，对于边缘组分修正，所述有效支持域为和的重叠区域；然后由有效支持域内a和b的数据组成数据对(ai,bi)，并对数据对(ai,bi)进行进一步筛选，由于噪声的存在，|bi|与|ai|的比值分布范围较大，但大部分数值集中在一定区域内，所以选择比值|bi|/|ai|分布在其最大概率附近的数据组作为有效数据。最后对所述有效数据对进行复数线性拟合，得到参数s的估计，即

本申请实施例采取的另一技术方案为：一种图像重构系统，包括：

频谱分离模块：用于对图像进行频谱分离，得到所述图像的频谱分量；

照明光参数计算模块：用于对所述频谱分量进行峰值点搜索，获取照明光的方向和空间频率；

参数修正模块：用于通过线性拟合算法计算所述照明光的调制度和初始相位，通过所述照明光的调制度和初始相位对所述频谱分量进行修正，将所述频谱分量的系数进行归一化；

频谱平移模块：用于对所述归一化的频谱分量进行去卷积和平移处理；

频谱组合模块：用于将所述去卷积和平移处理后的频谱分量进行频谱组合，得到所述图像的重构结果。

本申请实施例采取的技术方案还包括图像预处理模块，所述图像预处理模块用于对图像进行预处理；所述预处理包括：对所述图像进行光强修正，并对所述图像的边缘进行模糊化处理。

本申请实施例采取的技术方案还包括：所述照明光参数计算模块获取照明光的方向和空间频率具体为：根据照明光的空间频率初始值设定滤波函数，对所述频谱分量进行滤波，对所述滤波后的频谱分量进行峰值点搜索，确定滤波后的频谱分量的最大值，并在所述最大值附近±1个像素的范围内，采用质心法确定所述峰值点位置，得到所述照明光的方向和空间频率。

本申请实施例采取的技术方案还包括：所述参数修正模块通过照明光的调制度和初始相位对频谱分量进行修正具体包括：

中心组分的匹配：中心组分的表达式为：

式中，是光学传递函数，d＝1,2,3；n＝0,1,2；m＝-1,0,1,是频谱分量，假设b＝s*a，则s＝c1,0/cd,0，通过对表达式a和b进行线性拟合得到参数s，并将所有其他组分用参数s进行修正，即：其中表示将d1用d2替换；

边缘组分与中心组分进行匹配：

将边缘频谱变换到实域，乘以相移因子再变换至频域，得到移位之后的频谱：

获得平移之后的otf，即得到表达式a和b：

假设b＝s*a，通过复数线性拟合得到参数边缘组分修正为：

本申请实施例采取的技术方案还包括：所述参数修正模块通过线性拟合算法计算照明光的调制度和初始相位具体为：计算表达式a和b的有效支持域，取|otf|>0.1的区域，所述有效支持域为和的重叠区域；由有效支持域内a和b的数据组成数据对(ai,bi)，并对数据对(ai,bi)进行数据筛选，得到有效数据对，所述有效数据对通过复数线性拟合进行参数s的估计，即

相对于现有技术，本申请实施例产生的有益效果在于：本申请实施例的图像重构方法及系统利用理想相位对预处理图像进行频谱分离，根据照明光空间频率初始值设定滤波函数，对分离频谱进行滤波，再进行峰值点搜索，并通过质心法将峰值点的定位精度提高至亚像素级，以确定照明光的空间频率。该方法尽可能地排除了其他干扰信息，在此基础上进行峰值点查找，可以提高空间频率的求解精度。采用简单而有效的线性拟合算法求解照明光的调制度和初始相位，并在此基础上增加数据筛选功能，避免了繁琐而耗时的互相关算法，同时又提高了求解精度；并利用计算得到的调制度和初始相位对频谱分量进行参数归一化，从而有效提高重构图像的质量。

附图说明

图1是本申请实施例的图像重构方法的流程图；

图2为sim超分辨原理图；

图3为本申请实施例的频谱分量滤波示意图；

图4是本申请实施例的图像重构系统的结构示意图；

图5为本申请实施例的图像重构结果示意图；其中，图5(a)表示无参数修正部分的图像重构结果；图5(b)表示有参数修正部分的图像重构结果；图5(c)表示宽场成像结果。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本申请，并不用于限定本申请。

请参阅图1，是本申请实施例的图像重构方法的流程图。本申请实施例的图像重构方法包括以下步骤：

步骤10：对sim成像系统采集到的原始图像进行预处理；

在步骤10中，sim成像系统的成像原理为：sim是一个线性移不变光学系统，图像是样品发射信息与系统psf(pointspreadfunction，点扩散函数)的卷积：

在公式(1)中，是被测的样品信息，是照明光，是系统psf。在线性sim中，照明光通常为正弦条纹，可表示为：

在公式(2)中，为照明光的空间频率矢量，表示了照明光的周期和方向，a表示照明光的调制度。结合公式(1)和(2)，图像的频域表达式为：

在公式(3)中，m＝-1,0,1，符号(～)表示时域变量的频域表达式，是光学传递函数(opticaltransferfunction，otf)。可见，每幅采集得到的原始图像都是由0和±1级频谱分量组合而成。

为了将三个频谱分量分离，需要改变照明光的相位，照明光的相位值为初始相位和平移相位之和，即初始相位是未知的，与照明光的绝对位置有关，平移相位一般取0～2π上等间隔分布的相位，即为了均匀的提高分辨率，一般采用相隔60度的三个照明光方向，如此，sim至少需要采集9幅原始图像。由于系统误差，每个方向上的照明光的调制度、初始相位和空间频率是不同的，公式(3)又可以表示为：

在公式(4)中，d表示照明光的方向，n表示平移相位，m表示频谱分量，d＝1,2,3；n＝0,1,2；m＝-1,0,1,是待分离的频谱分量。

sim超分辨原理如图2所示。图2(a)中，13表示传统宽场显微镜可观测到的样品频域信息，由于系统otf的支持域有限，支持域外的高频信息无法获取，因此分辨率受限。采用正弦结构光对样品进行照明，如公式(4)所示，sim获取的样品信息中含有高频组分，图2(b)中，14和15则表示某方向上被平移至正确位置的高频分量，可见，在此方向上，分辨率有所提高。为了均匀提高系统分辨率，照明光采用三个方向，将所有方向上的频谱分量叠加，就得到图2(c)中16所示的区域。相比于13，区域16所示的系统分辨率最大可提高2倍。

在图像采集过程中，由于相机曝光时间和照明光强的浮动，需要对测量序列中的所有原始图像进行光强修正，保证每幅图像的平均值相同。另一方面，重构算法主要在傅里叶频域进行，为了防止图像边缘突变引入的频域误差，需要对原始图像的边缘进行模糊，本申请实施例中，取原始图像边缘10个像素进行正弦模糊化处理。

步骤20：对预处理后的图像进行频谱分离，得到图像的频谱分量；

在步骤20中，频谱分离具体为：根据平移相位，得到分离矩阵，根据公式(5)分别对每个方向上的3幅原始图像进行分离，最终得到三个方向上的9个频谱分量

步骤30：根据频谱分量计算照明光的方向和空间频率；

在步骤30中，由频谱分量的表达式可知，样品的零频位于处，根据傅里叶变换性质，样品的零频值是时域图像上所有像素的光强之和，即对于sim系统，照明光的空间频率漂移较小，因此根据其空间频率初始值和角度初始值设定滤波函数，对频谱分量进行滤波，通过搜索确定滤波之后的数据的最大值，然后在最大值附近±1个像素的范围内，采用质心法将峰值点位置精确到亚像素，进而得到照明光的方向和空间频率。具体如图3所示，为本申请实施例的频谱分量滤波示意图。由于频谱的共轭性，只需要对m≥0的组分进行峰值点查找。该寻峰算法尽可能地排除了其他干扰信息，在此基础上进行峰值点查找，可以以提高空间频率的求解精度。

步骤40：获取sim成像系统的otf参数；

在步骤40中，otf在sim重构算法中是一个重要的参数，获取otf的方法包括两种：一种是采用荧光小球对sim成像系统进行标定，通过数据处理，得到psf，进而得到otf，这种方法可以反映真实的实验系统，但是由于荧光信号较弱，测量数据受噪声干扰较大，且荧光小球并非理想点光源，所以很难得到效果较好的otf。另一种方法是采用psf理论模型，目前常用的psf理论模型有born&wolf模型、richard&wolf模型、gibson&lanni模型等。考虑到商用高精度物镜的像差可以控制在很小的范围内，本申请实施例中使用的psf理论模型是richard&wolf模型。模型中使用到的光学参数是根据实测的psf进行确定，以保证模型与实际系统更为接近。

步骤50：根据otf对频谱分量进行参数修正，通过线性拟合算法计算照明光的调制度和初始相位，并通过照明光的调制度和初始相位将所有频谱分量的系数归一化为相同值；

在步骤50中，由公式(4)可知，每个频谱分量的系数是不同的，因此需要对其进行归一化，这样才能通过各个频谱的线性组合，获得正确的样品信息。本申请实施例对频谱分量的系数进行归一化的方式具体为：

(a)首先进行中心组分的匹配，由上述推导可知，中心组分(m＝0)的表达式为：

假设b＝s*a，则s＝c1,0/cd,0，参数s可以通过对表达式a和b进行线性拟合得到，之后将所有其他组分用参数s进行修正，即：

其中表示将d1用d2替换。

(b)边缘组分(m＝±1)与中心组分进行匹配，根据傅里叶频移定理：

将边缘频谱变换到实域，乘以相移因子然后再变换至频域，得到移位之后的频谱：

同样地，可以获得平移之后的otf，即这样得到表达式a和b：

这两个表达式在otf的重合区域应有相同的数值所以假设b＝s*a，通过复数线性拟合得到参数如此，边缘组分可以修正为：

对于m＝1和m＝-1，参数s是共轭的，因此只需要计算m＝1时的参数。此时所有频谱分量的参数都归一化为c1,0。

考虑到噪声和计算误差，需要对表达式a和b中的数据进行筛选，然后再进行复数线性拟合计算。首先，傅里叶相移定理会引入不可避免的计算误差，通过仿真分析证明该计算误差对绝对值较小的数据的相位影响是显著的，因此数据筛选的第一步是去掉绝对值过小的数据，即根据otf确定有效支持域。本申请实施例中，对于中心组分修正，所述支持域取|otf|>0.1的区域，对于边缘组分修正，所述有效支持域为和的重叠区域。然后由有效支持域内a和b的数据组成数据对(ai,bi)，虽然理论上数据bi与数据ai成正比，但是由于噪声的存在，|bi|与|ai|的比值分布范围较大，但大部分数值集中在一定区域内，所以选择比值|bi|/|ai|分布在其最大概率附近的数据组作为有效数据。最后经过筛选的有效数据对(ai,bi)通过复数线性拟合用于参数s的估计，即

本申请实施例中通过将初始相位和调制度作为一个整体参数进行求解，在线性拟合算法的基础上增加数据筛选功能，避免了繁琐而耗时的自相关算法，同时可以降低误差，提高计算精度。利用求解得到的参数对频谱分量进行修正，最终将所有频谱分量的系数归一化为相同值。

步骤60：对频谱分量进行去卷积和频谱平移处理；

在步骤60中，对于光学成像系统，去卷积算法是提高分辨率的重要手段。将每个未移位的频谱分量分别乘以其对应的维纳去卷积函数，得到滤波之后的频谱组分为：

在公式(12)中，w是根据经验调节的常数，一般与数据信噪比成反比。可见，此时样品的零频位于根据相移定理，将样品的零频移至正确位置，即得到：

需要注意的是，频谱平移之后很可能会超过图像所能容纳的最大频率(2×px)^-1，其中px为像素大小。所以需要在频域对图像进行扩展。首先将n×n的频域图像周围用0填充至2n×2n，这样图像能容纳的最大频率增加至px^-1，足够容纳平移之后的频谱。

步骤70：将所有平移至正确位置的频谱分量进行加权组合，得到图像的超分辨重构结果：

在公式(14)中，a(k)是切趾滤波函数，用于滤除高频噪声。对进行傅里叶逆变换得到时域的超分辨重构结果。

请参阅图4，是本申请实施例的图像重构系统的结构示意图。本申请实施例的图像重构系统包括图像预处理模块、频谱分离模块、照明光参数计算模块、otf参数计算模块、参数修正模块、频谱平移模块和频谱组合模块；具体地：

图像预处理模块：用于对sim成像系统采集到的原始图像进行预处理；在图像采集过程中，由于相机曝光时间和照明光强的浮动，需要对测量序列中的所有原始图像进行光强修正，保证每幅图像的平均值相同。另一方面，重构算法主要在频域进行，需要用到傅里叶变换，为了防止图像边缘突变引入的频域误差，需要对原始图像的边缘进行模糊，本申请实施例中，取10个像素对原始图像的边缘进行模糊化处理。

频谱分离模块：用于对预处理后的图像进行频谱分离，得到图像的频谱分量；其中，频谱分离具体为：根据平移相位，得到分离矩阵，根据公式(5)分别对每个方向上的3幅原始图像进行分离，最终得到三个方向上的9个频谱分量

照明光参数计算模块：用于根据频谱分量计算照明光的方向和空间频率；由频谱分量的表达式可知，样品的零频位于处，根据傅里叶变换性质，样品的零频值是时域图像上所有像素的光强之和，即对于sim系统，照明光的空间频率漂移较小，因此根据其空间频率初始值和角度初始值设定滤波函数，对频谱分量进行滤波，通过搜索确定滤波之后的数据的最大值，然后在最大值附近±1个像素的范围内，采用质心法将峰值点位置精确到亚像素，进而得到照明光的方向和空间频率。由于频谱的共轭性，只需要对m≥0的组分进行峰值点查找。可以尽可能地排除其他干扰，提取出有用的信息量，在此基础上进行峰值点查找，以提高空间频率的求解精度。

otf参数计算模块：用于获取sim成像系统的otf参数；otf在sim重构算法中是一个重要的参数，获取otf参数的方法包括两种：一种是采用荧光小球对sim成像系统进行标定，通过数据处理，得到psf，进而得到otf，这种方法可以反映真实的实验系统，但是由于荧光信号较弱，测量数据受噪声干扰较大，且荧光小球并非理想点光源，所以很难得到效果较好的otf。另一种方法是采用psf理论模型，目前常用的psf理论模型有born&wolf模型、richard&wolf模型、gibson&lanni模型等。考虑到商用高精度物镜的像差可以控制在很小的范围内，本申请实施例中使用的psf理论模型是richard&wolf模型。模型中使用到的参数是根据实测的psf进行确定，以保证模型与实际系统更为接近。

参数修正模块：用于通过线性拟合算法计算照明光的调制度和初始相位，通过照明光的调制度和初始相位对频谱分量进行修正，将所有频谱分量的系数归一化为相同值；其中，由公式(4)可知，每个频谱分量的系数是不同的，因此需要对其进行归一化，这样才能通过各个频谱的线性组合，获得正确的样品信息。本申请实施例对频谱分量的系数进行归一化的方式具体为：

(a)首先进行中心组分的匹配，由上述推倒可知，中心组分(m＝0)的表达式为：

假设b＝s*a，则s＝c1,0/cd,0，参数s可以通过对表达式a和b进行线性拟合得到，之后将所有其他组分用参数s进行修正，即：

其中表示将d1用d2替换。

(b)边缘组分(m＝±1)与中心组分进行匹配，根据傅里叶频移定理：

将边缘频谱变换到实域，乘以相移因子然后再变换至频域，得到移位之后的频谱：

同样地，可以获得平移之后的otf，即这样得到表达式a和b：

这两个表达式在otf的重合区域应有相同的数值所以假设b＝s*a，通过复数线性拟合得到参数如此，边缘组分可以修正为：

对于m＝1和m＝-1，参数s是共轭的，因此只需要计算m＝1时的参数。此时所有频谱分量的参数都归一化为c1,0。

参数修正中使用到的线性拟合算法的步骤为：首先计算表达式a和b的有效支持域，一般取|otf|>0.1的区域，如公式(10)所述，有效支持域为和的重叠区域；然后由有效支持域内a和b的数据组成数据对(ai,bi)，由于噪声和计算误差的存在，需要对数据对(ai,bi)进行数据筛选，得到有效数据对，以提高计算精度，数据筛选过程主要包括：

(i)傅里叶相移定理会引入不可避免的计算误差，通过仿真证明该计算误差对绝对值较小的数据的相位影响是显著的，因此数据筛选的第一步是去掉绝对值过小的数据；

(ii)虽然理论上数据bi与数据ai成正比，但是由于噪声的存在，|bi|与|ai|的比值分布范围较大，但大部分数值集中在一定区域内，所以选择比值|bi|/|ai|分布在其最大概率附近的数据组作为有效数据。经过筛选的有效数据对(ai,bi)通过复数线性拟合用于参数s的估计，即

本申请实施例中通过将初始相位和调制度作为一个整体参数进行求解，在线性拟合算法的基础上增加数据筛选功能，避免了繁琐而耗时的自相关算法，同时可以降低误差，提高计算精度。利用求解得到的参数对频谱分量进行修正，最终将所有频谱分量的系数归一化为相同值。

频谱平移模块：用于对频谱分量进行去卷积和频谱平移处理；对于光学成像系统，去卷积算法是提高分辨率的重要手段。将每个未移位的频谱分量分别乘以其对应的维纳去卷积函数，得到滤波之后的频谱组分为：

在公式(12)中，w是根据经验调节的常数，一般与数据信噪比成反比。可见，此时样品的零频位于根据相移定理，将样品的零频移至正确位置，即得到：

需要注意的是，频谱平移之后很可能会超过图像所能容纳的最大频率(2×px)^-1，其中px为像素大小。所以需要在频域对图像进行扩展。首先将n×n的频域图像周围用0填充至2n×2n，这样图像能容纳的最大频率增加至px^-1，足够容纳平移之后的频谱。

频谱组合模块：用于将所有平移至正确位置的频谱分量进行加权组合，得到图像的超分辨重构结果：

在公式(14)中，a(k)是切趾滤波函数，用于滤除高频噪声。对进行傅里叶逆变换得到时域的超分辨重构结果。

实例说明：

采用本申请实施例的图像重构算法对sim原始图像(来自geomx)进行重构，得到的图像重构结果如图5所示，为本申请实施例的图像重构结果示意图。图示的左下角图像是黑色方框内的放大图，标尺为4um。其中，图5(a)表示无参数修正部分的图像重构结果；图5(b)表示有参数修正部分的图像重构结果；图5(c)表示宽场成像结果。对比图5(a)和图5(b)可知，参数修正部分是算法中必不可少的一项，对重构图像的质量有重要影响。对比图5(b)和图5(c)可知，本申请实施例的图像重构算法可有效提高系统分辨率。

本申请实施例的图像重构方法及系统利用理想相位对预处理图像进行频谱分离，根据照明光空间频率初始值设定滤波函数，对分离频谱进行滤波，再进行峰值点搜索，并通过质心法将峰值点的定位精度提高至亚像素级，以确定照明光的空间频率。采用简单而有效的线性拟合算法求解照明光的调制度和初始相位，并在此基础上增加数据筛选功能，避免了繁琐而耗时的互相关算法，同时又提高了求解精度；并利用计算得到的调制度和初始相位对频谱分量进行参数归一化，从而有效提高重构图像的质量。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本申请。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本申请的精神或范围的情况下，在其它实施例中实现。因此，本申请将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。