CN107180411B - 一种图像重构方法及系统 - Google Patents

一种图像重构方法及系统 Download PDF

Info

Publication number
CN107180411B
CN107180411B CN201710362081.3A CN201710362081A CN107180411B CN 107180411 B CN107180411 B CN 107180411B CN 201710362081 A CN201710362081 A CN 201710362081A CN 107180411 B CN107180411 B CN 107180411B
Authority
CN
China
Prior art keywords
image
component
illumination light
spectrum
frequency spectrum
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710362081.3A
Other languages
English (en)
Other versions
CN107180411A (zh
Inventor
唐玉国
张艳微
巩岩
骆聪
郎松
王宏伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Suzhou Institute of Biomedical Engineering and Technology of CAS
Original Assignee
Suzhou Institute of Biomedical Engineering and Technology of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Suzhou Institute of Biomedical Engineering and Technology of CAS filed Critical Suzhou Institute of Biomedical Engineering and Technology of CAS
Priority to CN201710362081.3A priority Critical patent/CN107180411B/zh
Publication of CN107180411A publication Critical patent/CN107180411A/zh
Application granted granted Critical
Publication of CN107180411B publication Critical patent/CN107180411B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4053—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Investigating, Analyzing Materials By Fluorescence Or Luminescence (AREA)
  • Image Analysis (AREA)

Abstract

本申请涉及图像重构技术领域，特别涉及一种图像重构方法及系统。所述图像重构方法包括：步骤a：对图像进行频谱分离，得到所述图像的频谱分量；步骤b：对所述频谱分量进行峰值点搜索，获取照明光的方向和空间频率；步骤c：通过线性拟合算法计算所述照明光的调制度和初始相位，通过所述照明光的调制度和初始相位对所述频谱分量进行修正，将所述频谱分量的系数进行归一化；步骤d：对所述归一化的频谱分量进行去卷积和平移处理，并将所述去卷积和平移处理后的频谱分量进行频谱组合，得到所述图像的重构结果。本申请根据照明光空间频率初始值设定滤波函数，对分离频谱进行滤波，再进行峰值点搜索，提高了求解精度，有效提高重构图像的质量。

Description

一种图像重构方法及系统
技术领域
本申请涉及图像重构技术领域，特别涉及一种图像重构方法及系统。
背景技术
作为光学显微镜的一种，荧光显微镜是现代生命科学研究中必不可少的研究工具，然而传统光学显微镜的分辨率受到衍射极限的限制。近年来，一系列突破衍射极限的超分辨显微成像技术相继出现，例如随机光学重构显微镜、光激活定位显微镜、荧光光激活定位显微镜、受激发射损耗显微镜、结构光照明显微镜等。结构光照明显微镜(StructuredIllumination Microscope,SIM)是目前主流的超分辨成像技术之一。该技术采用加载了一定频率信息的结构光对样品进行照明(以下统称为照明光)，将传统显微镜无法获取的高频信息编码到采集图像中，再通过图像处理算法对样品信息进行重构，其横向分辨率可以达到传统分辨率的两倍。
超分辨重构算法是SIM的核心技术之一，其算法主要步骤包括：图像预处理、频谱分离、计算照明光空间频率、参数修正、去卷积、频谱平移以及频谱融合。影响重构图像质量的因素有很多，例如：平移相位不准确而导致的频谱分离不彻底、照明光参数的求解误差、去卷积过程引入的噪声放大、点扩散函数误差、离焦光线等。实验系统中照明光参数主要有空间频率、方向、调制度和初始相位等参数，由于实验系统的不稳定，这些参数是无法预先知道的，只能通过一系列原始测量数据对这些参数进行求解，其求解精度对最终重构图像的质量有重要影响。
目前，照明光参数的求解算法主要有以下几种方法：文献 1[Phase-shiftestimation in sinusoidally illuminated images for lateral superresolution]通过分析原始数据在频域中的峰值点的位置及其相位来确定照明光的空间频率和相位，这种方法对于高频率的照明光和信噪较大的数据并不适用；文献2[Surpassing the lateralresolution limit by a factor of two using structured illumination microscopy]通过计算平移之后与之前的组分进行互相关，进而求解得到照明光参数，但此方法计算量大，耗时较长，不适用于SIM数据的快速处理；文献3[Phase optimisation for sturcturedillumination microscopy]在文献2 的基础上增加了迭代优化算法，加快了算法的计算效率，但算法较为复杂，计算精度与初始值有关，初始值不当会引起算法无法收敛至准确值。
发明内容
本申请提供了一种图像重构方法及系统，旨在至少在一定程度上解决现有技术中的上述技术问题之一。
为了解决上述问题，本申请提供了如下技术方案：
一种图像重构方法，包括：
步骤a：对图像进行频谱分离，得到所述图像的频谱分量；
步骤b：对所述频谱分量进行峰值点搜索，获取照明光的方向和空间频率；
步骤c：通过线性拟合算法计算所述照明光的调制度和初始相位，通过所述照明光的调制度和初始相位对所述频谱分量进行修正，将所述频谱分量的系数进行归一化；
步骤d：对所述归一化的频谱分量进行去卷积和平移处理，并将所述去卷积和平移处理后的频谱分量进行频谱组合，得到所述图像的重构结果。
本申请实施例采取的技术方案还包括：所述步骤a还包括：对所述图像进行预处理；所述预处理包括：对所述图像进行光强修正，并对所述图像的边缘进行模糊化处理。
本申请实施例采取的技术方案还包括：在所述步骤b中，所述获取照明光的方向和空间频率具体包括：根据照明光的空间频率初始值设定滤波函数，对所述频谱分量进行滤波，对所述滤波后的频谱分量进行峰值点搜索，确定滤波后的频谱分量的最大值，并在所述最大值附近±1个像素的范围内，采用质心法确定所述峰值点位置，并精确到亚像素，进而得到所述照明光的方向和空间频率。
本申请实施例采取的技术方案还包括：在所述步骤c中，所述通过照明光的调制度和初始相位对频谱分量进行修正具体包括：
步骤c1：中心组分的匹配：中心组分的表达式为：
Figure DEST_PATH_GDA0001354945190000041
Figure DEST_PATH_GDA0001354945190000042
上述公式中，
Figure DEST_PATH_GDA0001354945190000043
是光学传递函数， d＝1,2,3；n＝0,1,2；m＝-1,0,1，
Figure DEST_PATH_GDA0001354945190000044
是频谱分量，假设B＝s*A，则 s＝c1,0/cd,0，通过对表达式A和B进行线性拟合得到参数s，并将所有其他组分用参数s进行修正，即：
Figure DEST_PATH_GDA0001354945190000045
其中
Figure DEST_PATH_GDA0001354945190000046
表示将D1用D2替换；
步骤c2：边缘组分与中心组分进行匹配：根据傅里叶频移定理：
Figure DEST_PATH_GDA0001354945190000047
Figure DEST_PATH_GDA0001354945190000048
首先将边缘频谱变换到实域，乘以相移因子
Figure DEST_PATH_GDA0001354945190000049
再变换至频域，得到移位之后的频谱：
Figure DEST_PATH_GDA00013549451900000410
同理，获得平移之后的OTF，即
Figure DEST_PATH_GDA00013549451900000411
得到表达式A和B：
Figure DEST_PATH_GDA00013549451900000412
Figure DEST_PATH_GDA00013549451900000413
假设B＝s*A，通过复数线性拟合得到参数
Figure DEST_PATH_GDA00013549451900000414
边缘组分修正为：
Figure DEST_PATH_GDA00013549451900000415
本申请实施例采取的技术方案还包括：在所述步骤c中，所述通过线性拟合算法计算照明光的调制度和初始相位具体为：考虑到噪声和计算误差的存在，需要对表达式A和B内的数据进行筛选。首先计算表达式A和B的有效支持域，对于中心组分修正，所述支持域取|OTF|>0.1的区域，对于边缘组分修正，所述有效支持域为
Figure DEST_PATH_GDA0001354945190000051
和
Figure DEST_PATH_GDA0001354945190000052
的重叠区域；然后由有效支持域内A和B的数据组成数据对(ai,bi)，并对数据对(ai,bi)进行进一步筛选，由于噪声的存在，|bi|与|ai|的比值分布范围较大，但大部分数值集中在一定区域内，所以选择比值|bi|/|ai|分布在其最大概率附近的数据组作为有效数据。最后对所述有效数据对进行复数线性拟合，得到参数s的估计，即
Figure DEST_PATH_GDA0001354945190000053
本申请实施例采取的另一技术方案为：一种图像重构系统，包括：
频谱分离模块：用于对图像进行频谱分离，得到所述图像的频谱分量；
照明光参数计算模块：用于对所述频谱分量进行峰值点搜索，获取照明光的方向和空间频率；
参数修正模块：用于通过线性拟合算法计算所述照明光的调制度和初始相位，通过所述照明光的调制度和初始相位对所述频谱分量进行修正，将所述频谱分量的系数进行归一化；
频谱平移模块：用于对所述归一化的频谱分量进行去卷积和平移处理；
频谱组合模块：用于将所述去卷积和平移处理后的频谱分量进行频谱组合，得到所述图像的重构结果。
本申请实施例采取的技术方案还包括图像预处理模块，所述图像预处理模块用于对图像进行预处理；所述预处理包括：对所述图像进行光强修正，并对所述图像的边缘进行模糊化处理。
本申请实施例采取的技术方案还包括：所述照明光参数计算模块获取照明光的方向和空间频率具体为：根据照明光的空间频率初始值设定滤波函数，对所述频谱分量进行滤波，对所述滤波后的频谱分量进行峰值点搜索，确定滤波后的频谱分量的最大值，并在所述最大值附近±1个像素的范围内，采用质心法确定所述峰值点位置，得到所述照明光的方向和空间频率。
本申请实施例采取的技术方案还包括：所述参数修正模块通过照明光的调制度和初始相位对频谱分量进行修正具体包括：
中心组分的匹配：中心组分的表达式为：
Figure DEST_PATH_GDA0001354945190000061
Figure DEST_PATH_GDA0001354945190000062
式中，
Figure DEST_PATH_GDA0001354945190000063
是光学传递函数，d＝1,2,3；n＝0,1,2；m＝-1,0,1,
Figure DEST_PATH_GDA0001354945190000064
是频谱分量，假设B＝s*A，则s＝c1,0/cd,0，通过对表达式A和B进行线性拟合得到参数s，并将所有其他组分用参数s进行修正，即：
Figure DEST_PATH_GDA0001354945190000065
其中
Figure DEST_PATH_GDA0001354945190000066
表示将D1用D2替换；
边缘组分与中心组分进行匹配：
Figure DEST_PATH_GDA0001354945190000067
Figure DEST_PATH_GDA0001354945190000068
将边缘频谱变换到实域，乘以相移因子
Figure DEST_PATH_GDA0001354945190000071
再变换至频域，得到移位之后的频谱：
Figure DEST_PATH_GDA0001354945190000072
获得平移之后的OTF，即
Figure DEST_PATH_GDA0001354945190000073
得到表达式A和B：
Figure DEST_PATH_GDA0001354945190000074
Figure DEST_PATH_GDA0001354945190000075
假设B＝s*A，通过复数线性拟合得到参数
Figure DEST_PATH_GDA0001354945190000076
边缘组分修正为：
Figure DEST_PATH_GDA0001354945190000077
本申请实施例采取的技术方案还包括：所述参数修正模块通过线性拟合算法计算照明光的调制度和初始相位具体为：计算表达式A 和B的有效支持域，取|OTF|>0.1的区域，所述有效支持域为
Figure DEST_PATH_GDA0001354945190000078
和
Figure DEST_PATH_GDA0001354945190000079
的重叠区域；由有效支持域内A和B的数据组成数据对(ai,bi)，并对数据对(ai,bi)进行数据筛选，得到有效数据对，所述有效数据对通过复数线性拟合进行参数s的估计，即
Figure DEST_PATH_GDA00013549451900000710
相对于现有技术，本申请实施例产生的有益效果在于：本申请实施例的图像重构方法及系统利用理想相位对预处理图像进行频谱分离，根据照明光空间频率初始值设定滤波函数，对分离频谱进行滤波，再进行峰值点搜索，并通过质心法将峰值点的定位精度提高至亚像素级，以确定照明光的空间频率。该方法尽可能地排除了其他干扰信息，在此基础上进行峰值点查找，可以提高空间频率的求解精度。采用简单而有效的线性拟合算法求解照明光的调制度和初始相位，并在此基础上增加数据筛选功能，避免了繁琐而耗时的互相关算法，同时又提高了求解精度；并利用计算得到的调制度和初始相位对频谱分量进行参数归一化，从而有效提高重构图像的质量。
附图说明
图1是本申请实施例的图像重构方法的流程图；
图2为SIM超分辨原理图；
图3为本申请实施例的频谱分量滤波示意图；
图4是本申请实施例的图像重构系统的结构示意图；
图5为本申请实施例的图像重构结果示意图；其中，图5(a)表示无参数修正部分的图像重构结果；图5(b)表示有参数修正部分的图像重构结果；图5(c)表示宽场成像结果。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本申请，并不用于限定本申请。
请参阅图1，是本申请实施例的图像重构方法的流程图。本申请实施例的图像重构方法包括以下步骤：
步骤10：对SIM成像系统采集到的原始图像进行预处理；
在步骤10中，SIM成像系统的成像原理为：SIM是一个线性移不变光学系统，图像是样品发射信息与系统PSF(point spread function，点扩散函数)的卷积：
Figure DEST_PATH_GDA0001354945190000091
在公式(1)中，
Figure DEST_PATH_GDA0001354945190000092
是被测的样品信息，
Figure DEST_PATH_GDA0001354945190000093
是照明光，
Figure DEST_PATH_GDA0001354945190000094
是系统PSF。在线性SIM中，照明光通常为正弦条纹，可表示为：
Figure DEST_PATH_GDA0001354945190000095
在公式(2)中，
Figure DEST_PATH_GDA0001354945190000096
为照明光的空间频率矢量，表示了照明光的周期和方向，a表示照明光的调制度。结合公式(1)和(2)，图像的频域表达式为：
Figure DEST_PATH_GDA0001354945190000097
在公式(3)中，m＝-1,0,1，符号(～)表示时域变量的频域表达式，
Figure DEST_PATH_GDA0001354945190000098
是光学传递函数(optical transfer function，OTF)。可见，每幅采集得到的原始图像都是由0和±1级频谱分量组合而成。
为了将三个频谱分量分离，需要改变照明光的相位，照明光的相位值为初始相位和平移相位之和，即
Figure DEST_PATH_GDA0001354945190000099
初始相位是未知的，与照明光的绝对位置有关，平移相位一般取0～2π上等间隔分布的相位，即
Figure DEST_PATH_GDA00013549451900000910
为了均匀的提高分辨率，一般采用相隔60度的三个照明光方向，如此，SIM至少需要采集9幅原始图像。由于系统误差，每个方向上的照明光的调制度、初始相位和空间频率是不同的，公式(3)又可以表示为：
Figure DEST_PATH_GDA0001354945190000101
在公式(4)中，d表示照明光的方向，n表示平移相位，m表示频谱分量，d＝1,2,3；n＝0,1,2；m＝-1,0,1,
Figure DEST_PATH_GDA0001354945190000102
是待分离的频谱分量。
SIM超分辨原理如图2所示。图2(a)中，13表示传统宽场显微镜可观测到的样品频域信息，由于系统OTF的支持域有限，支持域外的高频信息无法获取，因此分辨率受限。采用正弦结构光对样品进行照明，如公式(4)所示，SIM获取的样品信息中含有高频组分，图2(b)中，14和15则表示某方向上被平移至正确位置的高频分量，可见，在此方向上，分辨率有所提高。为了均匀提高系统分辨率，照明光采用三个方向，将所有方向上的频谱分量叠加，就得到图2 (c)中16所示的区域。相比于13，区域16所示的系统分辨率最大可提高2倍。
在图像采集过程中，由于相机曝光时间和照明光强的浮动，需要对测量序列中的所有原始图像进行光强修正，保证每幅图像的平均值相同。另一方面，重构算法主要在傅里叶频域进行，为了防止图像边缘突变引入的频域误差，需要对原始图像的边缘进行模糊，本申请实施例中，取原始图像边缘10个像素进行正弦模糊化处理。
步骤20：对预处理后的图像进行频谱分离，得到图像的频谱分量；
在步骤20中，频谱分离具体为：根据平移相位，得到分离矩阵，根据公式(5)分别对每个方向上的3幅原始图像进行分离，最终得到三个方向上的9个频谱分量
Figure DEST_PATH_GDA0001354945190000111
Figure DEST_PATH_GDA0001354945190000112
步骤30：根据频谱分量计算照明光的方向和空间频率；
在步骤30中，由频谱分量
Figure DEST_PATH_GDA0001354945190000113
的表达式可知，样品的零频位于
Figure DEST_PATH_GDA0001354945190000114
处，根据傅里叶变换性质，样品的零频值是时域图像上所有像素的光强之和，即
Figure DEST_PATH_GDA0001354945190000115
对于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)的表达式为：
Figure DEST_PATH_GDA0001354945190000121
Figure DEST_PATH_GDA0001354945190000122
假设B＝s*A，则s＝c1,0/cd,0，参数s可以通过对表达式A和B进行线性拟合得到，之后将所有其他组分用参数s进行修正，即：
Figure DEST_PATH_GDA0001354945190000123
其中
Figure DEST_PATH_GDA0001354945190000124
表示将D1用D2替换。
(b)边缘组分(m＝±1)与中心组分进行匹配，根据傅里叶频移定理：
Figure DEST_PATH_GDA0001354945190000131
Figure DEST_PATH_GDA0001354945190000132
将边缘频谱变换到实域，乘以相移因子
Figure DEST_PATH_GDA0001354945190000133
然后再变换至频域，得到移位之后的频谱：
Figure DEST_PATH_GDA0001354945190000134
同样地，可以获得平移之后的OTF，即
Figure DEST_PATH_GDA0001354945190000135
这样得到表达式A 和B：
Figure DEST_PATH_GDA0001354945190000136
Figure DEST_PATH_GDA0001354945190000137
这两个表达式在OTF的重合区域应有相同的数值
Figure DEST_PATH_GDA0001354945190000138
所以假设B＝s*A，通过复数线性拟合得到参数
Figure DEST_PATH_GDA0001354945190000139
如此，边缘组分可以修正为：
Figure DEST_PATH_GDA00013549451900001310
对于m＝1和m＝-1，参数s是共轭的，因此只需要计算m＝1时的参数。此时所有频谱分量的参数都归一化为c1,0。
考虑到噪声和计算误差，需要对表达式A和B中的数据进行筛选，然后再进行复数线性拟合计算。首先，傅里叶相移定理会引入不可避免的计算误差，通过仿真分析证明该计算误差对绝对值较小的数据的相位影响是显著的，因此数据筛选的第一步是去掉绝对值过小的数据，即根据OTF确定有效支持域。本申请实施例中，对于中心组分修正，所述支持域取|OTF|>0.1的区域，对于边缘组分修正，所述有效支持域为
Figure DEST_PATH_GDA0001354945190000141
和
Figure DEST_PATH_GDA0001354945190000142
的重叠区域。然后由有效支持域内A和B的数据组成数据对(ai,bi)，虽然理论上数据bi与数据ai成正比，但是由于噪声的存在，|bi|与|ai|的比值分布范围较大，但大部分数值集中在一定区域内，所以选择比值|bi|/|ai|分布在其最大概率附近的数据组作为有效数据。最后经过筛选的有效数据对(ai,bi)通过复数线性拟合用于参数s的估计，即
Figure DEST_PATH_GDA0001354945190000143
本申请实施例中通过将初始相位和调制度作为一个整体参数进行求解，在线性拟合算法的基础上增加数据筛选功能，避免了繁琐而耗时的自相关算法，同时可以降低误差，提高计算精度。利用求解得到的参数对频谱分量进行修正，最终将所有频谱分量的系数归一化为相同值。
步骤60：对频谱分量进行去卷积和频谱平移处理；
在步骤60中，对于光学成像系统，去卷积算法是提高分辨率的重要手段。将每个未移位的频谱分量
Figure DEST_PATH_GDA0001354945190000144
分别乘以其对应的维纳去卷积函数，得到滤波之后的频谱组分为：
Figure DEST_PATH_GDA0001354945190000145
在公式(12)中，w是根据经验调节的常数，一般与数据信噪比成反比。可见，此时样品的零频位于
Figure DEST_PATH_GDA0001354945190000146
根据相移定理，将样品的零频移至正确位置，即得到：
Figure DEST_PATH_GDA0001354945190000147
需要注意的是，频谱平移之后很可能会超过图像所能容纳的最大频率(2×px)-1，其中px为像素大小。所以需要在频域对图像进行扩展。首先将N×N的频域图像周围用0填充至2N×2N，这样图像能容纳的最大频率增加至px-1，足够容纳平移之后的频谱。
步骤70：将所有平移至正确位置的频谱分量进行加权组合，得到图像的超分辨重构结果：
Figure DEST_PATH_GDA0001354945190000151
在公式(14)中，A(k)是切趾滤波函数，用于滤除高频噪声。对
Figure DEST_PATH_GDA0001354945190000152
进行傅里叶逆变换得到时域的超分辨重构结果。
请参阅图4，是本申请实施例的图像重构系统的结构示意图。本申请实施例的图像重构系统包括图像预处理模块、频谱分离模块、照明光参数计算模块、OTF参数计算模块、参数修正模块、频谱平移模块和频谱组合模块；具体地：
图像预处理模块：用于对SIM成像系统采集到的原始图像进行预处理；在图像采集过程中，由于相机曝光时间和照明光强的浮动，需要对测量序列中的所有原始图像进行光强修正，保证每幅图像的平均值相同。另一方面，重构算法主要在频域进行，需要用到傅里叶变换，为了防止图像边缘突变引入的频域误差，需要对原始图像的边缘进行模糊，本申请实施例中，取10个像素对原始图像的边缘进行模糊化处理。
频谱分离模块：用于对预处理后的图像进行频谱分离，得到图像的频谱分量；其中，频谱分离具体为：根据平移相位，得到分离矩阵，根据公式(5)分别对每个方向上的3幅原始图像进行分离，最终得到三个方向上的9个频谱分量
Figure DEST_PATH_GDA0001354945190000161
Figure DEST_PATH_GDA0001354945190000162
照明光参数计算模块：用于根据频谱分量计算照明光的方向和空间频率；由频谱分量
Figure DEST_PATH_GDA0001354945190000163
的表达式可知，样品的零频位于
Figure DEST_PATH_GDA0001354945190000164
处，根据傅里叶变换性质，样品的零频值是时域图像上所有像素的光强之和，即
Figure DEST_PATH_GDA0001354945190000165
对于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)的表达式为：
Figure DEST_PATH_GDA0001354945190000171
Figure DEST_PATH_GDA0001354945190000172
假设B＝s*A，则s＝c1,0/cd,0，参数s可以通过对表达式A和B进行线性拟合得到，之后将所有其他组分用参数s进行修正，即：
Figure DEST_PATH_GDA0001354945190000173
其中
Figure DEST_PATH_GDA0001354945190000174
表示将D1用D2替换。
(b)边缘组分(m＝±1)与中心组分进行匹配，根据傅里叶频移定理：
Figure DEST_PATH_GDA0001354945190000181
Figure DEST_PATH_GDA0001354945190000182
将边缘频谱变换到实域，乘以相移因子
Figure DEST_PATH_GDA0001354945190000183
然后再变换至频域，得到移位之后的频谱：
Figure DEST_PATH_GDA0001354945190000184
同样地，可以获得平移之后的OTF，即
Figure DEST_PATH_GDA0001354945190000185
这样得到表达式A 和B：
Figure DEST_PATH_GDA0001354945190000186
Figure DEST_PATH_GDA0001354945190000187
这两个表达式在OTF的重合区域应有相同的数值
Figure DEST_PATH_GDA0001354945190000188
所以假设B＝s*A，通过复数线性拟合得到参数
Figure DEST_PATH_GDA0001354945190000189
如此，边缘组分可以修正为：
Figure DEST_PATH_GDA00013549451900001810
对于m＝1和m＝-1，参数s是共轭的，因此只需要计算m＝1时的参数。此时所有频谱分量的参数都归一化为c1,0。
参数修正中使用到的线性拟合算法的步骤为：首先计算表达式A 和B的有效支持域，一般取|OTF|>0.1的区域，如公式(10)所述，有效支持域为
Figure DEST_PATH_GDA00013549451900001811
和
Figure DEST_PATH_GDA00013549451900001812
的重叠区域；然后由有效支持域内A 和B的数据组成数据对(ai,bi)，由于噪声和计算误差的存在，需要对数据对(ai,bi)进行数据筛选，得到有效数据对，以提高计算精度，数据筛选过程主要包括：
(I)傅里叶相移定理会引入不可避免的计算误差，通过仿真证明该计算误差对绝对值较小的数据的相位影响是显著的，因此数据筛选的第一步是去掉绝对值过小的数据；
(II)虽然理论上数据bi与数据ai成正比，但是由于噪声的存在， |bi|与|ai|的比值分布范围较大，但大部分数值集中在一定区域内，所以选择比值|bi|/|ai|分布在其最大概率附近的数据组作为有效数据。经过筛选的有效数据对(ai,bi)通过复数线性拟合用于参数s的估计，即
Figure DEST_PATH_GDA0001354945190000191
本申请实施例中通过将初始相位和调制度作为一个整体参数进行求解，在线性拟合算法的基础上增加数据筛选功能，避免了繁琐而耗时的自相关算法，同时可以降低误差，提高计算精度。利用求解得到的参数对频谱分量进行修正，最终将所有频谱分量的系数归一化为相同值。
频谱平移模块：用于对频谱分量进行去卷积和频谱平移处理；对于光学成像系统，去卷积算法是提高分辨率的重要手段。将每个未移位的频谱分量
Figure DEST_PATH_GDA0001354945190000192
分别乘以其对应的维纳去卷积函数，得到滤波之后的频谱组分为：
Figure DEST_PATH_GDA0001354945190000193
在公式(12)中，w是根据经验调节的常数，一般与数据信噪比成反比。可见，此时样品的零频位于
Figure DEST_PATH_GDA0001354945190000194
根据相移定理，将样品的零频移至正确位置，即得到：
Figure DEST_PATH_GDA0001354945190000201
需要注意的是，频谱平移之后很可能会超过图像所能容纳的最大频率(2×px)-1，其中px为像素大小。所以需要在频域对图像进行扩展。首先将N×N的频域图像周围用0填充至2N×2N，这样图像能容纳的最大频率增加至px-1，足够容纳平移之后的频谱。
频谱组合模块：用于将所有平移至正确位置的频谱分量进行加权组合，得到图像的超分辨重构结果：
Figure DEST_PATH_GDA0001354945190000202
在公式(14)中，A(k)是切趾滤波函数，用于滤除高频噪声。对
Figure DEST_PATH_GDA0001354945190000203
进行傅里叶逆变换得到时域的超分辨重构结果。
实例说明：
采用本申请实施例的图像重构算法对SIM原始图像(来自GE OMX)进行重构，得到的图像重构结果如图5所示，为本申请实施例的图像重构结果示意图。图示的左下角图像是黑色方框内的放大图，标尺为4um。其中，图5(a)表示无参数修正部分的图像重构结果；图5(b)表示有参数修正部分的图像重构结果；图5(c)表示宽场成像结果。对比图5(a)和图5(b)可知，参数修正部分是算法中必不可少的一项，对重构图像的质量有重要影响。对比图5(b)和图5(c)可知，本申请实施例的图像重构算法可有效提高系统分辨率。
本申请实施例的图像重构方法及系统利用理想相位对预处理图像进行频谱分离，根据照明光空间频率初始值设定滤波函数，对分离频谱进行滤波，再进行峰值点搜索，并通过质心法将峰值点的定位精度提高至亚像素级，以确定照明光的空间频率。采用简单而有效的线性拟合算法求解照明光的调制度和初始相位，并在此基础上增加数据筛选功能，避免了繁琐而耗时的互相关算法，同时又提高了求解精度；并利用计算得到的调制度和初始相位对频谱分量进行参数归一化，从而有效提高重构图像的质量。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本申请。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本申请的精神或范围的情况下，在其它实施例中实现。因此，本申请将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (8)

1.一种图像重构方法，其特征在于，包括：
步骤a：对图像进行频谱分离，得到所述图像的频谱分量；
步骤b：对所述频谱分量进行峰值点搜索，获取照明光的方向和空间频率；
步骤c：通过线性拟合算法计算所述照明光的调制度和初始相位，通过所述照明光的调制度和初始相位对所述频谱分量进行修正，将所述频谱分量的系数进行归一化；
步骤d：对所述归一化的频谱分量进行去卷积和平移处理，并将所述去卷积和平移处理后的频谱分量进行频谱组合，得到所述图像的重构结果；
其中，在所述步骤c中，通过所述照明光的调制度和初始相位对所述频谱分量进行修正具体包括：
步骤c1：中心组分的匹配：中心组分的表达式为：
Figure FDA0002664553480000011
Figure FDA0002664553480000012
上述公式中，
Figure FDA0002664553480000013
是光学传递函数，d∈{1,2,3}；m∈{-1,0,1}，
Figure FDA0002664553480000014
是频谱分量，假设B＝s*A，则s＝c1,0/cd,0，通过对表达式A和B进行线性拟合得到参数s，并将所有其他组分用参数s进行修正，即：
Figure FDA0002664553480000015
其中
Figure FDA0002664553480000016
表示将D1用D2替换；
步骤c2：边缘组分与中心组分进行匹配：根据傅里叶频移定理：
Figure FDA0002664553480000017
Figure FDA0002664553480000018
首先将边缘频谱变换到实域，乘以相移因子
Figure FDA0002664553480000021
再变换至频域，得到移位之后的频谱：
Figure FDA0002664553480000022
同理，获得平移之后的OTF，即
Figure FDA0002664553480000023
得到表达式A和B：
Figure FDA0002664553480000024
Figure FDA0002664553480000025
假设B＝s*A，通过复数线性拟合得到参数
Figure FDA0002664553480000026
边缘组分修正为：
Figure FDA0002664553480000027
2.根据权利要求1所述的图像重构方法，其特征在于，所述步骤a还包括：对所述图像进行预处理；所述预处理包括：对所述图像进行光强修正，并对所述图像的边缘进行模糊化处理。
3.根据权利要求1或2所述的图像重构方法，其特征在于，在所述步骤b中，所述获取照明光的方向和空间频率具体包括：根据照明光的空间频率初始值设定滤波函数，对所述频谱分量进行滤波，对所述滤波后的频谱分量进行峰值点搜索，确定滤波后的频谱分量的最大值，并在所述最大值附近±1个像素的范围内，采用质心法确定所述峰值点位置，并精确到亚像素，得到所述照明光的方向和空间频率。
4.根据权利要求3所述的图像重构方法，其特征在于，在所述步骤c中，通过线性拟合算法计算所述照明光的调制度和初始相位具体为：计算表达式A和B的有效支持域，对于中心组分修正，所述支持域取|OTF|>0.1的区域，对于边缘组分修正，所述有效支持域为
Figure FDA0002664553480000028
和
Figure FDA0002664553480000029
的重叠区域；分别取计算表达式A和B在有效支持域内的数据点ai和bi组成数据对(ai,bi)，并对数据对(ai,bi)进行数据筛选，得到有效数据对，对所述有效数据对通过复数线性拟合进行参数s的估计，即
Figure FDA0002664553480000031
5.一种图像重构系统，其特征在于，包括：
频谱分离模块：用于对图像进行频谱分离，得到所述图像的频谱分量；
照明光参数计算模块：用于对所述频谱分量进行峰值点搜索，获取照明光的方向和空间频率；
参数修正模块：用于通过线性拟合算法计算所述照明光的调制度和初始相位，通过所述照明光的调制度和初始相位对所述频谱分量进行修正，将所述频谱分量的系数进行归一化；
频谱平移模块：用于对所述归一化的频谱分量进行去卷积和平移处理；
频谱组合模块：用于将所述去卷积和平移处理后的频谱分量进行频谱组合，得到所述图像的重构结果；
所述参数修正模块通过照明光的调制度和初始相位对频谱分量进行修正具体包括：
中心组分的匹配：中心组分的表达式为：
Figure FDA0002664553480000037
Figure FDA0002664553480000032
上述公式中，
Figure FDA0002664553480000033
是光学传递函数，d∈{1,2,3}；m∈{-1,0,1},
Figure FDA0002664553480000034
是频谱分量，假设B＝s*A，则s＝c1,0/cd,0，通过对表达式A和B进行线性拟合得到参数s，并将所有其他组分用参数s进行修正，即：
Figure FDA0002664553480000035
其中
Figure FDA0002664553480000036
表示将D1用D2替换；
边缘组分与中心组分进行匹配：
Figure FDA0002664553480000041
Figure FDA0002664553480000042
将边缘频谱变换到实域，乘以相移因子
Figure FDA0002664553480000043
再变换至频域，得到移位之后的频谱：
Figure FDA0002664553480000044
获得平移之后的OTF，即
Figure FDA0002664553480000045
得到表达式A和B：
Figure FDA0002664553480000046
Figure FDA0002664553480000047
假设B＝s*A，通过复数线性拟合得到参数
Figure FDA0002664553480000048
边缘组分修正为：
Figure FDA0002664553480000049
6.根据权利要求5所述的图像重构系统，其特征在于，还包括图像预处理模块，所述图像预处理模块用于对图像进行预处理；所述预处理包括：对所述图像进行光强修正，并对所述图像的边缘进行模糊化处理。
7.根据权利要求5或6所述的图像重构系统，其特征在于，所述照明光参数计算模块获取照明光的方向和空间频率具体为：根据照明光的空间频率初始值设定滤波函数，对所述频谱分量进行滤波，对所述滤波后的频谱分量进行峰值点搜索，确定滤波后的频谱分量的最大值，并在所述最大值附近±1个像素的范围内，采用质心法确定所述峰值点位置，并精确到亚像素，进而得到所述照明光的方向和空间频率。
8.根据权利要求7所述的图像重构系统，其特征在于，所述参数修正模块通过线性拟合算法计算照明光的调制度和初始相位具体为：计算表达式A和B的有效支持域，取|OTF|>0.1的区域，所述有效支持域为
Figure FDA0002664553480000051
和
Figure FDA0002664553480000052
的重叠区域；分别取计算表达式A和B在有效支持域内的数据点ai和bi组成数据对(ai,bi)，并对数据对(ai,bi)进行数据筛选，得到有效数据对，所述有效数据对通过复数线性拟合进行参数s的估计，即
Figure FDA0002664553480000053
CN201710362081.3A 2017-05-19 2017-05-19 一种图像重构方法及系统 Active CN107180411B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710362081.3A CN107180411B (zh) 2017-05-19 2017-05-19 一种图像重构方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710362081.3A CN107180411B (zh) 2017-05-19 2017-05-19 一种图像重构方法及系统

Publications (2)

Publication Number Publication Date
CN107180411A CN107180411A (zh) 2017-09-19
CN107180411B true CN107180411B (zh) 2021-05-18

Family

ID=59832365

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710362081.3A Active CN107180411B (zh) 2017-05-19 2017-05-19 一种图像重构方法及系统

Country Status (1)

Country Link
CN (1) CN107180411B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109118453B (zh) * 2018-08-28 2022-03-04 西北核技术研究所 一种背景抑制的图像处理方法
CN109886907B (zh) * 2019-01-29 2023-05-12 云南大学 一种基于权值矩阵算法的基谐波b超融合方法
CN111077121B (zh) * 2019-12-06 2020-11-17 中国科学院西安光学精密机械研究所 空域中直接重构结构光照明超分辨图像的快速方法及系统
CN113496463B (zh) * 2020-03-19 2024-05-31 深圳华大生命科学研究院 图像重构方法、系统、电子装置及存储介质
CN112508791B (zh) * 2020-12-18 2022-11-04 中国工程物理研究院激光聚变研究中心 一种初始相位提取方法、装置、电子设备及存储介质
CN114092325B (zh) * 2021-09-24 2022-08-16 熵智科技(深圳)有限公司 一种荧光图像超分辨重建方法、装置、计算机设备及介质
US20230170368A1 (en) * 2021-11-30 2023-06-01 Visera Technologies Company Limited Image sensor and method for detecting images
CN114219715B (zh) * 2021-12-17 2022-07-08 广州超视计生物科技有限公司 一种结构照明光超分辨显微镜图像的质量增强方法及系统
CN116402678B (zh) * 2022-12-19 2023-10-20 中国科学院苏州生物医学工程技术研究所 超分辨结构光照明显微镜的频谱优化直接重建方法
CN116628759B (zh) * 2023-07-26 2023-10-03 徐州医科大学 一种MNSS平台通信Cookie数据模糊方法及数据管理方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104515759A (zh) * 2014-12-16 2015-04-15 中国科学院苏州生物医学工程技术研究所 非线性结构光照明显微成像方法及系统
CN105068232A (zh) * 2015-08-31 2015-11-18 福建师范大学 一种双通道结构光照明超分辨成像方法与装置
CN105589188A (zh) * 2016-03-10 2016-05-18 清华大学 一种结构光照明显微镜的成像方法及装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104515759A (zh) * 2014-12-16 2015-04-15 中国科学院苏州生物医学工程技术研究所 非线性结构光照明显微成像方法及系统
CN105068232A (zh) * 2015-08-31 2015-11-18 福建师范大学 一种双通道结构光照明超分辨成像方法与装置
CN105589188A (zh) * 2016-03-10 2016-05-18 清华大学 一种结构光照明显微镜的成像方法及装置

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Lateral superresolution using a posteriori phase shift estimation for a moving object:experimental results;Sapna A. Shroff等;《Optical Society of America》;20100831;第27卷(第8期);第1-13页 *
Saturated patterned excitation microscopy—a concept for optical resolution improvement;Rainer Heintzmann等;《https://www.researchgate.net/publication/11229619》;20020831;第19卷(第8期);第1-12页 *
基于激光干涉的结构光照明超分辨荧光显微镜系统;文刚等;《光学学报》;20170310;第37卷(第3期);第1-11页 *
结构光照明显微中的超分辨图像重建研究;周兴等;《光学学报》;20170310;第37卷(第3期);第1-12页 *

Also Published As

Publication number Publication date
CN107180411A (zh) 2017-09-19

Similar Documents

Publication Publication Date Title
CN107180411B (zh) 一种图像重构方法及系统
WO2019024491A1 (zh) 基于环状可编程led照明的高效率定量相位显微成像方法
US9560338B2 (en) Methods and systems for three dimensional optical imaging, sensing, particle localization and manipulation
Cannell et al. Image enhancement by deconvolution
CN107966801A (zh) 一种基于环形照明的高速傅立叶叠层成像装置及重构方法
WO2020258434A1 (zh) 一种基于tie的相位成像方法、装置及可读存储介质
CN111062889B (zh) 一种用于傅里叶叠层显微成像技术的光强校正方法
CN109416464B (zh) 在选择角度地照明时的伪迹减少
WO2017181044A1 (en) Optical phase retrieval systems using color-multiplexed illumination
US20060186311A1 (en) Method for improving depth discrimination in optical reproduction systems
CN110058392A (zh) 一种基于光强传输方程的散斑定量相位成像系统及其方法
Konda et al. Multi-aperture Fourier ptychographic microscopy, theory and validation
AU2010238375A1 (en) Improvements in imaging
CN113759535B (zh) 一种基于多角度照明反卷积的高分辨率显微成像方法
Li et al. PURE-LET deconvolution of 3D fluorescence microscopy images
Konda et al. Parallelized aperture synthesis using multi-aperture Fourier ptychographic microscopy
WO2019152216A1 (en) Systems and methods for robust background correction and/or emitter localization for super-resolution localization microscopy
Maalouf Contribution to fluorescence microscopy, 3D thick samples deconvolution and depth-variant PSF
CN114926357B (zh) 计算显微成像系统的led阵列光源位姿自校正方法
CN116430571A (zh) 多模态超分辨定量相位显微成像方法
CN111815544B (zh) 一种数字全息频谱中心亚像素搜索方法
CN112233040B (zh) 一种自动离焦校正的傅里叶叠层显微成像方法
US20230073901A1 (en) Systems and methods for performing multiple-wavelength quantitative phase imaging (qpi)
Hasanzade et al. Wide field-of-view Fourier Ptychography microscopy based on Fresnel propagation scheme
WO2003034010A1 (en) Phase determination of a radiation wavefield

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant