CN113012198B - 颅内外血管吻合术中的血流自动定量分析方法 - Google Patents
颅内外血管吻合术中的血流自动定量分析方法 Download PDFInfo
- Publication number
- CN113012198B CN113012198B CN202110303224.XA CN202110303224A CN113012198B CN 113012198 B CN113012198 B CN 113012198B CN 202110303224 A CN202110303224 A CN 202110303224A CN 113012198 B CN113012198 B CN 113012198B
- Authority
- CN
- China
- Prior art keywords
- image
- layer
- gaussian
- images
- blood
- 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
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
- G06T7/246—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/22—Matching criteria, e.g. proximity measures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
- G06T7/269—Analysis of motion using gradient-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
- G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/20—Image preprocessing
- G06V10/25—Determination of region of interest [ROI] or a volume of interest [VOI]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/20—Image preprocessing
- G06V10/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
- G06V10/267—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion by performing operations on regions, e.g. growing, shrinking or watersheds
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/40—Extraction of image or video features
- G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
- G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]
-
- 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/10016—Video; Image sequence
-
- 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/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
-
- 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/20081—Training; Learning
-
- 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/20084—Artificial neural networks [ANN]
-
- 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
- G06T2207/30104—Vascular flow; Blood flow; Perfusion
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Multimedia (AREA)
- Data Mining & Analysis (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Artificial Intelligence (AREA)
- General Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Image Analysis (AREA)
Abstract
本发明提供一种颅内外血管吻合术中血流自动定量分析方法,包括步骤:S1、从若干患者术前术后ICG‑VA视频中各提取一帧视频图像;S2、训练多任务Unet网络,用于分割一帧ICG‑VA视频图像中的所有血管以及受体血管;S3、基于训练好的网络分割待检测患者术前术后一帧视频图像的所有血管以及受体血管;S4、基于sift算法,将患者术前术后的视频图像进行配准;S5、基于HS光流法分析患者术前术后受体血管的血流方向。本发明不仅能对受体血管进行血液灌注分析,还可以自动检测和匹配手术前后图像的血管,测量任意血管的血流方向和强度曲线,并给出校正后的手术前后的血流灌注变化图,使术者能够直观地评估吻合术引起的血流灌注的变化。
Description
技术领域
本发明涉及计算机辅助诊断技术领域,具体涉及一种颅内外血管吻合术中的血流自动定量分析方法。
背景技术
颞浅动脉-大脑中动脉吻合术(STA-MCA)自1986年由Spetzler和Martin首次报道,已广泛用于临床出血或缺血性脑血管病及复杂颅内动脉瘤的血流重建。Raabe等人首次证明了吲哚菁绿(ICG)动态流过吻合血管,可辅助判定术中血流的通畅性和血流方向。大量研究证明术中吲哚菁绿血管造影(ICG-VA)和已有商业软件可用于生成延时彩色图以获取局部脑血流量(CBF),并可评估术后吻合口周围血流灌注情况。
但是作为图像后处理的该商业软件,无法实现CBF的连续实时和动态分析,不能自动检测手术前后血管的血流方向。
由于ICG-VA图像中存在许多直径差异大且分辨率低的血管,现有的商业软件对血管识别率低,识别精度不高。
另一方面,手术前后图像难免会产生差异,造成这些差异的主要原因有:(1)显微镜镜头位置的变化(显微镜镜头沿X轴,Y轴或Z轴移动);(2)相机内部参数的变化,如焦距和分辨率;(3)吻合术引起的血流灌注的改变。这些差异对于从视觉上直观比较搭桥手术前后的血流灌注情况产生了一定影响。
因此,亟需开发一种有助于临床研究的自动血流定量分析方法,该方法能够准确提取ICG-VA视频图像中的血管信息,且能够通过配准手术前后ICG-VA视频中的图像来直观地比较血流灌注差异,并基于手术前后ICG-VA视频进行血流向的实时测量,以协助术者根据手术效果随时调整后续治疗方案。
发明内容
本发明的目的是提供一种颅内外血管吻合术中的血流自动定量分析方法,能够自动检测和匹配患者手术前后ICG-VA视频图像中的血管,分析患者术前术后图像中受体血管的血流方向,以辅助医师分析手术前后血流灌注的变化情况,提高手术的成功率。
为了达到上述目的,本发明提供一种颅内外血管吻合术中的血流自动定量分析方法,包含步骤:
S1、采集若干患者术前术后的ICG-VA视频,从每个ICG-VA视频提取一帧视频图像;手工标注所提取视频图像中的所有血管及受体血管;基于标注的视频图像建立训练集和测试集;
S2、建立多任务Unet网络,通过所述训练集和测试集训练该多任务Unet网络;所述多任务Unet网络用于分割一帧ICG-VA视频图像中所有血管以及受体血管;
S3、从待检测患者的术前/术后ICG-VA视频中采集一帧术前/术后待测视频图像,通过训练好的多任务Unet网络,分割该术前/术后待测视频图像的所有血管得到对应的第一图像/第二图像,并分割该术前/术后待测视频图像中的受体血管得到对应的第三图像/第四图像;
S4、基于SIFT图像配准算法,将待检测患者的第一图像向第二图像进行配准;
S5、通过HS光流法的血流分析算法追踪待检测患者第三图像/第四图像中的术前/术后受体血管的血流方向。
优选的,所述多任务Unet网络包含一个编码器路径和结构相同的两个解码器路径;编码器路径的输出作为两个解码器路径的输入;其中一个解码器路径用于从输入的视频图像中分割所有血管,另一个解码器路径用于从输入的视频图像中分割受体血管;
所述编码器路径包括依序连接的四个编码器层,分别为第一至第四编码器层;所述编码器层包含依序连接的卷积层、非线性激活层和最大池化层;
所述解码器路径包含依序连接的四个解码器层,输出卷积层;所述四个解码器层分别为第一至第四解码器层;其中第i编码器层与第5-i解码器层的特征级联,i∈[1,4]。
优选的,多任务UNet网络的优化器为Adam函数;输出卷积层的激活函数采用softmax函数;多任务Unet网络的损失函数为:
Ltotal=αLall+βLreceip (1)
其中:Ltotal表示总损失,α和β是比例因子;
Lall=∑x∈Ω logpall(x;lall(x)),Lreceip=∑x∈Ω logpreceip(x;lreceip(x));
Lall和Lreceip表示分割所有血管和分割受体血管的分类错误函数;x是图像空间Ω中的像素位置;pall表示经过softmax激活函数后,对于真标签lall的预测概率;preceip表示经过softmax激活函数后,对于真实标签lreceip的预测概率。
优选的,步骤S4包含:
S41、通过高斯卷积,分别为待检测患者的第一、第二图像生成对应的第一高斯金字塔;
从第一高斯金字塔的底部至顶部,该第一高斯金字塔具有第一至第N组高斯图像;且从第一高斯金字塔的底部至顶部,每组高斯图像具有第一层至第L层高斯图像;第一图像/第二图像作为其对应的第一高斯金字塔第一组第一层图像;
第一高斯金字塔同组的第i+1层高斯图像由同组的第i层高斯图像经高斯卷积之后生成;
第一高斯金字塔中,第j+1组的第一层图像是由第j组的第4层图像经采样得到。
S42、基于所述第一高斯金字塔生成对应的差分金字塔;
从差分金字塔的底部至顶部,差分金字塔具有第一至第N组差分图像;差分金字塔的第O组第l层差分图像是有第一高斯金字塔的第O组第l+1层高斯图像减第O组第l层高斯图像得到;其中O∈[1,N],l∈[1,L-1]。
S43、检测差分图像的局部极大值或极小值,获取差分图像的关键点;
S44、为每个关键点生成多维度的特征描述符;
S45、使用最近邻搜索,将第一图像的关键点与第二图像的关键点进行匹配;
S46、得到匹配上的特征点对后,利用特征点对计算第一图像到第二图像的单应性矩阵;将第一图像乘以该单应性矩阵,得到调整后的第一图像。
优选的,步骤S41中由同组的第i层高斯图像通过高斯卷积生成同组第i+1层高斯图像的方法为:
令J(x,y)表示待高斯卷积的图像,L(x,y,σ′)为J(x,y)经过高斯卷积生成的图像,其中
L(x,y,σ′)=G(x,y,σ′)*J(x,y) (2);
其中(x,y)为坐标,σ′为J(x,y)的平滑度,*为卷积操作,G(x,y,σ′)是尺度变化的高斯函数,
优选的,步骤S5包含:
S51、通过高斯卷积,分别为待检测患者的第三、第四图像生成对应的第二高斯金字塔;从所述第二高斯金字塔的底部至顶部,第二高斯金字塔包含第一至第L层图像;
第三图像/第四图像作为其对应的第二高斯金字塔的第一层图像;
第二高斯金字塔的第i+1层图像由其第i层图像经高斯卷积之后生成;i∈[1,L-1];
S52、通过HS光流法计算所述第二高斯金字塔中第一层图像像素的光流场(u,v);
S53、计算第二高斯金字塔中第L层图像像素的光流场(uL,vL):
S54、基于所述光流场(uL,vL),统计术前/术后待测视频图像中受体血管的所有像素点在X轴、Y轴方向的流向,得出待检测患者受体血管血流的最终流向。
优选的,步骤S52包含:
S521、令第三图像/第四图像作为待测光流图像;第t时刻和第t+Δt时刻的待测光流图像中,受体血管对应的像素点的亮度一致,建立公式(4):
I(x,y,t)=I(x+Δx,y+Δy,t+Δt) (4)
(x,y)表示t时刻的待测光流图像中受体血管的像素点,I(x,y,t)表示该像素点(x,y)在t时刻的亮度值;I(x+Δx,y+Δy,t+Δt)表示该待测光流图像中的的像素点(x+Δx,y+Δy)在t+Δt时刻的亮度值;其中Δx表示像素点(x,y)在X轴方向移动的距离,Δy表示像素点(x,y)在Y轴方向移动的距离,Δt表示相邻帧图像的时间差;
S522、对公式(4)右侧使用泰勒公式展开得到,
Ixu+Iyv+It=0; (5)
其中Ix表示像素点(x,y)在X轴方向的偏导,Iy表示像素点(x,y)在Y轴方向的偏导,It表示像素点(x,y)在相邻时刻的偏导,u表示像素点(x,y)在X轴方向的光流场,v表示像素点(x,y)在Y轴方向的光流场;通过引入平滑能量函数Es对像素点的运动距离进行约束,即,
Es=∫∫[(ux)2+(uy)2+(vx)2+(vy)2]dxdy; (6)
其中ux表示u在X轴方向的偏导,uy表示u在Y轴方向的偏导,vx表示v在X轴方向的偏导,vy表示v在Y轴方向的偏导;
S523、得出求解光流场的能量方程式,
minE=∫∫{(Ixu+Iyv+It)2+λ[(ux)2+(uy)2+(vx)2+(vy)2]}dxdy; (7)
其中λ是平滑项系数,λ越大平滑度越高,则对光流场进行估计的精度越高;
S524、通过欧拉-拉格朗日方程求解式(7),得到式(8)所示的两个迭代公式:
uk+1、vk+1表示第k+1次迭代得到的u、v值;是第k次迭代在像素点(x,y)邻近范围内所有像素点的光流场平均值,初始值均为0;依次迭代计算u和v的值;直到满足条件,当连续两次迭代结果对应的能量函数值E之差小于给定阈值0.1,退出迭代,得到u,v最终值。
优选的,令像素点(x,y)为待测光流图像中第i行的第j个像素点,(x,y)邻近范围内所有像素点包含该待测光流图像中第i行的第j-1、第j+1个像素点,以及第i-1、第i+1的第j个像素点。
与现有技术相比,本发明的有益效果在于:
本发明不仅可以实现血液灌注分析,还可以自动检测和匹配患者手术前后ICG-VA视频图像的血管,测量任意血管的血流方向和强度曲线,并给出校正后的手术前后的血流灌注变化图,使术者能够直观地评估吻合术引起的血流灌注的变化。
附图说明
为了更清楚地说明本发明技术方案,下面将对描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图是本发明的一个实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图:
图1为本发明的颅内外血管吻合术中血流自动定量分析方法的流程图;
图2为本发明的多任务UNet网络结构框图;
图3为八例ICG血管造影图像的分割结果图;
图4为四例手术前后造影图像配准结果图;
图5为本发明基于HS光流法的血流定量分析算法流程图;
图6为一实施例的局部光流场图;
图7为一实施例的受体血管的时间-强度曲线图;
图8为一实施例的血流灌注图;
图9为SIFT图像配准算法中,生成的两组第一高斯金字塔图像和两组差分图像示意图;
图10为对差分图像中的一个像素点搜索局部极大值、局部极小值示意图;
图11为一实施例中对差分图像的关键点生成128维的特征描述符示意图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
本发明的自动血流自动定量分析主要由三部分组成:血管分割,手术前后图像配准,血流分析。血流分析是主要任务,血管分割和图像配准是前提任务。
本发明中将ICG-VA视频分解成多帧图像,由于图像中存在许多直径差异大且分辨率低的血管,我们使用多任务UNet网络进行像素级别的分类。为了从视觉上直观比较搭桥手术前后血流灌注的差异,将目标(术前)图像向参考(术后)图像的方向进行配准。血管及其分岔点的特征在常见的基于特征配准方法中被广泛使用。由于尺度不变特征变换(SIFT)具有旋转尺度不变性,对视角变化、仿射变换和噪声也保持一定的稳定性,因此本发明通过SIFT算法提取术前术后图像中的特征用于图像配准。
光流是空间移动物体在观察成像平面上的像素移动的瞬时速度。它利用图像序列中时域中像素的变化以及相邻帧之间的相关性来计算相邻帧之间物体的运动信息。光流法可以在不知道场景信息的情况下准确地识别出运动物体的位置,且适用于视频轻微晃动的情况。实验表明,对于ICG-VA中的任何相邻帧,血流对应的像素点的亮度和运动距离的变化都很小,满足使用光流法的条件。因此本发明中使用Horn-Schunck(HS)光流方法来确定ICG-VA图像中血流的流向,以评估本发明在血流重建中的实际应用价值。
本发明提供一种颅内外血管吻合术中的血流自动定量分析方法,如图1所示,包括步骤:
S1、采集若干患者术前术后的ICG-VA视频,从每个ICG-VA视频提取一帧视频图像(在本实施例中,从视频中选取亮度稳定的最后一帧视频图像);手工标注所提取视频图像中的所有血管及受体血管;基于标注的视频图像建立训练集和测试集;
S2、建立多任务Unet网络,通过所述训练集和测试集训练该多任务Unet网络;所述多任务Unet网络用于分割ICG-VA视频图像中的所有血管以及受体血管;
多任务UNet网络是一个编码器-解码器结构的网络。如图2所示,图2左侧是一个全卷积的编码器路径,图2右侧包含两个结构相同的反卷积的解码器路径。编码器路径的输出作为两个解码器路径的输入。两个解码器路径具有不同的训练任务:图2上方的解码器路径用于从输入的视频图像中分割所有血管(包含受体血管),图2下方的解码器路径用于从输入的视频图像中分割受体血管。两个解码器路径都是像素级的分类任务。
所述编码器路径包括依序连接的四个编码器层(分别为第一至第四编码器层),所述编码器层包含依序连接的卷积层(具有3×3的卷积核)、非线性激活层和最大池化层(用于下采样)。通过最大池化层使得卷积之后生成的特征图的大小减半。
所述反卷积的解码器路径包含依序连接的四个解码器层(分别为第一至第四解码器层)、输出卷积层(具有1×1的卷积核)。每个解码器层包含一个反卷积层(用于上采样)。其中第i编码器层对应第5-i解码器层,(i∈[1,4]);每个解码器层与其相对应的编码器层的特征级联,从而实现将解码器生成的每层特征图都有效使用到后续解码器层的计算中,这有利于保留多尺度特征。
本发明的实施例中,多任务UNet网络的优化器为Adam函数;输出卷积层的激活函数采用softmax函数。
本发明的实施例中,多任务Unet网络的损失函数如公式(1)所示:
Ltotal=αLall+βLreceip (1)
其中:Ltotal表示总损失,α和β是比例因子;
且Lall=∑x∈Ω logpall(x;lall(x)),Lreceip=∑x∈Ω logpreceip(x;lreceip(x))。
Lall和Lreceip表示分割所有血管和分割受体血管的分类错误函数。x是图像空间Ω中的像素位置。pall表示经过softmax激活函数后,对于真标签lall的预测概率。preceip表示经过softmax激活函数后,对于真实标签lreceip的预测概率。
本发明的实施例中首先对采集到的120例患者的220个ICG血管造影视频的最后一帧图像进行手工标注。为了提高模型的鲁棒性和准确率,采用公开数据集DRIVE和140张ICG血管造影图像作为训练集,80张ICG血管造影图像作为测试集。
由于训练图像的数量有限,因此增加了数据集以避免过度拟合。首先,通过水平,垂直和对角线翻转来扩充数据,以便将原始数据集中的每个图像扩充为4张图像。为了提高图像分割的鲁棒性,对测试图像进行了相同的数据扩充,这意味着每个图像被预测了5次。然后,我们将5个预测的平均值作为最终预测图。
所有模型分别使用Python和MATLAB进行实验,并基于PyTorch框架实现分割。电脑配置为NVIDIA GeForce GTX 1080GPU 6GB RAM。我们使用Adam优化器并将批处理大小设置为32。我们使用了降低学习率的平稳策略。在第0、20和100次迭代中,学习率分别设置为0.01、0.001、0.0001,总迭代次数设置为150。
我们使用两个指标来评估分割性能:Jaccard指数和Dice系数。Jaccard指数定义为预测结果与真实结果的交集除以预测结果与真实结果的并集,其计算公式如下,A和B分别表示本发明预测结果和真实结果:
我们获得了Dice系数为0.80,Jaccard系数为0.73的分割结果。图3显示了ICG造影图像中部分区域的分割结果。
S3、从待检测患者的术前/术后ICG-VA视频中依时序采集若干帧术前/术后待测视频图像,通过训练好的多任务Unet网络,分割所述术前/术后待测视频图像中的所有血管得到第一图像/第二图像,并分割所述术前/术后待测视频图像中的受体血管得到第三图像/第四图像;(需要进一步强调的是,本发明中的第一图像至第四图像是一帧视频图像,而不是一张静态的jpg图片。在第一图像至第四图像中血管是不动的,血液在血管中流动)
S4、基于SIFT图像配准算法,将待检测患者的第一图像向第二图像进行配准;
步骤S4包含:
S41、分别以患者的第一、第二图像作为原始图像,通过高斯卷积生成对应的第一高斯金字塔(图像尺度空间);
第一高斯金字塔具有N组高斯图像(从第一高斯金字塔的底部至顶部,分别记为第一组至第N组高斯图像),其中每组图像包含L层图像(从第一高斯金字塔的底部至顶部,分别记为第一层至第L层高斯图像)。
第一、第二图像分别作为对应第一高斯金字塔的第一组第一层高斯图像。
如图9所示,在第一高斯金字塔的同组图像中,每一层高斯图像的尺寸都是一样的,只是平滑系数不一样。一组高斯图像也叫作一个子八度(octave)。第i层的平滑系数记为从下至上,令同组内各层高斯图像对应的平滑系数分别为:0,σ,kσ,k2σ,......,k(L-2)σ。k为比例系数。σ为平滑度因子,在Sift算法中σ的取值通常为1.6。
第一高斯金字塔同组的第i+1层高斯图像是由同组的第i层高斯图像经高斯卷积(也称为高斯平滑或称高斯滤波)之后生成的。
令J(x,y)表示待高斯卷积的图像,L(x,y,σ)为该图像的一个尺度空间表示(也即经高斯卷积生成的图像),
L(x,y,σ′)=G(x,y,σ′)*J(x,y)) (2)
其中(x,y)为坐标,σ′为待高斯卷积的图像的平滑度,*为卷积操作,G(x,y,σ′)是尺度变化的高斯函数,
第一高斯金字塔中,第j+1组的第一层图像是由第j组的第4层图像经采样得到。
S42、基于所述第一高斯金字塔生成对应的差分(DOG Difference of Gaussian)金字塔;
差分金字塔也具有N组差分图像(从差分金字塔的底部至顶部,分别记为第一组至第N组),差分金字塔的一组差分图像对应第一高斯金字塔的一组图像。差分金字塔的第1组第1层是由第一高斯金字塔的第1组第2层减第1组第1层得到的…以此类推,逐组逐层生成每一个差分图像,所有差分图像构成差分金字塔。如图9所示,也即,差分金字塔的第O组第l层图像是有第一高斯金字塔的第O组第l+1层减第O组第l层得到的。每一组在层数上,差分金字塔比第一高斯金字塔少一层。后续Sift特征点的提取都是在DOG金字塔上进行的。O∈[1,N],l∈[1,L-1]。
S43、检测差分图像的局部极大值或极小值,获取差分图像的关键点;
如图10所示,在本发明的实施例中,对于差分图像中的一个像素点(图10中的“×”)而言,它需要与自己周围的8邻域,以及尺度空间中上下两层中的相邻的18个点(图10中的“〇”)比较(共与27个点相比)。如果该像素点的像素值是局部最大值或局部最小值,它就可能是一个关键点。
S44、为每个关键点生成特征描述符.在本发明的实施例中,所述特征描述符具有128维。
图11中左部分的中央黑点为一个关键点的位置,每个小格代表该关键点邻域所在尺度空间的一个像素。
如图11所示,将坐标移至关键点主方向,在本发明的实施例中,计算以关键点为中心的16×16邻域内所有像素的梯度方向并绘制直方图,图11中每个小格(窗口)代表关键点邻域所在尺度空间的一个像素。所述梯度方向是在0~360°范围内,将梯度方向归一化到8个方向内,每个方向包含相同大小的角度范围。所述直方图中显示每个梯度方向像素点的个数。
如图11所示,关键点的16×16邻域可以划分16个4×4的块,一个块对应一个种子点,也即一个关键点的邻域由16个种子点来描述。然后在每个块上计算8个方向的梯度方向直方图(每个种子点的8个梯度方向信息,在一些实施例中,还会对梯度方向进行加权,此为现有技术)。因此关键点的特征描述符为一个包含位置、尺度、梯度方向信息的128维(16个块×8个方向)特征向量。
S45、基于最小欧式距离使用最近邻搜索,将所述第一图像向第二图像进行配准。
具体来说,若第一图像的一个关键点A的特征描述符与第二图像的一个关键点B的特征描述符具有最小欧式距离,则关键点A为关键点B的最近邻,关键点A与B匹配。匹配验证的有效度量方式是计算最近邻与次近邻的距离之比。根据多次实验,将这个比率设为0.8。
S46、得到匹配上的特征点对后,利用特征点对计算第一图像到第二图像的单应性矩阵(此为现有技术);将第一图像乘以该单应性矩阵,得到调整后的第一图像。
我们客观地评估了配准算法的整体性能。对于至少35%的血管重叠率,和具有至少10个SIFT特征对应的图像对,则认为匹配成功。成功率是成功匹配的图像对的数量与所有图像对的数量之比。在120对图像中,本发明成功匹配了97对,成功率为0.81。
图3由6组2*2的图像块组成,每个图像块中包含四个图像。图像块中这四个图像的意义分别是:图像块左上角的图像是术前ICG图像,图像块右上角是本发明对该术前ICG图像的血管分割结果,左下角是术后ICG图像,右下角是本发明对该术后ICG图像的血管分割结果。图4是对部分术前术后图像对的血管分割图匹配了关键点后进行单应性变化后的结果。图4中箭头指向的黑圆点表示匹配成功的关键点。将进行单应性变化后的术前图像与原始术后图像(即图4的A组中的两个图像)的进行重叠操作后,得到图4中的B组图像。通过图4中的B组图像可以更直观地评估配准结果。
S5、如图5所示,基于所述第三/第四图像,通过HS(Horn-Schunck)光流法的血流分析算法追踪待检测患者受体血管的血流方向。
步骤S5包含:
S51、通过高斯卷积,分别为待检测患者的第三、第四图像生成对应的第二高斯金字塔,本发明的实施例中与第三、第四图像对应的第二高斯金字塔各包含四层图像;第三图像/第四图像作为其对应的第二高斯金字塔的第一层图像;第二高斯金字塔的第i+1层图像由其第i层图像经高斯卷积之后生成(图像的卷积生成方式与步骤S41中相同,此为现有技术);i∈[1,3]。
S52、通过HS光流计算法对第二高斯金字塔中的第一层图像像素的光流场(u,v);
步骤S52包含:
S521、令第三图像/第四图像作为待测光流图像;第t时刻和第t+Δt时刻的待测光流图像中,受体血管对应的像素点的亮度一致,建立公式(4):
I(x,y,t)=I(x+Δx,y+Δy,t+Δt) (4)
(x,y)表示t时刻的待测光流图像中受体血管的像素点,I(x,y,t)表示该像素点(x,y)在t时刻的亮度值;I(x+Δx,y+Δy,t+Δt)表示该待测光流图像中的的像素点(x+Δx,y+Δy)在t+Δt时刻的亮度值;其中Δx表示像素点(x,y)在X轴方向移动的距离,Δy表示像素点(x,y)在Y轴方向移动的距离,Δt表示相邻帧图像的时间差;
S522、对公式(4)右侧使用泰勒公式展开得到,
Ixu+Iyv+It=0; (5)
其中Ix表示像素点(x,y)在X轴方向的偏导,Iy表示像素点(x,y)在Y轴方向的偏导,It表示像素点(x,y)在相邻时刻的偏导,u表示像素点(x,y)在X轴方向的光流场,v表示像素点(x,y)在Y轴方向的光流场;通过引入平滑能量函数Es对像素点的运动距离进行约束,即,
Es=∫∫[(ux)2+(uy)2+(vx)2+(vy)2]dxdy (6)
其中ux表示u在X轴方向的偏导,uy表示u在Y轴方向的偏导,vx表示v在X轴方向的偏导,vy表示v在Y轴方向的偏导;
S523、得出求解光流场的能量方程式,
min E=∫∫{(Ixu+Iyv+It)2+λ[(ux)2+(uy)2+(vx)2+(vy)2]}dxdy;
(7)
其中λ是平滑项系数,λ越大平滑度越高,则对光流场进行估计的精度越高;
S524、使用迭代计算法,基于公式(7)计算得到u和v;u可以表示血液从左往右或从右往左流,v可以表示血液从上往下流,或从下往上流。
求解式(7)是一个泛函极值问题,可用欧拉-拉格朗日方程求解,得到式(8)迭代公式:
uk+1、vk+1表示第k+1次迭代得到的u、v值;是第k次迭代在像素点(x,y)邻近范围内所有像素点的光流场平均值;令像素点(x,y)为待测光流图像中第i行的第j个像素点,(x,y)邻近范围内所有像素点包含该待测光流图像中第i行的第j-1、第j+1个像素点,以及第i-1、第i+1的第j个像素点。
S53、计算第二高斯金字塔中第L层图像像素的光流场(uL,vL);由于第二高斯金字塔图像每向上压缩一次,压缩后的图像大小为下一层图像的1/4,图像中同一像素点在相同时间水平和垂直方向的位移为原来的1/2。因此,根据第二高斯金字塔的层数L,每一层的光流估算方程为:
S54、基于所述光流场,统计每帧术前/术后待测视频图像中受体血管的所有像素点在X轴、Y轴方向的流向,得出待检测患者受体血管血流的最终流向。
下面对本实施例的血流定量分析的具体实现过程进行说明。
120例患者术前及术后共240个视频,对其分别用光流法和已有商业软件进行结果的判读,若光流法与Flow800软件对术前视频中血流方向(正向、反向)、术后视频血流方向(正向、反向、中间向两侧流)判读一致,则认为正确。对每例患者经手术后是否改变血流方向进行判读,若光流法与商用软件均显示改变了血流方向,则认为判读一致。
采用SPSS 18.0软件对实验结果数据进行统计分析。为了验证本说明方法的有效性,将HS光流法的分析结果与商业软件产生的结果进行了比较。Kappa值≤0.2为是极低一致性,>0.2~0.4为一般的一致性,>0.4~0.6为中等的一致性,>0.6~0.8为高度的一致性,>0.8~1.0为几乎完全一致。以P<0.05为差异有统计学意义。
HS光流法获得107例术前和104例术后视频与商业软件结果一致的结论。HS光流法和商业软件在判读手术前后血管流向方面有高度一致性,术前Kappa值为0.775,术后Kappa值为0.768,均P<0.01,见表1。
图6显示了某个病例的ICG视频受体血管的4帧光流场图,矢量方向反映了像素点的移动方向。图6右边的箭头指向受体血管,虚线框表示吻合口。
HS光流法判读的经过手术改变受体血管血流方向的患者有55例,商业软件发现58例,Kappa值为0.716,P=0.002,见表1。
绘制商业软件选择的感兴趣区域(ROI)的时间-亮度曲线,并基于亮度值绘制血流灌注图。本发明获得的曲线与商用软件有相似的趋势(见图7),图7的横坐标为时间,纵坐标为荧光剂反应的亮度。每一条曲线为术者选择的每一个感兴趣区域ROI的时间-亮度曲线。且本发明绘制的血流灌注图反映的血液灌注信息也与软件相似(见图8)。图8中,A列为术前由商业软件得到的血流灌注图,B列为术前通过本发明的方法得到的血流灌注图,C列为术后由商业软件得到的血流灌注图,D列为术后由本发明的方法得到的血流灌注图。
将本发明计算出的ROI达到最大亮度值所需时间与商业软件计算的延迟时间进行比较。120例术前视频中有127个ROI,120例术后视频中有178个ROI。首先,使用正态Q-Q图检验来验证两种方法获得的数据差异服从正态分布,满足T检验的前提条件。然后,使用配对样本T检验判断这两种方法得到的时间差是否具有统计学意义。术前配对T检验结果见表2,P=0.745,术后P=0.671,P>0.05,本发明的方法相比于现有商业软件检测结果更准确。
表2
以上所述,仅为本发明的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到各种等效的修改或替换,这些修改或替换都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应以权利要求的保护范围为准。
Claims (7)
1.一种颅内外血管吻合术中的血流自动定量分析方法,其特征在于,包含步骤:
S1、采集若干患者术前术后的ICG-VA视频,从每个ICG-VA视频提取一帧视频图像;手工标注所提取视频图像中的所有血管及受体血管;基于标注的视频图像建立训练集和测试集;
S2、建立多任务Unet网络,通过所述训练集和测试集训练该多任务Unet网络;所述多任务Unet网络用于分割一帧ICG-VA视频图像中所有血管以及受体血管;
所述多任务Unet网络包含一个编码器路径和结构相同的两个解码器路径;编码器路径的输出作为两个解码器路径的输入;其中一个解码器路径用于从输入的视频图像中分割所有血管,另一个解码器路径用于从输入的视频图像中分割受体血管;
所述编码器路径包括依序连接的四个编码器层,分别为第一至第四编码器层;所述编码器层包含依序连接的卷积层、非线性激活层和最大池化层;
所述解码器路径包含依序连接的四个解码器层,输出卷积层;所述四个解码器层分别为第一至第四解码器层;其中第i编码器层与第5-i解码器层的特征级联,i∈[1,4];
S3、从待检测患者的术前/术后ICG-VA视频中采集一帧术前/术后待测视频图像,通过训练好的多任务Unet网络,分割该术前/术后待测视频图像的所有血管得到对应的第一图像/第二图像,并分割该术前/术后待测视频图像中的受体血管得到对应的第三图像/第四图像;
S4、基于SIFT图像配准算法,将待检测患者的第一图像向第二图像进行配准;
S5、通过HS光流法的血流分析算法追踪待检测患者第三图像/第四图像中的术前/术后受体血管的血流方向。
2.如权利要求1所述的颅内外血管吻合术中的血流自动定量分析方法,其特征在于,多任务UNet网络的优化器为Adam函数;输出卷积层的激活函数采用softmax函数;多任务Unet网络的损失函数为:
Ltotal=αLall+βLreceip (1)
其中:Ltotal表示总损失,α和β是比例因子;
Lall=∑x∈Ωlogpall(x;lall(x)),Lreceip=∑x∈Ωlogpreceip(x;lreceip(x));
Lall和Lreceip表示分割所有血管和分割受体血管的分类错误函数;x是图像空间Ω中的像素位置;pall表示经过softmax激活函数后,对于真标签lall的预测概率;preceip表示经过softmax激活函数后,对于真实标签lreceip的预测概率。
3.如权利要求1所述的颅内外血管吻合术中的血流自动定量分析方法,其特征在于,步骤S4包含:
S41、通过高斯卷积,分别为待检测患者的第一、第二图像生成对应的第一高斯金字塔;
从第一高斯金字塔的底部至顶部,该第一高斯金字塔具有第一至第N组高斯图像;且从第一高斯金字塔的底部至顶部,每组高斯图像具有第一层至第L层高斯图像;第一图像/第二图像作为其对应的第一高斯金字塔第一组第一层图像;
第一高斯金字塔同组的第i+1层高斯图像由同组的第i层高斯图像经高斯卷积之后生成;
第一高斯金字塔中,第j+1组的第一层图像是由第j组的第4层图像经采样得到;
S42、基于所述第一高斯金字塔生成对应的差分金字塔;
从差分金字塔的底部至顶部,差分金字塔具有第一至第N组差分图像;差分金字塔的第O组第l层差分图像是有第一高斯金字塔的第O组第l+1层高斯图像减第O组第l层高斯图像得到;其中O∈[1,N],l∈[1,L-1];
S43、检测差分图像的局部极大值或极小值,获取差分图像的关键点;
S44、为每个关键点生成多维度的特征描述符;
S45、使用最近邻搜索,将第一图像的关键点与第二图像的关键点进行匹配;
S46、得到匹配上的特征点对后,利用特征点对计算第一图像到第二图像的单应性矩阵;将第一图像乘以该单应性矩阵,得到调整后的第一图像。
5.如权利要求1所述的颅内外血管吻合术中的血流自动定量分析方法,其特征在于,步骤S5包含:
S51、通过高斯卷积,分别为待检测患者的第三、第四图像生成对应的第二高斯金字塔;从所述第二高斯金字塔的底部至顶部,第二高斯金字塔包含第一至第L层图像;
第三图像/第四图像作为其对应的第二高斯金字塔的第一层图像;
第二高斯金字塔的第i+1层图像由其第i层图像经高斯卷积之后生成;i∈[1,L-1];
S52、通过HS光流法计算所述第二高斯金字塔中第一层图像像素的光流场(u,v);
S53、计算第二高斯金字塔中第L层图像像素的光流场(uL,vL):
S54、基于所述光流场(uL,vL),统计术前/术后待测视频图像中受体血管的所有像素点在X轴、Y轴方向的流向,得出待检测患者受体血管血流的最终流向。
6.如权利要求5所述的颅内外血管吻合术中的血流自动定量分析方法,其特征在于,步骤S52包含:
S521、令第三图像/第四图像作为待测光流图像;第t时刻和第t+Δt时刻的待测光流图像中,受体血管对应的像素点的亮度一致,建立公式(4):
I(x,y,t)=I(x+Δx,y+Δy,t+Δt) (4)
(x,y)表示t时刻的待测光流图像中受体血管的像素点,I(x,y,t)表示该像素点(x,y)在t时刻的亮度值;I(x+Δx,y+Δy,t+Δt)表示该待测光流图像中的的像素点(x+Δx,y+Δy)在t+Δt时刻的亮度值;其中Δx表示像素点(x,y)在X轴方向移动的距离,Δy表示像素点(x,y)在Y轴方向移动的距离,Δt表示相邻帧图像的时间差;
S522、对公式(4)右侧使用泰勒公式展开得到,
Ixu+Iyv+It=0; (5)
其中Ix表示像素点(x,y)在X轴方向的偏导,Iy表示像素点(x,y)在Y轴方向的偏导,It表示像素点(x,y)在相邻时刻的偏导,u表示像素点(x,y)在X轴方向的光流场,v表示像素点(x,y)在Y轴方向的光流场;通过引入平滑能量函数Es对像素点的运动距离进行约束,即,
Es=∫∫[(ux)2+(uy)2+(vx)2+(vy)2]dxdy; (6)
其中ux表示u在X轴方向的偏导,uy表示u在Y轴方向的偏导,vx表示v在X轴方向的偏导,vy表示v在Y轴方向的偏导;
S523、得出求解光流场的能量方程式,
min E=∫∫{(Ixu+Iyv+It)2+λ[(ux)2+(uy)2+(vx)2+(vy)2]}dxdy; (7)
其中λ是平滑项系数,λ越大平滑度越高,则对光流场进行估计的精度越高;
S524、通过欧拉-拉格朗日方程求解式(7),得到式(8)所示的两个迭代公式:
7.如权利要求6所述的颅内外血管吻合术中的血流自动定量分析方法,其特征在于,令像素点(x,y)为待测光流图像中第i行的第j个像素点,(x,y)邻近范围内所有像素点包含该待测光流图像中第i行的第j-1、第j+1个像素点,以及第i-1、第i+1的第j个像素点。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110303224.XA CN113012198B (zh) | 2021-03-22 | 2021-03-22 | 颅内外血管吻合术中的血流自动定量分析方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110303224.XA CN113012198B (zh) | 2021-03-22 | 2021-03-22 | 颅内外血管吻合术中的血流自动定量分析方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113012198A CN113012198A (zh) | 2021-06-22 |
CN113012198B true CN113012198B (zh) | 2022-04-01 |
Family
ID=76404429
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110303224.XA Active CN113012198B (zh) | 2021-03-22 | 2021-03-22 | 颅内外血管吻合术中的血流自动定量分析方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113012198B (zh) |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108122236A (zh) * | 2017-12-18 | 2018-06-05 | 上海交通大学 | 基于距离调制损失的迭代式眼底图像血管分割方法 |
CN109146872A (zh) * | 2018-09-03 | 2019-01-04 | 北京邮电大学 | 基于深度学习和光流法的心脏冠状动脉影像分割识别方法 |
Family Cites Families (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106780477A (zh) * | 2016-12-30 | 2017-05-31 | 上海联影医疗科技有限公司 | 一种血流分析方法和系统 |
US20180333120A1 (en) * | 2017-05-19 | 2018-11-22 | The Chinese University Of Hong Kong | Methods and apparatuses for quantifying vascular fluid motions from dsa |
CN109523536B (zh) * | 2018-11-19 | 2020-10-16 | 复旦大学附属华山医院 | 颅内外血流重建术的搭桥血管确定方法及系统 |
CN110473188B (zh) * | 2019-08-08 | 2022-03-11 | 福州大学 | 一种基于Frangi增强和注意力机制UNet的眼底图像血管分割方法 |
CN112164046B (zh) * | 2020-09-24 | 2022-12-16 | 中山大学中山眼科中心 | 一种眼微血管血流动力学参数自动分析方法 |
-
2021
- 2021-03-22 CN CN202110303224.XA patent/CN113012198B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108122236A (zh) * | 2017-12-18 | 2018-06-05 | 上海交通大学 | 基于距离调制损失的迭代式眼底图像血管分割方法 |
CN109146872A (zh) * | 2018-09-03 | 2019-01-04 | 北京邮电大学 | 基于深度学习和光流法的心脏冠状动脉影像分割识别方法 |
Also Published As
Publication number | Publication date |
---|---|
CN113012198A (zh) | 2021-06-22 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Wang et al. | Patch-based output space adversarial learning for joint optic disc and cup segmentation | |
Lin et al. | Automatic retinal vessel segmentation via deeply supervised and smoothly regularized network | |
Farag et al. | A bottom-up approach for pancreas segmentation using cascaded superpixels and (deep) image patch labeling | |
Guo et al. | Semi-supervised WCE image classification with adaptive aggregated attention | |
WO2018120942A1 (zh) | 一种多模型融合自动检测医学图像中病变的系统及方法 | |
CN111369528B (zh) | 基于深度卷积网络的冠状动脉血管造影图像狭窄区域标示方法 | |
Zhang et al. | Short-term lesion change detection for melanoma screening with novel siamese neural network | |
Hatamizadeh et al. | Deep dilated convolutional nets for the automatic segmentation of retinal vessels | |
Wang et al. | An automatic approach for retinal vessel segmentation by multi-scale morphology and seed point tracking | |
Zhang et al. | TUnet-LBF: Retinal fundus image fine segmentation model based on transformer Unet network and LBF | |
CN113298830A (zh) | 一种基于自监督的急性颅内ich区域图像分割方法 | |
Noh et al. | Multimodal registration of fundus images with fluorescein angiography for fine-scale vessel segmentation | |
Zhang et al. | Artifact detection in endoscopic video with deep convolutional neural networks | |
CN112419246B (zh) | 量化食管粘膜IPCLs血管形态分布的深度检测网络 | |
Mahapatra | Registration of histopathogy images using structural information from fine grained feature maps | |
Lu et al. | Adaboost-based detection and segmentation of bioresorbable vascular scaffolds struts in IVOCT images | |
CN113012198B (zh) | 颅内外血管吻合术中的血流自动定量分析方法 | |
CN116958679A (zh) | 一种基于弱监督的目标检测方法及相关设备 | |
Liu et al. | A Transformer-based pyramid network for coronary calcified plaque segmentation in intravascular optical coherence tomography images | |
Garcia-Peraza-Herrera et al. | Interpretable fully convolutional classification of intrapapillary capillary loops for real-time detection of early squamous neoplasia | |
CN111047582A (zh) | 基于深度学习的小肠镜下克罗恩病辅助诊断系统 | |
Wang et al. | An efficient hierarchical optic disc and cup segmentation network combined with multi-task learning and adversarial learning | |
Chen et al. | Spatio-temporal multi-task network cascade for accurate assessment of cardiac CT perfusion | |
Wang et al. | PCRTAM-Net: A novel pre-activated convolution residual and triple attention mechanism network for retinal vessel segmentation | |
CN114330484A (zh) | 弱监督学习糖尿病视网膜病变分级与病灶识别方法及系统 |
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 |