CN115035124B - 基于Harris角点检测的随动定位系统导针计算方法 - Google Patents

基于Harris角点检测的随动定位系统导针计算方法 Download PDF

Info

Publication number
CN115035124B
CN115035124B CN202210971975.3A CN202210971975A CN115035124B CN 115035124 B CN115035124 B CN 115035124B CN 202210971975 A CN202210971975 A CN 202210971975A CN 115035124 B CN115035124 B CN 115035124B
Authority
CN
China
Prior art keywords
image
dmax
point
points
plane
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
CN202210971975.3A
Other languages
English (en)
Other versions
CN115035124A (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.)
Nanjing Vishee Medical Technology Co Ltd
Original Assignee
Nanjing Vishee Medical Technology Co Ltd
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 Nanjing Vishee Medical Technology Co Ltd filed Critical Nanjing Vishee Medical Technology Co Ltd
Priority to CN202210971975.3A priority Critical patent/CN115035124B/zh
Publication of CN115035124A publication Critical patent/CN115035124A/zh
Application granted granted Critical
Publication of CN115035124B publication Critical patent/CN115035124B/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
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10088—Magnetic resonance imaging [MRI]
    • 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/20112—Image segmentation details
    • G06T2207/20164—Salient point detection; Corner 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/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30016—Brain

Landscapes

  • Engineering & Computer Science (AREA)
  • Quality & Reliability (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Health & Medical Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)
  • Image Analysis (AREA)

Abstract

本发明公布了一种基于Harris角点检测的随动定位系统导针计算方法，该方法仅采用结构像制作DLPFC靶区就可以实现快速、精准的导航定位，得到一种操作简便、快捷、自动化、定位精准、成本低廉的经颅磁定位用导航系统，成功解决了磁刺激定位技术问题，降低患者拍摄磁共振的成本，节约治疗师的治疗时间。相比其他定位导航技术，本发明不需要附加视觉系统和机械臂系统即可达到同样的定位效果，将精度误差控制在毫米级甚至更高，且有效避免了机械设定的误差，系统自动进行靶点计算，整个治疗过程5min内可以输出靶点坐标。

Description

基于Harris角点检测的随动定位系统导针计算方法
技术领域
本发明涉及一种导针计算方法，具体涉及一种基于Harris角点检测的随动定位导航系统的导针计算方法，属于图像处理技术领域。
背景技术
经颅磁刺激可用于治疗抑郁症、偏瘫、失语等精神类疾病，在进行脑刺激治疗时，如何定位并找到头部正确的靶点或病灶是至关重要的，一般需要经验丰富的临床医生根据经验进行人工定位。具体操作手法如下：临床医生根据经验用拍头刺激患者头部，然后在设备上查看脉冲波，如波形不正确，则重复刺激，直到找到正确的靶点位置，但该方法操作十分繁琐，耗时费力，耽误治疗时间，且定位不精确，会增加患者的痛苦，在临床应用中十分受限。
随着技术的发展，除传统的人工定位方法外，用于治疗的经颅磁刺激系统能够依赖于视觉相机和机械臂实现定位导航，但是该导航架结构复杂、价格昂贵、不易普及。开发出的新型随动定位系统，对于功能磁共振成像技术而言，只需计算体素时间序列间的相关性即可确定磁刺激靶点，但功能磁共振成像采集方式复杂、对患者要求较高，所以其应用推广受到了一定程度的限制。结构磁共振成像采集方式具有简单的优点，但结构磁共振成像无时间序列，因此需要手动制作靶区来生成靶点，耗时费力且过度依赖于医生的技术。
鉴于上述原因，如何优化经颅定位用导航系统，解决结构磁共振成像技术中的磁刺激定位技术问题，是亟待解决的行业难题。
发明内容
为解决现有技术的不足，本发明的目的在于提供一种适用于随动定位导航系统的导针计算方法，以实现简便、快捷、准确、低成本地解决结构磁共振成像技术中的磁刺激定位问题。
为了实现上述目标，本发明采用如下的技术方案：
基于Harris角点检测的随动定位系统导针计算方法，包括如下步骤：
S1、构建DLPFC靶区图像
Figure 289032DEST_PATH_IMAGE001
；
根据磁共振模板空间下的坐标生成二值球体图像
Figure 923276DEST_PATH_IMAGE002
，将
Figure 286124DEST_PATH_IMAGE003
与磁共振模板空间下的脑区图像掩膜
Figure 728738DEST_PATH_IMAGE004
相乘，获取两者重合的区域，即为所需靶区图像
Figure 89312DEST_PATH_IMAGE001
；
S2、获取
Figure 742010DEST_PATH_IMAGE001
中的所有角点；
逐层取DLPFC三维靶区图像
Figure 518336DEST_PATH_IMAGE001
的每一层进行Harris检测，得到三维的角点分布图；
S3、求取任意两个角点之间的欧几里得距离；
S4、取距离排序在前两位的角点坐标，输出距离排序在前两位的dmax1、dmax2及两个距离对应的坐标
Figure 471249DEST_PATH_IMAGE005
、
Figure 2724DEST_PATH_IMAGE006
、
Figure 782199DEST_PATH_IMAGE007
和
Figure 424533DEST_PATH_IMAGE008
作为平面拟合的点，用于对图像
Figure 231952DEST_PATH_IMAGE001
的剖面拟合平面方程；
S5、拟合图像
Figure 934329DEST_PATH_IMAGE009
的剖面所在平面：用最小二乘法拟合步骤S4中获取的四个角点
Figure 702565DEST_PATH_IMAGE005
、
Figure 148589DEST_PATH_IMAGE010
、
Figure 76094DEST_PATH_IMAGE007
和
Figure 949372DEST_PATH_IMAGE008
坐标所在平面；
S6、计算MNI空间下的DLPFC中心点
Figure 939325DEST_PATH_IMAGE011
在步骤S5得到的平面上的投影点
Figure 189041DEST_PATH_IMAGE012
；
S7、计算得到导针：连接DLPFC中心点
Figure 971052DEST_PATH_IMAGE013
与投影点
Figure 15231DEST_PATH_IMAGE014
，即得导针。
优选地，前述步骤S1中，二值球体图像
Figure 492480DEST_PATH_IMAGE002
的半径为20mm。
更优选地，前述步骤S1中，得到靶区图像
Figure 280308DEST_PATH_IMAGE001
后，计算靶区图像
Figure 916825DEST_PATH_IMAGE001
的灰度值
Figure 131906DEST_PATH_IMAGE015
：
Figure 597916DEST_PATH_IMAGE016
，
其中，
Figure 986172DEST_PATH_IMAGE017
表示图像
Figure 680458DEST_PATH_IMAGE018
中的像素在
Figure 4123DEST_PATH_IMAGE019
位置上的灰度值，
Figure 518281DEST_PATH_IMAGE020
表示图像
Figure 710228DEST_PATH_IMAGE004
中的像素在
Figure 259021DEST_PATH_IMAGE021
位置上的灰度值，
Figure 753588DEST_PATH_IMAGE022
表示
Figure 551779DEST_PATH_IMAGE003
与
Figure 157204DEST_PATH_IMAGE023
对应像素灰度值的乘法运算。
更优选地，前述步骤S2包括如下子步骤：
S2.1、选择三维球体图像
Figure 560504DEST_PATH_IMAGE024
的任一
Figure 350605DEST_PATH_IMAGE025
层记为
Figure 839355DEST_PATH_IMAGE026
，则三维图像
Figure 12585DEST_PATH_IMAGE027
变成二维图像
Figure 270391DEST_PATH_IMAGE028
，其中
Figure 231394DEST_PATH_IMAGE029
表示图像中像素的位置；
利用水平、垂直差分算子对其图像灰度
Figure 207440DEST_PATH_IMAGE030
的每个像素进行滤波，在水平和垂直方向上计算梯度，分别求得
Figure 420247DEST_PATH_IMAGE031
和
Figure 532560DEST_PATH_IMAGE032
，
Figure 664464DEST_PATH_IMAGE031
和
Figure 799910DEST_PATH_IMAGE033
为图像灰度
Figure 878724DEST_PATH_IMAGE034
的偏导；
设滑动窗口的滑动变量为
Figure 907860DEST_PATH_IMAGE035
，则以目标像素为中心的矩形窗口
Figure 413928DEST_PATH_IMAGE036
在任意方向上的灰度变化量为
Figure 36670DEST_PATH_IMAGE037
，
Figure 802818DEST_PATH_IMAGE038
代表水平、垂直方向上的偏移量，定义如下：
Figure 479787DEST_PATH_IMAGE039
S2.2、获得互相关矩阵M：
Figure 825711DEST_PATH_IMAGE040
S2.3、对互相关矩阵M中的四个元素做高斯平滑滤波
Figure 308645DEST_PATH_IMAGE041
，对图像中的每一点构建新的互相关矩阵M’：
离散二维高斯滤波函数为：
Figure 922160DEST_PATH_IMAGE042
新的互相关矩阵为：
Figure 770030DEST_PATH_IMAGE043
，
S2.4、利用新的互相关矩阵构造角点响应函数R：
Figure 491999DEST_PATH_IMAGE044
，
其中，|-|表示矩阵的行列式，tr表示矩阵M’的迹，参数k一般取值0.04-0.06，利用R的取值来确定
Figure 450727DEST_PATH_IMAGE045
目标像素点的特征；
Figure 449907DEST_PATH_IMAGE046
和
Figure 203100DEST_PATH_IMAGE047
分别代表矩阵M’的最大特征值和最小特征值;
S2.5、对R进行非极大值抑制，同时满足R大于一阈值threshold且为某邻域内的局部极大值的点记为角点。
再优选地，前述步骤S2.5中，通过非极大值抑制将不是极大值的点去掉。
再优选地，前述步骤S2.5中，设定阈值R>0，则为角点；R<0则为边缘。
进一步优选地，前述步骤S3中，任意两角点的欧几里得距离
Figure 412364DEST_PATH_IMAGE048
的计算方法为：
Figure 174784DEST_PATH_IMAGE049
，
Figure 28470DEST_PATH_IMAGE050
和
Figure 14881DEST_PATH_IMAGE051
分别为两个角点在
Figure 383545DEST_PATH_IMAGE052
中的位置坐标；
包括如下子步骤：
S3.1设角点总数为N，未计算欧几里得距离的角点总数为N0，未计算过欧几里得距离的角点坐标数组为N[i]。
S3.2若N0≤0，则结束计算；若N0>0，则任取N[i]中一角点
Figure 120295DEST_PATH_IMAGE053
并将该点从N[i]中移除。
S3.3循环迭代N[i]数组，取N[i]中的点
Figure 156384DEST_PATH_IMAGE054
，其中n=i；
S3.4计算
Figure 579275DEST_PATH_IMAGE055
和
Figure 435236DEST_PATH_IMAGE054
之间的欧几里得距离；第一次计算的距离记为dmax，迭代计算的距离记为d[i]。
S3.5若d[i]≤dmax，则循环继续，重复步骤S3.3；若d[i]>dmax，则把d[i]赋值给dmax；
S3.6判断循环是否结束；是，则将dmax和坐标
Figure 477141DEST_PATH_IMAGE056
、
Figure 164474DEST_PATH_IMAGE057
存储到数组dmax[M-1]；否，则重复步骤S3.2-S3.5。
更进一步优选地，前述步骤S4包括如下子步骤：
S4.1设角点间的距离总数为Md，角点坐标和距离数组为dmax[Md]；
S4.2使用选择排序对最大距离进行排序，迭代数组dmax[Md]，取dmax[i]；
S4.3初始化dmax1的值为dmax[0]，dmax2的值为0；
S4.4若dmax[i] ＞dmax2，则把dmax[i]赋值给dmax2；若dmax[i]≤dmax2，则循环继续；
S4.5若dmax2＞dmax1，则交换这两个数的值；若dmax2≤dmax1，则循环继续。
S4.6判断循环是否结束；是，则输出距离排序在前两位的dmax1和dmax2及两个距离对应的坐标
Figure 695950DEST_PATH_IMAGE005
、
Figure 976890DEST_PATH_IMAGE006
、
Figure 619224DEST_PATH_IMAGE007
和
Figure 426643DEST_PATH_IMAGE008
作为平面拟合的点；否，则重复步骤S4.2-S4.5。
更进一步优选地，前述步骤S5采用最小二乘法拟合平面，具体拟合过程为：
设平面方程为：
Figure 129019DEST_PATH_IMAGE058
，
其中，x，y，z是图像
Figure 897255DEST_PATH_IMAGE059
的角点坐标，A、B、C是描述平面空间特征的常数，
使该平面到四个角点的距离最近，根据最小二乘法，其响应函数S为：
Figure 343280DEST_PATH_IMAGE060
计算出一组A、B、C，使得对四个角点来说S最小，对A、B、C分别求导，得到S最小时的A、B、C的值：
Figure 5206DEST_PATH_IMAGE061
求解方程组，即可得到四个角点所在平面方程的系数A、B、C，得到拟合平面。
更进一步优选地，前述步骤S6中，计算MNI空间下的DLPFC中心点
Figure 878484DEST_PATH_IMAGE011
在步骤S5平面上的投影点
Figure 635480DEST_PATH_IMAGE062
，两点形成的直线为
Figure 885196DEST_PATH_IMAGE063
，则
Figure 401628DEST_PATH_IMAGE063
的参数方程为：
Figure 445808DEST_PATH_IMAGE064
则投影点
Figure 188636DEST_PATH_IMAGE065
。
本发明的有益之处在于：
本发明的基于Harris角点检测的导航系统导针计算方法，采用结构磁共振成像技术，结合图像处理技术和数学模型准确地计算出导针位置和角度，计算过程中，利用二值图像运算、Harris角点检测、排序算法、平面拟合算法和投影点计算等多种方法，能够实现利用结构磁共振成像完成靶点计算，并将该方法与经颅磁刺激系统结合。最终得到了一种操作简便、快捷、自动化、定位精准、成本低廉的经颅磁定位用导航系统，成功解决了磁刺激定位技术问题，降低患者拍摄磁共振的成本，节约治疗师的治疗时间。
现有技术中，利用磁共振的导航技术需要结合结构磁共振成像和功能磁共振成像共同实现导航。通过功能像判断DLPFC的位置，结构像（更加清晰）进行查看，功能像可以计算DLPFC和SGACC间的负相关性，两者相辅相成。但，本发明的方法通过创新后，仅采用结构像制作DLPFC靶区就可以实现快速、精准的导航定位。在计算过程中，为了拟合半球剖面的方程，首先提取角点，然后利用欧几里得距离获取位于剖面的角点坐标来进行平面拟合，确保平面拟合不会局限在半球中间位置，从而达到理想的拟合效果。
相比其他定位导航技术，本发明不需要附加视觉系统和机械臂系统即可达到同样的定位效果，将精度误差控制在毫米级甚至更高，且有效避免了机械设定的误差，系统自动进行靶点计算，整个治疗过程5min内可以输出靶点坐标。
附图说明
图1是本发明的基于Harris角点检测的随动定位导航系统导针生成流程图；
图2是DLPFC球体与脑掩膜叠加显示图；
图3是DLPFC球体与脑掩膜相乘结果显示图；
图4是本发明的实施例中所得到的DLPFC所有角点分布图；
图5是本发明的实施例中计算任意两个角点之间的欧几里得距离的流程图；
图6是本发明的实施例中剖面拟合的算法流程图；
图7是DLPFC剖面拟合示意图。
具体实施方式
以下结合附图和具体实施例对本发明作具体的介绍。
本发明的基于Harris角点检测的随动定位系统导针计算方法流程如图1所示，利用图像处理技术和数学模型准确地计算导针位置和角度。计算过程中，结合DLPFC（背外侧前额皮层）制作方法、二值图像运算、Harris角点检测、排序算法、平面拟合算法和投影点计算等具体过程，能够实现利用结构磁共振成像快速、准确、低成本地完成靶点定位和计算。
具体包括如下步骤：
第一步（S1）、构建DLPFC靶区图像。
利用既有的医学图像处理软件（包括但不限于FSL、SPM等），根据磁共振模板（MNI）空间下的坐标，生成二值球体图像
Figure 976463DEST_PATH_IMAGE066
，将
Figure 612981DEST_PATH_IMAGE067
与磁共振模板空间下的脑区图像掩膜
Figure 93641DEST_PATH_IMAGE004
相乘获取两者重合的区域，即为所需靶区图像
Figure 386082DEST_PATH_IMAGE001
。
在本实施方式中，二值图像的半径为20mm，实际应用中可根据需求灵活调整，不限于此固定值。磁共振模板（MNI）空间下的坐标（靶点坐标）可以任意选择，一般选择经验性的坐标或者论文里给出的坐标，如：(-46,45,38)、(38.1,59.66,34.15)等。
如上所述的二值图像是像素值只有0和1数值的图像，0代表黑色，1代表白色，图2和图3中所示的白色即代表图像区域
Figure 649704DEST_PATH_IMAGE001
，图3即为两幅二值图像
Figure 78411DEST_PATH_IMAGE067
和
Figure 792289DEST_PATH_IMAGE004
相乘后的结果，不重叠的区域是0x1=0，重叠区域是1x1=1。
接着，计算靶区图像
Figure 306447DEST_PATH_IMAGE001
的灰度值
Figure 373760DEST_PATH_IMAGE015
：
Figure 922553DEST_PATH_IMAGE068
，
其中，
Figure 541754DEST_PATH_IMAGE069
表示图像
Figure 543208DEST_PATH_IMAGE070
中的像素在
Figure 912747DEST_PATH_IMAGE071
位置上的灰度值，
Figure 316046DEST_PATH_IMAGE020
表示图像
Figure 106148DEST_PATH_IMAGE072
中的像素在
Figure 594898DEST_PATH_IMAGE071
位置上的灰度值，“
Figure 4014DEST_PATH_IMAGE073
” 表示
Figure 261820DEST_PATH_IMAGE074
与
Figure 488402DEST_PATH_IMAGE075
对应像素灰度值的乘法运算。
第二步（S2）、采用Harris角点检测方法对靶区图像
Figure 464448DEST_PATH_IMAGE001
进行检测，获取
Figure 411675DEST_PATH_IMAGE001
中的所有角点。
需要说明的是，DLPFC实质上是一个三维半球，在本实施例的角点检测过程中，逐层取DLPFC三维图像的每一层进行Harris检测，最后得到的是如图4所示的三维角点分布图，DLPFC的维度与MNI模板的维度相同。
该步骤具体包括如下子步骤：
S2.1、对于三维球体图像
Figure 789567DEST_PATH_IMAGE024
，选择任一
Figure 921471DEST_PATH_IMAGE076
层记为
Figure 384814DEST_PATH_IMAGE026
，则三维图像
Figure 135732DEST_PATH_IMAGE077
变成二维图像
Figure 368130DEST_PATH_IMAGE028
。
利用水平、垂直差分算子对其图像灰度
Figure 139777DEST_PATH_IMAGE078
的每个像素进行滤波，分别计算水平和垂直方向上的梯度，求得
Figure 887153DEST_PATH_IMAGE031
和
Figure 504079DEST_PATH_IMAGE032
，其中，
Figure 41851DEST_PATH_IMAGE031
和
Figure 718820DEST_PATH_IMAGE079
分别为图像灰度
Figure 953492DEST_PATH_IMAGE034
在水平和垂直方向的偏导。
具体计算过程为：设滑动窗口的滑动变量为
Figure 108530DEST_PATH_IMAGE080
，则以目标像素为中心的矩形窗口
Figure 253203DEST_PATH_IMAGE036
在任意方向上的灰度变化量为
Figure 101074DEST_PATH_IMAGE037
，
Figure 823042DEST_PATH_IMAGE038
代表水平、垂直方向上的偏移量，定义如下：
Figure 781771DEST_PATH_IMAGE039
。
S2.2、获得互相关矩阵M：
Figure 515372DEST_PATH_IMAGE040
S2.3、对互相关矩阵M中的四个元素做高斯平滑滤波
Figure 534143DEST_PATH_IMAGE041
，对图像中的每一点构建新的互相关矩阵M’。
离散二维高斯滤波函数为：
Figure 743408DEST_PATH_IMAGE081
其中，
Figure 505827DEST_PATH_IMAGE082
决定了高斯函数的宽度，其取值决定了高斯函数窗口的大小，在实际运算中常取值0.8或者1。
新的互相关矩阵为：
Figure 359514DEST_PATH_IMAGE083
每个像素求偏导后都有
Figure 283607DEST_PATH_IMAGE031
和
Figure 980168DEST_PATH_IMAGE079
，
Figure 280699DEST_PATH_IMAGE084
则表示一张图像中多个像素进行高斯滤波后，求和的结果。
S2.4、利用新的互相关矩阵构造角点响应函数R：
Figure 753007DEST_PATH_IMAGE085
其中，|-|表示矩阵的行列式，tr表示矩阵M’的迹，参数k一般取值0.04-0.06，利用R的取值来确定
Figure 848002DEST_PATH_IMAGE086
目标像素点的特征。
Figure 31858DEST_PATH_IMAGE046
和
Figure 136081DEST_PATH_IMAGE087
分别代表矩阵M’的最大特征值和最小特征值。
S2.5、对R进行非极大值抑制，将不是极大值的点去掉，即将不是角点的值去掉。本实施例中，R<0，则认为不是角点，同时满足R大于一阈值且R为某邻域内的局部极大值的点记为角点。
具体到本实施例中，我们设定阈值R>0，则为角点；R<0则为边缘。
第三步（S3）、计算任意两个角点之间的欧几里得距离，任意两角点的欧几里得距离
Figure 698780DEST_PATH_IMAGE048
的计算方法如下：
Figure 230255DEST_PATH_IMAGE088
。
参见图5，包括如下子步骤：
S3.1设角点总数为N，未计算欧几里得距离的角点总数为N0，未计算过距离的角点坐标数组为N[i]。
S3.2若N0≤0，则结束计算；若N0>0，则任取N[i]中一角点
Figure 635829DEST_PATH_IMAGE089
并将该点从N[i]中移除。
S3.3循环迭代N[i]数组，取N[i]中的点
Figure 543742DEST_PATH_IMAGE090
，其中n=i，
Figure 226527DEST_PATH_IMAGE050
和
Figure 928904DEST_PATH_IMAGE091
为两个角点在
Figure 821774DEST_PATH_IMAGE052
中的位置坐标。
S3.4计算
Figure 533378DEST_PATH_IMAGE089
和
Figure 132986DEST_PATH_IMAGE090
之间的欧几里得距离。第一次计算的距离记为dmax，迭代计算的距离记为d[i]。
S3.5若d[i]≤dmax，则循环继续，重复步骤S3.3；若d[i]>dmax，则将d[i]赋值给dmax。
S3.6判断循环是否结束：是，则将dmax和坐标
Figure 943948DEST_PATH_IMAGE056
、
Figure 261797DEST_PATH_IMAGE092
存储到数组dmax[M-1]；否，则重复步骤S3.2-S3.5。
第四步（S4）、采用排序算法取距离排序在前两位的角点坐标，用于对图像
Figure 308250DEST_PATH_IMAGE093
的剖面拟合平面方程。
参见图6，具体的算法过程如下：
S4.1设角点间的距离总数为Md，角点坐标和距离数组为dmax[Md]。
S4.2使用选择排序对最大距离进行排序，迭代数组dmax[Md]，取dmax[i]。
S4.3初始化dmax1的值为dmax[0]，dmax2的值为0。
S4.4若dmax[i]>dmax2，则将dmax[i]赋值给dmax2；若dmax[i]≤dmax2，则循环继续。
S4.5若dmax2> dmax1，则交换这两个数的值；若dmax2≤ dmax1，则循环继续。
S4.6判断循环是否结束：是，则输出距离排序在前两位的dmax1和dmax2及两个距离对应的坐标
Figure 27944DEST_PATH_IMAGE005
、
Figure 511272DEST_PATH_IMAGE006
、
Figure 316417DEST_PATH_IMAGE007
和
Figure 166561DEST_PATH_IMAGE008
作为平面拟合的点；否，则重复步骤S4.2-S4.5。
第五步（S5）、采用最小二乘法拟合平面。
用最小二乘法拟合步骤S4中获取的四个角点
Figure 740762DEST_PATH_IMAGE005
、
Figure 893525DEST_PATH_IMAGE006
、
Figure 920387DEST_PATH_IMAGE007
和
Figure 574223DEST_PATH_IMAGE008
坐标所在平面，即为图像
Figure 2930DEST_PATH_IMAGE009
的剖面，设平面方程为：
Figure 654491DEST_PATH_IMAGE058
其中，x，y，z是图像
Figure 371911DEST_PATH_IMAGE059
的角点坐标，A、B、C是描述平面空间特征的常数。
使该平面到四个角点的距离最近，根据最小二乘法，其响应函数S为：
Figure 298279DEST_PATH_IMAGE060
计算出一组A、B、C，使得对四个角点来说S最小，对A、B、C分别求导，得到S最小时的A、B、C的值：
Figure 581493DEST_PATH_IMAGE061
其中，
Figure 341638DEST_PATH_IMAGE094
、
Figure 343093DEST_PATH_IMAGE095
、
Figure 338730DEST_PATH_IMAGE096
为平面中任一点的三维坐标，求解方程组，即可得到四个角点所在平面方程的系数A、B、C，得到的DLPFC剖面拟合示意图如图7所示。
第六步（S6）、计算MNI空间下的DLPFC中心点
Figure 476451DEST_PATH_IMAGE097
在步骤S5平面上的投影点
Figure 906033DEST_PATH_IMAGE062
，两点形成的直线为
Figure 394783DEST_PATH_IMAGE063
，则
Figure 928532DEST_PATH_IMAGE063
的参数方程为：
Figure 186338DEST_PATH_IMAGE098
则
Figure 288287DEST_PATH_IMAGE099
。
第七步（S7）、连接DLPFC中心点
Figure 998754DEST_PATH_IMAGE100
与投影点
Figure 336194DEST_PATH_IMAGE101
即可得到导针。
综上，本发明的导航系统导针计算方法是基于Harris角点检测而实现的，利用二值图像运算、Harris角点检测、排序算法、平面拟合算法和投影点计算等多种方法，实现利用结构磁共振成像完成靶点计算，并将该方法与经颅磁刺激系统结合。相比其他定位导航技术，本发明不需要附加视觉系统和机械臂系统即可达到同样的定位效果，将精度误差控制在毫米级甚至更高，且有效避免了机械设定的误差，系统自动进行靶点计算，整个治疗过程5min内可以输出靶点坐标。
以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，上述实施例不以任何形式限制本发明，凡采用等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。

Claims (9)

1.基于Harris角点检测的随动定位系统导针计算方法，其特征在于，包括如下步骤：
S1、构建DLPFC靶区图像
Figure DEST_PATH_IMAGE001
；
根据磁共振模板空间下的坐标生成二值球体图像
Figure DEST_PATH_IMAGE002
，将
Figure DEST_PATH_IMAGE003
与磁共振模板空间下的脑区图像掩膜
Figure DEST_PATH_IMAGE004
相乘，获取两者重合的区域，即为所需靶区图像
Figure 676083DEST_PATH_IMAGE001
；
S2、获取
Figure 327644DEST_PATH_IMAGE001
中的所有角点；
逐层取DLPFC三维靶区图像
Figure 904119DEST_PATH_IMAGE001
的每一层进行Harris检测，得到三维的角点分布图；
S3、求取任意两个角点之间的欧几里得距离；
S4、取距离排序在前两位的角点坐标，输出距离排序在前两位的dmax1、dmax2及两个距离对应的坐标
Figure DEST_PATH_IMAGE005
、
Figure DEST_PATH_IMAGE006
、
Figure DEST_PATH_IMAGE007
和
Figure DEST_PATH_IMAGE008
作为平面拟合的点，用于对图像
Figure 96066DEST_PATH_IMAGE001
的剖面拟合平面方程；
S5、拟合图像
Figure DEST_PATH_IMAGE009
的剖面所在平面：用最小二乘法拟合步骤S4中获取的四个角点
Figure 208641DEST_PATH_IMAGE005
、
Figure 765524DEST_PATH_IMAGE006
、
Figure 829295DEST_PATH_IMAGE007
和
Figure 762616DEST_PATH_IMAGE008
坐标所在平面；
具体拟合过程为：
设平面方程为：
Figure DEST_PATH_IMAGE010
其中，x，y，z是图像
Figure DEST_PATH_IMAGE011
的角点坐标，A、B、C是描述平面空间特征的常数，
使该平面到四个角点的距离最近，根据最小二乘法，其响应函数S为：
Figure DEST_PATH_IMAGE012
计算出一组A、B、C，使得对四个角点来说S最小，对A、B、C分别求导，得到S最小时的A、B、C的值：
Figure DEST_PATH_IMAGE013
求解方程组，即可得到四个角点所在平面方程的系数A、B、C，得到拟合平面；
S6、计算MNI空间下的DLPFC中心点
Figure DEST_PATH_IMAGE014
在步骤S5得到的平面上的投影点
Figure DEST_PATH_IMAGE015
；
S7、计算得到导针：连接DLPFC中心点
Figure DEST_PATH_IMAGE016
与投影点
Figure DEST_PATH_IMAGE017
，即得导针。
2.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S1中，二值球体图像
Figure 556128DEST_PATH_IMAGE002
的半径为20mm。
3.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S1中，得到靶区图像
Figure 346230DEST_PATH_IMAGE001
后，计算靶区图像
Figure 834980DEST_PATH_IMAGE001
的灰度值
Figure DEST_PATH_IMAGE018
：
Figure DEST_PATH_IMAGE019
，
其中，
Figure DEST_PATH_IMAGE020
表示图像
Figure DEST_PATH_IMAGE021
中的像素在
Figure DEST_PATH_IMAGE022
位置上的灰度值，
Figure DEST_PATH_IMAGE023
表示图像
Figure DEST_PATH_IMAGE024
中的像素在
Figure DEST_PATH_IMAGE025
位置上的灰度值，
Figure DEST_PATH_IMAGE026
表示
Figure 824906DEST_PATH_IMAGE021
与
Figure DEST_PATH_IMAGE027
对应像素灰度值的乘法运算。
4.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S2包括如下子步骤：
S2.1、选择三维球体图像
Figure DEST_PATH_IMAGE028
的任一
Figure DEST_PATH_IMAGE029
层记为
Figure DEST_PATH_IMAGE030
，则三维图像
Figure DEST_PATH_IMAGE031
变成二维图像
Figure DEST_PATH_IMAGE032
，其中
Figure DEST_PATH_IMAGE033
表示图像中像素的位置；
利用水平、垂直差分算子对其图像灰度
Figure DEST_PATH_IMAGE034
的每个像素进行滤波，在水平和垂直方向上计算梯度，求得
Figure DEST_PATH_IMAGE035
和
Figure DEST_PATH_IMAGE036
，
Figure 207346DEST_PATH_IMAGE035
和
Figure DEST_PATH_IMAGE037
为图像灰度
Figure DEST_PATH_IMAGE038
的偏导；
设滑动窗口的滑动变量为
Figure DEST_PATH_IMAGE039
，则以目标像素为中心的矩形窗口
Figure DEST_PATH_IMAGE040
在任意方向上的灰度变化量为
Figure DEST_PATH_IMAGE041
，
Figure DEST_PATH_IMAGE042
代表水平、垂直方向上的偏移量，定义如下：
Figure DEST_PATH_IMAGE043
S2.2、获得互相关矩阵M：
Figure DEST_PATH_IMAGE044
S2.3、对互相关矩阵M中的四个元素做高斯平滑滤波
Figure DEST_PATH_IMAGE045
，对图像中的每一点构建新的互相关矩阵M’：
离散二维高斯滤波函数为：
Figure DEST_PATH_IMAGE046
，其中，
Figure 591185DEST_PATH_IMAGE047
决定高斯函数的宽度，取值0.8或者1；
新的互相关矩阵为：
Figure DEST_PATH_IMAGE048
，
S2.4、利用新的互相关矩阵构造角点响应函数R：
Figure DEST_PATH_IMAGE049
，
其中，|-|表示矩阵的行列式，tr表示矩阵M’的迹，参数k取值0.04-0.06，利用R的取值来确定
Figure DEST_PATH_IMAGE050
目标像素点的特征；
Figure DEST_PATH_IMAGE051
和
Figure DEST_PATH_IMAGE052
分别代表矩阵M’的最大特征值和最小特征值;
S2.5、对R进行非极大值抑制，同时满足R大于一阈值threshold且为某邻域内的局部极大值的点记为角点。
5.根据权利要求4所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S2.5中，通过非极大值抑制将不是极大值的点去掉。
6.根据权利要求4所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S2.5中，设定阈值R>0，则为角点；R<0则为边缘。
7.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S3中，任意两角点的欧几里得距离
Figure DEST_PATH_IMAGE053
的计算方法为：
Figure DEST_PATH_IMAGE054
，
Figure DEST_PATH_IMAGE055
和
Figure DEST_PATH_IMAGE056
分别为两个角点在
Figure DEST_PATH_IMAGE057
中的位置坐标；
包括如下子步骤：
S3.1设角点总数为N，未计算欧几里得距离的角点总数为N0，未计算过欧几里得距离的角点坐标数组为N[i]；
S3.2若N0≤0，则结束计算；若N0>0，则任取N[i]中一角点
Figure DEST_PATH_IMAGE058
并将该点从N[i]中移除；
S3.3循环迭代N[i]数组，取N[i]中的点
Figure DEST_PATH_IMAGE059
，其中n=i；
S3.4计算
Figure 724488DEST_PATH_IMAGE058
和
Figure 796350DEST_PATH_IMAGE059
之间的欧几里得距离；第一次计算的距离记为dmax，迭代计算的距离记为d[i]；
S3.5若d[i]≤dmax，则循环继续，重复步骤S3.3；若d[i]>dmax，则把d[i]赋值给dmax；
S3.6判断循环是否结束；是，则将dmax和坐标
Figure DEST_PATH_IMAGE060
、
Figure DEST_PATH_IMAGE061
存储到数组dmax[M-1]；否，则重复步骤S3.2-S3.5。
8.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S4包括如下子步骤：
S4.1设角点间的距离总数为Md，角点坐标和距离数组为dmax[Md]；
S4.2使用选择排序对最大距离进行排序，迭代数组dmax[Md]，取dmax[i]；
S4.3初始化dmax1的值为dmax[0]，dmax2的值为0；
S4.4若dmax[i] ＞dmax2，则把dmax[i]赋值给dmax2；若dmax[i]≤dmax2，则循环继续；
S4.5若dmax2＞dmax1，则交换这两个数的值；若dmax2≤dmax1，则循环继续；
S4.6判断循环是否结束；是，则输出距离排序在前两位的dmax1和dmax2及两个距离对应的坐标
Figure 705400DEST_PATH_IMAGE005
、
Figure 837304DEST_PATH_IMAGE006
、
Figure 566225DEST_PATH_IMAGE007
和
Figure 176198DEST_PATH_IMAGE008
作为平面拟合的点；否，则重复步骤S4.2-S4.5。
9.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S6中，计算MNI空间下的DLPFC中心点
Figure 674176DEST_PATH_IMAGE014
在步骤S5平面上的投影点
Figure DEST_PATH_IMAGE062
，两点形成的直线为
Figure DEST_PATH_IMAGE063
，则
Figure 711402DEST_PATH_IMAGE063
的参数方程为：
Figure DEST_PATH_IMAGE064
，
则投影点
Figure DEST_PATH_IMAGE065
。
CN202210971975.3A 2022-08-15 2022-08-15 基于Harris角点检测的随动定位系统导针计算方法 Active CN115035124B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210971975.3A CN115035124B (zh) 2022-08-15 2022-08-15 基于Harris角点检测的随动定位系统导针计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210971975.3A CN115035124B (zh) 2022-08-15 2022-08-15 基于Harris角点检测的随动定位系统导针计算方法

Publications (2)

Publication Number Publication Date
CN115035124A CN115035124A (zh) 2022-09-09
CN115035124B true CN115035124B (zh) 2022-11-11

Family

ID=83129980

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210971975.3A Active CN115035124B (zh) 2022-08-15 2022-08-15 基于Harris角点检测的随动定位系统导针计算方法

Country Status (1)

Country Link
CN (1) CN115035124B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116543042B (zh) * 2023-05-10 2024-01-23 中国科学院心理研究所 基于组水平平均统计图的抑郁症tms个体化靶点定位方法及系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002078712A (ja) * 2000-09-05 2002-03-19 Koseki Ika Kk 胸骨縫合用ワイヤーの製造方法、材質、拘束帯、包装体及び錐
CN110857319A (zh) * 2018-08-24 2020-03-03 杭州康万达医药科技有限公司 一种分离的t细胞受体、其修饰的细胞、编码核酸及其应用
CN210250096U (zh) * 2019-05-05 2020-04-07 钦州市第二人民医院 一种半自动型活检针简易定位尺
CN113893036A (zh) * 2021-09-09 2022-01-07 上海交通大学 一种磁共振环境下的介入机器人装置

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8811699B2 (en) * 2010-09-22 2014-08-19 Siemens Aktiengesellschaft Detection of landmarks and key-frames in cardiac perfusion MRI using a joint spatial-temporal context model
US11071596B2 (en) * 2016-08-16 2021-07-27 Insight Medical Systems, Inc. Systems and methods for sensory augmentation in medical procedures
CN113538384B (zh) * 2021-07-19 2024-03-26 凌云光技术股份有限公司 一种特征的定位方法及装置

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002078712A (ja) * 2000-09-05 2002-03-19 Koseki Ika Kk 胸骨縫合用ワイヤーの製造方法、材質、拘束帯、包装体及び錐
CN110857319A (zh) * 2018-08-24 2020-03-03 杭州康万达医药科技有限公司 一种分离的t细胞受体、其修饰的细胞、编码核酸及其应用
CN210250096U (zh) * 2019-05-05 2020-04-07 钦州市第二人民医院 一种半自动型活检针简易定位尺
CN113893036A (zh) * 2021-09-09 2022-01-07 上海交通大学 一种磁共振环境下的介入机器人装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
P5C-5 Design and Validation of an Ultrasound Array Optimised for Epidural Needle Guidance;S. Cochran et al.;《2007 IEEE Ultrasonics Symposium Proceedings》;20071226;2255-2258 *
基于视觉的小型电机定子绕线机器人轨迹规划仿真研究;王林强;《中国优秀硕士学位论文全文数据库 (工程科技Ⅱ辑)》;20210915;C042-83 *

Also Published As

Publication number Publication date
CN115035124A (zh) 2022-09-09

Similar Documents

Publication Publication Date Title
CN110517202B (zh) 一种车身摄像头标定方法及其标定装置
CN106920234B (zh) 一种复合式自动放疗计划的方法
CN109272443B (zh) 一种基于全卷积神经网络的pet与ct图像配准方法
CN109598762A (zh) 一种高精度双目相机标定方法
CN106485695B (zh) 基于统计形状模型的医学图像Graph Cut分割方法
Mattes et al. Nonrigid multimodality image registration
Chanwimaluang et al. Hybrid retinal image registration
EP2919194B1 (en) Image data processing device and transcranial magnetic stimulation apparatus
US9466093B2 (en) Anterior commissure and posterior commissure segmentation system and method
US20160203604A1 (en) Method for automatic detection of anatomical landmarks in volumetric data
US20070053491A1 (en) Adaptive radiation therapy method with target detection
WO2014144019A1 (en) Methods, systems, and computer readable media for real-time 2d/3d deformable registration using metric learning
CN101258525A (zh) 用于对心脏右心室进行可靠的3d评价的超声系统及其方法
CN115035124B (zh) 基于Harris角点检测的随动定位系统导针计算方法
Tong et al. X-ray2Shape: reconstruction of 3D liver shape from a single 2D projection image
CN112348890B (zh) 一种空间定位方法、装置及计算机可读存储介质
CN113610887A (zh) 胶囊内窥镜运动拍摄路径的确定方法、存储介质和设备
CN106778660A (zh) 一种人脸姿态校正方法及装置
CN108549906A (zh) 放疗勾靶图像配准方法及装置
CN117094917A (zh) 一种心血管3d打印数据处理方法
US9147250B2 (en) System and method for automatic magnetic resonance volume composition and normalization
CN105389476B (zh) 基于梯度特征的调强放射治疗计划剂量数据的插值算法
Yao et al. Registrating oblique SAR images based on complementary integrated filtering and multilevel matching
CN112837225B (zh) 一种立位全脊柱图像自动无缝拼接的方法及装置
CN109785340A (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