CN114881907B - 光学显微图像多景深焦点合成方法、系统及图像处理方法 - Google Patents
光学显微图像多景深焦点合成方法、系统及图像处理方法 Download PDFInfo
- Publication number
- CN114881907B CN114881907B CN202210755701.0A CN202210755701A CN114881907B CN 114881907 B CN114881907 B CN 114881907B CN 202210755701 A CN202210755701 A CN 202210755701A CN 114881907 B CN114881907 B CN 114881907B
- Authority
- CN
- China
- Prior art keywords
- image
- pixel
- channel
- filtering
- image set
- 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
Links
- 230000003287 optical effect Effects 0.000 title claims abstract description 34
- 238000003672 processing method Methods 0.000 title claims abstract description 9
- 238000001308 synthesis method Methods 0.000 title claims description 16
- 238000001914 filtration Methods 0.000 claims abstract description 107
- 238000005259 measurement Methods 0.000 claims abstract description 92
- 238000000034 method Methods 0.000 claims abstract description 61
- 230000004927 fusion Effects 0.000 claims abstract description 47
- 238000004364 calculation method Methods 0.000 claims abstract description 41
- 238000012545 processing Methods 0.000 claims abstract description 38
- 230000002194 synthesizing effect Effects 0.000 claims abstract description 5
- 238000010586 diagram Methods 0.000 claims description 41
- 238000000879 optical micrograph Methods 0.000 claims description 29
- 230000009466 transformation Effects 0.000 claims description 13
- 238000004891 communication Methods 0.000 claims description 10
- 230000008569 process Effects 0.000 claims description 9
- 230000015572 biosynthetic process Effects 0.000 claims description 7
- 230000006870 function Effects 0.000 claims description 7
- 238000003786 synthesis reaction Methods 0.000 claims description 7
- 238000011524 similarity measure Methods 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 claims description 4
- 238000000399 optical microscopy Methods 0.000 claims description 4
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 description 6
- 230000000694 effects Effects 0.000 description 6
- 238000013519 translation Methods 0.000 description 6
- 230000002776 aggregation Effects 0.000 description 5
- 238000004220 aggregation Methods 0.000 description 5
- 238000003384 imaging method Methods 0.000 description 3
- 230000000007 visual effect Effects 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 238000012952 Resampling Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 230000007812 deficiency Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000009499 grossing Methods 0.000 description 2
- 238000007500 overflow downdraw method Methods 0.000 description 2
- YBJHBAHKTGYVGT-ZKWXMUAHSA-N (+)-Biotin Chemical compound N1C(=O)N[C@@H]2[C@H](CCCCC(=O)O)SC[C@@H]21 YBJHBAHKTGYVGT-ZKWXMUAHSA-N 0.000 description 1
- 206010013647 Drowning Diseases 0.000 description 1
- VYPSYNLAJGMNEJ-UHFFFAOYSA-N Silicium dioxide Chemical compound O=[Si]=O VYPSYNLAJGMNEJ-UHFFFAOYSA-N 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
- FEPMHVLSLDOMQC-UHFFFAOYSA-N virginiamycin-S1 Natural products CC1OC(=O)C(C=2C=CC=CC=2)NC(=O)C2CC(=O)CCN2C(=O)C(CC=2C=CC=CC=2)N(C)C(=O)C2CCCN2C(=O)C(CC)NC(=O)C1NC(=O)C1=NC=CC=C1O FEPMHVLSLDOMQC-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T1/00—General purpose image data processing
- G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration using local operators
-
- 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/10—Image acquisition modality
- G06T2207/10024—Color image
-
- 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/10—Image acquisition modality
- G06T2207/10052—Images from lightfield camera
-
- 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/20172—Image enhancement details
- G06T2207/20192—Edge enhancement; Edge preservation
-
- 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/20212—Image combination
- G06T2207/20221—Image fusion; Image merging
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
本发明涉及一种光学显微图像多景深焦点合成方法、系统及图像处理方法。该方法主要是将图像集合由RGB空间转换成LAB空间,获得L通道图像集合的每张图像的对焦测量图，对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波处理，对未滤波的对焦测量图集合进行级联滤波处理，对滤波后对焦测量图集合进行最大密度投影计算，获得像素索引图，基于像素索引图获得L通道融合图像、A通道融合图像和B通道融合图像，基于L通道融合图像、A通道融合图像和B通道融合图像，将LAB颜色空间转换成RGB颜色空间，获得融合后的光学显微彩色图像。本申请获得的融合图像更加平滑聚合，融合质量得到提升。
Description
技术领域
本发明涉及图像处理领域，特别是涉及光学显微图像多景深焦点合成方法、系统、计算机装置、存储介质及图像处理方法。
背景技术
许多检测鉴定行业需要借助光学显微成像技术在数十倍到数百倍的放大尺度下观测目标对象的微观结构特征。由于观测对象具有三维形态，而普通光学显微镜在单一固定位置的景深对焦能力有限，超出范围会由于衍射现象造成模糊。通常在实际应用中需要手动调节镜头在不同轴向位置实现局部对焦。该方式依赖手动操作，耗费时力，且无法一次获得当前视野下所有目标的清晰形貌。
目前，一些中高端显微镜成像系统会采用多景深图像融合方法来解决上述问题。该方法是通过获取显微镜根据预设采集顺序采集到的不同高度载物平台上同一样本的一组光学显微图像，将该组光学显微图像通过数字图像处理技术获得一幅全对焦图像。传统的多景深图像融合方法的缺点是：噪声干扰较多，图像平滑度较差，尤其是图像边缘纹理的融合效果较差。
发明内容
基于上述问题，提供一种光学显微图像多景深焦点合成方法。该方法可有效提高图像的融合效果。
一种光学显微图像多景深焦点合成方法，包括：
将光学显微图像集合由RGB颜色空间转换成LAB颜色空间，进而获得L通道图像集合，A通道图像集合和B通道图像集合，
计算获得L通道图像集合的每张图像的对焦测量图，形成未滤波的对焦测量图集合，
对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波处理，获得滤波后L通道图像集合、滤波后A通道图像集合、滤波后B通道图像集合，
对未滤波的对焦测量图集合进行级联滤波处理，获得滤波后对焦测量图集合，对滤波后对焦测量图集合进行最大密度投影计算，获得用于融合的像素索引图，
基于像素索引图、滤波后L通道图像集合、滤波后A通道图像集合和滤波后B通道图像集合，获得L通道融合图像、A通道融合图像和B通道融合图像，
基于L通道融合图像、A通道融合图像和B通道融合图像，将LAB颜色空间转换成RGB颜色空间，获得融合后的光学显微彩色图像。
本申请的上述方法通过将RGB颜色空间转换成LAB颜色空间，减少了多景深图像集合颜色差异的干扰，同时，对未滤波的对焦测量图集合进行级联滤波处理，以及对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波处理，可降低噪声的干扰，使得最终获得的光学显微彩色图像更加平滑聚合，尤其是图像边缘纹理更加清晰，融合质量得到提升。上述光学显微彩色图像也就是全对焦图像。
在其中一个实施例中，所述计算获得L通道图像集合的每张图像的对焦测量图，具体包括：通过改进的拉普拉斯能量交叉和算法来计算获得L通道图像集合的每张图像的对焦测量图，所述改进的拉普拉斯能量交叉和算法应用了以下公式：
其中，L(x, y)是像素(x, y)位置处的对焦测量值，(ε, η)是以像素(x, y)为中心的邻域U内一点，
为像素(ε, η)的水平方向二阶导的绝对值，
为像素(ε, η)的垂直方向二阶导的绝对值，
为像素(ε, η)的正45度方向二阶导的绝对值，
为像素(ε, η)的负45度方向二阶导的绝对值。r是拉普拉斯能量交叉和算子函数的窗口大小。
在其中一个实施例中，所述对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波处理，具体包括：基于对未滤波的对焦测量图集合进行三维加权中值滤波而获得的每个像素以其为中心的邻域内中值像素pK的位置，对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波，也就是将所述邻域内中值像素pK的值代替L通道图像集合、A通道图像集合和B通道图像集合上对应的邻域内中心像素的值，
所述三维加权中值滤波的计算公式为：
，该公式表达的意思是用邻域Ω内的中值像素pK的值来替换中心像素的值，中心像素也就是需要通过三维加权中值滤波处理的像素，其中的中值像素pK的确定方法为：
，邻域Ω内的每个像素pk相对中心像素的相似性测度值即为pk的权重，记为
，所述邻域Ω为三维区域。P*是经过滤波后的中心像素值。
在其中一个实施例中，所述对未滤波的对焦测量图集合进行级联滤波处理，获得滤波后对焦测量图集合，对滤波后对焦测量图集合进行最大密度投影计算，获得用于融合的像素索引图，具体包括：
先对未滤波的对焦测量图集合进行三维加权中值滤波处理，获得一次滤波后对焦测量图集合，对一次滤波后对焦测量图集合进行最大密度投影计算，获得导向滤波引导图，
再基于滤波引导图对一次滤波后对焦测量图集合进行高斯卷积导向滤波处理，进而获得二次滤波后对焦测量图集合，对二次滤波后对焦测量图集合进行最大密度投影计算，获得用于融合的像素索引图，
所述三维加权中值滤波的计算公式为：
，该公式表达的意思是用邻域Ω内的中值像素pK的值来替换中心像素的值，中心像素也就是需要通过三维加权中值滤波处理的像素，其中的中值像素pK的确定方法为：
，邻域Ω内的每个像素pk相对中心像素的相似性测度值即为pk的权重，记为
，所述邻域Ω为三维区域，
所述高斯卷积导向滤波的计算公式为：
其中，Ii和Pi分别为滤波引导图和一次滤波后对焦测量图集合的图像中的邻域
内的像素值，
为高斯卷积核，ϵ是为防止ak过大的正则化参数。ak和bk分别是通过高斯卷积导向滤波的计算公式计算得到的每个像素的线性变换参数。
在其中一个实施例中，所述高斯卷积导向滤波处理的计算公式还包括：
，其中
和
是以当前像素i为中心，对高斯卷积核的半径内所有像素j的线性变换参数
和
分别求平均。
在其中一个实施例中，在将光学显微图像集合由RGB颜色空间转换成LAB颜色空间之前，先对光学显微图像集合进行空间对齐。
在其中一个实施例中，所述光学显微图像集合的获取方法包括：按步进改变镜头距离样本高度依次采集得到若干张光学显微图像，所述若干张光学显微图像形成光学显微图像集合。
在其中一个实施例中，所述按步进改变镜头距离样本高度依次采集得到若干张光学显微图像是在单一视野下，沿轴向，以固定步进距离移动镜头来对样本进行拍摄。
一种光学显微图像多景深焦点合成系统，包括：
数据获取模块，所述数据获取模块用于获取光学显微图像集合，
数据处理模块，所述数据处理模块用于根据所述的方法来处理光学显微图像集合。
一种计算机装置，包括：处理器、存储器、通信接口和通信总线，所述处理器、存储器和通信接口通过所述通信总线完成相互间的通信，所述存储器用于存放至少一个可执行指令，所述可执行指令使所述处理器执行所述的光学显微图像多景深焦点合成方法对应的操作。
一种计算机存储介质，所述计算机存储介质中存储有至少一个可执行指令，所述可执行指令使处理器执行所述的光学显微图像多景深焦点合成方法对应的操作。
一种基于CPU/GPU异构计算的图像处理方法，根据所述的光学显微图像多景深焦点合成方法设定计算任务，通过CPU/GPU异构计算方法处理所述计算任务。
附图说明
图1为本申请的一个实施例的光学显微图像多景深焦点合成方法的流程图。
图2为根据最大密度投影得到的像素索引图来对多景深拍摄图像集合进行融合的示意图。其中，沿着轴向（Z轴）对所有对焦测量图做最大密度投影得到像素索引图，该索引图中每个像素值为该像素沿轴向对焦程度最好的图像层索引值，最后基于索引图从原始图像集合提取对应像素灰度信息形成融合图像。
图3为以像素(x, y)为中心的三维邻域Ω的示意图。该图中，三维邻域Ω为一个5×5×3的邻域。
图4为利用各种滤波方法以及本申请提出的级联滤波方法得到的像素索引图。其中，（a）图是只利用二维加权中值滤波得到的一组对焦测量图投影得到的像素索引图；（b）图是利用二维加权中值滤波和常规导向滤波级联进行滤波得到的一组对焦测量图投影得到的像素索引图；（c）图是只利用三维加权中值滤波得到的一组对焦测量图投影得到的像素索引图；（d）图是本申请提出的利用三维加权中值滤波和高斯卷积导向滤波处理得到的一组对焦测量图投影得到的像素索引图。通过对比可知，（d）图相比其它图在平滑与聚合效果上更好。
图5为本申请提出的级联滤波方法在各阶段得到的像素索引图。从左往右：不做滤波处理的对焦测量图集合投影得到的像素索引图；只做三维加权中值滤波的对焦测量图集合投影得到的像素索引图；采用三维加权中值滤波以及高斯卷积导向滤波处理后的对焦测量图集合投影得到的像素索引图。
图6为一组多景深硅藻光镜图像的融合结果图。其中，前两排为该组多景深光学显微图像集合，第三排为不同阶段的融合结果图（从左往右依次是：未做滤波处理直接生成的像素索引图用来融合得到的结果图；三维加权中值滤波后生成的像素索引图用来融合得到的结果图；采用三维加权中值滤波与高斯卷积导向滤波级联滤波后生成的像素索引图用来融合得到的结果图），第四排为对应结果图的局部放大展示图。
图7为本申请的实施例的CPU/GPU异构计算执行流程图。
图8为本申请的实施例的CUDA并行计算框架中的线程Block和Grid概念的示意图。其中，图8中左边的图表示一副48×64的图像划分为12个GPU线程Block来处理，即Grid大小为3×4。图8中右边的图表示一个Block包含16×16个线程，每个线程处理一个像素的相关计算。
具体实施方式
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。在下面的描述中阐述了很多具体细节以便于充分理解本发明。但是本发明能够以很多不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似改进，因此本发明不受下面公开的具体实施例的限制。
需要说明的是，当元件被称为“固定于”另一个元件，它可以直接在另一个元件上或者也可以存在居中的元件。当一个元件被认为是“连接”另一个元件，它可以是直接连接到另一个元件或者可能同时存在居中元件。
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。
如图1所示，本申请的实施例提供了一种光学显微图像多景深焦点合成方法，该方法包括：
将光学显微图像集合由RGB颜色空间转换成LAB颜色空间，进而获得L通道图像集合，A通道图像集合和B通道图像集合。
计算获得L通道图像集合的每张图像的对焦测量图，形成未滤波的对焦测量图集合。
对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波处理，获得滤波后L通道图像集合、滤波后A通道图像集合、滤波后B通道图像集合。
对未滤波的对焦测量图集合进行级联滤波处理，获得滤波后对焦测量图集合，对滤波后对焦测量图集合进行最大密度投影计算，获得用于融合的像素索引图。
基于像素索引图、滤波后L通道图像集合、滤波后A通道图像集合、滤波后B通道图像集合获得L通道融合图像、A通道融合图像和B通道融合图像。
基于L通道融合图像、A通道融合图像和B通道融合图像，将LAB颜色空间转换成RGB颜色空间，获得融合后的光学显微彩色图像。上述光学显微彩色图像也就是全对焦图像。
本申请的上述方法中，对于RGB三通道彩色图像，首先会转换到LAB图像空间，从L（亮度）通道图像计算对应的对焦测量图，L、A（颜色）、B（颜色）三通道图像则根据对焦测量图集合进行最大密度投影得到的像素索引图分别进行融合，最后再从LAB图像空间转换回RGB图像空间。
在其中一个实施例中，所述计算获得L通道图像集合的每张图像的对焦测量图，具体包括：通过改进的拉普拉斯能量交叉和算法来计算获得L通道图像集合的每张图像的对焦测量图，所述改进的拉普拉斯能量交叉和算法应用了以下公式：
其中，L(x, y)是像素(x, y)位置处的对焦测量值，(ε, η)是以像素(x, y)为中心的邻域U内一点，
为像素(ε, η)的水平方向二阶导的绝对值，
为像素(ε, η)的垂直方向二阶导的绝对值，
为像素(ε, η)的正45度方向二阶导的绝对值，
为像素(ε, η)的负45度方向二阶导的绝对值。r是拉普拉斯能量交叉和算子函数的窗口大小。
具体的，图像的对焦测量图用来衡量该图像每个像素的对焦程度，常用的计算方法包含图像的一阶导（梯度算子）或二阶导（拉普拉斯算子）。本申请的改进的拉普拉斯能量交叉和方法是对采用图像二阶导求取对焦测量图方法的一种改进。具体应用了以下公式：
需要说明的是，上述(ε, η)是以像素(x, y)为中心的邻域U内一点，通过上述公式计算该点水平方向、垂直方向以及±45°方向的二阶导绝对值和，并对邻域U内的所有点进行累加。本申请的上述方法具有更鲁棒的特性。
在其中一个实施例中，所述对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波处理，具体包括：基于对未滤波的对焦测量图集合进行三维加权中值滤波而获得的每个像素以其为中心的邻域内中值像素pK的位置，对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波，也就是将所述邻域内中值像素pK的值代替L通道图像集合、A通道图像集合和B通道图像集合上对应的邻域内中心像素的值。
由于未滤波的对焦测量图集合的各个像素坐标与L通道图像集合、A通道图像集合和B通道图像集合的各个像素坐标是一一对应的，所以根据未滤波的对焦测量图集合的各个邻域中心像素坐标可以确定L通道图像集合、A通道图像集合和B通道图像集合的对应的邻域中心像素坐标，而根据计算得到的对焦测量图集合的每个邻域中心像素对应的中值像素位置，也可以确定L通道图像集合、A通道图像集合和B通道图像集合的对应的中值像素的坐标，进而利用中值像素的值来替代中心像素的值。
例如，未滤波的对焦测量图集合的第k层图像里某个像素(x, y)，在以这个像素(x, y)为中心像素的三维邻域内做三维加权中值滤波计算，得到这个邻域内的中值像素pK，其位置是(ε, η)，然后就可根据这个信息对L通道图像集合，A通道图像集合，B通道图像集合中对应的第k层图像的这个像素(x, y)进行滤波，即用L通道图像集合，A通道图像集合，B通道图像集合的对应的第k层图像的(ε, η)位置处的中值像素pK的值分别代替L通道图像集合，A通道图像集合，B通道图像集合的对应的第k层图像的像素(x, y)的值。
所述三维加权中值滤波的计算公式为：
，该公式表达的意思是用邻域Ω内的中值像素pK的值来替换中心像素的值，中心像素也就是需要通过三维加权中值滤波处理的像素，其中的中值像素pK的确定方法为：
，邻域Ω内的每个像素pk相对中心像素的相似性测度值即为pk的权重，记为
，所述邻域Ω为三维区域。P*是经过滤波后的中心像素值。
以下举例介绍上述三维加权中值滤波方法。
对当前图像的某个像素(x, y)做三维加权中值滤波的过程为：以该像素(x, y)为中心像素，计算该像素(x, y)所在层以及上下n层的(2∙
+1)×(2∙
+1)×(2∙κ+1)大小邻域Ω内的每个像素(ε, η)∈Ω的权重，
为三维加权中值滤波算子在层内（XY向）的邻域半径，κ为层间（Z向）邻域半径, κ等于n。邻域Ω内的每个像素(ε, η)的权重为这个像素的对焦测量值F(ε, η)与中心像素(x, y)的对焦测量值F(x, y)的相似性测度
，可采用如高斯函数、L1距离等来度量：
a)高斯函数：
,e是高斯函数中的自然常数，σ为高斯函数中的标准差；
b)L1距离：
。
图3给出了一个以像素(x, y)为中心的三维邻域Ω的例子，该例子中，三维邻域Ω为一个5×5×3的邻域。也就是该三维邻域Ω包括像素(x, y)所在层以及上下各1层，每层的邻域范围是5×5。
这里，三维邻域Ω内的每个像素p=(ε, η)相对中心像素(x, y)的相似性测度即为该像素的权重，记为
。而对中心像素(x, y)的三维加权中值滤波处理可按如下步骤来计算：
1)首先对
进行累加得到权重累加和
；
2)对
进行排序，得到排序后的权重系数
；
3)从最小权重（k=1）开始累加直至大于等于
的
时，即为中值（此时的k=K），该步骤的数学描述（s.t.为subject to的意思）为：
；
找到该中值对应的三维邻域Ω内的中值像素pK，用该像素pK的值来替换中心像素(x, y)的值。
虽然本申请的上述改进的拉普拉斯能量交叉和方法可以表征图像对焦的程度并一定程度上抑制噪声，但也存在对焦点误响应的情况，需要进一步聚合平滑才能更好地描述一组图像中每一幅的对焦情况。因此，本申请采用级联方式对未滤波的对焦测量图集合进行滤波。
具体的，所述对未滤波的对焦测量图集合进行级联滤波处理包括：
先对未滤波的对焦测量图集合进行三维加权中值滤波处理，获得一次滤波后对焦测量图集合，对一次滤波后对焦测量图集合进行最大密度投影计算，获得导向滤波引导图；再基于滤波引导图对一次滤波后对焦测量图集合进行高斯卷积导向滤波处理，进而获得二次滤波后对焦测量图集合，对二次滤波后对焦测量图集合进行最大密度投影计算，获得用于融合的像素索引图。
具体的，上述对未滤波的对焦测量图集合进行三维加权中值滤波处理中应用的三维加权中值滤波的具体方法可参考上述对该方法的描述。
以下介绍本申请的高斯卷积导向滤波处理方法。
导向滤波是一种图像滤波技术，通过一张引导图，对输入图像进行滤波处理，使得最后的输出图像大体上与初始图像相似，但是纹理部分与引导图相似。该方法的基础为以像素k为中心的图像邻域Ωk内的线性拟合，其数学模型为：
本申请的引导图为上述滤波引导图，输入图像为一次滤波后对焦测量图集合的图像，对该组图实施最大密度投影得到的融合图作为引导图。上述高斯卷积导向滤波方法是在上述模型的基础上，施加高斯卷积
来求解相应的系数。具体应用了以下公式：
其中，Ii和Pi分别为滤波引导图和一次滤波后对焦测量图集合的图像中的邻域
内的像素值，
为高斯卷积核，ϵ是为防止ak过大的正则化参数。这里，ak和bk分别是通过上述公式计算得到的图像空间中每个像素（位置下标记为k）的线性变换参数，对每个像素都有两个参数ak和bk，因此对整个图像空间来说，是一组参数{ak，bk}。
本申请的上述高斯卷积导向滤波处理方法在邻域内能更好地起到平滑聚合的效果。
由上述高斯卷积导向滤波处理方法的对应公式可以看出，对于输出的每一个像素，都需要对输入图和引导图在窗口内进行高斯卷积计算得到(ak, bk)，这其中存在很多冗余和重复计算。为了克服上述不足，可以通过对经过所有覆盖当前像素的窗口求平均来计算单个像素的输出，也就是：
。
具体的，上述
和
可通过以下公式获得：
上述
表示的是哪些窗口覆盖了当前像素i，其中
和
是以当前像素i为中心，对高斯卷积核的半径R内所有像素j的线性变换参数
和
分别求平均，上述参数
和
也可用以下公式计算获得，也就是下标k换成j：
在其中一个实施例中，在将光学显微图像集合由RGB颜色空间转换成LAB颜色空间之前，先对光学显微图像集合进行空间对齐。
具体的，在拍摄一组光学显微图像过程中，由于系统机械抖动可能会存在非共轴的情况，即图像间会存在微小的平移。此外，由于镜头位置的改变，实际的放大倍率会有微小的变化，即图像间会存在微小的像素尺度差异。为了克服上述不足，可对光学显微图像集合的图像进行空间对齐。图像对齐计算得到的可以是单一的平移或缩放变换，也可以是平移与缩放变换的级联，亦或是结合了平移与缩放的复合图像空间变换Tμ(X)，如以下公式所示：
这里，C为图像中心位置，T为平移量参数，S为缩放因子参数。X代表的是参考图像空间的任一点像素位置。公式的意思是：经过变换Tμ（由缩放因子参数S、平移量参数T、参考图像中心位置C确定，上述参数可以通过相应的配准算法获得）找到待配准图形空间一点，再将该点的像素值填充到参考图像空间的X位置处（该过程即通过图像空间变换进行重采样对齐过程）。
本申请中，图像平移变换采用的方法包括但不限于基于图像灰度匹配的方法，如“归一化互相关”，及频域的方法，如“相位校正”等。
本申请中，图像缩放变换的计算可通过对图像（以中心为原点）进行缩放，再将得到的重采样图像与参考图像进行灰度匹配搜索最优的缩放系数这一策略；也可以基于已知的光学显微镜的焦距与物距信息，利用薄透镜成像公式（Thin Lens Equation）计算出像距，即公式：
。这里，u和v为物距和像距，f为焦距。计算每张图像拍摄时的像距和参考图像对应的像距的比值v/vRef即为缩放系数S。由计算得到的缩放系数对图像进行缩放重采样，本申请中采用常规的双线性插值和B样条插值，但不限于此。
在其中一个实施例中，所述光学显微图像集合的获取方法包括：按步进改变镜头距离样本高度依次采集得到若干张光学显微图像，所述若干张光学显微图像形成光学显微图像集合。
在其中一个实施例中，所述按步进改变镜头距离样本高度依次采集得到若干张光学显微图像是在单一视野下，沿轴向，以固定步进距离移动镜头来对样本进行拍摄。
本申请的实施例还提供了一种光学显微图像多景深焦点合成系统，该系统包括：
数据获取模块，所述数据获取模块用于获取光学显微图像集合，
数据处理模块，所述数据处理模块用于根据所述的方法来处理光学显微图像集合。
本申请的实施例还提供了一种计算机装置，包括：处理器、存储器、通信接口和通信总线，所述处理器、存储器和通信接口通过所述通信总线完成相互间的通信，所述存储器用于存放至少一个可执行指令，所述可执行指令使所述处理器执行所述的光学显微图像多景深焦点合成方法对应的操作。
本申请的实施例还提供了一种计算机存储介质，所述计算机存储介质中存储有至少一个可执行指令，所述可执行指令使处理器执行所述的光学显微图像多景深焦点合成方法对应的操作。
本申请的上述方法涉及到大量计算，为了提高计算效率，本申请的实施例还提供了一种基于CPU/GPU异构计算的图像处理实现方法。该方法具体包括：根据所述的光学显微图像多景深焦点合成方法设定计算任务，通过CPU/GPU异构计算方法处理所述计算任务。
本发明采用了结合CPU与GPU的异构计算来提升效率，CPU代码采用C++语言开发，而GPU代码采用CUDA语言开发，异构执行流程如图7所示。
具体的，首先，所述光学显微图像集合包括K张多景深拍摄图像，开辟N个CPU线程，N<K，每个CPU线程处理大约K/N张图像， 因为所有GPU的计算任务默认都会放在CUDA的第一个流上顺序执行，所以为每个CPU线程显式地创建一个对应的GPU流处理队列，即CUDAStream，这样每个流只会负责该线程内的GPU图像计算任务。
在不同流任务之间需要做线程同步，因为有些后续步骤需要所有线程的前置任务都完成才可执行，比如三维加权中值滤波需要上下多层对焦测量图、高斯卷积导向滤波需要从所有对焦测量图计算出引导图等。以下举一个实施例来详细说明。
计算一组多景深拍摄图像的融合需要做四次线程同步，具体为：
线程同步1：在该步操作前，每个线程首先会载入平均分配的大约K/N张图像（用下标1, 2, ..., s表示）到内存，该线程对应的GPU流处理队列包含如下按顺序执行的分步任务：
a)将第1幅RGB图像转换到LAB图像空间、将LAB三通道图像数据从内存传输到GPU显存空间、计算L通道图像对应的对焦测量图；
b)将第2幅RGB图像转换到LAB图像空间、将LAB三通道图像数据从内存传输到GPU显存空间、计算L通道图像对应的对焦测量图；
c)按照处理第1幅图的步骤计算第3到第s幅图像的对焦测量图。
线程同步2：在该步操作前，需要对每层未滤波的对焦测量图进行三维加权中值滤波，并基于获得的各个中心像素对应的邻域内的中值像素pK的位置，对L通道图像集合的图像、A通道图像集合的图像和B通道图像集合的图像分别进行滤波。
这里以未滤波的对焦测量图集合为例进行说明，每个线程处理对应的多张对焦测量图（同样用下标1, 2, ..., s表示），该线程对应的GPU流处理队列按顺序执行任务：
a)取该线程中第1幅对焦测量图和它上下相邻一层的对焦测量图总共三幅进行加权中值滤波，得到第1幅图的滤波后的对焦测量图；
b)取该线程中第2幅对焦测量图和它上下相邻一层的对焦测量图总共三幅进行加权中值滤波，得到第1幅图的滤波后的对焦测量图；
c)按照处理第1幅图的步骤对第3到第s幅对焦测量图进行同样的滤波。
在此线程同步操作后，会对所有对焦测量图实施高斯卷积导向滤波，这里首先需要从所有对焦测量图计算一幅引导图用于滤波。
线程同步3：在该步操作前，使用线程同步2后计算得到的引导图对每幅对焦测量图进一步进行滤波，同样是将总共K幅对焦测量图的滤波任务分配到N个线程/GPU流处理队列中；做完同步后从所有对焦测量图计算像素索引图。
线程同步4：在该步操作前，利用线程同步3得到的像素索引图对所有景深的L通道、A通道、B通道图像进行融合操作，这里每个通道图像的融合各用一个线程/GPU 流负责处理，各通道完成后进行同步，并完成从LAB图像空间到RGB图像空间的转换，得到最终的融合图，也就是全对焦图像。
可以理解，根据本申请的上述光学显微图像多景深焦点合成方法的具体步骤的不同，可以对上述异构计算方法的计算任务进行适应性调整。
对于计算对焦测量图及后续滤波处理等操作，这里采用CUDA的并行计算模式实现，针对图像中每个像素的具体计算任务（如求某像素的拉普拉斯交叉能量和），都会由GPU线程束（通常以32个为一组）中的一个线程来单独负责执行，而一幅图像可以划分为多个如16×16或32×32大小的图像块，每个图像块对应了CUDA并行架构中的一个包含若干组线程束的Block，而若干个Block又构成了CUDA并行计算架构中的Grid结构，如图8所示。
基于以上异构计算处理方式，在搭载i7处理器和一块英伟达RTX 3060显卡的工作站上能够在1到2秒完成10张2048×2448大小的图像的加载与融合（不包含图像对齐）。
以下举例说明本申请的上述方法。
以法医溺亡检验中的硅藻检测识别应用为例，采用Motic公司开发的相衬光学显微镜（40×物镜，10×目镜）拍摄包含硅藻的尸体组织样品，得到多组多景深局部对焦图像，每张图拍摄时的轴向步进为1um，每组包含10张图。
因为步进相比物距很小，图像间的缩放比例很小可不考虑。采用本申请的上述光学显微图像多景深焦点合成方法以及计算架构进行测试，不同阶段得到的像素索引图如图5所示，为便于展示灰度级拉伸到0~255。可以看到不引入滤波平滑聚合时，结果的噪声水平较大，随着三维加权中值滤波（基于L1距离测度）和高斯卷积导向滤波的加入，聚合效果逐渐明显。从图6所示的融合效果图来看，滤波的加入在不破坏清晰度的情况下能够起到对噪声的抑制，特别是在边缘。
以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。
Claims (10)
1.一种光学显微图像多景深焦点合成方法，其特征在于，包括：
将光学显微图像集合由RGB颜色空间转换成LAB颜色空间，进而获得L通道图像集合，A通道图像集合和B通道图像集合，
计算获得L通道图像集合的每张图像的对焦测量图，形成未滤波的对焦测量图集合，
对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波处理，获得滤波后L通道图像集合、滤波后A通道图像集合、滤波后B通道图像集合，
对未滤波的对焦测量图集合进行级联滤波处理，获得滤波后对焦测量图集合，对滤波后对焦测量图集合进行最大密度投影计算，获得用于融合的像素索引图，
基于像素索引图、滤波后L通道图像集合、滤波后A通道图像集合和滤波后B通道图像集合，获得L通道融合图像、A通道融合图像和B通道融合图像，
基于L通道融合图像、A通道融合图像和B通道融合图像，将LAB颜色空间转换成RGB颜色空间，获得融合后的光学显微彩色图像。
2.根据权利要求1所述的光学显微图像多景深焦点合成方法，其特征在于，所述计算获得L通道图像集合的每张图像的对焦测量图，具体包括：通过改进的拉普拉斯能量交叉和算法来计算获得L通道图像集合的每张图像的对焦测量图，所述改进的拉普拉斯能量交叉和算法应用了以下公式：
其中，L(x, y)是像素(x, y)位置处的对焦测量值，(ε, η)是以像素(x, y)为中心的邻域U内一点，
为像素(ε, η)的水平方向二阶导的绝对值，
为像素(ε, η)的垂直方向二阶导的绝对值，
为像素(ε, η)的正45度方向二阶导的绝对值，
为像素(ε, η)的负45度方向二阶导的绝对值, r是拉普拉斯能量交叉和算子函数的窗口大小。
3.根据权利要求1所述的光学显微图像多景深焦点合成方法，其特征在于，所述对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波处理，具体包括：基于对未滤波的对焦测量图集合进行三维加权中值滤波而获得的每个像素以其为中心的邻域内中值像素pK的位置，对L通道图像集合、A通道图像集合和B通道图像集合分别进行滤波，也就是将所述邻域内中值像素pK的值代替L通道图像集合、A通道图像集合和B通道图像集合上对应的邻域内中心像素的值，
所述三维加权中值滤波的计算公式为：
，该公式表达的意思是用邻域Ω内的中值像素pK的值来替换中心像素的值，中心像素也就是需要通过三维加权中值滤波处理的像素，其中的中值像素pK的确定方法为：
，邻域Ω内的每个像素pk相对中心像素的相似性测度值即为pk的权重，记为
，所述邻域Ω为三维区域，P*是经过滤波后的中心像素值。
4.根据权利要求1所述的光学显微图像多景深焦点合成方法，其特征在于，所述对未滤波的对焦测量图集合进行级联滤波处理，获得滤波后对焦测量图集合，对滤波后对焦测量图集合进行最大密度投影计算，获得用于融合的像素索引图，具体包括：
先对未滤波的对焦测量图集合进行三维加权中值滤波处理，获得一次滤波后对焦测量图集合，对一次滤波后对焦测量图集合进行最大密度投影计算，获得导向滤波引导图，
再基于滤波引导图对一次滤波后对焦测量图集合进行高斯卷积导向滤波处理，进而获得二次滤波后对焦测量图集合，对二次滤波后对焦测量图集合进行最大密度投影计算，获得用于融合的像素索引图，
所述三维加权中值滤波的计算公式为：
，该公式表达的意思是用邻域Ω内的中值像素pK的值来替换中心像素的值，中心像素也就是需要通过三维加权中值滤波处理的像素，其中的中值像素pK的确定方法为：
，邻域Ω内的每个像素pk相对中心像素的相似性测度值即为pk的权重，记为
，所述邻域Ω为三维区域，
所述高斯卷积导向滤波的计算公式为：
其中，Ii和Pi分别为滤波引导图和一次滤波后对焦测量图集合的图像中的邻域
内的像素值，
为高斯卷积核，ϵ是为防止ak过大的正则化参数，ak和bk分别是通过高斯卷积导向滤波的计算公式计算得到的每个像素的线性变换参数。
5.根据权利要求4所述的光学显微图像多景深焦点合成方法，其特征在于，所述高斯卷积导向滤波处理的计算公式还包括：
，其中
和
是以当前像素i为中心，对高斯卷积核的半径R内所有像素j的线性变换参数
和
分别求平均。
6.根据权利要求1所述的光学显微图像多景深焦点合成方法，其特征在于，在将光学显微图像集合由RGB颜色空间转换成LAB颜色空间之前，先对光学显微图像集合进行空间对齐。
7.一种光学显微图像多景深焦点合成系统，其特征在于，包括：
数据获取模块，所述数据获取模块用于获取光学显微图像集合，
数据处理模块，所述数据处理模块用于根据权利要求1至6中任意一项所述的方法来处理光学显微图像集合。
8.一种计算机装置，其特征在于，包括：处理器、存储器、通信接口和通信总线，所述处理器、存储器和通信接口通过所述通信总线完成相互间的通信，所述存储器用于存放至少一个可执行指令，所述可执行指令使所述处理器执行如权利要求1至6中任意一项所述的光学显微图像多景深焦点合成方法对应的操作。
9.一种计算机存储介质，其特征在于，所述计算机存储介质中存储有至少一个可执行指令，所述可执行指令使处理器执行如权利要求1至6中任意一项所述的光学显微图像多景深焦点合成方法对应的操作。
10.一种基于CPU/GPU异构计算的图像处理方法，其特征在于，根据权利要求1至6中任意一项所述的光学显微图像多景深焦点合成方法设定计算任务，通过CPU/GPU异构计算方法处理所述计算任务。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210755701.0A CN114881907B (zh) | 2022-06-30 | 2022-06-30 | 光学显微图像多景深焦点合成方法、系统及图像处理方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210755701.0A CN114881907B (zh) | 2022-06-30 | 2022-06-30 | 光学显微图像多景深焦点合成方法、系统及图像处理方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114881907A CN114881907A (zh) | 2022-08-09 |
| CN114881907B true CN114881907B (zh) | 2022-09-23 |
Family
ID=82682975
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210755701.0A Active CN114881907B (zh) | 2022-06-30 | 2022-06-30 | 光学显微图像多景深焦点合成方法、系统及图像处理方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114881907B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116309216B (zh) * | 2023-02-27 | 2024-01-09 | 南京博视医疗科技有限公司 | 基于多波段的伪彩色图像融合方法及图像融合系统 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103308452A (zh) * | 2013-05-27 | 2013-09-18 | 中国科学院自动化研究所 | 一种基于景深融合的光学投影断层成像图像获取方法 |
| CN104182952A (zh) * | 2014-08-19 | 2014-12-03 | 中国科学院西安光学精密机械研究所 | 多聚焦序列图像融合方法 |
| CN107578418A (zh) * | 2017-09-08 | 2018-01-12 | 华中科技大学 | 一种融合色彩和深度信息的室内场景轮廓检测方法 |
| CN109360235A (zh) * | 2018-09-29 | 2019-02-19 | 中国航空工业集团公司上海航空测控技术研究所 | 一种基于光场数据的混合深度估计方法 |
-
2022
- 2022-06-30 CN CN202210755701.0A patent/CN114881907B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103308452A (zh) * | 2013-05-27 | 2013-09-18 | 中国科学院自动化研究所 | 一种基于景深融合的光学投影断层成像图像获取方法 |
| CN104182952A (zh) * | 2014-08-19 | 2014-12-03 | 中国科学院西安光学精密机械研究所 | 多聚焦序列图像融合方法 |
| CN107578418A (zh) * | 2017-09-08 | 2018-01-12 | 华中科技大学 | 一种融合色彩和深度信息的室内场景轮廓检测方法 |
| CN109360235A (zh) * | 2018-09-29 | 2019-02-19 | 中国航空工业集团公司上海航空测控技术研究所 | 一种基于光场数据的混合深度估计方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114881907A (zh) | 2022-08-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111402146B (zh) | 图像处理方法以及图像处理装置 | |
| Li et al. | Nonnegative mixed-norm preconditioning for microscopy image segmentation | |
| US10168526B2 (en) | Cell contour formation apparatus and method of the same, and non-transitory computer readable storage medium storing a cell contour formation program | |
| CN112733950A (zh) | 一种基于图像融合与目标检测结合的电力设备故障诊断方法 | |
| CN103426148A (zh) | 生成低分辨率输入数据结构的超分辨率版本的方法和设备 | |
| CN104137143A (zh) | 生成低分辨率输入数据结构的超分辨率版本的方法和设备 | |
| JP2013117848A (ja) | 画像処理装置及び画像処理方法 | |
| CN111626936A (zh) | 一种显微图像的快速全景拼接方法及系统 | |
| CN111179170B (zh) | 一种显微血细胞图像快速全景拼接方法 | |
| CN114092325B (zh) | 一种荧光图像超分辨重建方法、装置、计算机设备及介质 | |
| CN114881907B (zh) | 光学显微图像多景深焦点合成方法、系统及图像处理方法 | |
| Ben Hadj et al. | Space variant blind image restoration | |
| JP2015108837A (ja) | 画像処理装置及び画像処理方法 | |
| CN113269672B (zh) | 一种超分辨率的细胞图像构建方法和系统 | |
| CN111062895A (zh) | 一种基于多视场分割的显微图像复原方法 | |
| Zhang et al. | Group-based sparse representation for Fourier ptychography microscopy | |
| Yoo et al. | 3D image reconstruction from multi-focus microscope: axial super-resolution and multiple-frame processing | |
| Salih et al. | Adaptive local exposure based region determination for non-uniform illumination and low contrast images | |
| CN112801913A (zh) | 一种解决显微镜景深限制的方法 | |
| CN111899166A (zh) | 一种基于深度学习的医学高光谱显微图像超分辨重构方法 | |
| CN112163996A (zh) | 一种基于图像处理的平角视频融合方法 | |
| CN115760622A (zh) | 一种拼接显微图像的无监督自适应条纹校正方法 | |
| CN113962904B (zh) | 一种高光谱图像滤波降噪的方法 | |
| CN113191949B (zh) | 多尺度超分辨率病理图像数字化方法、系统及存储介质 | |
| Michelin et al. | Embryo cell membranes reconstruction by tensor voting |
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 | ||
| CP03 | Change of name, title or address |
Address after: Room 601, Building D, Zhonghe (Suzhou) Science and Technology Innovation Port, No. 588 Xiangrong Road, High Speed Rail New City, Xiangcheng District, Suzhou City, Jiangsu Province, 215000 (6th and 7th floors) Patentee after: Jiangsu Jicui sukesi Technology Co.,Ltd. Country or region after: China Address before: 215000 18th floor, Ziguang building (Qidi building), No. 99, nantiancheng Road, Xiangcheng District, Suzhou City, Jiangsu Province Patentee before: Jiangsu Jicui sukesi Technology Co.,Ltd. Country or region before: China |
|
| CP03 | Change of name, title or address |