在STED(受激发射耗尽，stimulated emission depletion)荧光显微镜中，如例如在EP 0801 759 B1中描述，通过在显微镜物镜的焦点区域中通过光脉冲激发荧光分子或荧光染料来获得被荧光分子或荧光染料染色的样本的至少一个图像。此外，对样本应用去激发光以使焦点区域的中心外部的区域或焦点区域去激发。因此，这种显微镜可以生成空间分辨率低于通常的光衍射极限的图像。去激发可以源自脉冲光源或连续波(CW)光源。这两种光源通常都是激光，称为激发和耗尽激光。
US 9,551,658 B2描述了一种用于借助时间选通(time-gated)检测以更佳的空间分辨率提取源自样本的光脉冲的方法和设备。去激发光导致荧光寿命缩短。因此，荧光信号在焦点区域的中心外部的时间衰减短于中心处的时间衰减。在这种情况下，寿命τ由给出，其中knr是非辐射去激发率，kr是辐射去激发率，kSTED是由耗尽激光引入的额外的非辐射去激发率。时间门(time gate)排除掉在激发脉冲之后的短时期内的信号。剩余的信号被收集。在这种时间选通方法中，源自焦点区域的中心的光的收集也是不完整的。来自于在选通关闭周期内到达的光子的信号丢失，这导致比分辨率较差的非选通图像更低的信噪比。
出版物：Lanzanò等人(2015)，Encoding and decoding spatio-temporal information for super-resolution microscopy(《超分辨率显微镜的时空信息编码和解码》)，Nature Communications,6:6701,10.1038/ncomms7701，描述了一种SPLIT方法，该方法可以通过应用相量方法保存早期光子的信号。在这里，每个像素的寿命衰减被转换为极坐标图，称为相量图。在此图中，每个像素的坐标，称为相量，是高分辨率信号、低分辨率信号和背景三个分量坐标的线性组合。高分辨率分量的权重用作所得到图像的强度，并包含来自所有检测到的光子的信息。尽管使用了所有光子的信号，但该方法具有必须获取大量光子的局限性。在低光子计数条件下，相量图中的分布非常宽，以至于无法实现信噪比的提高。
为了实现更窄的分布，必须对相量坐标进行滤波。本书章节：Digman等人(2012年)，The phasor approach to fluorescence lifetime imaging:exploiting phasor linear properties,in Fluorescence Lifetime Spectroscopy and Imaging:Principles and Applications in Biomedical Diagnostics(《荧光寿命光谱和成像中荧光寿命成像的相量方法：利用相量线性特性：生物医学诊断中的原理和应用》)(编辑L Marcu,PMW French和DS Elson)，CRC出版社，描述了一种滤波相量图的方法，目的是通过使用中值滤波器来保持空间信息。然而，在超分辨率的情况下，即使是这种滤波方法也会导致空间分辨率的不必要的下降。
在数字图像处理中，更先进的滤波方法是已知的。例如，出版物：Sendur,L.&Selesnick,I.(2003)，Bivariate Shrinkage with Local Variance Estimation(《局部方差估计的双变量收缩》)，Signal Processing Letters，IEEE.9,438–441,10.1109/LSP.2002.806054，描述了一种基于小波变换的高级去噪滤波器，该滤波器在应用于图像时可以更好地保持空间分辨率。由于相量分量图像中噪声的强烈变化，以与中值滤波器相同的方式在相量分量上应用此滤波器或类似滤波器的尝试失败了。
技术实现要素：
本发明的目的是提供一种在低光子计数条件下增加STED图像中的信噪比同时保持选通STED方法的空间分辨率的方法。此目标是通过包括至少一些以下步骤的方法来实现的：
获取样本的STED图像，该STED图像包括像素；
针对图像的像素计算光子到达时间的傅立叶系数或光子到达时间数据，得到实部系数和虚部系数；
实部系数表示第一图像，且虚部系数表示第二图像；
从STED图像得出强度图像；
将空间滤波器应用于第一图像、第二图像和强度图像，得到各自的经滤波的图像；
基于经滤波的第一图像和经滤波的强度图像计算图像G；
基于经滤波的第二图像和经滤波的强度图像计算图像S；以及
基于两个图像G和S计算结果图像。
优选地，获取STED图像包括以下步骤：
用脉冲激发光照射所述样本，其中脉冲激发光在照射所述样本时包括焦点；
用围绕焦点的中心区域的去激发光照射所述样本。
通常，使用具有STED功能的扫描显微镜获取样本或样本部分的二维(也包括三维)STED图像。例如，所获取的二维STED图像包括x和y方向的像素(例如512x 512个像素)，取决于扫描显微镜的获取参数。针对每个像素，记录相对于激光脉冲的光子到达时间。典型荧光染料的荧光寿命在1至4ns的范围内。许多STED染料具有多种寿命组成。在焦点区域的中心没有或弱去激发光，并且在那里观察到自然寿命。在焦点区域的中心外部，寿命会减少/较低，但例如在应用去激发光的区域中，由于去激发过程，通常低于1ns。焦点区域的中心外部的寿命减少与来自耗尽激光器的耗尽能量成反比。对于典型的STED图像，在明亮的像素中记录大约100个光子计数。光子到达时间的傅立叶系数的计算是基于像素执行的，优选地针对图像的所有像素。因为每个像素都具有实部傅里叶系数和虚部傅里叶系数，所以可以将实部(傅里叶)系数视为第一图像，将虚部(傅里叶)系数视为第二图像。然而，第一图像或第二图像是否实际形成例如二维图像并不重要。第一图像和第二图像可以被表示或视为具有例如只是线性或其他合适的数据存储布置的数据结构。类似的情形适用于至少以下之一：结果图像、经滤波的图像、强度图像、图像G和图像S。
在优选实施例中，其中获取STED图像进一步包括分别地针对每个像素登记相对于激发脉冲的所检测到的荧光光子到达时间。
替代地或附加地，通过至少两个时间门登记到达时间。
优选地，傅立叶系数的计算包括第n阶傅立叶系数的计算，其中n是预选值(自然数)。
优选地，强度图像是从STED图像获取得出的，或者是从强度值的计算得出的，所述强度值来自于针对像素登记的光子数量。
在优选实施例中，图像G是基于将经滤波的第一图像除以经滤波的强度图像来计算的。
在另一个优选实施例中，图像S是基于将经滤波的第二图像除以经滤波的强度图像来计算的。
在进一步优选实施例中，空间滤波器是各向异性扩散、总变差去噪、最大熵去噪、非局部均值与小的波去噪中的至少一种。
上述目的进一步通过适用于执行上述方法的设备来解决。
此外，上述目的还通过执行上述方法的计算机程序来解决。本申请意义上的第一和/或第二图像可以以相应的傅立叶系数的二维阵列的形式来表示。
在本发明的特别优选的实施例中，在将傅立叶系数除以经滤波的强度之前应用第一图像、第二图像和强度图像的滤波器/空间滤波。先前，是在对图像执行除法的方法步骤之后应用滤波器。(除法后的)相邻像素的相量坐标具有相似的值，而原始傅里叶系数取决于强度而变化，尤其是在低光子计数条件下。在本上下文中，术语相量坐标的使用方式与Lanzanò等人描述的类似。对于较少数量的光子计数，相邻像素的傅立叶系数显示出比高光子计数更高的变化。相似平均值的不同变化与高级滤波方法的先决条件相矛盾。在现有技术中选择了中值滤波器，当相邻像素中的光子计数几乎没有变化时，该中值滤波器执行是合理的。然而，STED图像显示了必须保留的高度变化。当在对图像执行除法之前对所述图像进行滤波时，根据本发明，给出了高级滤波器的先决条件。此外，当对经滤波的图像执行除法时，存在一种偏见，即经滤波的图像可能包含来自具有不同权重的不同区域的信息，并且因此除法步骤可能会导致相量坐标的更广泛分布。
然而，通过根据特别优选实施例的过程所获得的结果出人意料地优于经典/先前的方法，以至于使得目标实际上可以实现。已发现，高级滤波算法对于相邻像素中的统计数据的强烈变化更加敏感。
附图说明
图1示出STED去激发对寿命衰减的影响，以及来自图1a中的焦点的中心区域、图1b中的中心外部区域的数据的典型寿命衰减的相量图。图1c示出组合的寿命衰减和时间门的典型选择。这些分量在相量图中的位置如图1.d所示。
图2示出根据本发明优选实施例的当应用不同滤波方法时来自二十个图像像素的相量分布：图2a没有滤波器，图2b是中值滤波器5x5，图2c是中值滤波器7x7，以及图2d是小波滤波器。使用的数据是针对每像素的十个光子计数的平均值进行模拟的。
图3示出了针对图3b中的最佳滤波器、图3c中的没有滤波器、图3d中的中值滤波器7x7以及图3e中的根据本发明的优选实施例的滤波器的关于非期望平滑边缘的滤波器效果。
具体实施方式
除了通常在共焦激光扫描显微镜中实施的标准STED显微镜之外，用于检测相对于激发光的激发脉冲的光子到达时间的装置对于根据本发明的方法是强制性的。所述装置可以例如包含根据时间相关的单光子计数原理(TCSPC原理)的记录设备或基于时间门的记录设备。
所记录的数据可以表示为对于像素i在到达时间t处检测到的光子的数量N：Ni(t)。对于每个像素，计算傅立叶系数Freal和Fimag：
整数值n是傅立叶系数的预选阶。通常，选择第一阶(1)，即n＝1。f是激发光的脉冲重复率。此外，对于每个像素，计算强度I：
在下一步骤中，使用空间滤波器，特别是使用去噪滤波器，对三个图像Freal、Fimag和I进行滤波。优选地，多种先进的去噪算法可以保持空间分辨率，并且因此可以在此步骤中使用。例如，可以应用以下算法：各向异性扩散(Anisotropic Diffusion)、总变差(Total Variation)去噪、最大熵(Maximum Entropy)去噪、非局部均值与小波(Non-local-means and Wavelet)去噪。
在优选实施例中，选择了小波滤波器。基本原理类似于用于去噪的傅里叶滤波器。高斯分布噪声对所有傅立叶系数的幅值造成额外的偏移。图像被变换到傅立叶空间。减去所述偏移，并对图像进行反变换(back transforme)。此过程也称为收缩。小波空间与傅立叶空间的区别在于小波系数表示用于图像中的区域的频率信息。对于小波去噪，使用系数收缩的相同原理。最简单的方法是以相同的偏移量收缩所有系数。已经开发了可以应用于本发明的更先进的方法，例如：Sendur,L.&Selesnick,I.，(2003)，Bivariate Shrinkage with Local Variance Estimation(《局部方差估计的双变量收缩》)，Signal Processing Letters,IEEE.Vol.9.438–441,10.1109/LSP.2002.806054。
收缩方法假设噪声的高斯分布。对于强度图像I，已知噪声分布是泊松分布。一种使用利用泊松分布噪声为高斯分布噪声设计的滤波器的方法是例如通过Anscombe变换(安斯科姆变换)来变换数据：
然后应用该滤波器，并对数据进行反变换。Anscombe变换后的噪声的标准差为1.0。
该情况对于两个图像Freal和Fimag有所不同。在积分运算后，噪声不是完全泊松分布的，而是取决于荧光寿命。
然而，该情况并不像对相量分量进行滤波的尝试那样严重。在第一阶中可以假设分布几乎是泊松，并且对于这些图像的噪声方差可以通过像中位数绝对偏差(MAD)这样的统计健壮方法来估计：
σ＝meduan(|Xi-median(Xi)|)/0.67449 (4)
其中σ是标准差，并且中位数运算应用于所有像素i的强度X。在σ和最佳收缩偏移之间存在依赖性，特别是线性依赖性。因数取决于小波变换的实现细节，并且可以针对特别使用的收缩算法推导出。例如。在Sendur,L.&Selesnick,I.中，所述依赖性在公式5中给出，部分II包含所述推导。
在优选实施例中，通过Anscombe变换(3)对所有三个图像Freal、Fimag和I进行变换，然后对于所有三个图像，通过中位数绝对偏差方法(4)估计噪声。经Anscombe变换的图像被变换到小波空间。通过从噪声的经估计的标准差得出的值对系数进行收缩。然后对经收缩的图像从小波空间进行反变换，并进行Anscombe反变换。结果得到经滤波的图像Freal_filtered,Fimag_filtered和I_filtered。
在下一步骤中，计算对于所有图像像素的相量系数：
Gi＝Freal_filteredi/I_filteredi,Si＝Fimag_filteredi/I_filteredi. (5)
在最后步骤中，根据所述相量系数计算结果图像。该计算可以类似于Lanzano等人(2015)，Encoding and decoding spatio-temporal information for super-resolution microscopy(《编码和解码超分辨率显微镜的时空信息》)，Nature Communications,10.1038/ncomms7701中的SPLIT方法来实现。
下面解释了根据本发明的方法表现得优于对相量G和S分量进行滤波的现有技术的原因。图1a示出沿纵坐标1的强度I和沿横坐标2的到达时间t的衰减图，其中图表3表示源自激发光的焦点中心区域的光子。图1b示出对于源自焦点的中心的外部区域的光子的图表6，在该区域应用去激发光。源自焦点中心的外部区域的光子的平均到达时间低于源自焦点的光子的平均到达时间。图1c示出对于所有检测到的光子的图表7。在时间选通的STED中，选择(例如通过选通)到达时间间隔8，其中来自中心区域的光子占主导地位。在此间隔之前检测到的光子被忽略。因此，具有关于该中心区域的信息的信号4(参见图1a)丢失。只有信号的部分5(参见图1a)被用于生成图像。与经典STED方法相比，时间选通方法得到空间分辨率更好但噪声较大的图像。
具有仅来自于中心区域的衰减的像素位于例如相量图中的位置9处(参见图1d)。来自外部的光子位于相量图中的不同位置10处。SPLIT方法(现有技术)背后的想法是取决于各个像素在相量图中的位置来将权重应用于所述像素的强度。这得到了与时间选通STED方法相当的空间分辨率，但由于没有丢失计数的光子，因此噪声更少。
图2a可视化了对于图像像素中10个光子计数的平均值的模拟数据。该模拟假设具有4n(4纳秒)寿命11的来自中心区域的50％光子和具有200ps(200皮秒)寿命12的来自外部的50％光子的混合。绘制了20个像素的相量坐标13。对于高光子计数，相量的位置必须在两个分量11和12的位置之间的一半处，即在14处。对于由15表示的没有滤波的低光子计数的情况是使得由于在没有滤波器的情况下的广泛分布，与选通方法相比，可能没有任何改进。
图2b示出当根据现有技术应用由16表示的5x5中值滤波器时的相同模拟数据。对于图2c，使用了由17表示的7x7中值滤波器。最后，图2d示出当使用根据本发明实施例的小波滤波器时的结果(由18表示)。为了实现相当的信噪比，需要7x7中值滤波器。
滤波器必须额外保留空间分辨率以便在STED应用中有用。具有两种混合的模拟示出沿边缘的滤波器的行为。在根据图3a的模拟中，在左半部分20中，图像19由2/3的低分辨率分量和1/3的高分辨率分量的混合创建。在右半部分21中，模拟了1/3的低分辨率分量和2/3的高分辨率分量。中心区域22中的8个像素的7行已经被垂直累加，并且已经用沿着相量中的在4ns和200ps位置之间的线的强度梯度进行了加权。
图3b示出了针对非常高光子计数的最佳条件，对于8个平均像素的强度23。在图3c中针对没有滤波的情况，在图3d中针对中值滤波器7x7的情况，以及在图2e中针对根据本发明的小波滤波器的情况，示出了对于每像素10次计数的情形，对于8个平均像素的强度。
根据本发明的小波滤波器比根据现有技术的中值滤波器(图3d)更好地保持了空间分辨率(图3e)，并且具有相当的信噪比(图2c和图2d)。
已针对相邻像素中的相似噪声特性设计了高级滤波器。这种滤波器当根据现有技术直接应用于相量分量时会失败。由于相量公式中的除以强度的除法，噪声改变太大。
一些或所有方法步骤可以通过(或使用)硬件装置来执行，例如处理器、微处理器、可编程计算机或电子电路。在一些实施例中，一个或多个最重要的方法步骤可以由这样的装置执行。
取决于某些实施要求，本发明的实施例可以以硬件或软件来实施。可以使用非临时性存储介质来执行该实施，诸如数字存储介质，例如软盘、DVD、蓝光、CD、ROM、PROM和EPROM、EEPROM或FLASH存储器，其上存储有电子可读控制信号，其与可编程计算机系统协作(或能够协作)以执行相应的方法。因此，数字存储介质可以是计算机可读的。
根据本发明的一些实施例包括具有电子可读控制信号的数据载体，其能够与可编程计算机系统协作，从而执行本文描述的方法之一。
通常，本发明的实施例可以实现为具有程序代码的计算机程序产品，当计算机程序产品在计算机上运行时，该程序代码可操作用于执行方法之一。例如，程序代码可以存储在机器可读载体上。
其他实施例包括存储在机器可读载体上的用于执行本文描述的方法之一的计算机程序。
换言之，因此，本发明的实施例是具有程序代码的计算机程序，该程序代码用于当计算机程序在计算机上运行时执行本文描述的方法之一。
因此，本发明的进一步实施例是一种存储介质(或数据载体，或计算机可读介质)，其包括存储在其上的计算机程序，用于在由处理器执行时执行本文所述的方法之一。数据载体、数字存储介质或记录介质通常是有形的和/或非过渡的。本发明的进一步实施例是如本文所述的装置，包括处理器和存储介质。
因此，本发明的进一步实施例是表示用于执行本文所述方法之一的计算机程序的数据流或信号序列。数据流或信号序列可以例如被配置为经由数据通信连接、例如经由互联网传送。
进一步实施例包括处理装置，例如计算机或可编程逻辑器件，其被配置为或适用于执行本文描述的方法之一。
进一步实施例包括其上安装有用于执行本文所述方法之一的计算机程序的计算机。
根据本发明的进一步实施例包括一种装置或系统，其被配置为将用于执行本文描述的方法之一的计算机程序(例如，电子地或光学地)传送到接收器。例如，接收器可以是计算机、移动设备、存储器设备等。例如，该装置或系统可以包括用于将计算机程序传送到接收器的文件服务器。
在一些实施例中，可编程逻辑器件(例如，现场可编程门阵列)可用于执行本文描述的方法的一些或所有功能。在一些实施例中，现场可编程门阵列可以与微处理器协作以执行本文描述的方法之一。通常，这些方法优选地由任何硬件装置来执行。
