技术领域
[0001] 本
发明涉及颅骨手术技术领域，尤其涉及一种基于颅骨CBCT图像的修复导板生成方法。
背景技术
[0002] 在以对颅顶进行重建来实现对颅骨整形的过程中，最为关键的技术就是对原有颅骨进行重新设计，即分割、
切除和拼合，从而达到整形的效果。当前的颅顶设计依靠三维CT来判断颅顶畸形的状况，并凭着医生的经验在二维CT片上或纸上设计切割形状和路线，然后在手术过程中根据实际情况进行实地调整切割形状、路线及拼接方式。
[0003] 然而，上述常规的方法存在着很多的问题：1、有些颅顶重建不一定需要将所有的颅骨进行切割拼接，只对部分颅骨进行切割分离拼接，因此很难在术前做好方案的设计；2、颅顶重建涉及到三维空间的
定位问题，在二维平面上很难完整的考虑到颅顶重建过程中各骨瓣间的三维空间关系；3、由于缺乏有效手段，手术方案不能有效的转移到手术的实施过程中，对切割下来的骨瓣不能得到最佳的利用，很难得到骨瓣之间准确的拼接关系，因此很难达到术前方案所设计的预期整形结果；4、由于手术方案在手术过程中不能够充分实施，因此手术过程中无法实现术前方案中所规划的对重要
位置及解剖结构的保护；5、缺少有效的辅助工具，手术过程中不能够做到细致、精确，从而会对颅骨出现过度的切割，对骨瓣带来了很大的损伤不易于术后的愈合；6、手术时间长，对患者的血供、血压、代谢产生影响，增加了手术的不确定因素，加大了手术
风险。
[0004] 目前，数字化外科技术发展迅速并广泛应用于外科临床中。其中，能将术前设计精确转术中操作的技术有
手术导航和导板。手术导航需要通过大型的、价格高昂的导航仪将术前或术中的CT、MRI等信息在术中与患者实体数据进行匹配
整理，不仅操作复杂，也有匹配错误导致导航失败的可能；外科导板技术已经在骨科、整形外科、颌面外科有一定的应用，但所使用的患者数据多是来源于CT或
三维扫描，对设备要求较高，而且设计导板的过程也很复杂，需要专
门的工程设计人员，不利于临床推广普及。
[0005] 为此，
现有技术中通过CBCT重建后得到颅骨的重建体数据，以此生成三维空间中分析颅骨的形态学结构，并获取相关的信息，也可以由重建体数据生成
牙齿的全景图像，这样对像就可少做X线全景扫描，从而减少
辐射剂量。但是，现有技术中的大部分产品无法支持基于
口腔CBCT图像的
牙根查看的功能，更不能达到3D视图中数字化诊断的定位效果；同时，在通过CT等扫描方法进行的手术中存在图像与实体匹配程度低，而颅骨的手术对精准度的要求非常高，若手术进行时对颅骨中的相对位置无法判断，将无法进行手术，进而导致医生数字化术前规划过程中无法预判颅骨的情况，继而增加了手术的风险系数和手术的失败率。
发明内容
[0006] 本发明的目的在于提供一种基于颅骨CBCT图像的修复导板方法，提高了颅骨重建数据与患者实体数据的匹配一致性，提高了提取
精度，提高了手术的安全系数和成功率。
[0007] 本发明解决其技术问题所采用的技术方案是：一种基于颅骨CBCT图像的修复导板生成方法，其特征在于，所述方法包括以下步骤：
[0008] S100，光学扫描获得三维头颅模型，建立投影图像及CBCT文件序列；
[0009] S200，通过Marching Cube、区域增长
算法解析CBCT文件序列；
[0010] S300，通过OpenCV边缘提取方法检测颅骨缺损边缘点
云，并利用点云中的对称点生成镜像；
[0011] S400，计算镜像中缺损区域的
曲率，该曲率包括平均曲率、高斯曲率和主曲率；
[0012] S500，基于S400中的曲率执行区域增长算法形成首期修复导板；
[0013] S600，基于Poisson方程补洞平滑方法对首期修复导板边缘进行优化，形成最终的修复导板。
[0014] 优选地，所述Marching Cube的解析过程包括以下步骤：
[0015] S201，将原始数据经过预处理之后，读入
指定的数组中；
[0016] S202，从网格数据体中提取一个单元体，成为当前单元体，同时获取该单元体的所有信息，其信息包括边界的
顶点函数值、和单元体的点云坐标位置；
[0017] S203，将当前单元体中顶点的函数值与给定等值面值C进行比较，得出该单元体的状态表；
[0018] S204，根据当前单元体的状态表索引，找出与等值面相交的单元体棱边，并采用线性插值的方法，计算出各个交点的位置坐标；
[0019] S205，利用中心差分法，求出当前单元体中顶点的法向量，再采用线性插值的方法，得到三
角面片各个顶点的法向量；
[0020] S206，根据各个三角面片顶点的坐标，顶点法向量进行等值面图象的绘制。
[0021] 优选地，所述S200中执行区域增长算法时，包括将
阈值控制在770-25584区间，重建出3D模型查看CBCT主要的解剖结构。
[0022] 优选地，所述S300中的OpenCV边缘提取方法包括以下过程：
[0023] 通过离散化的高斯函数产生一组归一化的高斯核，而后基于高斯核函数对图像数据中灰度矩阵的每一点进行加权求和；
[0024] 通过增强算法对图像数据中灰度点邻域强度值有显著变化的点进行凸显；
[0025] 通过阈值化方法对噪音边缘点进行剔除。
[0026] 优选地，所述S300中的对称点生成镜像包括特征点的提取，其提取公式如下：
[0027]
[0028] 其中，pi,pj是两个对称点，hi,hj是归一化后的结果，Ci,j是所得的向量，提取Ci,j最大值以获得特征点。
[0029] 优选地，所述S400中缺损区域的曲率均值通过以下公式取得：
[0030]
[0031] 其中，P,Q为特征点的曲率均值，C为形态距离，T为TPS再预估变换矩阵权重。
[0032] 优选地，所述S500包括区域增长算法，其步骤如下：
[0033] S501，选取
种子点p(x0，y0)，用堆栈表示种子区域，将种子点push到种子堆栈中；
[0034] S502，将种子堆栈中第一个种子点pop出堆栈，并以该点为中心，遍历该中心8邻域
像素；
[0035] S503，判断遍历像素点是否已经在种子区域中，如果否，则判断遍历像素点是否满足相邻种子点相似性，如果该像素点(x，y)满足相似性，则将(x，y)push到堆栈中；
[0036] S504，重复S502-S503，直至种子堆栈为空。
[0037] 优选地，所述步骤S600包括根据三角网格中孔洞边界生成一个初始化补洞网格，而后通过法向估算和Poisson方程来修正补洞网格中三角面片的几何形状，使其能够适应并与周围的原始网格融合，其过程包括以下步骤：
[0038] S601，检测孔洞边界并初始化补洞网格；
[0039] S602，计算补洞网格中顶点的期望法向，构建Laplace方程Δf＝0求解补洞网格内部顶点的法向分布，其Laplace算子为：
[0040]
[0041] 其中，N1(xi)表示顶点xi的1环邻域点，αij和βij为边eij对应的2个对角；
[0042] S603，基于期望法向旋转补洞网格中的三角面片；
[0043] S604，基于Poisson方程调整补洞网格顶点位置，并计算
三维网格的梯度场，公式为：
[0044] 其中，f为待求的调整后网格顶点位置，w为撕裂网格的梯度场；
[0045] 梯度算子为：
[0046] 其中，基函数梯度▽Φi的表达式为：
[0047] 其表示将向量逆
时针旋转90度，AT表示三角片T的面积；
[0048] 散度算子：
[0049] 其中，T1(xi)表示顶点xi的1环邻域三角片，AT表示三角片T的面积。
[0050] S605，基于S400中得到的缺损区域的曲率进行
迭代区域增长，得到3D点云面片。
[0051] 本发明的有益效果在于：本发明通过对颅骨形态的提取，形成模拟图像，帮助医生提前预测出术中及术后的颅骨趋势，极大的提高了手术的精确性，降低了手术风险。同时，本发明通过边缘提取方法检测颅骨缺损边缘点云，并计算出镜像中缺损区域的曲率，进而通过该曲率执行区域增长算法形成首期修复导板，最后还对首期修复导板进行了优化处理，以提高图像与实际形态的一致性。此模拟图像不仅可以用于数字化模拟术前规划，更可以通过3D打印产生实体用于术前模拟测试，拟真地预先展现出术中所出现的情况。本发明操作流程简单，医生在数字化术前规划和诊断过程中，对颅骨的修复结构进行提前的预判，使得术前规划和术中情况一致，同时使手术的规划精准程度得到量化，进而提高了医生数字化程度的意识，以及增强了用户体验。另外，经
发明人试验，这样的处理方法同样适用于其他医疗科室相关的CBCT图像、MRI图像中的精细解剖结构的提取过程。
附图说明
[0052] 为了更清楚地说明本发明
实施例或现有技术中的技术方案，下面将结合附图及实施例对本发明作进一步说明，下面描述中的附图仅仅是本发明的部分实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图：
[0053] 图1是本发明一种基于颅骨CBCT图像的修复导板生成方法中步骤S100的扫描模型示意图；
[0054] 图2是本发明一种基于颅骨CBCT图像的修复导板生成方法中步骤S100的带缺损区域的颅骨示意图；
[0055] 图3是本发明一种基于颅骨CBCT图像的修复导板生成方法中步骤S200利用点云中的对称点生成的镜像示意图；
[0056] 图4是本发明一种基于颅骨CBCT图像的修复导板生成方法中镜像中缺损区域的曲面示意图；
[0057] 图5是本发明一种基于颅骨CBCT图像的修复导板生成方法中首期修复导板示意图；
[0058] 图6是本发明一种基于颅骨CBCT图像的修复导板生成方法中镜像中最终的修复导板示意图；
[0059] 图7是本发明一种基于颅骨CBCT图像的修复导板生成方法中点云面片的示意图；
[0060] 图8是本发明一种基于颅骨CBCT图像的修复导板生成方法中高斯-博内定理结果的示意图。
具体实施方式
[0061] 为了使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明的部分实施例，而不是全部实施例。基于本发明的实施例，本领域普通技术人员在没有付出创造性劳动的前提下所获得的所有其他实施例，都属于本发明的保护范围。
[0062] 在实施例1中，如图1、2、3、4、5及6所示，一种基于颅骨CBCT图像的修复导板生成方法，其特征在于，所述方法包括以下步骤：首先，通过解析口腔CBCT的文件序列(dcm格式)，并执行区域增长算法，同时将CBCT图像
阀值控制在770-25584区间，进而重建出CBCT的主要解剖结构(3D视图)；而后通过OpenCV边缘提取方法检测颅骨缺损边缘点云，并利用点云中的对称点生成镜像；进而计算镜像中缺损区域的曲率，通过该曲率执行区域增长算法形成首期修复导板；最后基于Poisson方程补洞平滑方法对首期修复导板边缘进行优化，形成最终的修复导板。
[0063] 具体地，对于步骤S200的Marching Cube方法和区域增长算法，Marching Cube利用等值面，即三维立体空间中所有具有相同值的点的集合，其集合区域展现类比如地形图里的等高线。考虑到Marching Cube中的基本假设条件，即空间中存在立方体(六面体)，而沿六面体边的数据场是呈连续性变化，即当立方体一条边的两个顶点分别大于或小于等值面的值，则在该条边上有且仅有一点是这条边与等值面的交点。在空间中过程，即是用大量小正方体对空间进行切分，进而用小正方体内部的平面来近似表示当前的等值面，故此，小正方体的数量越多，逼近的效果越好，随之进行的是大量的计算。计算式子包括{(x,y,z)|f(x,y,z)＝c}，其中，,c为常数逐个处理数据场中的立方体(
体素)，其计算过程包括分离出与等值面相交的立方体后，采用双线性插值计算出等值面与立方体边的交点，根据立方体每一顶点与等值面的相对位置，将等值面与立方体边的交点按三维空间中垂直地方式连接生成等值面，从而作为等值面在该立方体内的一个逼近表示。
[0064] 进一步地，执行步骤S201，将原始数据经过预处理之后，读入指定的数组中；执行步骤S202，从网格数据体中提取一个单元体，成为当前单元体，同时获取该单元体的所有信息，其信息包括边界的顶点函数值、和单元体的点云坐标位置，在本实施例中，以8个顶点为例；执行步骤S203，将当前单元体8个顶点的函数值与给定等值面值C进行比较,得到该单元体的状态；(edgeTable、triTable)；执行步骤S204，根据当前单元体的状态表索引，找出与等值面相交的单元体棱边，并采用线性插值的方法，计算出各个交点的位置坐标；执行步骤S205，利用中心差分法，求出当前单元体8个顶点的法向量，在采用线性插值的方法，得到三角面片各个顶点的法向；执行步骤S206，根据各个三角面片顶点的坐标，顶点法向量进行等值面图象的绘制。
[0065] 具体地，在中心差分法计算过程中，例如采用等时间步长，Δ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)。
[0066] 其中，对于步骤S200中的区域增长算法，获取分割区域中的一个种子点(即特征点)，在种子点的周围搜索与该种子点有相似性质的像素点，并将其合并到种子区域中，而后将合并的像素作为新的种子点继续搜索，直到种子区域中所有像素周围没有相似的像素点，算法结束。其基本算法流程为：步骤1，选取种子点p(x0，y0)，用堆栈表示种子区域，将种子点push到种子堆栈中；步骤2，将种子堆栈中第一个种子点pop出堆栈，并以该点为中心，遍历该中心8邻域像素；判断遍历像素点是否已经在种子区域中，如果否，判断遍历像素点是否满足相邻种子点相似性，如果像素点(x，y)满足相似性，将(x，y)push到堆栈中；步骤4，重复步骤2至3，直至种子堆栈为空。由此可见，影响区域生长算法的要素有三个：种子点的选取；搜索路径的选择；像素相似性的判断。
[0067] 其中，对于种子点的选取，区域生长算法是半交互式的分割算法，需要用户选取种子点，也可以是通过其他算法计算出来的种子点；对于搜索路径的选择，其搜索路径是选择相邻的像素，以二维图像为例，一般为8邻域搜索，或者4邻域搜索，以三维图像为例，一般为26邻域搜索或者6邻域搜索；对于像素相似性的判断，其相似性以像素值的相近程度为判断标准，例如，设置一定灰度范围做为相似的标准，也可以通过计算满足某种形状或者性质作为判断标准。
[0068] 进一步地，所述步骤S300中的OpenCV边缘提取方法，其提取步骤如下：
[0069] 滤波：
[0070]
边缘检测的算法主要是基于图像强度的一阶和二阶导数，但导数通常对噪声很敏感，因此必须采用
滤波器来改善与噪声相关的边缘检测器的性能。例如采用高斯滤波方法进行滤波，即采用离散化的高斯函数产生一组归一化的高斯核，然后基于高斯核函数对图像灰度矩阵的每一点进行加权求和。
[0071] 增强：
[0072] 增强边缘的
基础是确定图像各点邻域强度的变化值。增强算法可以将图像灰度点邻域强度值有显著变化的点凸显出来。在具体编程实现时，通过计算相应的梯度幅值来确定。
[0073] 检测：
[0074] 经过增强的图像，其邻域中存在较多比较大的梯度值，而在特定的应用中，这些梯度值并不是要找的边缘点，所以需要对这些梯度值进行取舍。实际工程中，可通过阈值化方法对边缘点进行检测，以提取出边缘点。
[0075] 本实施例采用的是基于Laplacian算子–Laplacian函数的方法进行边缘提取，其函数调用形式为：
[0076] void Laplacian(inputArray,outputArray,int ddepth,int ksize＝1,double scale＝1,double delta＝0,int borderType＝BORDER_DEFAULT)。
[0077] 进一步地，所述S300中的对称点生成镜像包括特征点的提取，其提取公式如下：
[0078]
[0079] 其中，pi,pj是两个对称点，hi,hj是归一化后的结果，Ci,j是所得的向量，提取Ci,j最大值以获得特征点。
[0080] 对于对称点的选取，主要自动
选定为颅骨中的两个眼眶最远特征点和鼻梁垂直法线部分，同时自动复制点云矩阵，进行三角化生成镜像，在3D空间中取眼眶向量中Ci，j最大值即可获得特征点。
[0081] 进一步地，所述步骤S400中的曲率包括平均曲率、高斯曲率和主曲率步骤如下：
[0082] 平均曲率的计算原理如下：
[0083] 平均曲率与曲面的单位法向量相关，该法向量的选取影响曲率的正负号。曲率的符号取决于法向量的方向：如果曲面“远离”法向量则曲率是正的。上述定义对3D空间中任何方式定义的曲面均成立，其主要用于计算单位法向量的散度。
[0084] 对曲面是两个坐标的函数定义的曲面，例如z＝S(x,y)，使用向下的法向量平均曲率(的两倍)表示为
[0085]
[0086] 高斯曲率的计算原理如下：
[0087] 负曲率曲面上的三角形三角之和小于平面三角形的三角之和。
[0088] 测地三角形(即黎曼球面几何中的三角形)的总曲率等于它的内角和与的差。正曲率曲面上的三角形的内角和大于180°，而负曲率曲面上的三角形的内角和小于180°。零曲率曲面上(如欧几里得平面)，其内角和等于π。
[0089]
[0090] 如图1所示，更一般的结果是高斯-博内定理。由此可知，负曲率曲面上的三角形三角之和小于平面三角形的三角之和。
[0091] 主曲率的计算原理如下：
[0092] 对曲面S:r＝r(u,v)上一给定点P0(u0,v0)，法曲率kn是切方向du:dv的函数，称法曲率的每个临界值(critical value)为曲面在这一点的主曲率；对应的方向称为曲面在这一点的主方向。采用下列定理进行计算：
[0093] 曲面上一点由方程Pdu2+2Qdudv+Rdv2＝0所确定的两个切方向互相垂直的充要条件是ER-2FQ+GP＝0
[0094] 其中，E、F、G是曲面的第一类基本量。两个方向du:dv和δu:δv
正交的充要条件是Eduδu+F(duδv+dvδu)+Gdvδv＝0，即：
[0095]
[0096] 将已知的二次方程写为
[0097]
[0098] 其中， 是它的两个根，且均应满足上述方程，由根与系数的关系知
[0099]
[0100] 进一步地，所述S500包括区域增长算法，其步骤如下：
[0101] S501，选取种子点p(x0，y0)，用堆栈表示种子区域，将种子点push到种子堆栈中；
[0102] S502，将种子堆栈中第一个种子点pop出堆栈，并以该点为中心，遍历该中心8邻域像素；
[0103] S503，判断遍历像素点是否已经在种子区域中，如果否，则判断遍历像素点是否满足相邻种子点相似性，如果该像素点(x，y)满足相似性，则将(x，y)push到堆栈中；
[0104] S504，重复S502-S503，直至种子堆栈为空。
[0105] 具体地，影响区域生长算法的要素有三个：种子点的选取；搜索路径的选择；像素相似性的判断。
[0106] 对于种子点的选取：一般情况下，区域生长算法是半交互式的分割算法，需要用户选取种子点，也可以是通过其他算法计算出来的种子点。
[0107] 对于搜索路径的选择：搜索路径一般选择相邻的像素，以二维图像为例，一般为8邻域搜索，或者4邻域搜索；以三维图像为例，一般为26邻域搜索或者6邻域搜索。
[0108] 对于像素相似性的判断：相似性一般以像素值的相近程度为判断标准，例如，可以设置一定灰度范围做为相似的标准；也可以通过计算满足某种形状或者性质作为判断标准。
[0109] 进一步地，在步骤S500中，对于Poisson方程的三角网格补洞方法，该算法首先需要根据孔洞边界生成一个初始化补洞网格，然后通过法向估算和Poisson方程来修正补洞网格中三角面片的几何形状，使其能够适应并与周围的原始网格融合。算法的主要步骤如下：如图所示，检测孔洞边界并初始化补洞网格，由于初始化补洞网格无法与原始孔洞周围的网格有效融合，因此需要调整补洞网格的顶点位置使得补洞网格与原始网格之间光滑过渡；计算补洞网格中顶点的期望法向，由于已知原始网格孔洞边界的法向，将其作为补洞网格边界的法向，构建Laplace方程Δf＝0求解补洞网格内部顶点的法向分布，假设f表示在每个顶点上的标量，那么网格域上在顶点xi处的Laplace算子定义如下(不考虑面积影响)：
[0110]
[0111] 其中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内的梯度算子定义如下： 其中基函数梯度▽Φi的表达式是
±表示将向量逆时针旋转90度，AT表示三角片T的面积，散度算子：假设w表
示在每个三角片上的向量，那么网格域上向量场w在顶点xi处的散度算子定义如下：
[0112]
[0113] T1(xi)表示顶点xi的1环邻域三角片，AT表示三角片T的面积；
[0114] S605，基于S400中得到的缺损区域的曲率进行迭代区域增长，得到3D点云面片。具体的，如图1所示，基于步骤S400中得到的最佳轮廓曲率进行迭代区域增长，以得到最优的3D点云面片，点云面片生成公式为
[0115]
[0116] 其中，A(v)的面积视三角形顶角的大小而定，θ表示所计算区域的三角顶角角度。
