CN116109762A - 自适应低伪影的多层三维结构光超分辨显微成像重建方法 - Google Patents
自适应低伪影的多层三维结构光超分辨显微成像重建方法 Download PDFInfo
- Publication number
- CN116109762A CN116109762A CN202211605447.2A CN202211605447A CN116109762A CN 116109762 A CN116109762 A CN 116109762A CN 202211605447 A CN202211605447 A CN 202211605447A CN 116109762 A CN116109762 A CN 116109762A
- Authority
- CN
- China
- Prior art keywords
- frequency
- dimensional
- spectrum
- level
- modulation
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 35
- 238000003384 imaging method Methods 0.000 title claims abstract description 17
- 238000001228 spectrum Methods 0.000 claims abstract description 99
- 239000013598 vector Substances 0.000 claims abstract description 98
- 239000010410 layer Substances 0.000 claims description 27
- 238000011156 evaluation Methods 0.000 claims description 15
- 238000005286 illumination Methods 0.000 claims description 14
- 230000003595 spectral effect Effects 0.000 claims description 9
- 230000003287 optical effect Effects 0.000 claims description 8
- 238000012546 transfer Methods 0.000 claims description 8
- 239000011159 matrix material Substances 0.000 claims description 6
- 238000005191 phase separation Methods 0.000 claims description 6
- 238000005457 optimization Methods 0.000 claims description 5
- 238000000926 separation method Methods 0.000 claims description 5
- 239000002356 single layer Substances 0.000 claims description 5
- 230000005284 excitation Effects 0.000 claims description 4
- 238000012545 processing Methods 0.000 claims description 4
- 238000012935 Averaging Methods 0.000 claims description 3
- 230000007547 defect Effects 0.000 claims description 3
- 238000009499 grossing Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000013519 translation Methods 0.000 claims description 3
- 238000002474 experimental method Methods 0.000 claims description 2
- 230000000694 effects Effects 0.000 abstract description 6
- 230000008030 elimination Effects 0.000 abstract description 2
- 238000003379 elimination reaction Methods 0.000 abstract description 2
- 238000001914 filtration Methods 0.000 abstract description 2
- 230000001629 suppression Effects 0.000 abstract description 2
- 230000003044 adaptive effect Effects 0.000 description 5
- 102000002151 Microfilament Proteins Human genes 0.000 description 3
- 108010040897 Microfilament Proteins Proteins 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 210000003632 microfilament Anatomy 0.000 description 3
- 102000029749 Microtubule Human genes 0.000 description 2
- 108091022875 Microtubule Proteins 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 210000004688 microtubule Anatomy 0.000 description 2
- 206010034972 Photosensitivity reaction Diseases 0.000 description 1
- 210000004027 cell Anatomy 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000000799 fluorescence microscopy Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000005764 inhibitory process Effects 0.000 description 1
- 230000008611 intercellular interaction Effects 0.000 description 1
- 238000000386 microscopy Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 201000008968 osteosarcoma Diseases 0.000 description 1
- 208000007578 phototoxic dermatitis Diseases 0.000 description 1
- 231100000018 phototoxicity Toxicity 0.000 description 1
- 210000004895 subcellular structure Anatomy 0.000 description 1
- 238000010869 super-resolution microscopy Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformation in the plane of the image
- G06T3/40—Scaling the whole image or part thereof
- G06T3/4053—Super resolution, i.e. output image resolution higher than sensor resolution
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/10—Image enhancement or restoration by non-spatial domain filtering
-
- G06T5/90—
-
- 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]
Abstract
本发明公开了一种自适应低伪影的多层三维结构光超分辨显微成像重建方法。本发明首先通过多级频谱协同估计的方法对频率矢量进行估计,极大的提高了低信噪比、低调制度情况下参数估计的稳定性与准确性,避免了参数估计错误带来的伪影;同时通过构建三维陷波函数和两步空间频域滤波,极大的降低了由于频谱尖峰和高频噪声等因素带来的重建伪影,同时提高了图像的分辨率;本发明在保真度、伪影抑制和散焦信号消除等方面优于其他算法,尤其是在低信噪比情况下有突出的重建效果。
Description
技术领域
本发明涉及荧光显微成像技术,具体涉及一种自适应低伪影的多层三维结构光超分辨显微成像重建方法。
背景技术
由于光毒性低、成像速度快等优势,结构化照明显微镜(SIM)被广泛用于观察亚细胞结构和细胞间互作。SIM通过结构光编码样本的高频信息,通过频谱搬移进行图像融合,从而进行图像的超分辨重构。然而,现有的SIM软件通常只能重建2D或单层3DSIM,如Fair-SIM和HiFi-SIM,但由于光学传递函数中存在的漏锥问题,2D或单层3DSIM算法在针对厚样本的重建中极易产生离焦背景与伪影。
而现有的多层3DSIM算法中对低信噪比的参数估计往往依赖于初始参数的选择,稳定性较差,同时其重建结果中严重的伪影问题也极大地限制了3DSIM的进一步发展。因此,发展一种自适应的稳定参数估计方法,通过优化频谱以去除重建伪影的多层三维结构光照明显微成像重建方法成为一种迫切的需要。
发明内容
针对3DSIM算法伪影较多的问题,本发明提出了一种自适应低伪影的多层三维结构光超分辨显微成像重建方法。
本发明的自适应低伪影的多层三维结构光超分辨显微成像重建方法,包括以下步骤:
1)得到原始图像栈:
结构光对样品的照射具有M个调制方向,每一个调制方向具有N个相位,在每一个调制方向上的每一个相位均采集一张原始图像,从而在样本的一层得到M×N张原始图像,利用平移升降台驱动样本进行z方向的逐层扫描,得到一共T层图像序列,并在每一层的每个调制方向上取5个不同相位,构成包含M×5×T张原始图像的原始图像栈;
2)采用多级估计法对频率矢量进行估计:
a)计算每一层图像序列的调制对比度,并将所有的T层图像序列以调制对比度作为权重进行加权求和,得到一层调制对比度最高的图像序列,此单层图像序列包括M×5张原始图像;
b)将对比度最高的图像序列中的每长原始图像进行多次反卷积处理,从而突出结构光调制后样本的高频频率峰值点,接着设置初始相位为0,对每个调制方向的结构光调制后的5张原始图像利用二维相位分离矩阵进行分离频谱,得到±2,±1,0级共5个二维频谱分量C2D n(kx,y),n=0,±1,±2,kx,y为xoy平面的频谱的单位矢量;
c)估计结构光照明条纹的频率矢量pθ:
根据步骤b)得到的0,±1,±2级二维频谱分量,采用+2,+1级频谱对结构光照明条纹的频率矢量估计进行逐级估计;
i.通过+2级频谱进行频率矢量估计:
通过+2级频谱与0级频谱的互相关迭代操作,估计M个调制方向上的结构光照明条纹的频率矢量分别为通过反余弦函数定义求得频率矢量所对应的空间角度为θ=θ1,θ2,…,θM;
ii.建立评价指标:
由于理想的正确的频率矢量的模值应当相等,其空间角度应当相差2π/M,与该标准相差越大,说明频率矢量估计的准确度越低,因此建立第一和第二评价指标E1和E2,分别表征频率矢量pθ和空间角度的估计准确程度,频率矢量的模值越接近且角度之差越均匀,证明通过n级频谱进行频率矢量估计的结果越准确,其中第一评价指标第二评价指标E2=std{|θ2-θ1|,|θ3-θ2|,…,|θ1-θM|},std为方差符号;并且,设置第一和第二标准阈值E1_max和E2_max;
iii.估计判定:
如果满足判定条件E1≤E1_max且E2≤E2_max则证明通过+2级频谱估计的频率矢量正确,则得到估计的频率矢量并得到频率矢量所对应的空间角度,进入步骤3);
如果E1>E1_max或E2>E2_max则证明通过+2级频谱估计的频率矢量不正确,由于结构光照明条纹的+1级条纹的调制度的强度通常强于+2级条纹的调制度,因此若通过+2级频谱的频率估计不正确,则采用+1级频谱进行协同估计,通过+1级频谱与0级频谱的互相关迭代操作估计频率矢量,并建立相应的第一和第二评价指标,如果满足判定条件,则得到估计的频率矢量并得到频率矢量所对应的空间角度,进入步骤3);
如果通过+2和+1级频谱估计的频率矢量都无法满足判定条件E1≤E1_max且E2≤E2_max,则选择最小的级次估计的频率矢量,并得到频率矢量所对应的空间角度,进入步骤3);
对所有调制方向上的频率矢量的模求均值得到xoy平面的调制频率矢量的模值
3)通过级次之间的强度关系估计不同角度结构光的调制深度m;并且通过级次之间重叠部分的复数差异估计不同角度的结构光的初始相位
4)基于输入的物理参数仿真出三维光转换函数(OTF)或直接输入实验测得的三维OTF;
5)进行频域分离:
a)三维结构光强度的空间分布表示为
其中x,y,z分别为三维空间上xyz方向的单位矢量,r为xoy平面的极坐标的极径,I0表示光强大小,m表示调制深度,和分别表示xoy平面和yoz平面的调制频率矢量的模值,ψ是光束的入射角,λ为激发波长,n0为介质折射率,根据
求得光束的入射角ψ,再根据得到yoz平面的调制频率矢量pz,整理得到从而结构光表示为基频信号、一次谐波与二次谐波的组合,a0=1+2m2,a1z(z)=a1cos2πpzz,a1=4m,a2=2m2,a0、a1z(z)和a2分别为基频信号的权重、一次谐波的权重和二次谐波的权重,其中只有一次谐波的权重与z有关;
b)当结构光照射到样本上时,相机接收的荧光分布即原始图像表示为
其中,S(r,z)表示样本本身的荧光分布,H(r,z)表示三维点拓展函数;将其变换到三维频域得到
其中H(kx,y,kz)为三维频域点拓展函数,S(kx,y,kz)为频域荧光分布,kx,y和kz分别为xoy平面和yoz平面的频谱的单位矢量,由于a1z(z)在z方向的傅里叶变换表示为a1z(kz)=a1·
[δ(kz-pz)+δ(kz+pz)],因此中a1(kz)·S(kx,y±px,y,kz)的三维频域又表示为:a1·[S(kx,y±px,y,kz-pz)+S(kx,y±px,y,kz+pz)],从而的三维频域表示为
c)对原始图像栈中的原始图像进行边缘平滑操作,同时进行xoy平面的P倍上采样处理(其中P为≥2的正整数),得到同一调制方向上5个相位的原始图像的傅里叶频域为k为三维频域的频谱;利用每个调制方向不同相位的频率分量求解不同级次的频谱信息,列出求解5个三维频谱分量
的表达式表示为
其中相位分离矩阵
6)将步骤5)中得到的5个三维频谱分量n=0,±1,±2基于估计的频率矢量进行频谱搬移,每个三维频谱分量被搬移到正确的位置Cns(k),其中
F-1[·]为三维傅里叶逆变换的符号;接着通过各级三维频谱分量之间的相位进行初始相位补偿;
7)构建第一和第二三维陷波函数notch1(x,y,z,n)和notch2(x,y,z,n):
为了抑制高频频率峰值点带来的伪影,构建三维陷波函数对进行了初始相位补偿后的各级三维频谱分量进行陷波处理,第一三维陷波函数
和第二三维陷波函数
是根据估计的xoy平面的频率矢量px,y与xoz平面频率矢量pz来设计的针对各级三维频谱分量的三维陷波函数,其中d1和d2分别为第一和第二陷波深度,w1和w2分别为陷波范围;各级三维频谱分量的三维陷波函数经过第一三维陷波函数陷波后的超分辨图像频谱notch1(x,y,z,n)表示基于notch1(x,y,z,0)搬移到n级频谱相应位置的陷波函数;
8)基于第二三维陷波函数和OTF构建滤波器的中间函数 其中OTF(x,y,z,n)表示OTF搬移到n级频谱相应位置的光传递函数,m(n)表示各个级次之间的强度关系;根据中间函数构建第一和第二空间频域滤波器和其中第一空间频域滤波器参数ω1数值较大,第二空间频域滤波器参数ω2数值较小,Apo为三维切趾函数;第一空间频域滤波器Filter1抑制了对各级频谱外围的分量,从而很好的降低了高频噪声带来的伪影,第二空间频域滤波器Filter2提高了三维频域外围的信息权重,同时补充了由于陷波带来的高频尖峰的缺失,从而提高了超分辨图像的分辨率;
9)利用第一和第二空间频域滤波器先后对陷波后的陷波后的超分辨图像频谱CSR_1(k)进行频域优化处理:
先采用第一空间频域滤波器优化得到的频谱为然后采用第二空间频域滤波器优化得到的频谱为最终通过三维傅里叶逆变换得到最终的超分辨图像ISR=F-1[CSR_3(k)]。
其中,在步骤2)的c)中,第一标准阈值E1_max的范围为0.1~10,第二标准阈值E2_max的范围为0.1~10。
在步骤7)中,第一和第二陷波深度d1和d2为0.2~1,第一和第二陷波范围w1和w2为0~5。
在步骤8)中,当三维OTF和Apo归一化后,第一空间频域滤波器参数ω1数值大于0.4;第二空间频域滤波器参数ω2数值小于0.2;ω1越小,对伪影抑制效果越明显,ω2越小,对分辨率提升越大,在实际中需要考虑重建图像的质量,达到参数选择的平衡,但一般来说默认参数能够在大多数情况下达到较好的结果;三维切趾函数Apo根据估计的频率矢量确定。
本发明的优点:
本发明首先通过多级频谱协同估计的方法,极大的提高了低信噪比、低调制度情况下参数估计的稳定性与准确性,避免了参数估计错误带来的伪影;同时通过构建三维陷波函数和两步空间频域滤波,极大的降低了由于频谱尖峰、高频噪声等因素带来的重建伪影,同时提高了图像的分辨率;本发明在保真度、伪影抑制和散焦信号消除等方面优于其他算法,尤其是在低信噪比情况下有突出的重建效果。
附图说明
图1为根据本发明的自适应低伪影的多层三维结构光超分辨显微成像重建方法的一个实施例得到的参数估计与频谱优化的原理示意图,其中,(a)为低信噪比条件下的空域图像以及频谱分离后的+1级和+2级频谱,(b)为仿真的光学传递函数、陷波函数、中间变量函数、切趾函数以及第一和第二空间频域滤波器的xoy平面与yoz平面示意图,每个图像下方是沿着频率矢量方向的一维强度分布图;
图2为本发明的自适应低伪影的多层三维结构光超分辨显微成像重建方法的流程图;
图3为根据本发明自适应低伪影的多层三维结构光超分辨显微成像重建方法的一个实施例得到的进行两步骤频域优化的作用效果图;
图4为根据本发明自适应低伪影的多层三维结构光超分辨显微成像重建方法的一个实施例得到的对比图,其中,(a)为本发明的方法与其余算法对微丝样本在不同信噪比情况下的重建对比图,(b)为本发明的方法与其余算法(算法1和算法2)对微丝样本在不同信噪比情况下的平均绝对值误差与信噪比对比图,(c)为本发明的方法与算法1对微管样本在低信噪比情况下的重建对比图;(d)为本发明的方法与算法2对微丝样本在低信噪比情况下的重建对比图。
具体实施方式
下面结合附图,通过具体实施例,进一步阐述本发明。
在本实施例中,采用人骨肉瘤细胞的微丝结构作为样本,采用M=3,N=5的结构光对样本进行照明,本实施例的自适应低伪影的多层三维结构光超分辨显微成像重建方法,如图2所示,包括以下步骤:
1)得到原始图像栈:
结构光对样品的照射具有3个调制方向,即M=3,每一个调制方向具有5个相位,在每一个调制方向上的每一个相位均采集一张原始图像,从而在样本的一层得到3×5张原始图像,利用平移升降台驱动样本进行z方向的逐层扫描,得到一共7层图像序列,即T=7,构成包含3×5×7张原始图像的原始图像栈;
2)采用多级估计法对频率矢量进行估计:
a)计算每一层图像序列的调制对比度,并将所有的T层图像序列以调制对比度作为权重进行加权求和,得到一层调制对比度最高的图像序列,此单层图像序列包括3×5张原始图像;
b)将对比度最高的图像序列中的每长原始图像进行多次反卷积处理,从而突出结构光调制后样本的高频频率峰值点,接着设置初始相位为0,对每个调制方向的结构光调制后的5张原始图像利用二维相位分离矩阵进行分离频谱,得到±2,±1,0共5个二维频谱分量C2D n(kx,y),n=0,±1,±2,kx,y为xoy平面的频谱的单位矢量;
c)估计结构光照明条纹的频率矢量:
根据步骤b)得到的0,±1,±2级二维频谱分量,采用+2,+1级频谱对结构光照明条纹的频率矢量估计进行逐级估计;
i.通过+2级频谱进行频率矢量估计:
通过+2级频谱与0级频谱的互相关迭代操作,估计M个调制方向上的结构光照明条纹的频率矢量分别为通过反余弦函数定义求得频率矢量所对应的空间角度为θ=θ1,θ2,…,θM;
ii.建立评价指标:
由于理想的正确的频率矢量的模值应当相等,其空间角度应当相差2π/M,与该标准相差越大,说明频率矢量估计的准确度越低,因此建立第一和第二评价指标E1和E2,分别表征频率矢量pθ和空间角度的估计准确程度,频率矢量的模值越接近且角度之差越均匀,证明通过n级频谱进行频率矢量估计的结果越准确,其中第一评价指标第二评价指标E2=std{|θ2-θ1|,|θ3-θ2|,…,|θ1-θM|},std为方差符号;并且,设置第一和第二标准阈值E1_max=0.5和E2_max=5;
iii.估计判定:
如果满足判定条件E1≤E1_max且E2≤E2_max则证明通过+2级频谱估计的频率矢量正确,则得到估计的频率矢量并得到频率矢量所对应的空间角度,进入步骤3);
如果E1>E1_max或E2>E2_max则证明通过+2级频谱估计的频率矢量不正确,由于结构光照明条纹的+1级条纹的调制度的强度通常强于+2级条纹的调制度,因此若通过+2级频谱的频率估计不正确,则采用+1级频谱进行协同估计,通过
+1级频谱与0级频谱的互相关迭代操作估计频率矢量,并建立相应的第一和第二评价指标,如果满足判定条件,则得到估计的频率矢量并得到频率矢量所对应的空间角度,进入步骤3);如果通过+2和+1级频谱估计的频率矢量都无法满足判定条件E1≤E1_max且E2≤E2_max,则选择最小的级次估计的频率矢量,并得到频率矢量所对应的空间角度,进入步骤3);
对所有调制方向上的频率矢量的模求均值得到xoy平面的调制频率矢量的模值
px,y为xoy平面的调制频率矢量的模值;
3)通过级次之间的强度关系估计不同角度结构光的调制深度m;并且通过级次之间重叠部分的复数差异估计不同角度结构光的初始相位
4)基于输入的物理参数仿真出三维光转换函数(OTF)或直接输入实际的三维OTF;
5)进行频域分离:
a)三维结构光强度的空间分布表示为 其中x,y,z分别为三维空间上xyz方向的单位矢量,r为xoy平面的极坐标的极径,I0表示光强大小,m表示调制深度,和分别表示xoy平面和yoz平面的调制频率矢量的模值,ψ是光束的入射角,λ为激发波长,n0为介质折射率,根据
求得光束的入射角ψ,再根据得到yoz平面的调制频率矢量pz,整理得到 从而结构光表示为基频信号、一次谐波与二次谐波的组合,a0=1+2m2,a1z(z)=a1cos2πpzz,a1=4m,a2=2m2,a0、a1z(z)和a2分别为基频信号的权重、一次谐波的权重和二次谐波的权重,其中只有一次谐波的权重与z有关;
b)当结构光照射到样本上时,相机接收的荧光分布即原始图像表示为 其中,S(r,z)表示样本本身的荧光分布,H(r,z)表示三维点拓展函数;将其变换到三维频域得到
其中H(kx,y,kz)为三维频域点拓展函数,S(kx,y,kz)为频域荧光分布,kx,y和kz分别为xoy平面和yoz平面的频谱的单位矢量,由于a1z(z)在z方向的傅里叶变换表示为a1z(kz)=a1·
[δ(kz-pz)+δ(kz+pz)],因此中a1(kz)·S(kx,y±px,y,kz)的三维频域又表示为:a1·[S(kx,y±px,y,kz-pz)+S(kx,y±px,y,kz+pz)],从而
的三维频域表示为
a)对原始图像栈中的原始图像进行边缘平滑操作,同时进行xoy平面的两倍上采样处理,得到同一调制方向上5个相位的原始图像的傅里叶频域为
k为三维频域的频谱;利用每个调制方向不同相位的频率分量求解不同级次的频谱信息,列出求解5个三维频谱分量的表达式表示为
其中相位分离矩阵
6)将步骤5)中得到的N个三维频谱分量n=0,±1,±2基于估计的频率矢量进行频谱搬移,每个三维频谱分量被搬移到正确的位置Cns(k),其中
F-1[·]为三维傅里叶逆变换的符号;接着通过各级三维频谱分量之间的相位进行初始相位补偿;
7)构建三维陷波函数:为了抑制高频频率峰值点带来的伪影,根据估计的xoy平面的频率矢量px,y与xoz平面频率矢量pz来构建第一三维陷波函数 和第二三维陷波函数 对进行了初始相位补偿后的各级三维频谱分量进行陷波处理,其中d1和d2分别为第一和第二陷波深度,其中第一和第二陷波深度的取值互相独立分别为d1=0.98和d2=0.92,第一和第二陷波范围的取值互相独立分别为w1=0.5和w2=0.4;经过陷波后的超分辨图像频谱notch(x,y,z,n)表示基于notch(x,y,z,0)搬移到n级频谱相应位置的陷波函数;
8)基于第二三维陷波函数和OTF构建滤波器的中间函数 其中OTF(x,y,z,n)表示OTF搬移到n级频谱相应位置的光传递函数,m(n)表示各个级次之间的强度关系;根据中间函数构建第一和第二空间频域滤波器和其中第一空间频域滤波器参数ω1=0.5,第二空间频域滤波器参数ω2=0.1,Apo为三维切趾函数,根据和分别表示xoy平面和yoz平面的调制频率矢量px,y和pz确定,如图1(b)所示;第一空间频域滤波器Filter1抑制了对各个级次的频谱外围的分量,从而很好的降低了高频噪声带来的伪影,第二空间频域滤波器Filter2提高了三维频域外围的信息权重,同时补充了由于陷波带来的高频尖峰的缺失,从而提高了超分辨图像的分辨率;
9)利用第一和第二空间频域滤波器先后对陷波后的陷波后的超分辨图像频谱CSR_1(k)进行频域优化处理:
第一步采用第一空间频域滤波器优化得到的频谱为
第二步采用第二空间频域滤波器优化得到的频谱为最终通过三维傅里叶逆变换得到最终的超分辨图像ISR=F-1[CSR_3(k)],如图3所示。
图4给出了本发明的方法与其他方法的对比效果,图4中,算法1为SIMnoise,算法2为OMX。在图4(a)中,采样了四个梯度的光强激发微丝样本,能够看出本发明的方法在各种信噪比情况下的重建效果优于OMX与SIMnoise,有着更少的伪影与离焦背景。以各个算法的最高信噪比图像作为基准,计算平均绝对值误差,同时计算各个图像的信噪比数据,发现本发明的方法在各个信噪比的情况下优于其他算法,如图4(b)所示。进一步采集了低信噪比情况下的微管样本和微丝样本分别与SIMnoise和OMX比较,仍然发现本发明的方法优于其他算法,如图4(c)和(d)所示。
最后需要注意的是,公布实施例的目的在于帮助进一步理解本发明,但是本领域的技术人员可以理解:在不脱离本发明及所附的权利要求的精神和范围内,各种替换和修改都是可能的。因此,本发明不应局限于实施例所公开的内容,本发明要求保护的范围以权利要求书界定的范围为准。
Claims (4)
1.一种自适应低伪影的多层三维结构光超分辨显微成像重建方法,其特征在于,所述成像重建方法包括以下步骤:
1)得到原始图像栈:
结构光对样品的照射具有M个调制方向,每一个调制方向具有N个相位,在每一个调制方向上的每一个相位均采集一张原始图像,从而在样本的一层得到M×N张原始图像,利用平移升降台驱动样本进行z方向的逐层扫描,得到一共T层图像序列,并在每一层的每个调制方向上取5个不同相位,构成包含M×5×T张原始图像的原始图像栈;
2)采用多级估计法对频率矢量进行估计:
a)计算每一层图像序列的调制对比度,并将所有的T层图像序列以调制对比度作为权重进行加权求和,得到一层调制对比度最高的图像序列,此单层图像序列包括M×5张原始图像;
b)将对比度最高的图像序列中的每长原始图像进行多次反卷积处理,从而突出结构光调制后样本的高频频率峰值点,接着设置初始相位为0,对每个调制方向的结构光调制后的5张原始图像利用二维相位分离矩阵进行分离频谱,得到±2,±1,0级共5个二维频谱分量C2D n(kx,y),n=0,±1,±2,kx,y为xoy平面的频谱的单位矢量;
c)估计结构光照明条纹的频率矢量pθ:
根据步骤b)得到的0,±1,±2级二维频谱分量,采用+2和+1级频谱对结构光照明条纹的频率矢量估计进行逐级估计;
i.通过+2级频谱进行频率矢量估计:
通过+2级频谱与0级频谱的互相关迭代操作,估计M个调制方向上的结构光照明条纹的频率矢量分别为通过反余弦函数定义求得频率矢量所对应的空间角度为θ=θ1,θ2,…,θM;
ii.建立评价指标:
由于理想的正确的频率矢量的模值应当相等,其空间角度应当相差2π/M,与该标准相差越大,说明频率矢量估计的准确度越低,因此建立第一和第二评价指标E1和E2,分别表征频率矢量pθ和空间角度的估计准确程度,频率矢量的模值越接近且角度之差越均匀,证明通过n级频谱进行频率矢量估计的结果越准确,
其中第一评价指标第二评价指标E2=std{|θ2-θ1|,|θ3-θ2|,…,|θ1-θM|},std为方差符号;并且,设置第一和第二标准阈值E1_max和E2_max;
iii.估计判定:
如果满足判定条件E1≤E1_max且E2≤E2_max则证明通过+2级频谱估计的频率矢量正确,则得到估计的频率矢量并得到频率矢量所对应的空间角度,进入步骤3);
如果E1>E1_max或E2>E2_max则证明通过+2级频谱估计的频率矢量不正确,由于结构光照明条纹的+1级条纹的调制度的强度通常强于+2级条纹的调制度,
因此若通过+2级频谱的频率估计不正确,则采用+1级频谱进行协同估计,通过+1级频谱与0级频谱的互相关迭代操作估计频率矢量,并建立相应的第一和第二评价指标,如果满足判定条件,则得到估计的频率矢量并得到频率矢量所对应的空间角度,进入步骤3);
如果通过+2和+1级频谱估计的频率矢量都无法满足判定条件E1≤E1_max且E2≤E2_max,则选择最小的级次估计的频率矢量,并得到频率矢量所对应的空间角度,进入步骤3);
对所有调制方向上的频率矢量的模求均值得到xoy平面的调制频率矢量的模值
3)通过级次之间的强度关系估计不同角度结构光的调制深度m;并且通过级次之间重叠部分的复数差异估计不同角度的结构光的初始相位
4)基于输入的物理参数仿真出三维光转换函数OTF或直接输入实验测得的三维OTF;
5)进行频域分离:
a)三维结构光强度的空间分布表示为 其中x,y和z分别为三维空间上xyz方向的单位矢量,r为xoy平面的极坐标的极径,I0表示光强大小,m表示调制深度,和分别表示xoy平面和yoz平面的调制频率矢量的模值,ψ是光束的入射角,λ为激发波长,n0为介质折射率,根据求得光束的入射角ψ,再根据得到yoz平面的调制频率矢量pz,整理得到 从而结构光表示为基频信号、一次谐波与二次谐波的组合,a0=1+2m2,a1z(z)=a1cos2πpzz,a1=4m,a2=2m2,a0、a1z(z)和a2分别为基频信号的权重、一次谐波的权重和二次谐波的权重,其中只有一次谐波的权重与z有关;
b)当结构光照射到样本上时,相机接收的荧光分布即原始图像表示为 其中,S(r,z)表示样本本身的荧光分布,H(r,z)表示三维点拓展函数;将其变换到三维频域得到
其中H(kx,y,kz)为三维频域点拓展函数,S(kx,y,kz)为频域荧光分布,kx,y和kz分别为xoy平面和yoz平面的频谱的单位矢量,由于a1z(z)在z方向的傅里叶变换表示为a1z(kz)=a1·[δ(kz-pz)+δ(kz+pz)],因此中a1(kz)·S(kx,y±px,y,kz)的三维频域又表示为:a1·[S(kx,y±px,y,kz-pz)+S(kx,y±px,y,kz+pz)],从而的三维频域表示为
c)对原始图像栈中的原始图像进行边缘平滑操作,同时进行xoy平面的P倍上采样处理,P为≥2的正整数,得到同一调制方向上5个相位的原始图像的傅里叶频域为k为三维频域的频谱;利用每个调制方向不同相位的频率分量求解不同级次的频谱信息,列出求解5个三维频谱分量的表达式表示为
其中相位分离矩阵
6)将步骤5)中得到的5个三维频谱分量基于估计的频率矢量进行频谱搬移,每个三维频谱分量被搬移到正确的位置Cns(k),其中
F-1[·]为三维傅里叶逆变换的符号;接着通过各级三维频谱分量之间的相位进行初始相位补偿;
7)构建第一和第二三维陷波函数notch1(x,y,z,n)和notch2(x,y,z,n):
为了抑制高频频率峰值点带来的伪影,根据估计的xoy平面的频率矢量px,y与xoz平面频率矢量pz来构建第一三维陷波函数 和第二三维陷波函数 对进行了初始相位补偿后的各级三维频谱分量进行陷波处理,其中d1和d2分别为第一和第二陷波深度,w1和w2分别为陷波范围;各级三维频谱分量的三维陷波函数经过第一三维陷波函数陷波后的超分辨图像频谱 notch1(x,y,z,n)表示基于notch1(x,y,z,0)搬移到n级频谱相应位置的陷波函数;
8)基于第二三维陷波函数和OTF构建滤波器的中间函数 其中OTF(x,y,z,n)表示OTF搬移到n级频谱相应位置的光传递函数,m(n)表示各个级次之间的强度关系;根据中间函数构建第一和第二空间频域滤波器和其中第一空间频域滤波器参数ω1数值较大,第二空间频域滤波器参数ω2数值较小,Apo为三维切趾函数;第一空间频域滤波器Filter1抑制了对各级频谱外围的分量,从而很好的降低了高频噪声带来的伪影,第二空间频域滤波器Filter2提高了三维频域外围的信息权重,同时补充了由于陷波带来的高频尖峰的缺失,从而提高了超分辨图像的分辨率;
9)利用第一和第二空间频域滤波器先后对陷波后的陷波后的超分辨图像频谱CSR_1(k)进行频域优化处理:
先采用第一空间频域滤波器优化得到的频谱为然后采用第二空间频域滤波器优化得到的频谱为最终通过三维傅里叶逆变换得到最终的超分辨图像ISR=F-1[CSR_3(k)]。
2.如权利要求1所述的成像重建方法,其特征在于,在步骤2)的c)中,第一标准阈值E1_max的范围为0.1~10,第二标准阈值E2_max的范围为0.1~10。
3.如权利要求1所述的成像重建方法,其特征在于,在步骤7)中,第一和第二陷波深度d1和d2为0.2~1,第一和第二陷波范围w1和w2为0~5。
4.如权利要求1所述的成像重建方法,其特征在于,在步骤8)中,第一空间频域滤波器参数ω1数值大于0.4;第二空间频域滤波器参数ω2数值小于0.2。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211605447.2A CN116109762A (zh) | 2022-12-14 | 2022-12-14 | 自适应低伪影的多层三维结构光超分辨显微成像重建方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211605447.2A CN116109762A (zh) | 2022-12-14 | 2022-12-14 | 自适应低伪影的多层三维结构光超分辨显微成像重建方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN116109762A true CN116109762A (zh) | 2023-05-12 |
Family
ID=86258817
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202211605447.2A Pending CN116109762A (zh) | 2022-12-14 | 2022-12-14 | 自适应低伪影的多层三维结构光超分辨显微成像重建方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN116109762A (zh) |
-
2022
- 2022-12-14 CN CN202211605447.2A patent/CN116109762A/zh active Pending
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111145089B (zh) | 高保真图像重构方法、系统、计算机设备和存储介质 | |
Pertuz et al. | Analysis of focus measure operators for shape-from-focus | |
Forster et al. | Complex wavelets for extended depth‐of‐field: A new method for the fusion of multichannel microscopy images | |
Cannell et al. | Image enhancement by deconvolution | |
JP2022538406A (ja) | 位相分布を強調するための画像処理装置および方法 | |
Kulkarni et al. | Fringe denoising algorithms: a review | |
Roels et al. | An overview of state‐of‐the‐art image restoration in electron microscopy | |
EP3692360A1 (en) | Methods and systems for 3d structure estimation using non-uniform refinement | |
Forster et al. | Extended depth-of-focus for multi-channel microscopy images: a complex wavelet approach | |
CN110838093B (zh) | 一种fMOST或MOST显微图像条纹噪声去除方法 | |
Zhi-guo et al. | A wavelet based algorithm for multi-focus micro-image fusion | |
Li et al. | PURE-LET deconvolution of 3D fluorescence microscopy images | |
CN111458317B (zh) | 一种直接结构光照明超分辨显微重建方法 | |
Kim et al. | Blind deconvolution of 3D fluorescence microscopy using depth‐variant asymmetric PSF | |
CN116109762A (zh) | 自适应低伪影的多层三维结构光超分辨显微成像重建方法 | |
Šroubek et al. | Fusion of blurred images | |
CN115839936B (zh) | 一种基于锁相探测的结构光照明超分辨显微成像重构方法 | |
Sadaghiani et al. | Image interpolation based on 2D-DWT with novel regularity-preserving algorithm using RLS adaptive filters | |
Rahman et al. | Comparison of computational methods developed to address depth-variant imaging in fluorescence microscopy | |
Chacko et al. | Fast spatially variant deconvolution for optical microscopy via iterative shrinkage thresholding | |
Mahmood et al. | Measuring focus quality in vector valued images for shape from focus | |
CN116402678B (zh) | 超分辨结构光照明显微镜的频谱优化直接重建方法 | |
Zhang et al. | Simultaneous estimation of spatial frequency and phase based on an improved component cross-correlation algorithm for structured illumination microscopy | |
Xu et al. | Interferogram blind denoising using deep residual learning for phase-shifting interferometry | |
Homem et al. | Poisson noise reduction in deconvolution microscopy |
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 |