CN112215878B - 一种基于surf特征点的x光图像配准方法 - Google Patents
一种基于surf特征点的x光图像配准方法 Download PDFInfo
- Publication number
- CN112215878B CN112215878B CN202011216382.3A CN202011216382A CN112215878B CN 112215878 B CN112215878 B CN 112215878B CN 202011216382 A CN202011216382 A CN 202011216382A CN 112215878 B CN112215878 B CN 112215878B
- Authority
- CN
- China
- Prior art keywords
- image
- feature
- points
- point
- feature point
- 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
- 238000000034 method Methods 0.000 title claims abstract description 54
- 239000011159 matrix material Substances 0.000 claims abstract description 87
- 238000004364 calculation method Methods 0.000 claims abstract description 33
- 238000000605 extraction Methods 0.000 claims abstract description 10
- 238000007781 pre-processing Methods 0.000 claims abstract description 10
- 238000013507 mapping Methods 0.000 claims abstract description 7
- 239000013598 vector Substances 0.000 claims description 19
- 230000014509 gene expression Effects 0.000 claims description 16
- 238000012216 screening Methods 0.000 claims description 14
- 238000004590 computer program Methods 0.000 claims description 12
- 238000001914 filtration Methods 0.000 claims description 7
- 235000009508 confectionery Nutrition 0.000 claims description 6
- 238000009795 derivation Methods 0.000 claims description 6
- 230000017105 transposition Effects 0.000 claims description 5
- 238000012545 processing Methods 0.000 claims description 4
- 230000004044 response Effects 0.000 claims description 4
- 230000006870 function Effects 0.000 claims description 3
- 238000002156 mixing Methods 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 5
- 238000003745 diagnosis Methods 0.000 description 4
- 230000036544 posture Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 239000000284 extract Substances 0.000 description 3
- 230000009466 transformation Effects 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 230000000399 orthopedic effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 229960001948 caffeine Drugs 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 238000003708 edge detection Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 238000011897 real-time detection Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- RYYVLZVUVIJVGH-UHFFFAOYSA-N trimethylxanthine Natural products CN1C(=O)N(C)C(=O)C2=C1N=CN2C RYYVLZVUVIJVGH-UHFFFAOYSA-N 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- 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
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/13—Edge detection
-
- 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/10116—X-ray 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/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30008—Bone
- G06T2207/30012—Spine; Backbone
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
本发明涉及一种基于SURF特征点的X光图像配准方法、计算机设备及计算机可读存储介质,该方法包括如下步骤:选取一幅X光图像的感兴趣区域,进行预处理和线特征的提取;在用于配准的X光图像中选取相应的感兴趣区域,进行预处理和线特征的提取;基于SURF方法分别对感兴趣区域进行特征点提取;根据特征点的特征描述,组成匹配的特征点对;随机选取4组匹配的特征点对,进行投影矩阵的计算;根据投影矩阵映射,计算映射后的距离损失值;进行多轮计算,保留最小的距离损失值及对应的投影矩阵并输出,按照输出的投影矩阵进行图像配准。本发明能够快速、准确地实现多源图像配准。
Description
技术领域
本发明涉及医学图像处理技术领域,尤其涉及一种基于SURF特征点的X光图像配准方法、计算机设备及计算机可读存储介质。
背景技术
自十九世纪末伦琴发现X射线至今已有一百多年历史,医学成像技术经历了从模拟化到数字化的发展历程。最初,X射线仅应用于人体骨骼成像,至20世纪30年代新技术代旋转阳极X-线管的诞生,X光成像效果有了明显提高,开始在某些活动器官的诊断上展示出独特的优势,逐渐成为了临床研究不可缺少的检测方式。第十五届国际放射学会上,诸多专家学者对基于数字化X-光成像技术给予了肯定,随着成像技术的快速发展,数字化医学影像的应用有着广阔的发展前景。在数字信息化高速发展的今天,医疗影像的数字信息化发展为临床治疗提供了快速、便捷、精准的诊断信息,医疗影像的数字信息化也将逐渐成为新的临床研究基础应用平台。
目前,在骨科图像对病情的诊断任务中,利用X光成像检查颈椎、脊椎或其他受伤区域时,由于监测源包含多角度、多个不同姿势的X光骨科检查图像,难以清晰、完整地呈现医生所需要的所有细节信息,不利于对病人的进一步诊断和治疗,因此,迫切需要实现针对X光成像的多角度影像配准工作。
发明内容
本发明的目的是针对上述至少一部分缺陷,提供一种基于图像特征信息的不同角度、不同姿势的X光图像配准方法。
为了实现上述目的,本发明提供了一种基于SURF特征点的X光图像配准方法,该方法包括如下步骤:
S1、选取一幅X光图像的感兴趣区域,记为I1图像,并对I1图像进行预处理和线特征的提取;
S2、在用于配准的X光图像中选取相应的感兴趣区域,记为I2图像,并对I2图像进行预处理和线特征的提取;
S3、基于SURF方法分别对I1图像和I2图像进行特征点提取;
S4、根据特征点的特征描述,将取自I2图像的各特征点与取自I1图像的各特征点进行匹配,组成匹配的特征点对;其中,匹配的特征点对中两个特征点的特征描述相符度不超过设定阈值;
S5、随机选取4组匹配的特征点对,进行投影矩阵的计算;
S6、根据投影矩阵将I1图像映射到I2图像上,计算映射后I1图像和I2图像的距离损失值;其中,距离损失值为欧式距离下I1图像和I2图像同一线特征之间最近点的点均方误差;
S7、将距离损失值与当前的存储结果进行比较,保留更小的距离损失值及对应的投影矩阵作为新的存储结果;
S8、判断是否进行下一轮计算,若进行,则返回步骤S5,否则输出存储结果中的距离损失值和投影矩阵;
S9、将I1图像和I2图像按照输出的投影矩阵进行图像配准。
优选地,所述步骤S1对I1图像进行预处理时,对I1图像进行进行低通滤波,candy边缘检测;
所述步骤S2中,对I2图像进行预处理时,对I2图像进行进行低通滤波,candy边缘检测。
优选地,所述步骤S3中,基于SURF方法分别对I1图像和I2图像进行特征点提取时,包括如下步骤:
S3-1、构建Hessian矩阵,表达式为:
其中,x、y表示像素点的坐标,Dxx(x,σ)表示对x方向的二阶导数,Dyy(x,σ)表示对y方向的二阶导数,Dxy(x,σ)表示对x方向求一阶偏导后再对y方向求一阶偏导,σ表示滤波器的尺度;
S3-2、构造SURF方法的尺度空间,尺度空间由s组o层组成,滤波器为盒式滤波器,内部使用高斯二阶微分滤波,不同组间图像的尺寸一致,使用的滤波器的模板尺寸逐渐增大,同一组、不同层的图像使用相同尺寸的滤波器,但滤波器的尺度空间因子逐渐增大,关系式为:l=2o+1(s+1)+1、L=3*l、其中,l表示滤波器在微分方向上对正负斑点响应长度,L表示滤波器的模板尺寸,s为组数,o为层数;
S3-3、以Hessian矩阵处理需进行特征点提取的图像中每个像素点,计算当前像素点的Hessian矩阵判别式:
H(x,σ)=Dxx(x,σ)*Dyy(x,σ)-(Dxy(x,σ))2
当Hessian矩阵判别式取得局部极大值时,选择当前像素点为兴趣点;
对选择出来的兴趣点,在连续空间进行泰勒展开,表达式为:
其中,D(X)表示选择出来的兴趣点信息,X=(x,y,σ)为两点之间的偏移位置,XT表示X的转置,DT表示D的转置,D为高斯差分函数;
通过求导为0的方式,确定连续空间上的极值点,求导公式为:
S3-4、在筛选得到的特征点的圆形邻域内,以0.2rad为步长进行扫描,统计60°扇形范围内所有特征点的水平、垂直的Harr小波特征绝对值总和,将最大值对应的扇形的方向作为该特征点的主方向;
S3-5、对确定主方向的各特征点,以特征点为中心,沿着该特征点的主方向提取4×4个矩形区域块,单个矩形区域块大小为5*5,对每个矩形区域块,分别统计25个像素点在相对于该特征点的主方向的水平方向、垂直方向的Haar小波特征,最终汇总所有矩形块区域的统计结果,得到该特征点的特征描述,包括水平方向Haar小波特征值之和、垂直方向Haar小波特征值之和、水平方向Haar小波特征值绝对值之和以及垂直方向Haar小波特征值绝对值之和,特征描述以特征向量的形式表示。
优选地,所述步骤S4中,将取自I2图像的各特征点与取自I1图像的各特征点进行匹配时,计算各取自I2图像的特征点与一个取自I1图像的特征点之间的特征向量误差值,并将特征向量误差值最小的、取自I2图像的特征点与该取自I1图像的特征点构成一组初步匹配特征点对;
以特征向量误差值作为特征描述相符度,将所有组初步匹配特征点对的特征向量误差值从小到大排序,将特征向量误差值不超过设定阈值的各组特征点对确定为匹配的特征点对。
优选地,所述步骤S5中,随机选取4组匹配的特征点对,进行投影矩阵的计算时,包括如下步骤:
S5-1、随机选取4组匹配的特征点对,对每一组匹配的特征点对建立方程:
得到两个关于投影矩阵元素的显式表达式:
(h31xi+h32yi+h33)·xi'=h11xi+h12yi+h13
(h31xi+h32yi+h33)·yi'=h21xi+h22yi+h23
其中,(xi,yi)表示第i组匹配的特征点对中取自I1图像的特征点,(x'i,y'i)表示第i组匹配的特征点对中取自I2图像的特征点,hjk表示投影矩阵元素,j、k=1、2、3;
S5-2、将投影矩阵元素重写为向量形式,得到如下关系式:
根据矩阵系数的线性特性,设h33恒为1,将每一组匹配的特征点对得到的两个关于投影矩阵元素的显式表达式代入,得到投影矩阵。
优选地,所述步骤S8中,判断是否进行下一轮计算时,检查本轮计算的计算轮数,若计算轮数未达到预设轮数,则返回步骤S5进行下一轮计算,且计算轮数加1。
优选地,所述步骤S8中,判断是否进行下一轮计算时,判断当前计算过的匹配的特征点对,若任意4组匹配的特征点对均已进行过投影矩阵的计算,则输出存储结果中的距离损失值和投影矩阵。
本发明还提供了一种计算机设备,包括存储器和处理器,所述存储器存储有计算机程序,所述处理器执行所述计算机程序时实现上述任一项所述基于SURF特征点的X光图像配准方法的步骤。
本发明还提供了一种计算机可读存储介质,其上存储有计算机程序,所述计算机程序被处理器执行时实现上述任一项所述基于SURF特征点的X光图像配准方法的步骤。
本发明的上述技术方案具有如下优点:本发明提供了一种基于SURF特征点的X光图像配准方法、计算机设备及计算机可读存储介质,该方法通过分别提取同一场景、不同角度和/或姿势的X光图像的线特征和点特征信息,先通过点特征信息计算出用于配准的投影矩阵,然后通过线特征进行评价,实现对投影矩阵的评估与更新,最终确定最佳的投影矩阵,实现高精度、高匹配度的多源图像配准,可应用于对X光成像采集的颈椎、脊椎等多类型目标图像,能够进行准确、实时地检测与识别配准。
附图说明
图1是本发明实施例中一种基于SURF特征点的X光图像配准方法步骤示意图;
图2是本发明实施例中基于SURF方法进行特征点提取的步骤示意图;
图3是本发明实施例中基于SURF特征点的X光图像配准方法进行配准得到的结果图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明的一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都属于本发明保护的范围。
如图1至图3所示,本发明实施例提供的一种基于SURF特征点的X光图像配准方法,包括如下步骤:
S1、选取一幅X光图像的感兴趣区域,记为I1图像,并对I1图像进行预处理和线特征的提取。
对于一幅X光图像(通常为多源图像中的主要图像,可记为图像1),选取其感兴趣区域(ROI,region ofinterest)后,可得到更小幅的图像,即I1图像,对I1图像进行预处理能够提高图像清晰度,便于后续的线特征提取。
提取线特征时,可根据ROI中的具体目标的特点,选取一条或多条线特征。
S2、在用于配准的X光图像中选取相应的感兴趣区域,记为I2图像,并对I2图像进行预处理和线特征的提取。
用于配准的X光图像(可记为图像2),即与主要图像(即图像1)在同一场景、不同角度和/或姿势的X光图像,其感兴趣区域(即I2图像)通常为与I1图像相同或相近的区域,以便进行配准。对于I2图像进行的预处理和线特征提取,应与对于I1图像进行的预处理和线特征提取相同。
在一个优选的实施方式中,步骤S1对I1图像进行预处理时,对I1图像进行进行低通滤波,candy边缘检测。相应的,步骤S2中,对I2图像进行预处理时,也对I2图像进行进行相同的低通滤波,candy边缘检测。
S3、基于SURF方法分别对I1图像和I2图像进行特征点提取。
SURF(Speeded Up Robust Features,加速稳健特征)方法是对SIFT(Scale-invariant feature transform,尺度不变特征转换)方法的改进,特别是在执行效率上。通过SURF方法提取X光图像感兴趣区域的点特征信息,能够实现图像实时配准。
S4、根据特征点的特征描述,将取自I2图像的各特征点与取自I1图像的各特征点进行匹配,组成匹配的特征点对。其中,匹配的特征点对中,两个特征点的特征描述相符度不超过设定阈值。
两个特征点的特征描述相符度,即取自I2图像的特征点与取自I1图像的特征点之间对应的特征向量误差值。
S5、随机选取4组匹配的特征点对,进行投影矩阵的计算。
投影矩阵,即将I1图像映射到I2图像的变换矩阵。
S6、根据投影矩阵将I1图像映射到I2图像上,计算映射后I1图像和I2图像的距离损失值。
其中,距离损失值为欧式距离下,(映射到I2图像上的)I1图像和I2图像同一线特征之间最近点的点均方误差。
S7、将距离损失值与当前的存储结果进行比较,保留更小的距离损失值及对应的投影矩阵作为新的存储结果。
将距离损失值与当前的存储结果进行比较时,若当前的存储结果为空集,则直接保存本轮计算得到的距离损失值及对应的投影矩阵作为新的存储结果,若当前的存储结果并非空集,则将本轮计算得到的距离损失值与当前存储结果中的距离损失值比较,以便在多轮计算后,保存最小的距离损失值及对应的投影矩阵作为存储结果。
S8、判断是否进行下一轮计算,若进行,则返回步骤S5,否则输出存储结果中的距离损失值和投影矩阵。
通过多次重复计算投影矩阵,能够通过迭代筛选出配准效果最好的投影矩阵,从而提高多源图像配准的精度与匹配度。
S9、将I1图像和I2图像按照输出的投影矩阵进行图像配准。
本发明对两幅X光图像分别进行线特征和点特征信息的提取,通过点特征信息配对,计算出投影矩阵,通过线特征对投影矩阵进行评估,选取效果最好的投影矩阵,可实现对不同角度和/或姿态的多源X光图像的高质量配准。除此以外,在整个配准流程中,预处理去除噪声等干扰信息,后续计算只使用了X光图像本身的线特征和点特征信息,尤其是其中的点特征信息,对不同类型、不同尺度、不同亮度的目标是鲁棒的,能够实现对多种场景配准的鲁棒性。
优选地,如图2所示,步骤S3中,基于SURF方法分别对I1图像和I2图像进行特征点提取时,包括如下步骤:
S3-1、构建Hessian矩阵。
其中,x、y表示像素点的坐标,Dxx(x,σ)表示对x方向的二阶导数,Dyy(x,σ)表示对y方向的二阶导数,Dxy(x,σ)表示对x方向求一阶偏导后再对y方向求一阶偏导,σ表示滤波器的尺度。
当Hessian矩阵的判别式H(x,σ)=Dxx(x,σ)*Dyy(x,σ)-(Dxy(x,σ))2取得局部极大值时,可判定当前像素点是比周围邻域内其他点更亮或更暗的点。
S3-2、构造尺度空间。
SURF方法的尺度空间由s组o层组成,滤波器为盒式滤波器,内部使用高斯二阶微分滤波,不同组之间的图像的尺寸一致,使用的滤波器的模板尺寸逐渐增大,同一组但不同层的图像使用相同尺寸的滤波器,但滤波器的尺度空间因子逐渐增大,具体关系式为:l=2o+1(s+1)+1、L=3*l、其中,l表示滤波器在微分方向上对正负斑点响应长度,L表示滤波器的模板尺寸,s为组数,o为层数,σ为滤波器的尺度。
S3-3、定位特征点。
以Hessian矩阵处理需进行特征点提取的图像(即I1图像、I2图像)中的每个像素点,计算当前像素点的Hessian矩阵判别式:H(x,σ)=Dxx(x,σ)*Dyy(x,σ)-(Dxy(x,σ))2,当Hessian矩阵判别式取得局部极大值时,选择当前像素点为兴趣点(也即找到像素点的极值点)。
将每个像素点的Hessian矩阵判别式与其图像域(相同大小的图像)和尺度域(相邻的尺度空间)的所有相邻点进行比较。当一个像素点像素值大于(或者小于)所有相邻点时,该像素点为极值点。检测的像素点要和其所在图像的3×33×3邻域的8个像素点(图像域相邻点),以及其相邻的上下两层3×33×3邻域的18个像素点(尺度域相邻点),共26个像素点,进行比较。通过Hessian矩阵判别式能去除像素值比较结果中的边缘响应的点。
对各个选择出来的兴趣点,在连续空间上进行泰勒展开,表达式为:
其中,D(X)表示选择出来的兴趣点信息,X=(x,y,σ)为两点之间的偏移位置,XT表示X的转置,DT表示D的转置,D为高斯差分函数。
通过求导为0的方式,确定兴趣点中的极值点,也即连续空间上的极值点,求导公式为:
对通过求导为0确定的极值点中的所有极大值点求Hessian矩阵,并通过筛选得到特征点,其中r为指定筛选参数,Tr(H)表示Hessian矩阵的迹,Det(H)表示Hessian矩阵的行列式。指定筛选参数r的具体值可根据实际情况进行设定。
进一步地,由于Hessian矩阵判别式中使用的是原始图像的高斯卷积,由于高斯核实服从正态分布的,从中心点往外,系数越来越低,为了提高运算速度,SURF方法使用了盒式滤波器来近似替代高斯滤波器,所以,可在Dxy上乘一个加权系数0.9,目的是为了平衡因使用盒式滤波器近似所带来的误差,即步骤S3-3中,通过筛选得到特征点时,Det(H)=Dxx*Dyy-(0.9*Dxy)2。
S3-4、确定各特征点的主方向。
在筛选得到的各个特征点的圆形邻域内,以0.2rad为步长进行扫描,统计60°扇形范围内所有特征点的水平的Harr小波特征绝对值、垂直的Harr小波特征绝对值总和,将最大值对应的扇形的方向作为该特征点的主方向。
特征点的圆形邻域,即以其自身为中心,6σ为半径的圆形范围内,σ为(特征点所在图像所使用的)滤波器的尺度。
S3-5、生成各特征点的特征描述。
对确定主方向的各特征点,以特征点为中心,沿着该特征点的主方向提取4×4个矩形区域块,单个矩形区域块大小为5*5,对每个矩形区域块,分别统计25个像素点在相对于该特征点的主方向的水平方向和垂直方向的Haar小波特征,最终汇总所有矩形块区域的统计结果,得到该特征点的特征描述,一个特征点的特征描述包括水平方向Haar小波特征值之和、垂直方向Haar小波特征值之和、水平方向Haar小波特征值绝对值之和以及垂直方向Haar小波特征值绝对值之和,特征描述以特征向量的形式表示。
优选地,步骤S4中,将取自I2图像的各特征点与取自I1图像的各特征点进行匹配时,计算各取自I2图像的特征点与一个取自I1图像的特征点之间的特征向量误差值,并将特征向量误差值最小的、取自I2图像的特征点与该取自I1图像的特征点构成一组初步匹配特征点对。按照同样的方式处理其他取自I1图像的特征点,直到将所有特征点都构成初步匹配特征点对。
以(一组初步匹配特征点对中)取自I1图像的特征点与取自I2图像的特征点之间的特征向量误差值,作为(一组初步匹配特征点对的)特征描述相符度,将所有组初步匹配特征点对的特征向量误差值从小到大排序,将特征向量误差值不超过设定阈值的各组特征点对确定为匹配的特征点对。筛选出相符度最高的特征点对,有利于后续计算确定最佳的投影矩阵。
优选地,步骤S5中,随机选取4组匹配的特征点对,进行投影矩阵的计算时,包括如下步骤:
S5-1、随机选取4组匹配的特征点对,对每一组匹配的特征点对建立方程:
得到两个关于投影矩阵元素的显式表达式:
(h31xi+h32yi+h33)·xi'=h11xi+h12yi+h13
(h31xi+h32yi+h33)·yi'=h21xi+h22yi+h23
其中,(xi,yi)表示第i组匹配的特征点对中取自I1图像的特征点,(x'i,y'i)表示第i组匹配的特征点对中取自I2图像的特征点,hjk表示投影矩阵元素,j、k=1、2、3;
S5-2、将投影矩阵元素重写为向量形式,得到如下关系式:
根据矩阵系数的线性特性,设h33恒为1,得到唯一变换矩阵,至此,变换矩阵的自由度为8,将每一组匹配的特征点对得到的两个关于投影矩阵元素的显式表达式(共8个显式表达式)代入,计算得到投影矩阵。
在一些实施方式中,步骤S8中,判断是否进行下一轮计算时,可检查本轮计算的计算轮数,若计算轮数未达到预设轮数,则返回步骤S5进行下一轮计算,且计算轮数加1。
在另一些实施方式中,步骤S8中,判断是否进行下一轮计算时,可判断当前计算过的匹配的特征点对,若任意4组匹配的特征点对均已进行过投影矩阵的计算,则不再返回步骤S5进行下一轮计算,输出存储结果中的距离损失值和投影矩阵。
特别地,在本发明一些优选的实施方式中,还提供了一种计算机设备,包括存储器和处理器,所述存储器存储有计算机程序,所述处理器执行所述计算机程序时实现上述任一实施方式中所述基于SURF特征点的X光图像配准方法的步骤。
在本发明另一些优选的实施方式中,还提供了一种计算机可读存储介质,其上存储有计算机程序,计算机程序被处理器执行时实现上述任一实施方式中所述基于SURF特征点的X光图像配准方法的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程,是可以通过计算机程序来指令相关的硬件来完成,计算机程序可存储于一非易失性计算机可读取存储介质中,该计算机程序在执行时,可包括如上述基于SURF特征点的X光图像配准方法实施例的流程,在此不再重复说明。
最后应说明的是:以上实施例仅用以说明本发明的技术方案,而非对其限制;尽管参照前述实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换;而这些修改或者替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
Claims (10)
1.一种基于SURF特征点的X光图像配准方法,其特征在于,包括如下步骤:
S1、选取一幅X光图像的感兴趣区域,记为I1图像,并对I1图像进行预处理和线特征的提取;
S2、在用于配准的X光图像中选取相应的感兴趣区域,记为I2图像,并对I2图像进行预处理和线特征的提取;
S3、基于SURF方法分别对I1图像和I2图像进行特征点提取;
S4、根据特征点的特征描述,将取自I2图像的各特征点与取自I1图像的各特征点进行匹配,组成匹配的特征点对;其中,匹配的特征点对中两个特征点的特征描述相符度不超过设定阈值;
S5、随机选取4组匹配的特征点对,进行投影矩阵的计算;
S6、根据投影矩阵将I1图像映射到I2图像上,计算映射后I1图像和I2图像的距离损失值;其中,距离损失值为欧式距离下I1图像和I2图像同一线特征之间最近点的点均方误差;
S7、将距离损失值与当前的存储结果进行比较,保留更小的距离损失值及对应的投影矩阵作为新的存储结果;
S8、判断是否进行下一轮计算,若进行,则返回步骤S5,否则输出存储结果中的距离损失值和投影矩阵;
S9、将I1图像和I2图像按照输出的投影矩阵进行图像配准。
2.根据权利要求1所述的基于SURF特征点的X光图像配准方法,其特征在于:
所述步骤S1对I1图像进行预处理时,对I1图像进行低通滤波,candy边缘检测;
所述步骤S2中,对I2图像进行预处理时,对I2图像进行低通滤波,candy边缘检测。
3.根据权利要求1所述的基于SURF特征点的X光图像配准方法,其特征在于,
所述步骤S3中,基于SURF方法分别对I1图像和I2图像进行特征点提取时,包括如下步骤:
S3-1、构建Hessian矩阵,表达式为:
其中,x、y表示像素点的坐标,Dxx(x,σ)表示对x方向的二阶导数,Dyy(x,σ)表示对y方向的二阶导数,Dxy(x,σ)表示对x方向求一阶偏导后再对y方向求一阶偏导,σ表示滤波器的尺度;
S3-2、构造SURF方法的尺度空间,尺度空间由s组o层组成,滤波器为盒式滤波器,内部使用高斯二阶微分滤波,不同组间图像的尺寸一致,使用的滤波器的模板尺寸逐渐增大,同一组、不同层的图像使用相同尺寸的滤波器,但滤波器的尺度空间因子逐渐增大,关系式为:l=2o+1(s+1)+1、L=3*l、其中,l表示滤波器在微分方向上对正负斑点响应长度,L表示滤波器的模板尺寸,s为组数,o为层数;
S3-3、以Hessian矩阵处理需进行特征点提取的图像中每个像素点,计算当前像素点的Hessian矩阵判别式:
H(x,σ)=Dxx(x,σ)*Dyy(x,σ)-(Dxy(x,σ))2
当Hessian矩阵判别式取得局部极大值时,选择当前像素点为兴趣点;
对选择出来的兴趣点,在连续空间进行泰勒展开,表达式为:
其中,D(X)表示选择出来的兴趣点信息,X=(x,y,σ)为两点之间的偏移位置,XT表示X的转置,DT表示D的转置,D为高斯差分函数;
通过求导为0的方式,确定连续空间上的极值点,求导公式为:
S3-4、在筛选得到的特征点的圆形邻域内,以0.2rad为步长进行扫描,统计60°扇形范围内所有特征点的水平、垂直的Harr小波特征绝对值总和,将最大值对应的扇形的方向作为该特征点的主方向;
S3-5、对确定主方向的各特征点,以特征点为中心,沿着该特征点的主方向提取4×4个矩形区域块,单个矩形区域块大小为5*5,对每个矩形区域块,分别统计25个像素点在相对于该特征点的主方向的水平方向、垂直方向的Haar小波特征,最终汇总所有矩形块区域的统计结果,得到该特征点的特征描述,包括水平方向Haar小波特征值之和、垂直方向Haar小波特征值之和、水平方向Haar小波特征值绝对值之和以及垂直方向Haar小波特征值绝对值之和,特征描述以特征向量的形式表示。
5.根据权利要求3所述的基于SURF特征点的X光图像配准方法,其特征在于:
所述步骤S4中,将取自I2图像的各特征点与取自I1图像的各特征点进行匹配时,计算各取自I2图像的特征点与一个取自I1图像的特征点之间的特征向量误差值,并将特征向量误差值最小的、取自I2图像的特征点与该取自I1图像的特征点构成一组初步匹配特征点对;
以特征向量误差值作为特征描述相符度,将所有组初步匹配特征点对的特征向量误差值从小到大排序,将特征向量误差值不超过设定阈值的各组特征点对确定为匹配的特征点对。
6.根据权利要求1所述的基于SURF特征点的X光图像配准方法,其特征在于:
所述步骤S5中,随机选取4组匹配的特征点对,进行投影矩阵的计算时,包括如下步骤:
S5-1、随机选取4组匹配的特征点对,对每一组匹配的特征点对建立方程:
得到两个关于投影矩阵元素的显式表达式:
(h31xi+h32yi+h33)·xi'=h11xi+h12yi+h13
(h31xi+h32yi+h33)·yi'=h21xi+h22yi+h23
其中,(xi,yi)表示第i组匹配的特征点对中取自I1图像的特征点,(x'i,y'i)表示第i组匹配的特征点对中取自I2图像的特征点,hjk表示投影矩阵元素,j、k=1、2、3;
S5-2、将投影矩阵元素重写为向量形式,得到如下关系式:
根据矩阵系数的线性特性,设h33恒为1,将每一组匹配的特征点对得到的两个关于投影矩阵元素的显式表达式代入,得到投影矩阵。
7.根据权利要求1所述的基于SURF特征点的X光图像配准方法,其特征在于:
所述步骤S8中,判断是否进行下一轮计算时,检查本轮计算的计算轮数,若计算轮数未达到预设轮数,则返回步骤S5进行下一轮计算,且计算轮数加1。
8.根据权利要求1所述的基于SURF特征点的X光图像配准方法,其特征在于:
所述步骤S8中,判断是否进行下一轮计算时,判断当前计算过的匹配的特征点对,若任意4组匹配的特征点对均已进行过投影矩阵的计算,则输出存储结果中的距离损失值和投影矩阵。
9.一种计算机设备,包括存储器和处理器,所述存储器存储有计算机程序,其特征在于,所述处理器执行所述计算机程序时实现权利要求1至8中任一项所述基于SURF特征点的X光图像配准方法的步骤。
10.一种计算机可读存储介质,其上存储有计算机程序,其特征在于,所述计算机程序被处理器执行时实现权利要求1至8中任一项所述基于SURF特征点的X光图像配准方法的步骤。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011216382.3A CN112215878B (zh) | 2020-11-04 | 2020-11-04 | 一种基于surf特征点的x光图像配准方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011216382.3A CN112215878B (zh) | 2020-11-04 | 2020-11-04 | 一种基于surf特征点的x光图像配准方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112215878A CN112215878A (zh) | 2021-01-12 |
CN112215878B true CN112215878B (zh) | 2023-03-24 |
Family
ID=74058150
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202011216382.3A Active CN112215878B (zh) | 2020-11-04 | 2020-11-04 | 一种基于surf特征点的x光图像配准方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112215878B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115841506B (zh) * | 2023-02-20 | 2023-05-02 | 广东省人民医院 | 一种荧光分子图像处理方法及*** |
CN116728420B (zh) * | 2023-08-11 | 2023-11-03 | 苏州安博医疗科技有限公司 | 一种脊柱外科手术用机械臂调控方法及*** |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106683127A (zh) * | 2017-01-05 | 2017-05-17 | 南京觅踪电子科技有限公司 | 一种基于surf算法的多模态医学图像配准方法 |
CN109271996A (zh) * | 2018-08-21 | 2019-01-25 | 南京理工大学 | 基于surf特征和哈希感知算法的fpc图像自动配准方法 |
CN109727239A (zh) * | 2018-12-27 | 2019-05-07 | 北京航天福道高技术股份有限公司 | 基于surf特征对巡检图与基准图的配准方法 |
CN111784675A (zh) * | 2020-07-01 | 2020-10-16 | 云南易见纹语科技有限公司 | 物品纹理信息处理的方法、装置、存储介质及电子设备 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106356757B (zh) * | 2016-08-11 | 2018-03-20 | 河海大学常州校区 | 一种基于人眼视觉特性的电力线路无人机巡检方法 |
-
2020
- 2020-11-04 CN CN202011216382.3A patent/CN112215878B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106683127A (zh) * | 2017-01-05 | 2017-05-17 | 南京觅踪电子科技有限公司 | 一种基于surf算法的多模态医学图像配准方法 |
CN109271996A (zh) * | 2018-08-21 | 2019-01-25 | 南京理工大学 | 基于surf特征和哈希感知算法的fpc图像自动配准方法 |
CN109727239A (zh) * | 2018-12-27 | 2019-05-07 | 北京航天福道高技术股份有限公司 | 基于surf特征对巡检图与基准图的配准方法 |
CN111784675A (zh) * | 2020-07-01 | 2020-10-16 | 云南易见纹语科技有限公司 | 物品纹理信息处理的方法、装置、存储介质及电子设备 |
Non-Patent Citations (1)
Title |
---|
"红外与可见光遥感图像自动配准算法研究";廉蔺;《中国博士学位论文全文数据库 信息科技辑》;20141015;第二章至第六章 * |
Also Published As
Publication number | Publication date |
---|---|
CN112215878A (zh) | 2021-01-12 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP6077993B2 (ja) | 画像の異形を識別するための画像データの処理方法、システムおよびプログラム | |
JP4130661B2 (ja) | 時間的に連続する胸部画像間の経時変化を検出する装置 | |
CN106934821B (zh) | 一种基于icp算法和b样条的锥形束ct和ct图像配准方法 | |
CN112215878B (zh) | 一种基于surf特征点的x光图像配准方法 | |
JP4138371B2 (ja) | 解剖学的特徴位置検出装置並びに記録媒体および被写体構造計測装置並びに記録媒体 | |
CN110782428B (zh) | 一种用于构建临床脑部ct图像roi模板的方法及*** | |
CN106296613B (zh) | 一种基于dr机器的双能量减影方法 | |
Wang et al. | Automatic fundus images mosaic based on SIFT feature | |
EP1551296A1 (en) | Method, code, and system for assaying joint deformity | |
CN113298853A (zh) | 一种逐级递进的二阶段医学图像配准方法 | |
CN111062936B (zh) | 用于面部变形诊疗效果的量化指标评估方法 | |
CN111179173B (zh) | 一种基于离散小波变换和坡度融合算法的图像拼接方法 | |
CN114782715B (zh) | 一种基于统计信息的静脉识别方法 | |
CN115861656A (zh) | 用于自动处理医学图像以输出警报的方法、设备和*** | |
CN113077502B (zh) | 基于标志球内部多层球面生成点的牙颌空间配准方法 | |
CN112801141B (zh) | 基于模板匹配和孪生神经网络优化的异源图像匹配方法 | |
Kusuma | Determination of spinal curvature from scoliosis X-ray images using K-means and curve fitting for early detection of scoliosis disease | |
CN113610746A (zh) | 一种图像处理方法、装置、计算机设备及存储介质 | |
JP2004188202A (ja) | デジタル胸部放射線写真の自動分析方法 | |
Martín-Fernández et al. | Automatic articulated registration of hand radiographs | |
Malode | New approach of statistical analysis for lung disease diagnosis using microscopy images | |
Iakovidis et al. | Robust model-based detection of the lung field boundaries in portable chest radiographs supported by selective thresholding | |
CN114581340A (zh) | 一种图像校正方法及设备 | |
Cao et al. | DSA image registration based on multiscale Gabor filters and mutual information | |
CN112766332A (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 |