一种肋骨可视化的肋骨骨折辅助诊断方法与流程

本发明是关于图像处理技术领域，特别涉及一种肋骨可视化的肋骨骨折辅助诊断方法。

背景技术：

计算机断层扫描(ct)是诊断肋骨骨折的主要方法，但肋骨特殊的组织结构导致肋骨遮挡因而微小骨折及骨裂容易漏诊，精确的诊断需要经验的积累，常规的二维方法需要手动追踪多个切片，传统的三维可视化方法在诊断时需要对图像进行旋转观察从而避免肋骨遮挡的影响。常用重建方法有多层面重组法、最大密度投影法及表面阴影显示法。多层面重组法需要做一系列图像且非目标器官存在形变；最大密度投影法低密度组织器官易被遗漏且噪声影响大；表面阴影显示法容易产生伪影且看不出密度信息。

技术实现要素：

本发明的主要目的在于克服现有技术中的不足，提供一种通过把肋骨投影到二维平面上，消除在临床诊断时肋骨遮挡导致微小骨折及骨裂的骨折辅助诊断方法。为解决上述技术问题，本发明的解决方案是：

提供一种肋骨可视化的肋骨骨折辅助诊断方法，包括下述过程：

一、获取图片；

二、提取肋骨；

三、人工修正，去掉脊柱、胸骨；

四、展开肋骨；

所述过程一是指：取得原始的人体胸部ct扫描图像；

所述过程二是指：将人体胸部ct扫描图像转化为256阶灰度图像，然后用增强滤波器对肋骨进行增强处理，调整像素值大小取一个范围提取出整个胸部骨骼，整个胸部骨骼包括肋骨、脊柱及胸骨，按照肋骨空间位置大小从上向下标记左右肋骨分支；

所述过程三是指：创建三维窗口，用户通过鼠标右键点击划圈标示出不属于肋骨的部分，计算机根据用户的标示将圈内部分(脊柱、胸骨)去除；

所述过程四是指：泊松曲面重建肋骨轮廓，拉普拉斯网格收缩提取肋骨骨架；应用最近点迭代矫正骨架，对准骨架位置；应用曲面重建算法(curvedplanarreformation，cpr)对每一只肋骨进行展开，得到最终的肋骨可视化展开结果，用于辅助医生进行肋骨病变诊断。

在本发明中，所述过程二具体包括下述步骤：

步骤a：把人体胸部ct扫描图像转化为256色灰度图像

步骤b：然后使用一个增强滤波算子对肋骨区域进行强化，从图像中逐点提取平滑的黑塞矩阵：

其中，所述i是指每点的二阶导数矩阵；所述iij是指二阶导数，i、j分别取x、y或者z；

步骤c：黑塞矩阵i通过高斯核的二阶偏导数与原图像进行卷积求得：

其中，gσ(v)＝exp(-v^tv/σ²)，v∈r³为高斯核；所述δ为拉普拉斯算符；所述为步骤a中的256色灰度图像；所述σ是指当前高斯核的标准差；所述t是指向量的转置；所述exp是指自然指数函数；所述*指的是函数间的卷积；所述r³为三维欧氏空间；所述v为欧氏空间中的一个三维向量；

步骤d：求出黑塞矩阵i的特征值λ1，λ2，λ3，且|λ1|＜|λ2|＜|λ3|；

计算每点的增强值：

其中，所述v0(i)用于表示像素值类似血管(vesselness)的程度；所述α、β、c都分别大于0(均为可调整经验参数，可根据实际灰度图像自行调节)；

步骤e：通过调整像素值范围提取出胸部骨骼，包含肋骨、脊柱及胸骨；

步骤f：按照肋骨空间位置大小从上到下对左右肋骨进行标记。

在本发明中，所述过程四中，泊松曲面重建肋骨轮廓具体为：

输入经过过程三人工修正后提取的肋骨分割结果的表面点集，记为点集合s，该点集合满足任意元素s∈s，且任意元素s都具有两个已知参数：元素s的三维空间中对应的位置s.p、元素s的内法线方向同时，记vec为三维欧式空间上的向量场，该向量场满足对所有的s∈s，向量场的值在点s.p上的值为而向量场在其余点的值为0：

泊松重建是寻找三维空间上的一个标量函数χ，使得式子最小；其中，符号||·||代表欧氏距离范数；

泊松重建转而求解将作为上述χ的近似解；上式中，△代表拉普拉斯算符，为三维空间上的一个标量函数；

最后，泊松重建通过求取同属一个水平集的点集：获得表面重建结果的顶点集合vp，并对vp构造相应八叉树，即获得对应表面重建结果gp，输出的泊松重建的结果记为：gp＝(vp，ep)；

上式中，vp代表待求水平集上的点集，q代表三维欧氏空间中的点，其中|·|代表集合的元素个数；gp代表泊松重建表面网格结果，其中ep是表面网格中边的集合，ep中任意元素是vp中两点间的连线。

在本发明中，所述过程四中，拉普拉斯网格收缩提取肋骨骨架是指：在提取骨架过程中保持了细节形变，通过估计肋骨上点云邻域建立散乱点之间的连接关系，然后构建拉普拉斯加权矩阵，将所有点作为全局的位置约束，进行点云收缩，对收缩后的点云进行聚类，连接成线得到肋骨骨架；具体为：

对于任意一个表面重建结果g＝(v，e)；其中，v为该表面的顶点集合(泊松重建所提取的等值面中的点)；e是该表面重建结果中边的集合(表面网格中边的集合)，是v中点的连线；假设g中的顶点数，即集合v的元素个数为n，记v＝{v1，...，vn}，vi∈r³；定义并计算该表面重建结果的拉普拉斯矩阵为l，对于该表面重建结果g(v，e)，记αij和βij为边(i，j)∈e的对角，其中(i，j)∈e代表从vi到vj的连线包含在e中，i，j分别为v中点的索引；定义角度权重：

wi，j＝cotαi，j+cotβi，j

通过如下方式计算局部拉普拉斯矩阵l＝{li，j}n×n

上式中，(i，j)中的i和j分别指代当前边的两个顶点的索引；k指代在边(i，k)中不同于当前索引i的那个顶点的索引；li，j是指拉普拉斯矩阵第i行第j列的元素；

拉普拉斯网格收缩提取肋骨骨架利用一个迭代过程实现，记上一步骤中得到的肋骨泊松表面重建结果为gp＝g¹＝(v¹，e)，其中g¹，v¹和e分别代表算法初始时的表面重建结果、与之对应的初始顶点集合以及边的集合，v¹＝vp，e＝ep；设初始权重矩阵和

首先计算对应的拉普拉斯矩阵l¹，在第t次迭代过程中，通过已知的第t步的结果其变量分别为第t步的表面重建结果，权重，以及拉普拉斯矩阵，算法重复以下步骤，并更新结果到第t+1次迭代中，直到当前迭代生成的表面重建结果体积为0，迭代过程具体为：

1)更新第t+1步的表面重建结果g^t+1(v^t+1，e)，其中v^t+1是第t+1次更新后的顶点集合，通过求解最小二乘：[wl^tl^twh^t]^tv^t+1＝[0wh^tv^t]^t；

其中，大写的t代表向量的转置，v^t为第t次迭代生成的顶点集合，wl^t以及wh^t为第t次迭代的权重矩阵；

2)更新第t+1步的权重矩阵为：

其中，和分别是顶点i当前邻域的面积和最初邻域的面积，sl值可调(一般取2.0)；

3)通过表面重建结果g^t+1(v^t+1，e)，计算第t+1步的拉普拉斯矩阵l^t+1；

迭代过程结束后，最后生成的表面重建结果为g^∞(v^∞，ep)，并且获得该最后收敛的表面重建结果的骨架v^∞，v^∞为该表面的顶点集合。

在本发明中，所述过程四中，应用最近点迭代矫正骨架具体是指：对给定两组不同的三维数据点集，其中一个为源数据点集(待配准点集)p＝{p1，..，pn}，另一个为目标点集q＝{q1，...，qm}，其中m≥n，icp试图找出一种刚性变换，将点集p与目标q进行匹配，即icp通过迭代的方法来寻找这样的最佳刚性匹配；

在第k次迭代过程中，icp算法找到与源数据点集p中的每个点对应的目标点集q中的点，并且计算最佳的刚性变换，使得下述距离最小：

其中，分别为第k步迭代的点集p，q中的对应点，t＝(r；t)，r∈so(3)，t∈r³为刚性变换，γ代表所有刚性变换的集合；其中so(3)代表三阶正交群，||·||代表欧氏度量，符号代表函数复合；

具体来说icp主要包括下述步骤：

(1)对每一个pi^k，计算目标点集q中对应点使最小；

(2)计算旋转矩阵r^k与平移向量t^k，极小化距离：

(3)更新第k+1步源点云中点的位置：

(4)计算l²距离误差：

(5)如果d^k+1不小于给定的阈值δ，返回到(1)执行，直到d^k+1＜δ或迭代次数大于预设的最大迭代次数为止；

最终，得到的变换为：其中ti为第i次迭代中求得的刚性变换；

因源数据点集p为肋骨骨架内的骨骼点v^∞，而点集q为与之对应的肋骨二值分割结果(看作点集)，通过上述应用最近点迭代矫正骨架的迭代算法，即实现将肋骨骨架与分割结果对齐。

在本发明中，所述过程四中，cpr方法展开肋骨具体包括下述步骤：

(1)应用得到的肋骨骨架(即通过应用最近点迭代矫正后的骨架)，均匀选取上面的3d采样点；

(2)在每个采样点上计算基于相邻点和前一点的正交观测方向，在每个点附近采样一个平面；

(3)在每个点上重复上述过程，取一个切片，重复第二步；

(4)把采样的切片相互堆叠在一起形成一个直卷；

(5)重复每根肋骨从而把肋骨展开。

与现有技术相比，本发明的有益效果是：

本发明通过把肋骨展开投影到平而上，消除了肋骨弯曲遮挡对诊断的影响，相对与传统临床诊断方法，在平面上观察提高率诊断的直观性，本发明可以作为诊断微小骨折骨裂的辅助方法。

附图说明

图1为肋骨二值图。

图2为修正后去掉脊柱和胸骨。

图3为肋骨xz方向展开图。

图4为肋骨yz方向展开图。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述：

一种新的肋骨可视化肋骨骨折辅助诊断方法，包括下述过程：

一、获取图片；

二、提取肋骨；

三、人工修正，去掉脊柱、胸骨；

四、展开肋骨。

所述过程一具体包括：通过和医院联系得到人体胸部ct扫描图片。

所述过程二具体包括：

对腹部ct图像进行处理，是基于骨骼的密度比较大，在ct扫描中ct值比较大，表现为骨骼与其他组织器官相比亮度高，因而有明显区别，通过选取适当的窗宽窗位使用多尺度的增强滤波算子提取出肋骨，在保留骨骼区域的同时避免了其他组织和器官(如血管、钙化点、肾脏)的干扰。具体步骤是：

步骤a：把图像转化为256色灰度图像

步骤b：然后使用一个增强滤波算子对目标区域进项强化，从图像中逐点提取平滑的黑塞矩阵

其中，所述i是指每点的二阶导数矩阵；所述iij是指二阶导数，i、j分别取x、y或者z；

步骤c：黑塞矩阵i通过高斯核的二阶偏导数与原图像进行卷积求得：

其中，gσ(v)＝exp(-v^tv/σ²)，v∈r³为高斯核；所述△为拉普拉斯算符；所述为步骤a中的256色灰度图像；所述σ是指当前高斯核的标准差；所述t是指向量的转置；所述exp是指自然指数函数；所述*指的是函数间的卷积，所述r³为三维欧氏空间，所述v为欧氏空间中的一个三维向量；

步骤d：求出黑塞矩阵i的特征值λ1，λ2，λ3，且|λ1|＜|λ2|＜|λ3|；

计算每点的增强值：

其中，所述v0(i)用于表示像素值类似血管(vesselness)的程度；所述α、β、c都分别大于0(均为可调整经验参数，可根据实际灰度图像自行调节)；

步骤e：通过调整像素值范围提取出胸部骨骼，包含胸骨、脊柱；

步骤f：按照肋骨空间位置大小从上到下对左右肋骨进行标记。

所述过程三具体包括：

创建三维窗口，用户通过鼠标右键点击划圈标示出不属于肋骨的部分，计算机根据用户的标示将圈内部分(脊柱、胸骨)去除。

所述过程四具体包括：

泊松曲面重建肋骨轮廓是通过估计模型的指示函数和提取等值面，对表面重建一个无缝的三角逼近；向量场v为空间m(在模型表面具有高分辨率，远离表面地方具有低分辨率)中的函数线性总和，构造并解泊松方程，提取求取的指示函数的等值面(它是一种全局的方法，把有法向点集的表面重建转化为一个空间的泊松问题，从而把表面重建问题表示为泊松方程的解)。

输入经过过程三人工修正后提取的肋骨分割结果的表面点集，记为点集合s，该点集合满足任意元素s∈s都具有两个已知参数：该元素s的三维空间中对应的位置s.p，以及该元素的内法线方向同时，记vec为三维欧式空间上的向量场，该向量场满足对所有的s∈s，向量场的值在点s.p上的值为而向量场在其余点的值为0。泊松重建希望寻找三维空间上的一个标量函数χ，使得式了最小。其中符号||·||代表欧氏距离范数。泊松重建转而求解然后将作为上述χ的近似解。上式中△代表拉普拉斯算符，为三维空间上的一个标量函数。最后，泊松重建通过求取同属一个水平集的点集：获得表面重建结果的顶点集合vp，并对vp构造相应八叉树并获得对应表面重建结果gp。上式中vp代表待求水平集上的点集，q代表三维欧氏空间中的点，其中|·|代表集合的元素个数。gp代表泊松重建表面网格结果。记gp＝(vp，ep)，其中ep是表面网格中边的集合，ep中任意元素vp中两点间的连线，下同。

在提取骨架过程中保持了细节形变，通过估计肋骨上点云邻域建立散乱点之间的连接关系，然后构建拉普拉斯加权矩阵，将所有点作为全局的位置约束，进行点云收缩，对收缩后的点云进行聚类，连接成线得到肋骨骨架。

对任意一个表面重建结果g＝(v，e)，其中v为该表面的顶点集合。而e是该表面重建结果中边的集合，为v中点的连线。假设g中的顶点数，即集合v的元素个数为n。在此情况下，记v＝{v1，...，vn}，vi∈r³。在这样的定义下，我们定义并计算该表面重建结果的拉普拉斯矩阵为l。对于该表面g(v，e)，记αij和βij为边(i，j)∈e的对角，其中(i，j)∈e代表从vi到vj的连线包含在e中，i，j分别为v中点的索引。定义角度权重：

wi，j＝cotαi，j+cotβi，j

我们通过如下方式计算局部拉普拉斯矩阵l＝{li，j}n×n

上式中，(i，j)中的i和j分别指代当前边的两个顶点的索引。k指代在边(i，k)中不同于当前索引i的那个顶点的索引。li，j是指拉普拉斯矩阵第i行第j列的元素。

拉普拉斯网格收缩提取肋骨骨架利用一个迭代过程实现，我们记上一步骤中得到的肋骨泊松表面重建结果为gp＝g¹＝(v¹，e)，其中g¹，v¹和e分别代表算法初始时的表面重建结果、与之对应的初始顶点集合以及边的集合，v¹＝vp，e＝ep。设初始权重矩阵和算法首先计算对应的拉普拉斯矩阵l¹。在第t次迭代过程中，通过已知的第t步的结果其变量分别为第t步的表面重建结果，权重，以及拉普拉斯矩阵，算法重复以下步骤，并更新结果到第t+1次迭代中，直到当前迭代生成的表面重建结果体积为0：

1)更新第t+1步的表面重建结果g^t+1(v^t+1，e)，其中v^t+1是第t+1次更新后的顶点集合，通过求解最小二乘：[wl^tl^twh^t]^tv^t+1＝[0wh^tv^t]^t。其中大写的t代表向量的转置，v^t为第t次迭代生成的顶点集合，wl^t以及wh^t为第t次迭代的权重矩阵。

2)更新第t+1步的权重矩阵为：其中，和分别是顶点i当前邻域的面积和最初邻域的面积，sl值可调，一般取2.0。

3)通过表面重建结果g^t+1(v^t+1，e)计算第t+1步的图像拉普拉斯矩阵l^t+1。

记最后生成的表面重建结果为g^∞(v^∞，ep)，则该表面的骨架即为v^∞，其中v^∞为该表面的顶点集合。

在本专利中，算法将第三步生成的表面重建结果gp＝(vp，ep)作为初始值即g¹＝(v¹，e)带入上述算法，并且获得该表面重建结果的骨架。这个骨架是一系列点的集合，该集合即为最后收敛的表面重建结果g^∞(v^∞，ep)中的v^∞。

应用最近点迭代矫正骨架具体是指：对给定两组不同的三维数据点集，其中一个为源数据点集(待配准点集)p＝{p1，...，pn}，另一个为目标点集q＝{q1，...，qm}，其中m≥n，icp试图找出一种刚性变换，将点集p与目标q进行匹配，即icp通过迭代的方法来寻找这样的最佳刚性匹配)。

在第k次迭代过程中，icp算法找到与源点集p中的每个点对应的目标点集q中的点，并且计算最佳的刚性变换，使得下述距离最小：

其中，分别为第k步迭代的点集p，q中的对应点，t＝(r；t)，r∈so(3)，t∈r³为刚性变换，γ代表所有刚性变换的集合；其中so(3)代表三阶正交群，||·||代表欧氏度量，符号代表函数复合。(所有参数及符号定义下同)

具体来说icp主要包括下述步骤：

(1)对每一个pi^k，计算参考点集q中对应点使最小；

(2)计算旋转矩阵r^k与平移向量t^k，极小化距离：

(3)更新第k+1步源点云中点的位置：

(4)计算l²距离误差：

(5)如果d^k+1不小于给定的阈值δ，返回到(1)执行，直到d^k+1＜δ或迭代次数大于预设的最大迭代次数为止；

最终，得到的变换为：其中ti为第i次迭代中求得的刚性变换。

因点集p为肋骨骨架内的骨骼点v^∞，而点集q为与之对应的肋骨二值分割结果(看作点集)，通过该算法，即实现将肋骨骨架与分割结果对齐。

利用cpr方法展开肋骨：

(1)应用得到的肋骨骨架(即通过应用最近点迭代矫正后的骨架)，均匀选取上面的3d采样点；

(2)在每个采样点上计算基于相邻点和前一点的正交观测方向，在每个点附近采样一个平面；

(3)在每个点上重复上述过程，取一个切片，重复第二步；

(4)把采样的切片相互堆叠在一起形成一个直卷；

(5)重复每根肋骨从而把肋骨展开。图3和图4是过程四得到的肋骨xz和yz方向的展开图，可以从图上观察到骨断裂处。

最后，需要注意的是，以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有很多变形。本领域的普通技术人员能从本发明公开的内容中直接导出或联想到的所有变形，均应认为是本发明的保护范围。