本发明公布了一种基于Harris角点检测的随动定位系统导针计算方法，该方法仅采用结构像制作DLPFC靶区就可以实现快速、精准的导航定位，得到一种操作简便、快捷、自动化、定位精准、成本低廉的经颅磁定位用导航系统，成功解决了磁刺激定位技术问题，降低患者拍摄磁共振的成本，节约治疗师的治疗时间。相比其他定位导航技术，本发明不需要附加视觉系统和机械臂系统即可达到同样的定位效果，将精度误差控制在毫米级甚至更高，且有效避免了机械设定的误差，系统自动进行靶点计算，整个治疗过程5min内可以输出靶点坐标。
1.基于Harris角点检测的随动定位系统导针计算方法，其特征在于，包括如下步骤：S1、构建DLPFC靶区图像 ；
根据磁共振模板空间下的坐标生成二值球体图像 ，将 与磁共振模板空间下的脑区图像掩膜 相乘，获取两者重合的区域，即为所需靶区图像 ；
逐层取DLPFC三维靶区图像 的每一层进行Harris检测，得到三维的角点分布图；
S4、取距离排序在前两位的角点坐标，输出距离排序在前两位的dmax1、dmax2及两个距离对应的坐标 、 、 和 作为平面拟合的点，用于对图像 的剖面拟合平面方程；
S5、拟合图像 的剖面所在平面：用最小二乘法拟合步骤S4中获取的四个角点、 、 和 坐标所在平面；
其中，x，y，z是图像 的角点坐标，A、B、C是描述平面空间特征的常数，使该平面到四个角点的距离最近，根据最小二乘法，其响应函数S为：计算出一组A、B、C，使得对四个角点来说S最小，对A、B、C分别求导，得到S最小时的A、B、C的值：求解方程组，即可得到四个角点所在平面方程的系数A、B、C，得到拟合平面；
S6、计算MNI空间下的DLPFC中心点 在步骤S5得到的平面上的投影点；
S7、计算得到导针：连接DLPFC中心点 与投影点 ，即得导针。
2.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S1中，二值球体图像 的半径为20mm。
3.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S1中，得到靶区图像 后，计算靶区图像 的灰度值 ：，
其中， 表示图像 中的像素在 位置上的灰度值， 表示图像 中的像素在 位置上的灰度值， 表示 与 对应像素灰度值的乘法运算。
4.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S2包括如下子步骤：S2.1、选择三维球体图像 的任一 层记为 ，则三维图像 变成二维图像 ，其中 表示图像中像素的位置；
利用水平、垂直差分算子对其图像灰度 的每个像素进行滤波，在水平和垂直方向上计算梯度，求得 和 ， 和 为图像灰度 的偏导；
设滑动窗口的滑动变量为 ，则以目标像素为中心的矩形窗口 在任意方向上的灰度变化量为 ， 代表水平、垂直方向上的偏移量，定义如下：S2.2、获得互相关矩阵M：
S2.3、对互相关矩阵M中的四个元素做高斯平滑滤波 ，对图像中的每一点构建新的互相关矩阵M’：离散二维高斯滤波函数为： ，其中，决定高斯函数的宽度，取值0.8或者1；
S2.4、利用新的互相关矩阵构造角点响应函数R：，
其中，|‑|表示矩阵的行列式，tr表示矩阵M’的迹，参数k取值0.04‑0.06，利用R的取值来确定 目标像素点的特征； 和 分别代表矩阵M’的最大特征值和最小特征值;S2.5、对R进行非极大值抑制，同时满足R大于一阈值threshold且为某邻域内的局部极大值的点记为角点。
5.根据权利要求4所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S2.5中，通过非极大值抑制将不是极大值的点去掉。
6.根据权利要求4所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S2.5中，设定阈值R>0，则为角点；R<0则为边缘。
7.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S3中，任意两角点的欧几里得距离 的计算方法为：， 和 分别为两个角点在
S3.1设角点总数为N，未计算欧几里得距离的角点总数为N0，未计算过欧几里得距离的角点坐标数组为N[i]；
S3.2若N0≤0，则结束计算；若N0>0，则任取N[i]中一角点 并将该点从N[i]中移除；
S3.3循环迭代N[i]数组，取N[i]中的点 ，其中n=i；
S3.4计算 和 之间的欧几里得距离；第一次计算的距离记为dmax，迭代计算的距离记为d[i]；
S3.5若d[i]≤dmax，则循环继续，重复步骤S3.3；若d[i]>dmax，则把d[i]赋值给dmax；
S3.6判断循环是否结束；是，则将dmax和坐标 、 存储到数组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及两个距离对应的坐标 、 、 和 作为平面拟合的点；否，则重复步骤S4.2‑S4.5。
9.根据权利要求1所述的基于Harris角点检测的随动定位系统导针计算方法，其特征在于，所述步骤S6中，计算MNI空间下的DLPFC中心点 在步骤S5平面上的投影点 ，两点形成的直线为 ，则 的参数方程为：，
基于Harris角点检测的随动定位系统导针计算方法
技术领域
[0001] 本发明涉及一种导针计算方法，具体涉及一种基于Harris角点检测的随动定位导航系统的导针计算方法，属于图像处理技术领域。
背景技术
[0002] 经颅磁刺激可用于治疗抑郁症、偏瘫、失语等精神类疾病，在进行脑刺激治疗时，如何定位并找到头部正确的靶点或病灶是至关重要的，一般需要经验丰富的临床医生根据经验进行人工定位。具体操作手法如下：临床医生根据经验用拍头刺激患者头部，然后在设备上查看脉冲波，如波形不正确，则重复刺激，直到找到正确的靶点位置，但该方法操作十分繁琐，耗时费力，耽误治疗时间，且定位不精确，会增加患者的痛苦，在临床应用中十分受限。
[0003] 随着技术的发展，除传统的人工定位方法外，用于治疗的经颅磁刺激系统能够依赖于视觉相机和机械臂实现定位导航，但是该导航架结构复杂、价格昂贵、不易普及。开发出的新型随动定位系统，对于功能磁共振成像技术而言，只需计算体素时间序列间的相关性即可确定磁刺激靶点，但功能磁共振成像采集方式复杂、对患者要求较高，所以其应用推广受到了一定程度的限制。结构磁共振成像采集方式具有简单的优点，但结构磁共振成像无时间序列，因此需要手动制作靶区来生成靶点，耗时费力且过度依赖于医生的技术。
[0004] 鉴于上述原因，如何优化经颅定位用导航系统，解决结构磁共振成像技术中的磁刺激定位技术问题，是亟待解决的行业难题。
发明内容
[0005] 为解决现有技术的不足，本发明的目的在于提供一种适用于随动定位导航系统的导针计算方法，以实现简便、快捷、准确、低成本地解决结构磁共振成像技术中的磁刺激定位问题。
[0006] 为了实现上述目标，本发明采用如下的技术方案：
[0007] 基于Harris角点检测的随动定位系统导针计算方法，包括如下步骤：
[0008] S1、构建DLPFC靶区图像 ；
[0009] 根据磁共振模板空间下的坐标生成二值球体图像 ，将 与磁共振模板空间下的脑区图像掩膜 相乘，获取两者重合的区域，即为所需靶区图像 ；
[0010] S2、获取 中的所有角点；
[0011] 逐层取DLPFC三维靶区图像 的每一层进行Harris检测，得到三维的角点分布图；
[0012] S3、求取任意两个角点之间的欧几里得距离；
[0013] S4、取距离排序在前两位的角点坐标，输出距离排序在前两位的dmax1、dmax2及两个距离对应的坐标 、 、 和 作为平面拟合的点，用于对图像 的剖面拟合平面方程；
[0014] S5、拟合图像 的剖面所在平面：用最小二乘法拟合步骤S4中获取的四个角点、 、 和 坐标所在平面；
[0015] S6、计算MNI空间下的DLPFC中心点 在步骤S5得到的平面上的投影点 ；
[0016] S7、计算得到导针：连接DLPFC中心点 与投影点 ，即得导针。
[0017] 优选地，前述步骤S1中，二值球体图像 的半径为20mm。
[0018] 更优选地，前述步骤S1中，得到靶区图像 后，计算靶区图像 的灰度值：
[0019] ，
[0020] 其中， 表示图像 中的像素在 位置上的灰度值，表示图像 中的像素在 位置上的灰度值， 表示 与 对应像素灰度值的
乘法运算。
[0021] 更优选地，前述步骤S2包括如下子步骤：
[0022] S2.1、选择三维球体图像 的任一 层记为 ，则三维图像 变成二维图像 ，其中 表示图像中像素的位置；
[0023] 利用水平、垂直差分算子对其图像灰度 的每个像素进行滤波，在水平和垂直方向上计算梯度，分别求得 和 ， 和 为图像灰度 的偏导；
[0024] 设滑动窗口的滑动变量为 ，则以目标像素为中心的矩形窗口 在任意方向上的灰度变化量为 ， 代表水平、垂直方向上的偏移量，定义如下：
[0025]
[0026] S2.2、获得互相关矩阵M：
[0027]
[0028] S2.3、对互相关矩阵M中的四个元素做高斯平滑滤波 ，对图像中的每一点构建新的互相关矩阵M’：
[0029] 离散二维高斯滤波函数为：
[0030] 新的互相关矩阵为： ，
[0031] S2.4、利用新的互相关矩阵构造角点响应函数R：
[0032] ，
[0033] 其中，|‑|表示矩阵的行列式，tr表示矩阵M’的迹，参数k一般取值0.04‑0.06，利用R的取值来确定 目标像素点的特征； 和 分别代表矩阵M’的最大特征值和最小特征值;
[0034] S2.5、对R进行非极大值抑制，同时满足R大于一阈值threshold且为某邻域内的局部极大值的点记为角点。
[0035] 再优选地，前述步骤S2.5中，通过非极大值抑制将不是极大值的点去掉。
[0036] 再优选地，前述步骤S2.5中，设定阈值R>0，则为角点；R<0则为边缘。
[0037] 进一步优选地，前述步骤S3中，任意两角点的欧几里得距离 的计算方法为：
[0038] ， 和 分别为两个角点在 中的位置坐标；
[0039] 包括如下子步骤：
[0040] S3.1设角点总数为N，未计算欧几里得距离的角点总数为N0，未计算过欧几里得距离的角点坐标数组为N[i]。
[0041] S3.2若N0≤0，则结束计算；若N0>0，则任取N[i]中一角点 并将该点从N[i]中移除。
[0042] S3.3循环迭代N[i]数组，取N[i]中的点 ，其中n=i；
[0043] S3.4计算 和 之间的欧几里得距离；第一次计算的距离记为dmax，迭代计算的距离记为d[i]。
[0044] S3.5若d[i]≤dmax，则循环继续，重复步骤S3.3；若d[i]>dmax，则把d[i]赋值给dmax；
[0045] S3.6判断循环是否结束；是，则将dmax和坐标 、 存储到数组dmax[M‑1]；否，则重复步骤S3.2‑S3.5。
[0046] 更进一步优选地，前述步骤S4包括如下子步骤：
[0047] S4.1设角点间的距离总数为Md，角点坐标和距离数组为dmax[Md]；
[0048] S4.2使用选择排序对最大距离进行排序，迭代数组dmax[Md]，取dmax[i]；
[0049] S4.3初始化dmax1的值为dmax[0]，dmax2的值为0；
[0050] S4.4若dmax[i] ＞dmax2，则把dmax[i]赋值给dmax2；若dmax[i]≤dmax2，则循环继续；
[0051] S4.5若dmax2＞dmax1，则交换这两个数的值；若dmax2≤dmax1，则循环继续。
[0052] S4.6判断循环是否结束；是，则输出距离排序在前两位的dmax1和dmax2及两个距离对应的坐标 、 、 和 作为平面拟合的点；否，则重复步骤S4.2‑S4.5。
[0053] 更进一步优选地，前述步骤S5采用最小二乘法拟合平面，具体拟合过程为：
[0054] 设平面方程为：
[0055] ，
[0056] 其中，x，y，z是图像 的角点坐标，A、B、C是描述平面空间特征的常数，[0057] 使该平面到四个角点的距离最近，根据最小二乘法，其响应函数S为：
[0058]
[0059] 计算出一组A、B、C，使得对四个角点来说S最小，对A、B、C分别求导，得到S最小时的A、B、C的值：
[0060]
[0061] 求解方程组，即可得到四个角点所在平面方程的系数A、B、C，得到拟合平面。
[0062] 更进一步优选地，前述步骤S6中，计算MNI空间下的DLPFC中心点在步骤S5平面上的投影点 ，两点形成的直线为 ，则 的参数方程为：
[0063]
[0064] 则投影点 。
[0065] 本发明的有益之处在于：
[0066] 本发明的基于Harris角点检测的导航系统导针计算方法，采用结构磁共振成像技术，结合图像处理技术和数学模型准确地计算出导针位置和角度，计算过程中，利用二值图像运算、Harris角点检测、排序算法、平面拟合算法和投影点计算等多种方法，能够实现利用结构磁共振成像完成靶点计算，并将该方法与经颅磁刺激系统结合。最终得到了一种操作简便、快捷、自动化、定位精准、成本低廉的经颅磁定位用导航系统，成功解决了磁刺激定位技术问题，降低患者拍摄磁共振的成本，节约治疗师的治疗时间。
[0067] 现有技术中，利用磁共振的导航技术需要结合结构磁共振成像和功能磁共振成像共同实现导航。通过功能像判断DLPFC的位置，结构像（更加清晰）进行查看，功能像可以计算DLPFC和SGACC间的负相关性，两者相辅相成。但，本发明的方法通过创新后，仅采用结构像制作DLPFC靶区就可以实现快速、精准的导航定位。在计算过程中，为了拟合半球剖面的方程，首先提取角点，然后利用欧几里得距离获取位于剖面的角点坐标来进行平面拟合，确保平面拟合不会局限在半球中间位置，从而达到理想的拟合效果。
[0068] 相比其他定位导航技术，本发明不需要附加视觉系统和机械臂系统即可达到同样的定位效果，将精度误差控制在毫米级甚至更高，且有效避免了机械设定的误差，系统自动进行靶点计算，整个治疗过程5min内可以输出靶点坐标。
附图说明
[0069] 图1是本发明的基于Harris角点检测的随动定位导航系统导针生成流程图；
[0070] 图2是DLPFC球体与脑掩膜叠加显示图；
[0071] 图3是DLPFC球体与脑掩膜相乘结果显示图；
[0072] 图4是本发明的实施例中所得到的DLPFC所有角点分布图；
[0073] 图5是本发明的实施例中计算任意两个角点之间的欧几里得距离的流程图；
[0074] 图6是本发明的实施例中剖面拟合的算法流程图；
[0075] 图7是DLPFC剖面拟合示意图。
具体实施方式
[0076] 以下结合附图和具体实施例对本发明作具体的介绍。
[0077] 本发明的基于Harris角点检测的随动定位系统导针计算方法流程如图1所示，利用图像处理技术和数学模型准确地计算导针位置和角度。计算过程中，结合DLPFC（背外侧前额皮层）制作方法、二值图像运算、Harris角点检测、排序算法、平面拟合算法和投影点计算等具体过程，能够实现利用结构磁共振成像快速、准确、低成本地完成靶点定位和计算。
[0078] 具体包括如下步骤：
[0079] 第一步（S1）、构建DLPFC靶区图像。
[0080] 利用既有的医学图像处理软件（包括但不限于FSL、SPM等），根据磁共振模板（MNI）空间下的坐标，生成二值球体图像 ，将 与磁共振模板空间下的脑区图像掩膜 相乘获取两者重合的区域，即为所需靶区图像 。
[0081] 在本实施方式中，二值图像的半径为20mm，实际应用中可根据需求灵活调整，不限于此固定值。磁共振模板（MNI）空间下的坐标（靶点坐标）可以任意选择，一般选择经验性的坐标或者论文里给出的坐标，如：(‑46,45,38)、(38.1,59.66,34.15)等。
[0082] 如上所述的二值图像是像素值只有0和1数值的图像，0代表黑色，1代表白色，图2和图3中所示的白色即代表图像区域 ，图3即为两幅二值图像 和 相乘后的结果，不重叠的区域是0x1=0，重叠区域是1x1=1。
[0083] 接着，计算靶区图像 的灰度值 ：
[0084] ，
[0085] 其中， 表示图像 中的像素在 位置上的灰度值，表示图像 中的像素在 位置上的灰度值，“” 表示 与 对应像素灰度值
的乘法运算。
[0086] 第二步（S2）、采用Harris角点检测方法对靶区图像 进行检测，获取 中的所有角点。
[0087] 需要说明的是，DLPFC实质上是一个三维半球，在本实施例的角点检测过程中，逐层取DLPFC三维图像的每一层进行Harris检测，最后得到的是如图4所示的三维角点分布图，DLPFC的维度与MNI模板的维度相同。
[0088] 该步骤具体包括如下子步骤：
[0089] S2.1、对于三维球体图像 ，选择任一 层记为 ，则三维图像 变成二维图像 。
[0090] 利用水平、垂直差分算子对其图像灰度 的每个像素进行滤波，分别计算水平和垂直方向上的梯度，求得 和 ，其中， 和 分别为图像灰度 在水平和垂直方向的偏导。
[0091] 具体计算过程为：设滑动窗口的滑动变量为 ，则以目标像素为中心的矩形窗口 在任意方向上的灰度变化量为 ， 代表水平、垂直方向上的偏移量，定义如下：
[0092] 。
[0093] S2.2、获得互相关矩阵M：
[0094]
[0095] S2.3、对互相关矩阵M中的四个元素做高斯平滑滤波 ，对图像中的每一点构建新的互相关矩阵M’。
[0096] 离散二维高斯滤波函数为：
[0097] 其中， 决定了高斯函数的宽度，其取值决定了高斯函数窗口的大小，在实际运算中常取值0.8或者1。
[0098] 新的互相关矩阵为：
[0099] 每个像素求偏导后都有 和 ， 则表示一张图像中多个像素进行高斯滤波后，求和的结果。
[0100] S2.4、利用新的互相关矩阵构造角点响应函数R：
[0101]
[0102] 其中，|‑|表示矩阵的行列式，tr表示矩阵M’的迹，参数k一般取值0.04‑0.06，利用R的取值来确定 目标像素点的特征。 和 分别代表矩阵M’的最大特征值和最小特征值。
[0103] S2.5、对R进行非极大值抑制，将不是极大值的点去掉，即将不是角点的值去掉。本实施例中，R<0，则认为不是角点，同时满足R大于一阈值且R为某邻域内的局部极大值的点记为角点。
[0104] 具体到本实施例中，我们设定阈值R>0，则为角点；R<0则为边缘。
[0105] 第三步（S3）、计算任意两个角点之间的欧几里得距离，任意两角点的欧几里得距离 的计算方法如下：
[0106] 。
[0107] 参见图5，包括如下子步骤：
[0108] S3.1设角点总数为N，未计算欧几里得距离的角点总数为N0，未计算过距离的角点坐标数组为N[i]。
[0109] S3.2若N0≤0，则结束计算；若N0>0，则任取N[i]中一角点 并将该点从N[i]中移除。
[0110] S3.3循环迭代N[i]数组，取N[i]中的点 ，其中n=i， 和为两个角点在 中的位置坐标。
[0111] S3.4计算 和 之间的欧几里得距离。第一次计算的距离记为dmax，迭代计算的距离记为d[i]。
[0112] S3.5若d[i]≤dmax，则循环继续，重复步骤S3.3；若d[i]>dmax，则将d[i]赋值给dmax。
[0113] S3.6判断循环是否结束：是，则将dmax和坐标 、 存储到数组dmax[M‑1]；否，则重复步骤S3.2‑S3.5。
[0114] 第四步（S4）、采用排序算法取距离排序在前两位的角点坐标，用于对图像 的剖面拟合平面方程。
[0115] 参见图6，具体的算法过程如下：
[0116] S4.1设角点间的距离总数为Md，角点坐标和距离数组为dmax[Md]。
[0117] S4.2使用选择排序对最大距离进行排序，迭代数组dmax[Md]，取dmax[i]。
[0118] S4.3初始化dmax1的值为dmax[0]，dmax2的值为0。
[0119] S4.4若dmax[i]>dmax2，则将dmax[i]赋值给dmax2；若dmax[i]≤dmax2，则循环继续。
[0120] S4.5若dmax2> dmax1，则交换这两个数的值；若dmax2≤ dmax1，则循环继续。
[0121] S4.6判断循环是否结束：是，则输出距离排序在前两位的dmax1和dmax2及两个距离对应的坐标 、 、 和 作为平面拟合的点；否，则重复步骤S4.2‑S4.5。
[0122] 第五步（S5）、采用最小二乘法拟合平面。
[0123] 用最小二乘法拟合步骤S4中获取的四个角点 、 、 和坐标所在平面，即为图像 的剖面，设平面方程为：
[0124]
[0125] 其中，x，y，z是图像 的角点坐标，A、B、C是描述平面空间特征的常数。
[0126] 使该平面到四个角点的距离最近，根据最小二乘法，其响应函数S为：
[0127]
[0128] 计算出一组A、B、C，使得对四个角点来说S最小，对A、B、C分别求导，得到S最小时的A、B、C的值：
[0129]
[0130] 其中， 、 、 为平面中任一点的三维坐标，求解方程组，即可得到四个角点所在平面方程的系数A、B、C，得到的DLPFC剖面拟合示意图如图7所示。
[0131] 第六步（S6）、计算MNI空间下的DLPFC中心点 在步骤S5平面上的投影点 ，两点形成的直线为 ，则 的参数方程为：
[0132]
[0133] 则 。
[0134] 第七步（S7）、连接DLPFC中心点 与投影点 即可得到导针。
[0135] 综上，本发明的导航系统导针计算方法是基于Harris角点检测而实现的，利用二值图像运算、Harris角点检测、排序算法、平面拟合算法和投影点计算等多种方法，实现利用结构磁共振成像完成靶点计算，并将该方法与经颅磁刺激系统结合。相比其他定位导航技术，本发明不需要附加视觉系统和机械臂系统即可达到同样的定位效果，将精度误差控制在毫米级甚至更高，且有效避免了机械设定的误差，系统自动进行靶点计算，整个治疗过程5min内可以输出靶点坐标。
[0136] 以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，上述实施例不以任何形式限制本发明，凡采用等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。
