基于口腔CBCT图像的牙根提取方法.pdf

《基于口腔CBCT图像的牙根提取方法.pdf》由会员分享，可在线阅读，更多相关《基于口腔CBCT图像的牙根提取方法.pdf（23页完成版）》请在专利查询网上搜索。

1、(19)中华人民共和国国家知识产权局 (12)发明专利申请 (10)申请公布号 (43)申请公布日 (21)申请号 201910421368.8 (22)申请日 2019.05.21 (71)申请人 重庆布瑞斯科技有限公司 地址 400000 重庆市九龙坡区石桥铺渝州 路51号附8号 (72)发明人 刘宇 (74)专利代理机构 重庆中流知识产权代理事务 所(普通合伙) 50214 代理人 陈立荣 (51)Int.Cl. G06T 7/136(2017.01) G06T 7/11(2017.01) G06T 5/00(2006.01) G16H 30/20(2018.01) G16H 50/50。

2、(2018.01) (54)发明名称 一种基于口腔CBCT图像的牙根提取方法 (57)摘要 本发明公开了一种基于口腔CBCT图像的牙 根提取方法， 包括以下步骤： S100， 光学扫描获得 三维牙颌和牙列模型， 建立投影图像及CBCT文件 序列； S200， 通过MarchingCube、 区域增长算法 解析CBCT文件序列； S300， 基于CBCT文件序列， 对 上下牙颌进行区分， 其区分过程包括通过VTK定 位咬合面和矢状面的关键帧， 执行限制性区域增 长算法； S400， 优化上下牙颌的分割效果， 并对牙 颌模型进行手动裁剪； S500， 去除裁剪后牙颌模 型的噪音图像， 并对其进行平。

3、滑。 本发明降低了 手术风险， 提高了手术成功率。 权利要求书3页 说明书7页 附图12页 CN 110189352 A 2019.08.30 CN 110189352 A 1.一种基于口腔CBCT图像的牙根提取方法， 其特征在于， 所述方法包括以下步骤： S100， 光学扫描获得三维牙颌和牙列模型， 建立投影图像及CBCT文件序列； S200， 通过Marching Cube、 区域增长算法解析CBCT文件序列； S300， 基于CBCT文件序列， 对上下牙颌进行区分， 其区分过程包括通过VTK定位咬合面 和矢状面的关键帧， 执行限制性区域增长算法； S400， 优化上下牙颌的分割效果， 并。

4、对牙颌模型进行手动裁剪； S500， 去除裁剪后牙颌模型的噪音图像， 并对其进行平滑。 2.根据权利要求1所述的基于口腔CBCT图像的牙根提取方法， 其特征在于， 所述 Marching Cube的实现包括以下步骤： S201， 将原始数据经过预处理之后， 读入指定的数组中； S202， 从网格数据体中提取一个单元体， 成为当前单元体， 同时获取该单元体的所有信 息， 其信息包括边界的顶点函数值、 和单元体的点云坐标位置； S203， 将当前单元体中顶点的函数值与给定等值面值C进行比较， 得出该单元体的状态 表； S204， 根据当前单元体的状态表索引， 找出与等值面相交的单元体棱边， 并采用。

5、线性插 值的方法， 计算出各个交点的位置坐标； S205， 利用中心差分法， 求出当前单元体中顶点的法向量， 再采用线性插值的方法， 得 到三角面片各个顶点的法向量； S206， 根据各个三角面片顶点的坐标， 顶点法向量进行等值面图象的绘制。 3.根据权利要求1所述的基于口腔CBCT图像的牙根提取方法， 其特征在于， 所述步骤 S200执行区域增长算法时， 包括将阈值控制在770-25584区间， 重建出3D模型查看CBCT主要 的解剖结构。 4.根据权利要求1所述的基于口腔CBCT图像的牙根提取方法， 其特征在于， 所述步骤 S300中的区域增长算法包括其输入阈值和拼接向量的改变， 其拼接向。

6、量公式如下： O1(x,y)F(x,y)s3； O2(x,y)s3t3； OO1O2(x,y)F(x,y)s3(x,y)s3t3； B wf1,(x,y)O,0,(x,y)其他； 其中， (s3,t3)为二维阈值向量， s3为灰度阈值， t3为梯度阈值， 将阈值小于s3像点作为 目标像点集； 为变量参数， 其大小与图像模糊程度对应； B wf为CT中的灰度图像值。 5.根据权利要求1所述的基于口腔CBCT图像的牙根提取方法， 其特征在于， 所述步骤 S400包括建立算法和搜索算法， 构建算法的步骤如下： S401， 按XX轴的分割线对空间进行分割； S402， 计算所有点的xx坐标的平均值， 。

7、选择所有点中最接近平均值的点作为分割线， 把 空间中的点进行分割； S403， 对分隔后的空间按照yy轴的分割线进行分割； S404， 计算所有点的xx坐标的平均值， 选择所有点中最接近平均值的点作为分割线， 把 空间中的点进行分割； S405， 对分割后的空间继续按xx轴进行分割， 依次类推， 以XX轴和YY轴为目标对象循环 权利要求书 1/3 页 2 CN 110189352 A 2 分割， 每个空间中分割至一个点时， 结束分割。 6.根据权利要求5所述的基于口腔CBCT图像的牙根提取方法， 其特征在于， 所述步骤 S401至步骤S405的分割过程对应于一个二叉树， 其每条分割线对应于二叉。

8、树中的一个分 支， 其空间中每个点对应于一个叶子节点。 7.根据权利要求6所述的基于口腔CBCT图像的牙根提取方法， 其特征在于， 所述二叉树 中包括空间点(X， Y)， 查找空间点(X， Y)的近邻点， 其查找方法步骤如下： S1， 对二叉树进行遍历， 到达叶子节点(X ， Y )； S2， 计算(X， Y)与(X ， Y )的距离； S3， 进行回溯， 计算上一层节点(x” ,y” )与(X， Y)的距离； S4， 比较(X ， Y )和(x” ,y” )分别与(X， Y)的距离， 以(X， Y)为圆心， 以两者分别相比之下 最近距离的距离为半径画一个圆； S5， 若上述步骤中的圆与三维空。

9、间中及圆直径长度一致的立方体的分割线有交点， 则 对上述步骤中所比较最近距离的点进行遍历， 比较遍历点与其的距离， 继续回溯， 并以(X， Y)为圆心， 遍历点与其的距离为半径画圆， 当所画之圆与分割线没有交点时， 即该遍历点为 最近邻点， 查找结束。 8.根据权利要求1或5所述的基于口腔CBCT图像的牙根提取方法， 其特征在于， 所述步 骤S400还包括根据k-d树进行无序点云去噪， 去噪过程包括如下步骤： S6， 根据点云数据生成k-d树， 建立整个3D图形的点云三维坐标及其在三维空间中相互 关联的拓扑关系； S7， 查找任一点的的邻域； S8， 计算该点与邻域内各点的距离， 并取其平均值。

10、； S9， 判断该平均值是否超过阈值， 若超过则判定该点为噪点， 进行去除。 9.根据权利要求1所述的基于口腔CBCT图像的牙根提取方法， 其特征在于， 所述步骤 S500包括通过K-D树建立之后进行无序点云的去噪算法， 还包括基于Poisson方程的三角网 格补洞算法进行平滑。 10.根据权利要求1或9所述的基于口腔CBCT图像的牙根提取方法， 其特征在于， 所述步 骤S500包括根据三角网格中孔洞边界生成一个初始化补洞网格， 而后通过法向估算和 Poisson方程来修正补洞网格中三角面片的几何形状， 使其能够适应并与周围的原始网格 融合， 其过程包括以下步骤： S501， 检测孔洞边界并初。

11、始化补洞网格； S502， 计算补洞网格中顶点的期望法向， 构建Laplace方程f0求解补洞网格内部顶 点的法向分布， 其Laplace算子为： 其中， N1(xi)表示顶点xi的1环邻域点， ij和 ij为边eij对应的2个对角； S503， 基于期望法向旋转补洞网格中的三角面片； S504， 基于Poisson方程调整补洞网格顶点位置， 并计算三维网格的梯度场， 公式为： 权利要求书 2/3 页 3 CN 110189352 A 3 其中， f为待求的调整后网格顶点位置， w为撕裂网格的梯度场； 梯度算子为： 其中， 基函数梯度i的表达式为： 其表示将向量逆时针旋转90度， AT表示三角。

12、片T的面积； 散度算子： 其中， T1(xi)表示顶点xi的1环邻域三角片， AT表示三角片T的面积。 权利要求书 3/3 页 4 CN 110189352 A 4 一种基于口腔CBCT图像的牙根提取方法 技术领域 0001 本发明涉及口腔CBCT图像处理技术领域， 尤其涉及一种基于口腔CBCT图像的牙根 提取方法。 背景技术 0002 目前临床中主要有两种牙科全景图像生成方法： 牙科全景X射线摄影方法和口腔 CBCT重建后得到体数据的方法。 0003 牙科全景X射线摄影方法主要采用X摄像和平板探测器以一种特殊的轨迹进行运 动， 在运动过程中采集特定位置的图像， 并最终将采集到的序列图像组合成。

13、牙科全景图像， 该图像是有设备直接采集获取， 对人体具有一定的辐射， 该方法缺点是只能获取牙齿的二 维投影数据， 无法获取牙齿三维空间的形态学信息。 0004 口腔CBCT重建后得到体数据的方法， 能够方便医生在三维空间中分析牙齿的形态 学结构， 并获取相关的信息， 也可以由重建体数据生成牙齿的全景图像， 这样对像就可少做 X线全景扫描， 从而减少辐射剂量， 这种方式就需要由重建体数据从而生成牙科全景图像。 但是， 现有技术中的大部分产品无法支持基于口腔CBCT图像的牙根查看的功能， 更不能达 到3D视图中数字化诊断的定位效果， 这样导致医生数字化术前规划过程中无法预判牙根的 情况， 进而增加。

14、手术的风险和手术的失败率。 例如在口腔根管根尖定位过程中， 提前预知牙 根和器械的相对位置将决定手术的成败； 种植修复过程中， 邻牙牙根的持续监控和判断将 决定手术种植体的生命周期； 正畸过程中需要持续对牙根的走势提前预判并控制， 决定了 正畸治疗的时间和最终效果。 0005 近年来， 口腔修复、 种植及正畸诊断分析过程中对于牙根的查看越发受到重视， 现 有技术中又缺乏一种精准牙根提取、 并使其3D可视化以及建模测试的方法， 无法降低相关 手术的难度和风险， 尤其是在精度要求很高的复杂手术中， 牙根的相对位置无法预判的情 况下， 手术无法进行。 0006 综上所述， 亟待一种解决上述问题的基于。

15、口腔CBCT图像的牙根提取方法。 发明内容 0007 本发明提供了一种基于口腔CBCT图像的牙根提取方法， 降低了手术风险， 提高了 手术成功率。 0008 本发明解决其技术问题所采用的技术方案是： 一种基于口腔CBCT图像的牙根提取 方法， 包括以下步骤： 0009 S100， 光学扫描获得三维牙颌和牙列模型， 建立投影图像及CBCT文件序列； 0010 S200， 通过Marching Cube、 区域增长算法解析CBCT文件序列； 0011 S300， 基于CBCT文件序列， 对上下牙颌进行区分， 其区分过程包括通过VTK定位咬 合面和矢状面的关键帧， 执行限制性区域增长算法； 0012。

16、 S400， 优化上下牙颌的分割效果， 并对牙颌模型进行手动裁剪； 说明书 1/7 页 5 CN 110189352 A 5 0013 S500， 去除裁剪后牙颌模型的噪音图像， 并对其进行平滑。 0014 进一步地， 所述Marching Cube的实现包括以下步骤： 0015 S201， 将原始数据经过预处理之后， 读入指定的数组中； 0016 S202， 从网格数据体中提取一个单元体， 成为当前单元体， 同时获取该单元体的所 有信息， 其信息包括边界的顶点函数值、 和单元体的点云坐标位置； 0017 S203， 将当前单元体中顶点的函数值与给定等值面值C进行比较， 得出该单元体的 状态表。

17、； 0018 S204， 根据当前单元体的状态表索引， 找出与等值面相交的单元体棱边， 并采用线 性插值的方法， 计算出各个交点的位置坐标； 0019 S205， 利用中心差分法， 求出当前单元体中顶点的法向量， 再采用线性插值的方 法， 得到三角面片各个顶点的法向量； 0020 S206， 根据各个三角面片顶点的坐标， 顶点法向量进行等值面图象的绘制。 0021 进一步地， 所述步骤S200执行区域增长算法时， 包括将阈值控制在770-25584区 间， 重建出3D模型查看CBCT主要的解剖结构。 0022 进一步地， 所述步骤S300中的区域增长算法包括其输入阈值和拼接向量的改变， 其拼接。

18、向量公式如下： 0023 O1(x,y)F(x,y)s3； 0024 O2(x,y)s3t3； 0025 OO1O2(x,y)F(x,y)s3(x,y)s3t3； 0026 B wf1,(x,y)O,0,(x,y)其他； 0027 其中， (s3,t3)为二维阈值向量， s3为灰度阈值， t3为梯度阈值， 将阈值小于s3像点 作为目标像点集； 为变量参数， 其大小与图像模糊程度对应； B wf为CT中的灰度图像值。 0028 进一步地， 所述步骤S400包括建立算法和搜索算法， 构建算法的步骤如下： 0029 S401， 按XX轴的分割线对空间进行分割； 0030 S402， 计算所有点的xx。

19、坐标的平均值， 选择所有点中最接近平均值的点作为分割 线， 把空间中的点进行分割； 0031 S403， 对分隔后的空间按照yy轴的分割线进行分割； 0032 S404， 计算所有点的xx坐标的平均值， 选择所有点中最接近平均值的点作为分割 线， 把空间中的点进行分割； 0033 S405， 对分割后的空间继续按xx轴进行分割， 依次类推， 以XX轴和YY轴为目标对象 循环分割， 每个空间中分割至一个点时， 结束分割。 0034 进一步地， 所述步骤S401至步骤S405的分割过程对应于一个二叉树， 其每条分割 线对应于二叉树中的一个分支， 其空间中每个点对应于一个叶子节点。 0035 进一步。

20、地， 所述二叉树中包括空间点(X， Y)， 查找空间点(X， Y)的近邻点， 其查找方 法步骤如下： 0036 S1， 对二叉树进行遍历， 到达叶子节点(X ， Y )； 0037 S2， 计算(X， Y)与(X ， Y )的距离； 0038 S3， 进行回溯， 计算上一层节点(x” ,y” )与(X， Y)的距离； 0039 S4， 比较(X ， Y )和(x” ,y” )分别与(X， Y)的距离， 以(X， Y)为圆心， 以两者分别相比 说明书 2/7 页 6 CN 110189352 A 6 之下最近距离的距离为半径画一个圆； 0040 S5， 若上述步骤中的圆与三维空间中及圆直径长度一。

21、致的立方体的分割线有交 点， 则对上述步骤中所比较最近距离的点进行遍历， 比较遍历点与其的距离， 继续回溯， 并 以(X， Y)为圆心， 遍历点与其的距离为半径画圆， 当所画之圆与分割线没有交点时， 即该遍 历点为最近邻点， 查找结束。 0041 进一步地， 所述步骤S400还包括根据k-d树进行无序点云去噪， 去噪过程包括如下 步骤： 0042 S6， 根据点云数据生成k-d树， 建立整个3D图形的点云三维坐标及其在三维空间中 相互关联的拓扑关系； 0043 S7， 查找任一点的的邻域； 0044 S8， 计算该点与邻域内各点的距离， 并取其平均值； 0045 S9， 判断该平均值是否超过阈。

22、值， 若超过则判定该点为噪点， 进行去除。 0046 进一步地， 所述步骤S500包括通过K-D树建立之后进行无序点云的去噪算法， 还包 括基于Poisson方程的三角网格补洞算法进行平滑。 0047 进一步地， 所述步骤S500包括根据三角网格中孔洞边界生成一个初始化补洞网 格， 而后通过法向估算和Poisson方程来修正补洞网格中三角面片的几何形状， 使其能够适 应并与周围的原始网格融合， 其过程包括以下步骤： 0048 S501， 检测孔洞边界并初始化补洞网格； 0049 S502， 计算补洞网格中顶点的期望法向， 构建Laplace方程求解补洞网格内部顶点 的法向分布， 其Laplac。

23、e算子为： 0050 0051 其中， N1(xi)表示顶点xi的1环邻域点， ij和 ij为边eij对应的2个对角； 0052 S503， 基于期望法向旋转补洞网格中的三角面片； 0053 S504， 基于Poisson方程调整补洞网格顶点位置， 并计算撕裂网格的梯度场， 公式 为： 0054 其中， f为待求的调整后网格顶点位置， w为撕裂网格的梯度场； 0055梯度算子为： 0056其中基函数梯度的表达式为： 0057 其表示将向量逆时针旋转90度， AT表示三角片T的面积； 0058散度算子： 0059 其中， T1(xi)表示顶点xi的1环邻域三角片， AT表示三角片T的面积。 00。

24、60 本发明的有益效果在于： 本发明通过对牙根形态的提取， 形成模拟图像， 帮助医生 提前预测出术中甚至术后的牙根趋势， 极大的提高了手术的精确性， 降低了手术风险。 同 说明书 3/7 页 7 CN 110189352 A 7 时， 本发明通过优化、 去噪和平滑等图像处理， 提取出图像质量较高的牙根模型， 此模型不 仅可以用于数字化模拟术前规划， 更可以通过3D打印产生实体用于术前模拟测试， 拟真地 预先展现出术中所出现的情况。 本发明操作流程简单， 医生在数字化术前规划和诊断过程 中， 对牙根的解剖结构进行提前的预判， 使得术前规划和术中情况一致， 同时使手术的规划 精准程度得到量化， 进。

25、而提高了医生数字化程度的意识， 以及增强了用户体验。 另外， 经发 明人试验， 这样的处理方法同样适用于其他医疗科室相关的CBCT图像、 MRI图像中的精细解 剖结构的提取过程。 附图说明 0061 为了更清楚地说明本发明实施例或现有技术中的技术方案， 下面将结合附图及实 施例对本发明作进一步说明， 下面描述中的附图仅仅是本发明的部分实施例， 对于本领域 普通技术人员来讲， 在不付出创造性劳动的前提下， 还可以根据这些附图获得其他附图： 0062 图1为本发明实施例1中基于口腔CBCT图像的牙根提取方法的步骤流程图； 0063 图2为本发明实施例1中Marching Cube实现的步骤流程图；。

26、 0064 图3为本发明实施例1中构建建立算法和搜索算法的步骤流程图； 0065 图4为本发明实施例1中空间最邻近点查找方法的步骤流程图； 0066 图5为本发明实施例1中去噪步骤流程图； 0067 图6为本发明实施例1中空间切分示意图； 0068 图7为本发明实施例1中空间最邻近点搜索过程示意图； 0069 图8为本发明实施例1中二叉树示意图； 0070 图9为本发明实施例1中空间分割线交点示意图； 0071 图10为本发明实施例1中空间分割线交点示意图； 0072 图11为本发明实施例1中孔洞边界示意图； 0073 图12为本发明实施例1中补洞网络示意图； 0074 图13为本发明实施例1。

27、中三角面片示意图； 0075 图14为本发明实施例1中调整补洞网格的结构示意图； 0076 图15为本发明实施例1中调整补洞网格的结构示意图； 0077 图16为本发明实施例1中调整补洞网格的结构示意图； 0078 图17为本发明实施例1中CBCT的3D建模结构示意图； 0079 图18为本发明实施例1中关键帧的定位示意图； 0080 图19为本发明实施例1中限制性区域增长效果示意图； 0081 图20为本发明实施例1中meshlab手动裁剪示意图； 0082 图21为本发明实施例1中裁剪效果示意图； 0083 图22为本发明实施例1中最终裁剪效果示意图。 具体实施方式 0084 为了使本发明。

28、实施例的目的、 技术方案和优点更加清楚， 下面将结合本发明实施 例中的技术方案进行清楚、 完整的描述， 显然， 所描述的实施例是本发明的部分实施例， 而 不是全部实施例。 基于本发明的实施例， 本领域普通技术人员在没有付出创造性劳动的前 说明书 4/7 页 8 CN 110189352 A 8 提下所获得的所有其他实施例， 都属于本发明的保护范围。 0085 在实施例1中， 如图1、 2、 3、 4、 17、 18、 19、 20、 21及22所示， 首先， 通过解析口腔CBCT 的文件序列(dcm格式)， 并执行区域增长算法， 同时将CBCT图像阀值控制在770-25584区间， 进而重建出。

29、CBCT的主要解剖结构(3D视图)； 而后基于vtk(第三方库)中针对咬合面和矢状 面， 对关键帧进行交叉定位， 递归执行限制性区域增长算法， 同时调整CBCT图像阀值在 1565-9925区间， 进行K-D Tree算法做特征提取， 利用meshlab的手动裁剪功能进行粗略去 噪； 最后调用K-DTree算法进行精细化去噪， 并采用基于Poisson方程的三角网格补洞方法 对模型进行表面处理。 0086 具体地， 如图6所示， 对于步骤S200的Marching Cube方法和区域增长算法， Marching Cube利用等值面， 即三维立体空间中所有具有相同值的点的集合， 其集合区域展 现。

30、类比如地形图里的等高线。 考虑到Marching Cube中的基本假设条件， 即空间中存在立方 体(六面体)， 而沿六面体边的数据场是呈连续性变化， 即当立方体一条边的两个顶点分别 大于或小于等值面的值， 则在该条边上有且仅有一点是这条边与等值面的交点。 在空间中 过程， 即是用大量小正方体对空间进行切分， 进而用小正方体内部的平面来近似表示当前 的等值面， 故此， 小正方体的数量越多， 逼近的效果越好， 随之进行的是大量的计算。 计算式 子包括(x,y,z)|f(x,y,z)c， 其中， ,c为常数逐个处理数据场中的立方体(体素)， 其计 算过程包括分离出与等值面相交的立方体后， 采用双线性。

31、插值计算出等值面与立方体边的 交点， 根据立方体每一顶点与等值面的相对位置， 将等值面与立方体边的交点按三维空间 中垂直地方式连接生成等值面， 从而作为等值面在该立方体内的一个逼近表示。 0087 进一步地， 如图2所示， 执行步骤S201， 将原始数据经过预处理之后， 读入指定的数 组中； 执行步骤S202， 从网格数据体中提取一个单元体， 成为当前单元体， 同时获取该单元 体的所有信息， 其信息包括边界的顶点函数值、 和单元体的点云坐标位置， 在本实施例中， 以8个顶点为例； 执行步骤S203， 3.将当前单元体8个顶点的函数值与给定等值面值C进行比 较,得到该单元体的状态； (edgeT。

32、able、 triTable)； 执行步骤S204， 根据当前单元体的状态 表索引， 找出与等值面相交的单元体棱边， 并采用线性插值的方法， 计算出各个交点的位置 坐标； 执行步骤S205， 利用中心差分法， 求出当前单元体8个顶点的法向量， 在采用线性插值 的方法， 得到三角面片各个顶点的法向； 执行步骤S206， 根据各个三角面片顶点的坐标， 顶 点法向量进行等值面图象的绘制。 0088 具体的， 在中心差分法计算过程中， 例如采用等时间步长， t(i)t(t为常 数)， 用u表示位移， 那么速度和加速度的中心差分近似为： u(i)u(i+1)-u(i-1)/(2 t)； u” (i)u(。

33、i+1)-2u(i)+u(i-1)/(t*t)。 0089 其中， 对于步骤S200中的区域增长算法， 获取分割区域中的一个种子点(即特征 点)， 在种子点的周围搜索与该种子点有相似性质的像素点， 并将其合并到种子区域中， 而 后将合并的像素作为新的种子点继续搜索， 直到种子区域中所有像素周围没有相似的像素 点， 算法结束。 其基本算法流程为： 步骤1， 选取种子点p(x0， y0)， 用堆栈表示种子区域， 将种 子点push到种子堆栈中； 步骤2， 将种子堆栈中第一个种子点pop出堆栈， 并以该点为中心， 遍历该中心8邻域像素； 判断遍历像素点是否已经在种子区域中， 如果否， 判断遍历像素点。

34、 是否满足相邻种子点相似性， 如果像素点(x， y)满足相似性， 将(x， y)push到堆栈中； 步骤4， 重复步骤2至3， 直至种子堆栈为空。 由此可见， 影响区域生长算法的要素有三个： 种子点的 说明书 5/7 页 9 CN 110189352 A 9 选取； 搜索路径的选择； 像素相似性的判断。 0090 具体地， 对于种子点的选取， 区域生长算法是半交互式的分割算法， 需要用户选取 种子点， 也可以是通过其他算法计算出来的种子点； 对于搜索路径的选择， 其搜索路径是选 择相邻的像素， 以二维图像为例， 一般为8邻域搜索， 或者4邻域搜索， 以三维图像为例， 一般 为26邻域搜索或者6。

35、邻域搜索； 对于像素相似性的判断， 其相似性以像素值的相近程度为判 断标准， 例如， 设置一定灰度范围做为相似的标准， 也可以通过计算满足某种形状或者性质 作为判断标准。 0091 具体地， 对于拼接向量公式， 目标集合： O1(x,y)F(x,y)s3； 由于足迹图像中 噪声和模糊边缘像素灰度值要大于s3， 所以若用单阈值s3对图像分割， 噪声和边缘像点将 被错划为背景， 此外， 噪声和CT图像中模糊边缘像素点具有较大梯度值,所以用二维阈值向 量(s3,t3)限制这部分像素,即把集合O2(x,y)s3t3也作为 目标集合， 因此， 图像分割后目标像素集应为OO1O2(x,y)F(x,y)s3。

36、(x,y)s3t3。 其中， 参数 选取与图像模糊程度有关， 在足迹图像中 取5 或7即可满足要求。 令分割后二值图像， 即像素值为0-255的灰度图像为B wf,则Bwf1,(x, y)O,0,(x,y)其他。 0092 具体地， 如图7及8所示， 步骤S300无法完全性地区分上下颌的所有点云， 需要采取 K-DTree算法和手动裁剪的方法提取牙根， 在步骤S400中， 在包括K-DTree的建立算法和搜 索算法的同时， 其手动裁剪还使用了meshlab的固有功能。 在执行步骤S401至步骤S405的过 程中， 其搜索过程如图所示， 在本实施例中， 建立起的二叉树包括如下一组二维点集： (2。

37、, 3),(8,1),(9,6),(4,7),(7,2),(5,4)(2,3),(8,1),(9,6),(4,7),(7,2),(5,4)， 例如对 于点(2,4.5)(2,4.5)， 对其最近邻点的查找方法如下： 如图所示， 在所建立好的k-d树上进 行遍历， 因为(2,4.5)(2,4.5)的x坐标比7小， 因此进入(7,2)(7,2)的左子树， 由于(2,4.5) (2,4.5)的纵坐标比4大， 因此进入(5,4)(5,4)的右子树， 到达叶子节点(4,7)(4,7)； 计算 (2,4.5)(2,4.5)与(4,7)(4,7)的距离为3.202， 现在假定(4,7)(4,7)为最近邻点，。

38、 进行回 溯， 计算(2,4.5)(2,4.5)与(5,4)(5,4)的距离为3.041， 则假定最近邻点为(5,4)(5,4)； 以 (2,4.5)(2,4.5)为圆心， 以3.041为半径画一个圆； 如图9所示， 可见该圆与(5,4)(5,4)所在 的分割线有交点， 因此需要对(5,4)(5,4)的左子树进行遍历， 发现与(2,3)(2,3)的距离为 1.5， 继续回溯， 而以(2,4.5)(2,4.5)为圆心， 以1.5为半径画圆， 与(7,2)(7,2)所在的分割 线没有交点， 因此不需要对(7,2)(7,2)的右子树进行遍历， 如图10所示， 所以最终找的最近 邻点为(2,3)(2,。

39、3)， 查找结束。 0093 具体地， 如图11、 12、 13、 14、 15及16所示， 手动裁剪过后依然会有不少噪音， 执行K- D Tree算法去噪去噪完成后， 表面依然存在粗糙情况， 需要执行平滑算法进行优化， 在步骤 S500中， 对于Poisson方程的三角网格补洞方法， 该算法首先需要根据孔洞边界生成一个初 始化补洞网格， 然后通过法向估算和Poisson方程来修正补洞网格中三角面片的几何形状， 使其能够适应并与周围的原始网格融合。 算法的主要步骤如下： 如图所示， 检测孔洞边界并 初始化补洞网格， 由于初始化补洞网格无法与原始孔洞周围的网格有效融合， 因此需要调 整补洞网格的。

40、顶点位置使得补洞网格与原始网格之间光滑过渡； 计算补洞网格中顶点的期 望法向， 由于已知原始网格孔洞边界的法向， 将其作为补洞网格边界的法向， 构建Laplace 方程f0求解补洞网格内部顶点的法向分布， 假设f表示在每个顶点上的标量， 那么网格 说明书 6/7 页 10 CN 110189352 A 10 域上在顶点xi处的Laplace算子定义如下(不考虑面积影响)： 0094 0095 其中N1(xi)表示顶点xi的1环邻域点， ij和 ij为边eij对应的2个对角； 基于期望 法向旋转补洞网格中的三角面片， 计算得到补洞网格中顶点的期望法向之后， 可以进一步 求得三角面片的期望法向， 。

41、三角面片的期望法向是其三个顶点期望法向的平均值， 然后补 洞网格中所有的三角面片根据期望法向进行旋转， 旋转参数计算方法如下： 如图所示， 假设 ni、 ni 和ci为三角面片fi的原始法向、 期望法向和重心位置， ni与ni 的叉乘方向a为三角 面片fi的旋转轴方向， ni与ni 之间的夹角为三角面片fi的旋转角度， 那么三角面片fi将 以ci为旋转中心， 绕旋转轴a旋转角度到新的位置； 基于Poisson方程调整补洞网格顶点 位置， 旋转补洞网格的三角面片会撕裂补洞网格， 因此利用Poisson方程将其重构成连续的 网格曲面， 在建立Poisson方程时先计算撕裂网格的梯度场， 将其作为P。

42、oisson方程的引导 场， 从而进行网格顶点位置的调整，其中f为待求的调整后网格顶点位置， w 为撕裂网格的梯度场； 梯度算子， 假设f表示在每个顶点上的标量， 那么网格域上标量场f在 任意三角面片T内的梯度算子定义如下：其中基函数梯度的表达式是， 表示将向量逆时针旋转90度， AT表示三角片T的面积， 散度算子： 假设w表示在每个三角片上 的向量， 那么网格域上向量场w在顶点xi处的散度算子定义如下： 0096 0097 其中T1(xi)表示顶点xi的1环邻域三角片， AT表示三角片T的面积。 说明书 7/7 页 11 CN 110189352 A 11 图1 说明书附图 1/12 页 1。

43、2 CN 110189352 A 12 图2 说明书附图 2/12 页 13 CN 110189352 A 13 图3 说明书附图 3/12 页 14 CN 110189352 A 14 图4 图5 说明书附图 4/12 页 15 CN 110189352 A 15 图6 图7 说明书附图 5/12 页 16 CN 110189352 A 16 图8 图9 说明书附图 6/12 页 17 CN 110189352 A 17 图10 图11 图12 说明书附图 7/12 页 18 CN 110189352 A 18 图13 图14 图15 说明书附图 8/12 页 19 CN 110189352 A 19 图16 图17 说明书附图 9/12 页 20 CN 110189352 A 20 图18 图19 说明书附图 10/12 页 21 CN 110189352 A 21 图20 图21 说明书附图 11/12 页 22 CN 110189352 A 22 图22 说明书附图 12/12 页 23 CN 110189352 A 23 。