CN113989246B - 一种基于血液流动特性的透明型血管图像分割方法 - Google Patents

一种基于血液流动特性的透明型血管图像分割方法 Download PDF

Info

Publication number
CN113989246B
CN113989246B CN202111280338.3A CN202111280338A CN113989246B CN 113989246 B CN113989246 B CN 113989246B CN 202111280338 A CN202111280338 A CN 202111280338A CN 113989246 B CN113989246 B CN 113989246B
Authority
CN
China
Prior art keywords
optical flow
pixel point
blood vessel
image
current pixel
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
CN202111280338.3A
Other languages
English (en)
Other versions
CN113989246A (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.)
Nankai University
Original Assignee
Nankai University
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 Nankai University filed Critical Nankai University
Priority to CN202111280338.3A priority Critical patent/CN113989246B/zh
Publication of CN113989246A publication Critical patent/CN113989246A/zh
Application granted granted Critical
Publication of CN113989246B publication Critical patent/CN113989246B/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
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • 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/4038—Image mosaicing, e.g. composing plane images from plane sub-images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2200/00—Indexing scheme for image data processing or generation, in general
    • G06T2200/32—Indexing scheme for image data processing or generation, in general involving image mosaicing
    • 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/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30101—Blood vessel; Artery; Vein; Vascular

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Software Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • General Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Medical Informatics (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Molecular Biology (AREA)
  • Quality & Reliability (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Image Analysis (AREA)

Abstract

本发明涉及一种基于血液流动特性的透明型血管图像分割方法，包括如下步骤：分段拍摄活体模式动物血液流动视频；对标准血管结构二值图以及对应的平均光流场、光流可视化图像进行特征提取，训练机器学习分类模型；计算局部血液流动视频中图像序列平均光流场，并将平均光流场可视化；使用神经网络分割光流可视化图像，累加视频图像序列的背景减除结果，将两步处理结果相加，作为区域生长种子图像；对光流可视化图像、平均光流场、血管结构种子图像进行特征提取，使用训练好的机器学习分类模型进行区域生长运算，得到局部血管结构二值图；进行血管结构图像拼接，得到全局血管结构二值图像。该方法可以简便、快速地获得连通性强、完整性高的血管结构。

Description

一种基于血液流动特性的透明型血管图像分割方法
技术领域
本发明涉及显微图像处理领域，具体涉及一种基于活体模式动物血液流动特性，实现简便、快速获取连通性强、完整性高的透明型血管图像的分割方法。
背景技术
在活体模式动物中准确检测在体血管位置，提供显微视觉反馈，有利于进一步的显微操作和机理研究。血管检测的研究方式主要是视网膜显微图像和荧光血管图像，这类图像边缘明显、纹理相对清晰，主流的血管检测方法对这类图像有较好的效果。然而，部分活体模式动物的血管在显微镜下呈透明状、血管连通性差、边缘模糊、成像信噪比低，现有常规血管检测方法并不适用于这些活体模式动物的血管检测。
目前，人们多使用共聚焦显微镜拍摄模式动物透明型血管的荧光显微图像，以进一步研究血管自动分割方法。然而，荧光显微图像的拍摄过程极其复杂、耗时，图像质量容易出现不连续、模糊等问题；而且，由于共聚焦显微镜的激光强度高，活体模式动物在拍摄过程中很容易死亡，无法进行后续的显微操作。因此，针对显微镜下模式动物难以识别的透明型在体血管，如何在不伤害活体模式动物的基础上，获得边缘清晰、连通性强的完整血管拓扑结构，是血管图像分割领域亟需解决的技术问题。
发明内容
基于现有利用共聚焦显微镜拍摄方式获取模式动物透明型在体血管荧光显微图像技术，存在的上述缺陷和不足，本发明提供了一种基于血液流动特性的透明型血管图像分割方法，该方法可以在不伤害模式动物的基础上，同时免去使用荧光显微镜拍摄血管的繁琐工作，直接通过活体模式生物的血液流动视频，计算血细胞的运动连续性特征，继而简便、快速地获得连通性强、完整性高的血管结构。
本发明解决上述技术问题所采用的技术方案是：一种基于血液流动特性的透明型血管图像分割方法，包括如下步骤：
S1，通过分段显微拍摄的方式，获取具有透明型血管结构的活体模式动物全局血液流动视频；
S2，手动获取活体模式动物局部标准血管结构二值图真值，使用稠密光流算法，计算与所述标准血管结构二值图真值相对应的血液流动视频中任意相邻两帧图像之间的瞬时光流场，对视频图像序列中所有瞬时光流场取平均，获得平均光流场，并将平均光流场可视化，获得光流可视化图像；输入得到的平均光流场，提取每个像素点的光流连续性特征；输入得到的光流可视化图像，提取每个像素点的灰度特征；输入标准血管结构二值图真值，提取每个像素点的血管连续性特征；将得到的每个像素点的光流连续性特征、血管连续性特征以及灰度特征组成特征向量序列，作为训练数据集，训练机器学习分类模型；
S3，使用稠密光流算法，计算其中一段局部血液流动视频中任意相邻两帧图像之间的瞬时光流场，对视频图像序列中所有瞬时光流场取平均，获得平均光流场，并将平均光流场可视化，获得光流可视化图像；
S4，利用DenseUNet神经网络分割步骤S3处理得到的光流可视化图像，得到与光流可视化图像相对应的血管结构二值图初步结果；
S5，采用背景减除法处理步骤S3中的视频图像序列，并将结果累加，获得含部分毛细血管结构的血管结构二值图；
S6，将步骤S4、S5得到的二值图进行累加，作为血管区域生长的种子图像；
S7，利用步骤S2中的特征提取方法，对步骤S3处理得到的光流可视化图像、平均光流场，以及步骤S6处理得到的血管结构种子图像，进行特征提取，计算当前所有候选点的灰度特征、血管连续性特征以及光流连续性特征，组成特征向量序列，使用步骤S2训练好的机器学习分类模型进行区域生长运算，更新血管结构种子图像；其中：区域生长未结束时，重复该步骤直至区域生长结束，得到血管结构二值图像；
S8，采用步骤S3至S7分别处理步骤S1拍摄得到的全局血液流动视频中其它局部血液流动视频，分别得到与该段视频相对应的血管结构二值图；
S9，对步骤S7以及步骤S8处理得到的血管结构二值图，采用基于SIFT特征提取的全局血管结构图像拼接，最后得到活体模式动物的全局血管结构二值图像。
进一步地，步骤S2、步骤S3中，所述稠密光流算法为Fameback光流算法；光流可视化图像通过Munsell颜色系统处理得到。
进一步地，步骤S2中，平均光流场中每个像素点Ii，j的光流连续性特征分别包括直接光流连续性特征(fi)与径向光流连续性特征(fr)；血管结构二值图中，每个像素点Ii，j的血管连续性特征分别包括直接血管连续性特征(ci)与径向血管连续性特征(cr)；其中：
像素点Ii，j直接光流连续性特征(fi)提取方法为：
输入平均光流场，计算当前像素点所在的八邻域内，当前像素点对应的光流矢量和血管像素点对应的光流矢量的相似性，并映射到一个概率值；若在八邻域中，当前像素点的矢量与血管像素点的矢量的相似性Pi(Ii，j)大于阈值Ti，则fi(Ii，j)＝1，否则等于0；具体计算公式如下：
Figure BDA0003328081930000031
Figure BDA0003328081930000032
Figure BDA0003328081930000033
Figure BDA0003328081930000034
Figure BDA0003328081930000035
其中：Tcos表示光流矢量的夹角阈值；Tspeed表示光流矢量的速度阈值；k表示八邻域；Normi表示直接归一化常数；Ti表示直接光流相似性概率阈值；isVessel(Ii，j)判断像素点Ii，j是否是血管像素点，1表示是，0表示不是；isSimilar(a，b)判断a和b是否相似，1表示相似，0表示不相似；Mij表示当前像素点Ii，j所在的八邻域内，与当前像素点Ii，j所对应的光流矢量相似的血管像素点所对应的光流矢量的数量；Pi(Ii，j)表示当前像素点Ii，j所在八邻域中，当前像素点Ii，j的光流矢量与血管像素点的光流矢量的相似性；
像素点Ii，j径向光流连续性特征(fr)提取方法为：
输入平均光流场，计算当前像素点所在的圆形邻域内，当前像素点对应的光流矢量和血管像素点对应的光流矢量的相似性，并映射到一个概率值；若在圆形邻域中，当前像素点的矢量与血管像素点的矢量的相似性Pr(Ii，j)大于阈值Tr，则fr(Ii，j)＝1，否则等于0；具体计算公式如下：
Figure BDA0003328081930000041
Figure BDA0003328081930000042
Figure BDA0003328081930000043
Figure BDA0003328081930000044
Figure BDA0003328081930000045
其中：k表示圆形邻域；Normr表示径向归一化常数；Tr表示径向光流相似性概率阈值；Mij表示当前像素点Ii，j所在的圆形邻域内，与当前像素点Ii，j所对应的光流矢量相似的血管像素点所对应的光流矢量的数量；Pr(Ii，j)表示当前像素点Ii，j所在圆形邻域中，当前像素点Ii，j的光流矢量与血管像素点的光流矢量的相似性；
像素点Ii，j直接血管连续性特征(ci)提取方法为：
输入血管结构二值图真值，统计当前像素点所在的八邻域内血管像素点的数量，并映射到一个概率值；若在八邻域中，当前像素点Ii，j属于血管结构的概率Pi(Ii，j)大于阈值Ti，则ci(Ii，j)＝1，否则等于0；具体计算公式如下：
Figure BDA0003328081930000046
Figure BDA0003328081930000047
Figure BDA0003328081930000048
Figure BDA0003328081930000051
其中：k表示八邻域，Normi表示直接归一化常数，Ti表示直接血管相似性概率阈值；isVessel(Ii，j)判断像素点Ii，j是否是血管像素点，1表示是，0表示不是；Mij表示当前像素点Ii，j所在的八邻域内，血管像素点的数量；Pi(Ii，j)表示当前像素点Ii，j所在八邻域中，当前像素点Ii，j属于血管结构的概率；
像素点Ii，j径向血管连续性特征(cr)的提取方法为：
输入血管结构二值图真值，统计当前像素点所在的圆形邻域内血管像素点的数量，并映射到一个概率值；若在圆形邻域中，当前像素点Ii，j属于血管结构的概率Pr(Ii，j)大于阈值Tr，则cr(Ii，j)＝1，否则等于0，具体计算公式如下：
Figure BDA0003328081930000052
Figure BDA0003328081930000053
Figure BDA0003328081930000054
Figure BDA0003328081930000055
其中上述公式中：k表示圆形邻域，Normr表示径向归一化常数，Tr表示径向血管相似性概率阈值；Mij表示当前像素点Ii，j所在的圆形邻域内，血管像素点的数量；Pr(Ii，j)表示当前像素点Ii，j所在圆形邻域中，当前像素点Ii，j属于血管结构的概率；
光流可视化图像中每个像素点的灰度特征共有53个，分别为：Gaussians filter特征4个，Hessian matrix特征10个，Frangi filter特征3个，Laplacian or sharpenfilter特征3个，Mean max min and median values特征5个，Anisotropic diffusion特征2个，Morphological opening and closing特征6个，Gradient information特征3个，Kuwahara filter特征2个，gabor filter特征4个，high order derivatives特征4个，Entropy特征3个，Rodrigues sobel特征2个，shake effect特征1个，enhancementfiltering特征1个。
进一步地，步骤S4具体为：采用视网膜血管公开数据集DRIVE训练DenseUNet神经网络，使用训练好的DenseUNet神经网络对光流可视化图像进行分割。
进一步地，步骤S5中，使用基于混合高斯模型的背景减除法对视频图像序列进行处理。
进一步地，步骤S7具体为：将步骤S6处理得到的血管结构种子图像中血管像素特定邻域内未访问过的像素点作为候选点；计算当前所有候选点的灰度特征、血管连续性特征、光流连续性特征，组成多个特征向量序列，输入到训练好的机器学习模型中进行分类；当机器学习分类结果为1时，即将该像素点标记为血管像素点，当机器学习分类结果为0时，将该像素点标记为非血管像素点；更新血管结构种子图像，完成一次区域生长；迭代进行上述生长过程，当图像中所有像素点均被访问过时，区域生长结束，得到血管结构二值图像。
进一步地，步骤S9具体为：对于两条相邻视频，读入两张与该视频对应的血管结构二值图像，分别提取其SIFT特征，并将特征点控制在4000个之内；然后基于欧氏距离进行特征匹配，使用FSC算法筛选特征匹配点并计算变换矩阵，使用该变换矩阵完成两张血管结构二值图的图像拼接；依次采用该方法逐步完成全局血管结构图像拼接。
本发明同现有技术相比具有以下优点及效果：
1、本发明所述透明型血管图像分割方法中，通过对分段拍摄的活体模式动物局部血液流动视频分别进行步骤S3至S7的分割计算，可以分别得到与该段视频相对应的血管结构二值图像，再利用基于SIFT特征提取的全局血管结构图像拼接，可以简便、快速得到显微镜下模式动物难以识别的透明型在体血管拓扑结构。
2、本发明所述基于血液流动特性的透明型血管图像分割方法，直接利用活体模式动物血液流动视频，分析血细胞的运动信息，获得血管结构区域，与传统研究方法相比，可以在不伤害模式动物的基础上，同时免去使用荧光显微镜拍摄血管的繁琐工作，是一种简便、快速的血管分割方法。
3、本发明所述基于血液流动特性的透明型血管图像分割方法中，根据血细胞运动的连续性，计算血细胞的运动连续性特征，并在结合特定分割算法的基础上，可以有效提高血管结构分割结果的连通性和完整性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例所述基于血液流动特性的透明型血管图像分割方法流程图。
图2为本发明实施例中通过显微拍摄方式得到的斑马鱼幼鱼血液流动局部视频图像。
图3为本发明实施例步骤S3计算得到的局部血液流动视频中相邻两帧图像的瞬时光流场示意图。
图4为本发明实施例步骤S3中与平均光流场相对应的光流可视化图像。
图5为本发明实施例步骤S2中，机器学习分类模型训练数据集中标准血管结构二值图真值图像示意图。
图6为本发明实施例步骤S4中，利用DenseUNet神经网络对光流可视化图像进行分割的结果示意图。
图7为本发明实施例中步骤S5处理得到的背景减除累加结果示意图。
图8为本发明实施例中步骤S6处理得到的血管区域生长种子图像示意图。
图9为本发明实施例中步骤S7处理得到的，与局部血液流动视频相对应的血管结构二值图像。
图10为本发明实施例中步骤S9处理得到的斑马鱼全局血管结构二值图像。
具体实施方式
下面结合实施例对本发明做进一步的详细说明，以下实施例是对本发明的解释而本发明并不局限于以下实施例。
实施例1：如图1所示，一种基于血液流动特性的透明型血管图像分割方法，该方法中，活体模式动物为具有透明型血管结构的斑马鱼幼鱼，具体包括以下步骤：
S1，通过分段显微拍摄的方式，获取斑马鱼幼鱼全局血液流动视频；具体实验方法如下：使用受精后60小时水浴饲养的斑马鱼幼鱼进行实验，它含有较少的黑色素，可以提供更干净、少遮挡的视野；使用电动倒置显微镜(TiE，尼康)在十倍物镜(NA：0.25)下观察目标；将斑马鱼幼鱼侧卧在培养皿中，通过电动X-Y工作台(ProScan III，Prior)进入显微镜视野，工作台的运动范围为120-80mm，定位分辨率为0.05m；使用高速相机拍摄斑马鱼幼鱼血液流动视频，视频帧率为50fps；通过工作台的移动，共拍摄由7段局部血液流动视频(视频之间具有重叠部分)构成的全局血液流动视频，每段视频时长1～2秒，帧宽度512像素，帧高度320像素；其中一段局部血液流动视频图像如图2所示；
S2，获取活体模式动物局部标准血管结构二值图真值，使用Farneback光流算法，计算与所述标准血管结构二值图真值所对应的血液流动视频中任意相邻两帧图像之间的瞬时光流场，对视频图像序列中所有瞬时光流场取平均，获得平均光流场，并使用Munsell颜色系统将平均光流场可视化，获得光流可视化图像；输入得到的平均光流场，提取每个像素点的光流连续性特征；输入得到的光流可视化图像，提取每个像素点的灰度特征；输入标准血管结构二值图真值，提取每个像素点的血管连续性特征；将得到的每个像素点的光流连续性特征、血管连续性特征以及灰度特征组成特征向量序列，作为训练数据集，训练机器学习分类模型；
具体地，平均光流场中每个像素点Ii，j的光流连续性特征分别包括直接光流连续性特征(fi)与径向光流连续性特征(fr)；血管结构二值图中，每个像素点Ii，j的血管连续性特征分别包括直接血管连续性特征(ci)与径向血管连续性特征(cr)；其中：
(1)像素点Ii，j直接光流连续性特征(fi)提取方法为：
输入平均光流场，计算当前像素点所在的八邻域内，当前像素点对应的光流矢量和血管像素点对应的光流矢量的相似性，并映射到一个概率值；若在八邻域中，当前像素点的矢量与血管像素点的矢量的相似性Pi(Ii，j)大于阈值Ti，则fi(Ii，j)＝1，否则等于0；具体计算公式如下：
Figure BDA0003328081930000081
Figure BDA0003328081930000091
Figure BDA0003328081930000092
Figure BDA0003328081930000093
Figure BDA0003328081930000094
其中：Tcos表示光流矢量的夹角阈值；Tspeed表示光流矢量的速度阈值；k表示八邻域；Normi表示直接归一化常数；Ti表示直接光流相似性概率阈值；isVessel(Ii，j)判断像素点Ii，j是否是血管像素点，1表示是，0表示不是；isSimilar(a，b)判断a和b是否相似，1表示相似，0表示不相似；Mij表示当前像素点Ii，j所在的八邻域内，与当前像素点Ii，j所对应的光流矢量相似的血管像素点所对应的光流矢量的数量；Pi(Ii，j)表示当前像素点Ii，j所在八邻域中，当前像素点Ii，j的光流矢量与血管像素点的光流矢量的相似性；
(2)像素点Ii，j径向光流连续性特征(fr)提取方法为：
输入平均光流场，计算当前像素点所在的圆形邻域内，当前像素点对应的光流矢量和血管像素点对应的光流矢量的相似性，并映射到一个概率值；若在圆形邻域中，当前像素点的矢量与血管像素点的矢量的相似性Pr(Ii，j)大于阈值Tr，则fr(Ii，j)＝1，否则等于0；具体计算公式如下：
Figure BDA0003328081930000095
Figure BDA0003328081930000096
Figure BDA0003328081930000097
Figure BDA0003328081930000098
Figure BDA0003328081930000101
上述公式中：Tcos表示光流矢量的夹角阈值；Tspeed表示光流矢量的速度阈值；k表示圆形邻域；Normr表示径向归一化常数；Tr表示径向光流相似性概率阈值；isVessel(Ii，j)判断像素点Ii，j是否是血管像素点，1表示是，0表示不是；isSimilar(a，b)判断a和b是否相似，1表示相似，0表示不相似；Mij表示当前像素点Ii，j所在的圆形邻域内，与当前像素点Ii，j所对应的光流矢量相似的血管像素点所对应的光流矢量的数量；Pr(Ii，j)表示当前像素点Ii，j所在圆形邻域中，当前像素点Ii，j的光流矢量与血管像素点的光流矢量的相似性；
(3)像素点Ii，j直接血管连续性特征(ci)提取方法为：
输入血管结构二值图真值，统计当前像素点所在的八邻域内血管像素点的数量，并映射到一个概率值；若在八邻域中，当前像素点Ii，j属于血管结构的概率Pi(Ii，j)大于阈值Ti，则ci(Ii，j)＝1，否则等于0；具体计算公式如下：
Figure BDA0003328081930000102
Figure BDA0003328081930000103
Figure BDA0003328081930000104
Figure BDA0003328081930000105
其中：k表示八邻域，Normi表示直接归一化常数，Ti表示直接血管相似性概率阈值；isVessel(Ii，j)判断像素点Ii，j是否是血管像素点，1表示是，0表示不是；Mij表示当前像素点Ii，j所在的八邻域内，血管像素点的数量；Pi(Ii，j)表示当前像素点Ii，j所在八邻域中，当前像素点Ii，j属于血管结构的概率；
(4)像素点Ii，j径向血管连续性特征(cr)的提取方法为：
输入血管结构二值图真值，统计当前像素点所在的圆形邻域内血管像素点的数量，并映射到一个概率值；若在圆形邻域中，当前像素点Ii，j属于血管结构的概率Pr(Ii，j)大于阈值Tr，则cr(Ii，j)＝1，否则等于0，具体计算公式如下。
Figure BDA0003328081930000111
Figure BDA0003328081930000112
Figure BDA0003328081930000113
Figure BDA0003328081930000114
其中：k表示圆形邻域，Normr表示径向归一化常数，Tr表示径向血管相似性概率阈值；isVessel(Ii，j)判断像素点Ii，j是否是血管像素点，1表示是，0表示不是；Mij表示当前像素点Ii，j所在的圆形邻域内，血管像素点的数量；Pr(Ii，j)表示当前像素点Ii，j所在圆形邻域中，当前像素点Ii，j属于血管结构的概率；
(5)光流可视化图像中，每个像素点Ii，j的灰度特征包括：Gaussians filter特征、Hessian matrix特征、Frangi filter特征、Laplacian or sharpen filter特征、Mean maxmin and median values特征、Anisotropic diffusion特征、Morphological opening andclosing特征、Gradient information特征、Kuwahara filter特征、Gabor filter特征、High order derivatives特征、Entropy特征、Rodrigues sobel特征、Shake effect特征、Enhancement filtering特征；其中：
Gaussians filter特征：一种线性平滑滤波，去除高斯噪声，降低图像灰度的尖锐变化，使图像变得模糊；使用不同的参数组合，共生成4个特征值；
Hessian matrix特征：一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率，在图像处理中通常用于边缘检测；使用不同的参数组合，共生成10个特征值；
Frangi filter特征：一种基于Hessian matrix的滤波器，用于边缘检测，并实现血管图像增强；使用不同的参数组合，共生成3个特征值；
Laplacian or sharpen filter特征：一种二阶空间导数的各向同性微分算子，用于增强图像中灰度突变的区域，减弱灰度的缓慢变化区域；使用不同的参数组合，共生成3个特征值；
Mean max min and median values特征：以像素为中心，在其7像素×7像素的正方形窗口中计算像素灰度的算术平均值、几何平均值、最大值、最小值和中值，共生成5个特征值；
Anisotropic diffusion特征：一种各向异性扩散滤波器，用于平滑图像，保持图像细节特征，并减少噪声；使用不同的参数组合，共生成2个特征值；
Morphological opening and closing特征：形态学中的开运算和闭运算；使用不同的参数组合，共生成6个特征值；
Gradient information特征：利用线性核和高斯核进行梯度信息估计，计算图像的梯度特征；使用不同的参数组合，共生成3个特征值；
Kuwahara filter特征：一种边缘保留滤波器，能在平滑图像的同时保留图像有意义的边缘信息；使用不同的参数组合，共生成2个特征值；
Gabor filter特征：一种用于边缘提取的线性滤波器，其频率和方向表达与人类视觉系统类似，能够提供良好的方向选择和尺度选择特性，对于光照变化不敏感；使用不同的参数组合，共生成4个特征值；
High order derivatives特征：计算图像的高阶导数；使用不同的参数组合，共生成4个特征值；
Entropy特征：一种特征的统计形式，反映图像灰度分布的聚集特征和空间特征；使用不同的参数组合，共生成3个特征值；
Rodrigues sobel特征：一种离散微分边缘检测算子，用来计算图像灰度的近似梯度；使用不同的参数组合，共生成2个特征值；
Shake effect特征：利用高斯模糊增强图像，减少图像噪声；使用不同的参数组合，共生成1个特征值；
Enhancement filtering特征：一种多尺度血管图像增强滤波，实现在不同尺度下的图像增强；使用不同的参数组合，共生成1个特征值；
对一幅长m、宽n的图像，提取特征后，共得到m*n个57×1的特征向量序列，(本实施例中共得到512*320个57×1的特征向量序列)将其作为训练数据集；将血管结构二值图真值中血管的像素值标注为1，非血管的像素值标注为0；使用随机森林算法作为机器学习模型；其中，图5为机器学习模型训练数据集中标准血管结构二值图真值示意图，由于目前还没有关于斑马鱼血管图像的数据和标准参考图像，本实施例步骤S2中所述标准血管结构二值图真值为实验室采用手动分割方式得到的标准参考图像，具体手动分割的方法属于现有技术部分，本发明不再赘述；
S3，使用稠密光流算法，计算如图2所示的局部血液流动视频中任意相邻两帧图像之间的瞬时光流场；对视频图像序列中所有瞬时光流场取平均，获得平均光流场，并将平均光流场可视化，获得光流可视化图像；
其中：所述稠密光流算法采用Farneback光流算法，计算得到的瞬时光流场如图3所示，(为直观清晰地描述光流场，图3仅画出部分像素点的光流矢量)其中箭头方向表示速度方向，箭头长度表示速度大小；进一步地，使用Munsell颜色系统将平均光流场可视化，可视化图像中，不同的颜色表示不同的运动方向，颜色的深浅表示运动速度的快慢，光流可视化图像如图4所示；
S4，利用DenseUNet神经网络分割步骤S3处理得到的光流可视化图像，得到与光流可视化图像相对应的血管结构二值图初步结果，如图6所示；具体为：使用视网膜血管图像数据集DRIVE训练DenseUNet神经网络，然后使用训练好的DenseUNet网络分割光流可视化图像；
S5，采用基于混合高斯模型的背景减除法处理步骤S3中的视频图像序列，累加每两帧图像的背景减除结果，获得含部分毛细血管结构的血管结构二值图，图7(a)是两帧图像的背景减除结果，图7(b)是视频图像序列背景减除的累加结果；
S6，将步骤S4、S5得到的二值图进行累加，作为血管区域生长的种子图像；其中：将步骤S4、S5得到的二值图进行累加后，需要去除面积小于20像素点的小连通区域，一定程度上去除噪声干扰；
S7，利用步骤S2中的特征提取方法，对步骤S3处理得到的光流可视化图像、平均光流场，以及步骤S6处理得到的血管结构种子图像，进行特征提取，计算当前所有候选点的灰度特征、血管连续性特征以及光流连续性特征，组成特征向量序列，使用步骤S2训练好的机器学习分类模型进行区域生长运算，更新血管结构种子图像；其中：区域生长未结束时，重复该步骤直至区域生长结束，得到血管结构二值图像；
具体为：输入步骤S6处理得到的血管结构种子图像，每次迭代时，将血管结构种子图像中血管像素特定邻域内未访问过的像素点作为候选点，并提取每个像素点的血管连续性特征；输入步骤S3得到的平均光流场，提取每个候选点的光流连续性特征；输入步骤S3中得到的光流可视化图像，提取每个候选点的灰度特征；将当前所有候选点的灰度特征、血管连续性特征、光流连续性特征，组成多个57×1的特征向量序列，输入到训练好的随机森林中进行分类；当特征向量的随机森林分类结果为1时，即将该像素点标记为血管像素点，当特征向量的随机森林分类结果为0时，将该像素点标记为非血管像素点；更新血管结构种子图像，完成一次区域生长；迭代进行上述生长过程，当血管结构种子图像中所有像素点均被访问过时，区域生长结束，得到血管结构二值图像；
S8，采用步骤S3至S7分别处理步骤S1拍摄得到的全局血液流动视频中其它局部血液流动视频，分别得到与该段视频相对应的血管结构二值图；
S9，对步骤S7以及步骤S8处理得到的血管结构二值图，采用基于SIFT特征提取的全局血管结构图像拼接，最后得到活体模式动物的全局血管结构二值图像；
具体为：对于两条相邻视频，读入两张与该视频对应的血管结构二值图像，分别提取其SIFT特征，并将特征点控制在4000个之内；然后基于欧氏距离进行特征匹配，使用FSC算法筛选特征匹配点并计算变换矩阵，使用该变换矩阵完成两张血管结构二值图的图像拼接；依次采用该方法逐步完成全局血管结构图像拼接。图10为斑马鱼幼鱼的全局血管结构二值图像。
此外，需要说明的是，本说明书中所描述的具体实施例，其零、部件的形状、所取名称等可以不同。凡依本发明专利构思所述的构造、特征及原理所做的等效或简单变化，均包括于本发明专利的保护范围内。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，只要不偏离本发明的结构或者超越本权利要求书所定义的范围，均应属于本发明的保护范围。

Claims (7)

1.一种基于血液流动特性的透明型血管图像分割方法，其特征在于，包括如下步骤：
S1，通过分段显微拍摄的方式，获取具有透明型血管结构的活体模式动物全局血液流动视频；
S2，获取活体模式动物局部标准血管结构二值图真值，使用稠密光流算法，计算与所述标准血管结构二值图真值相对应的血液流动视频中任意相邻两帧图像之间的瞬时光流场，对视频图像序列中所有瞬时光流场取平均，获得平均光流场，并将平均光流场可视化，获得光流可视化图像；输入得到的平均光流场，提取每个像素点的光流连续性特征；输入得到的光流可视化图像，提取每个像素点的灰度特征；输入标准血管结构二值图真值，提取每个像素点的血管连续性特征；将得到的每个像素点的光流连续性特征、血管连续性特征以及灰度特征组成特征向量序列，作为训练数据集，训练机器学习分类模型；
S3，使用稠密光流算法，计算步骤S1其中一段局部血液流动视频中任意相邻两帧图像之间的瞬时光流场，对视频图像序列中所有瞬时光流场取平均，获得平均光流场，并将平均光流场可视化，获得光流可视化图像；
S4，利用DenseUNet神经网络分割步骤S3处理得到的光流可视化图像，得到与光流可视化图像相对应的血管结构二值图初步结果；
S5，采用背景减除法处理步骤S3中的视频图像序列，并将结果累加，获得含部分毛细血管结构的血管结构二值图；
S6，将步骤S4、S5得到的二值图进行累加，作为血管区域生长的种子图像；
S7，利用步骤S2中的特征提取方法，对步骤S3处理得到的光流可视化图像、平均光流场以及步骤S6处理得到的血管结构种子图像，进行特征提取，计算当前所有候选点的灰度特征、血管连续性特征以及光流连续性特征，组成特征向量序列，使用步骤S2训练好的机器学习分类模型进行区域生长运算，更新血管结构种子图像；其中：区域生长未结束时，重复该步骤直至区域生长结束，得到血管结构二值图像；
S8，采用步骤S3至S7分别处理步骤S1拍摄得到的全局血液流动视频中其它局部血液流动视频，分别得到与该段视频相对应的血管结构二值图；
S9，对步骤S7以及步骤S8处理得到的血管结构二值图，采用基于SIFT特征提取的全局血管结构图像拼接，最后得到活体模式动物的全局血管结构二值图像；步骤S9具体为：对于两条相邻视频，读入两张与该视频对应的血管结构二值图像，分别提取其SIFT特征，并将特征点控制在4000个之内；然后基于欧氏距离进行特征匹配，使用FSC算法筛选特征匹配点并计算变换矩阵，使用该变换矩阵完成两张血管结构二值图的图像拼接；依次采用该方法逐步完成全局血管结构图像拼接。
2.根据权利要求1所述的一种基于血液流动特性的透明型血管图像分割方法，其特征在于，步骤S2、步骤S3中，所述稠密光流算法为Farneback光流算法；光流可视化图像通过Munsell颜色系统处理得到。
3.根据权利要求1所述的基于血液流动特性的透明型血管图像分割方法，其特征在于，步骤S2中，平均光流场中每个像素点I i,j 的光流连续性特征分别包括直接光流连续性特征(
Figure DEST_PATH_IMAGE001
)与径向光流连续性特征 (
Figure DEST_PATH_IMAGE002
)；血管结构二值图中，每个像素点I i,j 的血管连续性特征分别包括直接血管连续性特征(
Figure DEST_PATH_IMAGE003
)与径向血管连续性特征(
Figure DEST_PATH_IMAGE004
)；其中：
像素点I i,j 直接光流连续性特征(
Figure 460110DEST_PATH_IMAGE001
)提取方法为：
输入平均光流场，计算当前像素点所在的八邻域内，当前像素点对应的光流矢量和血管像素点对应的光流矢量的相似性，并映射到一个概率值；若在八邻域中，当前像素点的矢量与血管像素点的矢量的相似性
Figure DEST_PATH_IMAGE005
大于阈值
Figure DEST_PATH_IMAGE006
，则
Figure DEST_PATH_IMAGE007
，否则等于0；具体计算公式如下：
Figure DEST_PATH_IMAGE008
；
Figure DEST_PATH_IMAGE009
；
Figure DEST_PATH_IMAGE010
；
Figure DEST_PATH_IMAGE011
；
Figure DEST_PATH_IMAGE012
；
其中：T cos 表示光流矢量的夹角阈值；T speed 表示光流矢量的速度阈值；k表示八邻域；Norm i 表示直接归一化常数；T i 表示直接光流相似性概率阈值；isVessel(I i,j )判断像素点I i,j 是否是血管像素点，1表示是，0表示不是；isSimilar(a,b)判断a和b是否相似，1表示相似，0表示不相似；M ij 表示当前像素点I i,j 所在的八邻域内，与当前像素点I i,j 所对应的光流矢量相似的血管像素点所对应的光流矢量的数量；P i (I i,j )表示当前像素点I i,j 所在八邻域中，当前像素点I i,j 的光流矢量与血管像素点的光流矢量的相似性；
像素点I i,j 径向光流连续性特征 (
Figure DEST_PATH_IMAGE013
)提取方法为：
输入平均光流场，计算当前像素点所在的圆形邻域内，当前像素点对应的光流矢量和血管像素点对应的光流矢量的相似性，并映射到一个概率值；若在圆形邻域中，当前像素点的矢量与血管像素点的矢量的相似性
Figure DEST_PATH_IMAGE014
大于阈值
Figure DEST_PATH_IMAGE015
，则
Figure DEST_PATH_IMAGE016
，否则等于0；具体计算公式如下：
Figure DEST_PATH_IMAGE017
；
Figure DEST_PATH_IMAGE018
；
Figure DEST_PATH_IMAGE019
；
Figure DEST_PATH_IMAGE020
；
Figure DEST_PATH_IMAGE021
；
其中： k表示圆形邻域；Norm r 表示径向归一化常数；T r 表示径向光流相似性概率阈值；M ij 表示当前像素点I i,j 所在的圆形邻域内，与当前像素点I i,j 所对应的光流矢量相似的血管像素点所对应的光流矢量的数量；P r (I i,j )表示当前像素点I i,j 所在圆形邻域中，当前像素点I i,j 的光流矢量与血管像素点的光流矢量的相似性；
像素点I i,j 直接血管连续性特征(
Figure DEST_PATH_IMAGE022
)提取方法为：
输入血管结构二值图真值，统计当前像素点所在的八邻域内血管像素点的数量，并映射到一个概率值；若在八邻域中，当前像素点
Figure DEST_PATH_IMAGE023
属于血管结构的概率
Figure DEST_PATH_IMAGE024
大于阈值
Figure DEST_PATH_IMAGE025
，则
Figure DEST_PATH_IMAGE026
，否则等于0；具体计算公式如下：
Figure DEST_PATH_IMAGE027
；
Figure DEST_PATH_IMAGE028
；
Figure DEST_PATH_IMAGE029
；
Figure DEST_PATH_IMAGE030
；
其中：k表示八邻域，Norm i 表示直接归一化常数，T i 表示直接血管相似性概率阈值；isVessel(I i,j )判断像素点I i,j 是否是血管像素点，1表示是，0表示不是；M ij 表示当前像素点I i,j 所在的八邻域内，血管像素点的数量；P i (I i,j )表示当前像素点I i,j 所在八邻域中，当前像素点
Figure DEST_PATH_IMAGE031
属于血管结构的概率；
像素点I i,j 径向血管连续性特征(
Figure DEST_PATH_IMAGE032
)的提取方法为：
输入血管结构二值图真值，统计当前像素点所在的圆形邻域内血管像素点的数量，并映射到一个概率值；若在圆形邻域中，当前像素点
Figure DEST_PATH_IMAGE033
属于血管结构的概率
Figure DEST_PATH_IMAGE034
大于阈值
Figure DEST_PATH_IMAGE035
，则
Figure DEST_PATH_IMAGE036
，否则等于0，具体计算公式如下：
Figure DEST_PATH_IMAGE037
；
Figure DEST_PATH_IMAGE038
；
Figure DEST_PATH_IMAGE039
；
Figure DEST_PATH_IMAGE040
；
其中上述公式中：k表示圆形邻域，Norm r 表示径向归一化常数，T r 表示径向血管相似性概率阈值；M ij 表示当前像素点I i,j 所在的圆形邻域内，血管像素点的数量；P r (I i,j )表示当前像素点I i,j 所在圆形邻域中，当前像素点
Figure 384773DEST_PATH_IMAGE033
属于血管结构的概率。
4.根据权利要求1所述的一种基于血液流动特性的透明型血管图像分割方法，其特征在于，步骤S2中，光流可视化图像中每个像素点的灰度特征共有53个，分别为：Gaussiansfilter特征4个，Hessian matrix特征 10个，Frangi filter 特征3个，Laplacian orsharpen filter 特征3个，Mean max min and median values 特征5个，Anisotropicdiffusion 特征2个，Morphological opening and closing 特征6个，Gradientinformation 特征3个，Kuwahara filter特征 2个，gabor filter 特征4个，high orderderivatives 特征4个，Entropy 特征3个，Rodrigues sobel 特征2个，shake effect特征1个，enhancement filtering特征1个。
5.根据权利要求1所述的一种基于血液流动特性的透明型血管图像分割方法，其特征在于，步骤S4具体为：采用视网膜血管公开数据集DRIVE训练DenseUNet神经网络，使用训练好的DenseUNet神经网络对光流可视化图像进行分割。
6.根据权利要求1所述的一种基于血液流动特性的透明型血管图像分割方法，其特征在于，步骤S5中，使用基于混合高斯模型的背景减除法对视频图像序列进行处理。
7.根据权利要求1所述的一种基于血液流动特性的透明型血管图像分割方法，其特征在于，步骤S7具体为：将步骤S6处理得到的血管结构种子图像中血管像素特定邻域内未访问过的像素点作为候选点；计算当前所有候选点的灰度特征、血管连续性特征、光流连续性特征，组成多个特征向量序列，输入到训练好的机器学习模型中进行分类；当机器学习分类结果为1时，即将该像素点标记为血管像素点，当机器学习分类结果为0时，将该像素点标记为非血管像素点；更新血管结构种子图像，完成一次区域生长；迭代进行上述生长过程，当图像中所有像素点均被访问过时，区域生长结束，得到血管结构二值图像。
CN202111280338.3A 2021-10-29 2021-10-29 一种基于血液流动特性的透明型血管图像分割方法 Active CN113989246B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111280338.3A CN113989246B (zh) 2021-10-29 2021-10-29 一种基于血液流动特性的透明型血管图像分割方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111280338.3A CN113989246B (zh) 2021-10-29 2021-10-29 一种基于血液流动特性的透明型血管图像分割方法

Publications (2)

Publication Number Publication Date
CN113989246A CN113989246A (zh) 2022-01-28
CN113989246B true CN113989246B (zh) 2023-01-24

Family

ID=79745161

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111280338.3A Active CN113989246B (zh) 2021-10-29 2021-10-29 一种基于血液流动特性的透明型血管图像分割方法

Country Status (1)

Country Link
CN (1) CN113989246B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106407917A (zh) * 2016-09-05 2017-02-15 山东大学 基于动态尺度分配的视网膜血管提取方法及系统
CN106408562A (zh) * 2016-09-22 2017-02-15 华南理工大学 基于深度学习的眼底图像视网膜血管分割方法及系统
CN106530283A (zh) * 2016-10-20 2017-03-22 北京工业大学 一种基于svm的医疗图像血管识别方法
CN107248161A (zh) * 2017-05-11 2017-10-13 江西理工大学 一种多特征融合的有监督视网膜血管提取方法
CN109003279A (zh) * 2018-07-06 2018-12-14 东北大学 一种基于K-Means聚类标注和朴素贝叶斯模型的眼底视网膜血管分割方法及系统

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106407917A (zh) * 2016-09-05 2017-02-15 山东大学 基于动态尺度分配的视网膜血管提取方法及系统
CN106408562A (zh) * 2016-09-22 2017-02-15 华南理工大学 基于深度学习的眼底图像视网膜血管分割方法及系统
CN106530283A (zh) * 2016-10-20 2017-03-22 北京工业大学 一种基于svm的医疗图像血管识别方法
CN107248161A (zh) * 2017-05-11 2017-10-13 江西理工大学 一种多特征融合的有监督视网膜血管提取方法
CN109003279A (zh) * 2018-07-06 2018-12-14 东北大学 一种基于K-Means聚类标注和朴素贝叶斯模型的眼底视网膜血管分割方法及系统

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
A Discriminatively Trained Fully Connected Conditional Random Field Model for Blood Vessel Segmentation in Fundus Images;Orlando Jose Ignacio;《IEEE Transactions on Biomedical Engineering》;20171231;第34卷(第1期);全文 *
An improved U-Net method with High-resolution Feature Maps for Retinal Blood Vessel Segmentation;Dongye Chang-lei;《Journal of Physics: Conference Series》;20210131;第1848卷(第1期);全文 *
MFCSNet: Multi-Scale Deep Features Fusion and Cost-Sensitive Loss Function Based Segmentation Network for Remote Sensing Images;Ende Wang;《Applied Sciences》;20191231;第9卷(第19期);全文 *
基于改进U-Net视网膜血管图像分割算法;李大湘;《光学学报》;20200531;第40卷(第10期);全文 *
基于融合相位特征的视网膜血管分割算法;李媛媛;《计算机应用》;20180710;第38卷(第7期);全文 *
多特征融合神经网络的眼底血管分割算法;宋姝洁;《计算机科学与探索》;20210723;第12卷(第15期);全文 *

Also Published As

Publication number Publication date
CN113989246A (zh) 2022-01-28

Similar Documents

Publication Publication Date Title
Bhalla et al. A fuzzy convolutional neural network for enhancing multi-focus image fusion
CN108596849B (zh) 一种基于天空区域分割的单幅图像去雾方法
Mukherjee et al. Tubularity flow field—a technique for automatic neuron segmentation
EP3271894B1 (en) System and method for image segmentation
Jurrus et al. Axon tracking in serial block-face scanning electron microscopy
CN110736747B (zh) 一种细胞液基涂片镜下定位的方法及系统
US11538261B2 (en) Systems and methods for automated cell segmentation and labeling in immunofluorescence microscopy
Plissiti et al. Automated segmentation of cell nuclei in PAP smear images
Shen et al. Segmenting multiple overlapping nuclei in H&E stained breast cancer histopathology images based on an improved watershed
CN110648312A (zh) 一种基于鳞片形态特征分析的羊毛与羊绒纤维识别方法
CN113989246B (zh) 一种基于血液流动特性的透明型血管图像分割方法
CN117058018A (zh) 面向水下结构检测的悬浮杂质视觉遮挡区域修复方法
CN113989170B (zh) 一种基于血液流动特性的透明型血管类型识别方法
Abdelsamea An enhancement neighborhood connected segmentation for 2D-cellular image
US11854281B2 (en) System, method, and computer-accessible medium for processing brain images and extracting neuronal structures
Chen et al. Interactive image segmentation based on object contour feature image
Lelandais et al. Deconvolution regularized using fuzzy c-means algorithm for biomedical image deblurring and segmentation
Quelhas et al. Optical flow based arabidopsis thaliana root meristem cell division detection
Marcuzzo et al. Automatic cell segmentation from confocal microscopy images of the Arabidopsis root
Meng et al. An extraction method of fish vector field for marine cage culture based on UAV vision
Yu et al. Automatic cell object extraction of red tide algae in microscopic images
Marcuzzo et al. A hybrid approach for Arabidopsis root cell image segmentation
Krämer et al. Comparison of segmentation algorithms for the zebrafish heart in fluorescent microscopy images
Jadhav et al. SEGMENTATION OF BRIGHT REGION OF THE OPTIC DISC FOR EYE DISEASE PREDICTION.
Chen et al. Glottis segmentation using dynamic programming

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