本发明涉及口腔cbct图像处理
技术领域：
：，尤其涉及一种基于口腔cbct图像的牙根提取方法。
背景技术：
：：目前临床中主要有两种牙科全景图像生成方法：牙科全景x射线摄影方法和口腔cbct重建后得到体数据的方法。牙科全景x射线摄影方法主要采用x摄像和平板探测器以一种特殊的轨迹进行运动，在运动过程中采集特定位置的图像，并最终将采集到的序列图像组合成牙科全景图像，该图像是有设备直接采集获取，对人体具有一定的辐射，该方法缺点是只能获取牙齿的二维投影数据，无法获取牙齿三维空间的形态学信息。口腔cbct重建后得到体数据的方法，能够方便医生在三维空间中分析牙齿的形态学结构，并获取相关的信息，也可以由重建体数据生成牙齿的全景图像，这样对像就可少做x线全景扫描，从而减少辐射剂量，这种方式就需要由重建体数据从而生成牙科全景图像。但是，现有技术中的大部分产品无法支持基于口腔cbct图像的牙根查看的功能，更不能达到3d视图中数字化诊断的定位效果，这样导致医生数字化术前规划过程中无法预判牙根的情况，进而增加手术的风险和手术的失败率。例如在口腔根管根尖定位过程中，提前预知牙根和器械的相对位置将决定手术的成败；种植修复过程中，邻牙牙根的持续监控和判断将决定手术种植体的生命周期；正畸过程中需要持续对牙根的走势提前预判并控制，决定了正畸治疗的时间和最终效果。近年来，口腔修复、种植及正畸诊断分析过程中对于牙根的查看越发受到重视，现有技术中又缺乏一种精准牙根提取、并使其3d可视化以及建模测试的方法，无法降低相关手术的难度和风险，尤其是在精度要求很高的复杂手术中，牙根的相对位置无法预判的情况下，手术无法进行。综上所述，亟待一种解决上述问题的基于口腔cbct图像的牙根提取方法。技术实现要素：本发明提供了一种基于口腔cbct图像的牙根提取方法，降低了手术风险，提高了手术成功率。本发明解决其技术问题所采用的技术方案是：一种基于口腔cbct图像的牙根提取方法，包括以下步骤：s100，光学扫描获得三维牙颌和牙列模型，建立投影图像及cbct文件序列；s200，通过marchingcube、区域增长算法解析cbct文件序列；s300，基于cbct文件序列，对上下牙颌进行区分，其区分过程包括通过vtk定位咬合面和矢状面的关键帧，执行限制性区域增长算法；s400，优化上下牙颌的分割效果，并对牙颌模型进行手动裁剪；s500，去除裁剪后牙颌模型的噪音图像，并对其进行平滑。进一步地，所述marchingcube的实现包括以下步骤：s201，将原始数据经过预处理之后，读入指定的数组中；s202，从网格数据体中提取一个单元体，成为当前单元体，同时获取该单元体的所有信息，其信息包括边界的顶点函数值、和单元体的点云坐标位置；s203，将当前单元体中顶点的函数值与给定等值面值c进行比较，得出该单元体的状态表；s204，根据当前单元体的状态表索引，找出与等值面相交的单元体棱边，并采用线性插值的方法，计算出各个交点的位置坐标；s205，利用中心差分法，求出当前单元体中顶点的法向量，再采用线性插值的方法，得到三角面片各个顶点的法向量；s206，根据各个三角面片顶点的坐标，顶点法向量进行等值面图象的绘制。进一步地，所述步骤s200执行区域增长算法时，包括将阈值控制在770-25584区间，重建出3d模型查看cbct主要的解剖结构。进一步地，所述步骤s300中的区域增长算法包括其输入阈值和拼接向量的改变，其拼接向量公式如下：o1＝{(x,y)f(x,y)≤s3}；o2＝{(x,y)s3<f(x,y)≤(s3+σ)∧g(x,y)>t3}；o＝o1∪o2＝{(x,y)f(x,y)≤s3}∪{(x,y)s3<f(x,y)≤(s3+σ)∧g(x,y)>t3}；bwf＝1,(x,y)∈o,0,(x,y)∈其他；其中，(s3,t3)为二维阈值向量，s3为灰度阈值，t3为梯度阈值，将阈值小于s3像点作为目标像点集；σ为变量参数，其大小与图像模糊程度对应；bwf为ct中的灰度图像值。进一步地，所述步骤s400包括建立算法和搜索算法，构建算法的步骤如下：s401，按xx轴的分割线对空间进行分割；s402，计算所有点的xx坐标的平均值，选择所有点中最接近平均值的点作为分割线，把空间中的点进行分割；s403，对分隔后的空间按照yy轴的分割线进行分割；s404，计算所有点的xx坐标的平均值，选择所有点中最接近平均值的点作为分割线，把空间中的点进行分割；s405，对分割后的空间继续按xx轴进行分割，依次类推，以xx轴和yy轴为目标对象循环分割，每个空间中分割至一个点时，结束分割。进一步地，所述步骤s401至步骤s405的分割过程对应于一个二叉树，其每条分割线对应于二叉树中的一个分支，其空间中每个点对应于一个叶子节点。进一步地，所述二叉树中包括空间点(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，若上述步骤中的圆与三维空间中及圆直径长度一致的立方体的分割线有交点，则对上述步骤中所比较最近距离的点进行遍历，比较遍历点与其的距离，继续回溯，并以(x，y)为圆心，遍历点与其的距离为半径画圆，当所画之圆与分割线没有交点时，即该遍历点为最近邻点，查找结束。进一步地，所述步骤s400还包括根据k-d树进行无序点云去噪，去噪过程包括如下步骤：s6，根据点云数据生成k-d树，建立整个3d图形的点云三维坐标及其在三维空间中相互关联的拓扑关系；s7，查找任一点的的邻域；s8，计算该点与邻域内各点的距离，并取其平均值；s9，判断该平均值是否超过阈值，若超过则判定该点为噪点，进行去除。进一步地，所述步骤s500包括通过k-d树建立之后进行无序点云的去噪算法，还包括基于poisson方程的三角网格补洞算法进行平滑。进一步地，所述步骤s500包括根据三角网格中孔洞边界生成一个初始化补洞网格，而后通过法向估算和poisson方程来修正补洞网格中三角面片的几何形状，使其能够适应并与周围的原始网格融合，其过程包括以下步骤：s501，检测孔洞边界并初始化补洞网格；s502，计算补洞网格中顶点的期望法向，构建laplace方程求解补洞网格内部顶点的法向分布，其laplace算子为：其中，n1(xi)表示顶点xi的1环邻域点，αij和βij为边eij对应的2个对角；s503，基于期望法向旋转补洞网格中的三角面片；s504，基于poisson方程调整补洞网格顶点位置，并计算撕裂网格的梯度场，公式为：其中，f为待求的调整后网格顶点位置，w为撕裂网格的梯度场；梯度算子为：其中基函数梯度的表达式为：其表示将向量逆时针旋转90度，at表示三角片t的面积；散度算子：其中，t1(xi)表示顶点xi的1环邻域三角片，at表示三角片t的面积。本发明的有益效果在于：本发明通过对牙根形态的提取，形成模拟图像，帮助医生提前预测出术中甚至术后的牙根趋势，极大的提高了手术的精确性，降低了手术风险。同时，本发明通过优化、去噪和平滑等图像处理，提取出图像质量较高的牙根模型，此模型不仅可以用于数字化模拟术前规划，更可以通过3d打印产生实体用于术前模拟测试，拟真地预先展现出术中所出现的情况。本发明操作流程简单，医生在数字化术前规划和诊断过程中，对牙根的解剖结构进行提前的预判，使得术前规划和术中情况一致，同时使手术的规划精准程度得到量化，进而提高了医生数字化程度的意识，以及增强了用户体验。另外，经发明人试验，这样的处理方法同样适用于其他医疗科室相关的cbct图像、mri图像中的精细解剖结构的提取过程。附图说明为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将结合附图及实施例对本发明作进一步说明，下面描述中的附图仅仅是本发明的部分实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图：图1为本发明实施例1中基于口腔cbct图像的牙根提取方法的步骤流程图；图2为本发明实施例1中marchingcube实现的步骤流程图；图3为本发明实施例1中构建建立算法和搜索算法的步骤流程图；图4为本发明实施例1中空间最邻近点查找方法的步骤流程图；图5为本发明实施例1中去噪步骤流程图；图6为本发明实施例1中空间切分示意图；图7为本发明实施例1中空间最邻近点搜索过程示意图；图8为本发明实施例1中二叉树示意图；图9为本发明实施例1中空间分割线交点示意图；图10为本发明实施例1中空间分割线交点示意图；图11为本发明实施例1中孔洞边界示意图；图12为本发明实施例1中补洞网络示意图；图13为本发明实施例1中三角面片示意图；图14为本发明实施例1中调整补洞网格的结构示意图；图15为本发明实施例1中调整补洞网格的结构示意图；图16为本发明实施例1中调整补洞网格的结构示意图；图17为本发明实施例1中cbct的3d建模结构示意图；图18为本发明实施例1中关键帧的定位示意图；图19为本发明实施例1中限制性区域增长效果示意图；图20为本发明实施例1中meshlab手动裁剪示意图；图21为本发明实施例1中裁剪效果示意图；图22为本发明实施例1中最终裁剪效果示意图。具体实施方式为了使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明的部分实施例，而不是全部实施例。基于本发明的实施例，本领域普通技术人员在没有付出创造性劳动的前提下所获得的所有其他实施例，都属于本发明的保护范围。在实施例1中，如图1、2、3、4、17、18、19、20、21及22所示，首先，通过解析口腔cbct的文件序列(dcm格式)，并执行区域增长算法，同时将cbct图像阀值控制在770-25584区间，进而重建出cbct的主要解剖结构(3d视图)；而后基于vtk(第三方库)中针对咬合面和矢状面，对关键帧进行交叉定位，递归执行限制性区域增长算法，同时调整cbct图像阀值在1565-9925区间，进行k-dtree算法做特征提取，利用meshlab的手动裁剪功能进行粗略去噪；最后调用k-dtree算法进行精细化去噪，并采用基于poisson方程的三角网格补洞方法对模型进行表面处理。具体地，如图6所示，对于步骤s200的marchingcube方法和区域增长算法，marchingcube利用等值面，即三维立体空间中所有具有相同值的点的集合，其集合区域展现类比如地形图里的等高线。考虑到marchingcube中的基本假设条件，即空间中存在立方体(六面体)，而沿六面体边的数据场是呈连续性变化，即当立方体一条边的两个顶点分别大于或小于等值面的值，则在该条边上有且仅有一点是这条边与等值面的交点。在空间中过程，即是用大量小正方体对空间进行切分，进而用小正方体内部的平面来近似表示当前的等值面，故此，小正方体的数量越多，逼近的效果越好，随之进行的是大量的计算。计算式子包括{(x,y,z)|f(x,y,z)＝c}，其中，,c为常数逐个处理数据场中的立方体(体素)，其计算过程包括分离出与等值面相交的立方体后，采用双线性插值计算出等值面与立方体边的交点，根据立方体每一顶点与等值面的相对位置，将等值面与立方体边的交点按三维空间中垂直地方式连接生成等值面，从而作为等值面在该立方体内的一个逼近表示。进一步地，如图2所示，执行步骤s201，将原始数据经过预处理之后，读入指定的数组中；执行步骤s202，从网格数据体中提取一个单元体，成为当前单元体，同时获取该单元体的所有信息，其信息包括边界的顶点函数值、和单元体的点云坐标位置，在本实施例中，以8个顶点为例；执行步骤s203，3.将当前单元体8个顶点的函数值与给定等值面值c进行比较,得到该单元体的状态；(edgetable、tritable)；执行步骤s204，根据当前单元体的状态表索引，找出与等值面相交的单元体棱边，并采用线性插值的方法，计算出各个交点的位置坐标；执行步骤s205，利用中心差分法，求出当前单元体8个顶点的法向量，在采用线性插值的方法，得到三角面片各个顶点的法向；执行步骤s206，根据各个三角面片顶点的坐标，顶点法向量进行等值面图象的绘制。具体的，在中心差分法计算过程中，例如采用等时间步长，δt(i)＝δt(δt为常数)，用u表示位移，那么速度和加速度的中心差分近似为：u'(i)＝[u(i+1)-u(i-1)]/(2δt)；u”(i)＝[u(i+1)-2u(i)+u(i-1)]/(δt*δt)。其中，对于步骤s200中的区域增长算法，获取分割区域中的一个种子点(即特征点)，在种子点的周围搜索与该种子点有相似性质的像素点，并将其合并到种子区域中，而后将合并的像素作为新的种子点继续搜索，直到种子区域中所有像素周围没有相似的像素点，算法结束。其基本算法流程为：步骤1，选取种子点p(x0，y0)，用堆栈表示种子区域，将种子点push到种子堆栈中；步骤2，将种子堆栈中第一个种子点pop出堆栈，并以该点为中心，遍历该中心8邻域像素；判断遍历像素点是否已经在种子区域中，如果否，判断遍历像素点是否满足相邻种子点相似性，如果像素点(x，y)满足相似性，将(x，y)push到堆栈中；步骤4，重复步骤2至3，直至种子堆栈为空。由此可见，影响区域生长算法的要素有三个：种子点的选取；搜索路径的选择；像素相似性的判断。具体地，对于种子点的选取，区域生长算法是半交互式的分割算法，需要用户选取种子点，也可以是通过其他算法计算出来的种子点；对于搜索路径的选择，其搜索路径是选择相邻的像素，以二维图像为例，一般为8邻域搜索，或者4邻域搜索，以三维图像为例，一般为26邻域搜索或者6邻域搜索；对于像素相似性的判断，其相似性以像素值的相近程度为判断标准，例如，设置一定灰度范围做为相似的标准，也可以通过计算满足某种形状或者性质作为判断标准。具体地，对于拼接向量公式，目标集合：o1＝{(x,y)f(x,y)≤s3}；由于足迹图像中噪声和模糊边缘像素灰度值要大于s3，所以若用单阈值s3对图像分割，噪声和边缘像点将被错划为背景，此外，噪声和ct图像中模糊边缘像素点具有较大梯度值,所以用二维阈值向量(s3,t3)限制这部分像素,即把集合o2＝{(x,y)s3<f(x,y)≤(s3+σ)∧g(x,y)>t3}也作为目标集合，因此，图像分割后目标像素集应为o＝o1∪o2＝{(x,y)f(x,y)≤s3}∪{(x,y)s3<f(x,y)≤(s3+σ)∧g(x,y)>t3}。其中，参数σ选取与图像模糊程度有关，在足迹图像中σ取5或7即可满足要求。令分割后二值图像，即像素值为0-255的灰度图像为bwf,则bwf＝1,(x,y)∈o,0,(x,y)∈其他。具体地，如图7及8所示，步骤s300无法完全性地区分上下颌的所有点云，需要采取k-dtree算法和手动裁剪的方法提取牙根，在步骤s400中，在包括k-dtree的建立算法和搜索算法的同时，其手动裁剪还使用了meshlab的固有功能。在执行步骤s401至步骤s405的过程中，其搜索过程如图所示，在本实施例中，建立起的二叉树包括如下一组二维点集：{(2,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)为最近邻点，进行回溯，计算(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,3)，查找结束。具体地，如图11、12、13、14、15及16所示，手动裁剪过后依然会有不少噪音，执行k-dtree算法去噪去噪完成后，表面依然存在粗糙情况，需要执行平滑算法进行优化，在步骤s500中，对于poisson方程的三角网格补洞方法，该算法首先需要根据孔洞边界生成一个初始化补洞网格，然后通过法向估算和poisson方程来修正补洞网格中三角面片的几何形状，使其能够适应并与周围的原始网格融合。算法的主要步骤如下：如图所示，检测孔洞边界并初始化补洞网格，由于初始化补洞网格无法与原始孔洞周围的网格有效融合，因此需要调整补洞网格的顶点位置使得补洞网格与原始网格之间光滑过渡；计算补洞网格中顶点的期望法向，由于已知原始网格孔洞边界的法向，将其作为补洞网格边界的法向，构建laplace方程δf＝0求解补洞网格内部顶点的法向分布，假设f表示在每个顶点上的标量，那么网格域上在顶点xi处的laplace算子定义如下(不考虑面积影响)：其中n1(xi)表示顶点xi的1环邻域点，αij和βij为边eij对应的2个对角；基于期望法向旋转补洞网格中的三角面片，计算得到补洞网格中顶点的期望法向之后，可以进一步求得三角面片的期望法向，三角面片的期望法向是其三个顶点期望法向的平均值，然后补洞网格中所有的三角面片根据期望法向进行旋转，旋转参数计算方法如下：如图所示，假设ni、ni’和ci为三角面片fi的原始法向、期望法向和重心位置，ni与ni’的叉乘方向a为三角面片fi的旋转轴方向，ni与ni’之间的夹角φ为三角面片fi的旋转角度，那么三角面片fi将以ci为旋转中心，绕旋转轴a旋转角度φ到新的位置；基于poisson方程调整补洞网格顶点位置，旋转补洞网格的三角面片会撕裂补洞网格，因此利用poisson方程将其重构成连续的网格曲面，在建立poisson方程时先计算撕裂网格的梯度场，将其作为poisson方程的引导场，从而进行网格顶点位置的调整，其中f为待求的调整后网格顶点位置，w为撕裂网格的梯度场；梯度算子，假设f表示在每个顶点上的标量，那么网格域上标量场f在任意三角面片t内的梯度算子定义如下：其中基函数梯度的表达式是，±表示将向量逆时针旋转90度，at表示三角片t的面积，散度算子：假设w表示在每个三角片上的向量，那么网格域上向量场w在顶点xi处的散度算子定义如下：其中t1(xi)表示顶点xi的1环邻域三角片，at表示三角片t的面积。当前第1页12当前第1页12