一种标记物识别及标记点定位方法及手术导航系统与流程

本发明属于医学图像处理技术领域，尤其涉及一种标记物识别及标记点定位方法及手术导航系统。

背景技术：

肺癌是世界上具有最高死亡率的疾病之一，活检是肺部肿瘤定性的金标准。常规的肺活检需要在CT的引导下进行肺穿刺手术取出活体组织进行病理检查。肺活检手术需要依靠医生经验进行临床操作，而CT不能在术中实时成像，因此若要成功取到活体组织，需要对患者穿刺多次，并进行多次CT扫描确认。肺穿刺手术导航系统可以避免这一问题。手术导航系统基于患者术前的CT图像进行路径规划，术中利用光学定位设备将术前CT图像空间与患者空间进行配准，从而实现由术前CT图像引导术中手术的目的，将会大大减少患者所受穿刺次数和放射性辐射，提高手术成功率。

早期研究表明基于人工标记物的点配准精度高于基于解剖特征的配准。目前，在图像空间和患者空间配准过程中，基于人工标记物的点配准是临床应用中最准确和广泛使用的配准方法。点配准需要在患者空间选定一组标记点，并在图像空间中选定相应的标记点，基于两组一一对应的标记点对两个空间进行配准。在患者空间中，标记点由刚性工具选择人工标记物的底部中心，并由光学定位仪和跟踪器获取其坐标。而在图像空间坐标系就要相应获取标记物的底部中心坐标作为图像空间中的标记点。标记点的定位精度是影响配准和导航精度的重要因素之一。

目前，对图像中的标记物识别以及标记点定位的方法已经有了多种算法的尝试。关于标记物识别研究人员提出过以下几类方法：一类方法是利用边缘检测分割皮肤及标记物表面，并基于模板进行标记物识别。然而，算法的精度对图像边缘过于敏感。一类方法是使用kmeans聚类和多项式曲线进行标记物识别。但是这种方法受到标记物粘贴位置的影响。还有一类方法是使用直方图驱动的3-D区域生长，并移除太接近的候选点从而识别标记物。但是该方法没有使用临床数据进行评估。关于标记点定位研究人员提出过以下几种方法：一类方法是使用阈值和几何滤波方法来识别候选标记物的标记点，然后通过计算得出的疑似标记物的连通量大小和形状，与已知的标记物的大小和半径进行比较，定位标记点。但在此方法中计算得出的加权中心点不是真正的标记点。一类方法是将投影图像与模板相匹配来定位标记点，另一类方法是基于标记物表面提取和几何先验知识来进行标记点定位。但这两种方法都使用了标记物的工厂制造参数，鲁棒性不高。还有一类方法是基于CT图像提取皮肤表面，根据形状指数和曲率来计算标记点。但这种方法需要大量点云计算，比较耗时。另外一类方法是基于配准进行标记点定位的，使用了标记物不同方向的互信息测度。但是这个方法是半自动的，必须在人为干预下才能实现。因而，目前已经开发出来的胸部CT图像的标记物识别及标记点定位的方法都存在着一些缺陷。

技术实现要素：

(一)要解决的技术问题

针对现有存在的技术问题，本发明提供一种标记物识别及标记点定位方法。

(二)技术方案

为了达到上述目的，本发明采用的主要技术方案包括：

一种标记物识别及标记点定位方法，其包括以下步骤：

步骤1：读取CT图像；

步骤2：对疑似标记物初分割，获取最终疑似标记物；

步骤3：基于序列特性的标记物识别，判断最终疑似标记物是否属于正确标记物；

步骤4：对正确标记物进行投影，根据投影图像计算正确标记物的法向量和厚度，从而实现对标记点的定位。

作为上述标记物识别及标记点定位方法的一种优选方案，

在步骤2中，进行疑似标记物初分割，获取最终疑似标记物包括以下步骤：

步骤2.1：种子点选取；

步骤2.2：对种子点进行区域生长，通过区域生长得到疑似标记物；

步骤2.3：去除所有的疑似标记物中噪声点，获取最终疑似标记物。

作为上述标记物识别及标记点定位方法的一种优选方案，

在步骤2.1中，种子点的选取方式为：在体素阈值上下限区间为[V₁，V₂]的范围中选取种子点，其中，V₁＝2800，V₂＝3025。

作为上述标记物识别及标记点定位方法的一种优选方案，

在步骤2.2中，进行区域生长的方式为：以初始种子点为中心，遍历初始种子点的邻域体素点，如果邻域体素点满足CT值大于设定值时，则将初始种子点与其进行合并，接着把它当作新的初始种子点，直到图像中的每一个点都有归属时，生长结束。

作为上述标记物识别及标记点定位方法的一种优选方案，

在步骤2.3中，去除噪声点的方法为：

将体素点个数大于3000或者小于500的疑似标记物去除。

作为上述标记物识别及标记点定位方法的一种优选方案，

在步骤3中，基于序列特性的标记物识别包括以下步骤：

步骤3.1：对每一个疑似标记物进行坐标转换；

步骤3.2：建立切割平面集合，对切割平面集合中的每一个平面做区域生长并计算切割平面特性；

步骤3.3：根据标记物各切割平面的特性得出标记物的序列特性，并判断该疑似标记物是否属于正确标记物。

作为上述标记物识别及标记点定位方法的一种优选方案，

在步骤3.2中，建立切割平面集合，对切割平面集合中的每一个平面做区域生长并计算切割平面特性具体方法为：

从疑似标记物中心向Z轴的正方向和反方向取n个切割平面；

对每一个切割平面做区域生长；

统计该平面区域生长的次数，定义该切割平面特性。

作为上述标记物识别及标记点定位方法的一种优选方案，

在步骤3.3中，确定是否属于正确标记物的具体方法为：

疑似标记物的序列特性是沿着Z轴的正方向将每个切割平面的特性组合到一起，并且将连续重复出现的切割平面特性合并为一个。根据不同的序列特性确定疑似标记物是否为真实标记物。

作为上述标记物识别及标记点定位方法的一种优选方案，

在步骤4中，对标记点的定位具体的方法为：

步骤4.1：遍历所有正确标记物的像素，然后将他们沿着Y轴方向投影到Z-X平面，得到投影图像；

步骤4.2：计算投影图像的长轴和短轴，求出投影图像的参数；

步骤4.3：根据公式得到标记物厚度d以及标记物与投影图像的夹角θ；

步骤4.4：以外环的长轴OA为旋转轴将Y轴旋转θ角得到标记物法向量，并将标记物体中心沿着法向量反方向向底部移动d/2距离，得到标记点。

一种手术导航系统，其采用以上所述的标记物识别及标记点定位方法。

(三)有益效果

本发明的有益效果是：本发明提供的标记物识别及标记点定位方法运用于肺穿刺手术导航系统的配准过程中，其基于标记物的几何形状特性，采用序列特性的方法对标记物进行识别，并采用投影图像的方法对标记点进行定位。该方法无需使用标记物的工厂制造参数，避免了由于长时间使用的老化等原因，使标记物的工厂制造参数与实际尺寸存在误差而对定位精度产生影响。并且本方法不受标记物粘贴的位置和状态的影响。因此，本算法的精度和鲁棒性都较好。

附图说明

附图1是本发明所涉及方法的总体流程图。

附图2是本发明所设计方法的步骤3中基于序列特性进行标记物识别的流程图。

附图3是本发明所设计方法的步骤4中基于投影图像进行标记点定位的流程图。

附图4为本发明设计方法的步骤3.1中新坐标系O-XYZ和原始坐标系O'-X'Y'Z'变换关系示意图。

附图5为本发明设计方法的步骤3.2中标记物切割平面序列示意图。

附图6为本发明设计方法的步骤4.1，4.2中标记物投影图像示意图。

附图7为本发明设计方法的步骤4.3，4.4中标记点计算示意图。

附图8为本发明设计方法的步骤4.4标记物法向量确定示意图。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。

在本实施方式中，提供了一种标记物识别及标记点定位方法及手术导航系统，适用于肺部手术导航系统配准过程中。在本实施方式中，以胸部CT图像的标记物识别及标记点定位为例进行详细的说明，具体的如以下所述。

一种基于胸部CT图像的标记物识别及标记点定位方法运行在Visual Studio平台，适用于肺穿刺手术导航系统的配准过程中，实现基于胸部CT图像的标记物识别及标记点定位方法。

步骤1：读取胸部CT图像：

按照DICOM图像的格式标准，读取胸部CT图像数据。

步骤2：疑似标记物初分割，获取最终疑似标记物：

进行疑似标记物初分割具体步骤如下：

步骤2.1：种子点选取。在体素阈值上下限区间为[V₁，V²]的范围中选取种子点，V₁＝2800，V₂＝3025。

步骤2.2：进行26邻域区域生长。以初始种子点为中心，遍历初始种子点的26邻域体素点，如果邻域体素点满足CT值大于1320，则将初始种子点与其进行合并，接着把它当作新的初始种子点，直到图像中的每一个点都有归属时，生长结束。通过区域生长得到疑似标记物。

步骤2.3：去除噪声点。将体素点个数大于3000或者小于500的疑似标记物去除，进一步去除噪声。获得最终的疑似标记物，并将每个疑似标记物进行编号，使该疑似标记物的所有体素点的CT值设置为编号值。

步骤3：基于序列特性的标记物识别，判断最终疑似标记物是否属于正确标记物；

进行标记物识别具体步骤如下：

步骤3.1：对每一个疑似标记物进行坐标转换：

如图4所示，以数据体的中心O’到标记物的中心O的延长线所在的矢量为新坐标系的Y轴，经过标记物的中心O，并与Y轴垂直的平面为X-Z平面。而新坐标系O-XYZ与原坐标系O’-X’Y’Z’的坐标转换关系如公式(1)所示。

其中，α,β,γ是Y轴与Z’轴，X’轴和Y’轴的夹角，O(x₀,y₀,z₀)是新坐标系的原点，即标记物的中心。T是从新坐标系O-XYZ到原坐标系O’-X’Y’Z’的坐标转换矩阵。

所以，原坐标系上的任一点P'(x',y',z')在新坐标系下的坐标P(x,y,z)为：

P(x,y,z)＝[x',y',z',1]·T^-1。

步骤3.2：建立切割平面集合，对切割平面集合中的每一个平面做8邻域2维区域生长并计算切割平面特性。

从疑似标记物中心向Z轴的正方向和反方向以1mm的间距各取与Z轴垂直的10个切割平面，建立切割平面集合S＝{Si，i＝-10,......,10},该集合共21个切割平面。

对每一个切割平面做8邻域2维区域生长。以属于该疑似标记物的第一个点为种子点进行区域生长，把相同标记物编号的点扩充进种子点所在区域。当不再有像素满足这个准则时，区域生长停止。

为了进一步去除噪声，将区域生长后像素点个数小于N₁的连通量去除，其中N₁＝20。

统计每层区域生长的次数，定义该切割平面特性：

切割平面特性分为三种情况，第一种情况是切割平面只区域生长一次，得到一块连通量的标记物切面，定义为A特性，第二种情况是区域生长次数为两次，得到两块连通量的标记物切面，定义为B特性，其他情况我们定义为C特性。如图5中(a),(b),(c),(f),(g),(h)为A特性，(d),(e)为B特性。

步骤3.3：根据标记物各切割平面的特性得出标记物的序列特性，并判断该疑似标记物是否属于正确标记物。

疑似标记物的序列特性是沿着Z轴的正方向将每个切割平面的特性组合到一起，并且将连续重复出现的切割平面特性合并为一个，比如“AAA”可以定义为“A”。例如，在图5中，标记物的序列特性为“ABA”。对序列特性为“ABA”的疑似标记物识别为真实标记物。

步骤4：对标记物进行投影，根据投影图像计算标记物的法向量和厚度，从而实现对标记点的定位。此步骤中的标记点指的是图像空间中的标记点。

步骤4.1：标记物投影图像获取：遍历所有标记物的体素，然后将他们沿着Y轴方向投影到X-Z平面，得到投影图像。在沿着Y轴方向投影过程中，不论遇到一个还是多个标记物体素，投影图像相应像素点均置1。

通常，投影图像是椭圆环，如图6所示。如果投影方向平行于标记的法向量，则投影图像的形状将是圆形。外环的长轴和短轴为OA和OB，内环的长轴和短轴OA'和OB'。

步骤4.2：投影图像长轴和短轴的计算：遍历投影图像中所有像素并计算点O到每个像素的距离。如图6所示，最短距离记为OB’。OB’的延长线和投影图像外环的交点记为点B。同理，A’和A分别是OB’的垂直线与投影图像内环和外环的交点。

步骤4.3：标记物的厚度以及标记物与投影平面的夹角θ计算：

如图7所示，经过投影图像的Y轴和短轴OB'的平面显示出了如何计算标记物与投影平面的夹角θ和厚度。点O既是标记物的中心又是投影图像的中心。标记物的厚度可以描述为d，θ是标记物和投影平面之间的夹角，理论上在0°和180°之间。

d和θ可以通过以下公式计算：

其中R和r分别定义为标记物外圆的半径和内圆的半径。标记物的半径R对应于投影图像的长轴OA，r对应于投影图像的OA'。可以从投影图像计算OA，OA'，OB，OB'。

步骤4.4：确定标记物法向量及标记点定位：

如图8所示，若Y轴的单位向量与CT扫描轴Z’轴上的单位向量的内积大于0，则将Y方向上的单位向量沿着OA顺时针旋转θ度得到标记物法向量的单位向量，若Y轴的单位向量与CT扫描轴Z’轴上的单位向量的内积小于0，则将Y方向上的单位向量沿着OA逆时针旋转θ度得到标记物法向量的单位向量。如图7所示，将标记物的体中心沿着单位法向量的反方向向标记物底部移动d/2距离，得到标记点。

综上所述，该方法根据人工标记物的序列特性对标记物识别，根据投影图像对标记点(即标记物的底部中心点)定位。具体步骤如下：对粘贴人工标记物的胸部CT图像预处理，得到疑似标记物；对疑似标记物进行坐标变换，用新的坐标轴作为切割轴得到标记物的切割平面，计算所有切割平面的序列特性，并根据序列特性进行标记物识别；对标记物进行投影，根据投影图像计算标记物的法向量和厚度，从而实现对标记点的计算。该方法不依赖于标记物的工厂制造参数，较少受标记物粘贴的位置和条件的影响，具有较好的准确率和鲁棒性。

以上结合具体实施例描述了本发明的技术原理，这些描述只是为了解释本发明的原理，不能以任何方式解释为对本发明保护范围的限制。基于此处解释，本领域的技术人员不需要付出创造性的劳动即可联想到本发明的其它具体实施方式，这些方式都将落入本发明的保护范围之内。