CN110782470B - 一种基于形状信息的腕骨区域分割方法 - Google Patents

Abstract

一种基于形状信息的腕骨区域分割方法，包括：采用星型中值滤波，灰度图像二值化方法，依次对X光图像进行预处理；根据腕骨区域分布特点，使用双向迭代扫描算法，确定腕骨区域上边界的纵坐标位置；根据手腕关节形状信息，计算手腕两侧边缘梯度变化，确定下边界左右候选点纵坐标位置；综合上边界、左右候选点和腕骨区域信息，对坐标位置进行优化矫正，并根据得到的左上角、右下角坐标进行腕骨区域裁剪。本发明准确分割出手腕图像中的腕骨区域；在确定上边界时，采用双向迭代扫描算法，极大提升了扫描速率；在确定下边界时，同时考虑了左右两个边界点的情况，并进行优化选择，减弱了手掌关节位置摆放不同对定位准确率的影响。

Description

一种基于形状信息的腕骨区域分割方法

技术领域

本发明涉及一种基于形状信息的腕骨区域分割方法。

背景技术

骨龄评估的方法有很多种，最早如20世纪初期的骨龄计数法，主要观察骨化中心在一定年龄出现的数目，主要分布在腕骨区域。后来，当X线被发现并被投入骨骼发育研究，许多专家相继提出并修订了较精确的骨骼成熟图谱。于是有了图谱法，评价时将待测X线片与图谱逐个对照，取最相近者为其骨龄。近期随着人工智能技术在医学领域上的广泛应用，以及深度卷积神经网络在图像处理技术上的重大突破，越来越多的学者从机器学习的角度上研究骨龄。

同时，大量的实验和研究结果表明，腕骨作为人体手腕部一个重要组成部位，在骨龄评测中是不可忽视的一个潜在信息来源。发育早期，它们在图像中呈现为密集的针尖形状，在发育过程中，它们的体型不断增大，最终达到成年时的体型和特征形状。这一特点在儿童1～8岁时期尤其明显。因此，不管是传统的骨龄评估方法，还是最先进的深度学习算法，对腕骨进行一个针对性的提取分割，以及进一步研究就显得极具重要意义。

腕骨，主要集中在人手腕X光图像中指骨及掌骨相对靠下位置。而鉴于其特殊的形状规则，往往可以从阈值、区域和边缘等方向实现定位分割。目前常见的图像处理算法，如纵向方波扫描，通过固定像素距离从上到下扫描二值化图像，根据相邻最小或最大连续横坐标序列宽度位置，来定位腕骨区域左上角或右下角坐标。也有利用神经网络对大量人工标注的腕骨区域进行学习，从而自动预测新输入图像的腕骨区域。

然而，由于常见的人手腕X光图像在拍摄体位上要求不严格，不同人手腕形状各异，手掌关节摆放不正等因素，给统一的图像处理方法带来了许多困难。例如根据腕关节最短横向宽度的特点定位腕骨尾部区域，会因部分群体不满足该特点而引入巨大误差。根据神经网络预测的腕骨区域，极大程度上取决于训练样本，对于新采集的色度、亮度等迥异的图像，预测精度会大大下降。同时，也因其自动定位的特点，导致了最终腕骨区域大小的不可控。

发明内容

本发明要克服现有技术的上述缺点，提供一种基于形状信息的腕骨区域分割方法。

本发明解决技术问题采用如下技术方案：

一种基于形状信息的腕骨区域分割方法，包括以下步骤：

步骤1：采用星型中值滤波，灰度图像二值化方法，依次对X光图像进行预处理；

步骤2：根据腕骨区域分布特点，使用双向迭代扫描算法，确定腕骨区域上边界的纵坐标位置；

步骤3：根据手腕关节形状信息，计算手腕两侧边缘梯度变化，确定下边界左右候选点纵坐标位置；

步骤4：综合上边界、左右候选点和腕骨区域信息，对坐标位置进行优化矫正，并根据得到的左上角、右下角坐标进行腕骨区域裁剪。

作为本发明的进一步优化，所述的步骤1具体包括：

1)将原始图像转换为灰度图，公式为：

grey＝(299*R+587*G+114*B+500)/1000 (1-1)

其中，R、G、B代表三通道原始图像具体位置的三原色像素值，grey表示由此转换成的灰度值；

2)采用星型中值滤波器，对灰度图像进行降噪处理，公式为：

Gray(i,j)＝median[x(i,j),x(i,j+1),x(i,j-1),x(i-1,j),x(i+1,j)] (1-2)

其中，x代表经过上一步处理后的灰度图像，x(i,j)代表图像中坐标(i,j)位置处的灰度值，median表示对列表内的数值进行中值运算操作，Gray(i,j)表示图像坐标(i,j)位置经过滤波后的灰度值；

3)对图像进行二值化处理，将手腕和背景区分出来，手腕部分使用白色像素点填充，背景部分使用黑色像素点填充，得到二值化后的图像，二值化的公式为：

其中，0代表黑色像素填充，255代表白色像素点填充，threshold表示像素区分的阈值，其取值方法为：

首先计算灰度图像Gray的平均灰度值ths，其中mean表示对列表内的数值进行均值运算操作：

ths＝mean[Gray(i,j)] (1-4)

分别挑选出灰度值大于等于和小于ths值的所有像素点c1，c2；

c1＝{Gray(i,j)|Gray(i,j)＞＝ths} (1-5)

c2＝{Gray(i,j)|Gray(i,j)＜ths} (1-6)

然后计算阈值threshold，公式为：

判断此时fabs(threshold-ths)值是否小于0.001，若是则结束运算，否则运算ths＝threshold，重新挑选c1、c2，并根据(1-7)计算新的threshold值，直到满足以上条件；

4)对二值化图像进行8连通域的标记，使得同一连通区域内的所有像素相同，不同连通区域像素不同，并统计所有连通区域内白色像素点S_n，n代表标记为n的连通域。寻找白色像素点最多的一个连通区域m，即为手掌区域，满足S_m≥S_n，

N为连通区域的数量。其余的连通区域当作背景区域，使用0像素值填充：

其中，X表示经过预处理后的图像。

作为本发明的更进一步优化，所述的步骤2具体包括：

1)求出连通区域S_m的最小外接矩阵，高和宽分别为BH，BW；预处理后的图像高和宽分别为XH，XW；

2)定义一个函数F(R_t)，接收图像中某一行t的所有像素值集合R_t作为参数，R_t＝{X(t,j)|j＝1,2,…,XW}，t∈{1,2,…,XH}，用于判断第t行是否满足以下条件：

p_i＝len(G(R_i)) (2-2)

其中，D值表示判定取样间隙，初始值为5，函数G(R_i)表示找出R_i集合中所有像素大于0的连续子序列段，p_i表示G(R_i)中序列段的数量；

3)从下往上，每隔50像素点对图像X进行扫描，起始位置为XH-BH/3，结束位置为XH-BH，当某行t，满足F(R_t)为TRUE时，XH-BH/3≤t≤XH-BH，跳转到4)，否则继续扫描；

4)从上往下，每隔5像素点对图像X进行扫描，起始位置为t，结束位置为t+50，当某行u，满足F(R_u)为TRUE时，t≤u≤t+50，结束扫描，否则继续扫描；

5)确定腕骨区域上边界的纵坐标位置top，公式为：

作为本发明的更进一步优化，所述的步骤3具体包括：

1)求出手腕图像下半部分左边缘集LC，右边缘集RC，公式为：

2)遍历左边缘集LC，并以step＝50为步长计算梯度grad_i，当某个点i，0<i<len(LC)-step，其左邻域内梯度呈现从一个极值持续衰减的趋势，且满足以下条件时停止遍历：

确定下边界左候选点loc，即为此时i的值；

3)遍历右边缘集RC，并以step＝50为步长计算梯度grad_i，当某个点i，0<i<len(RC)-step，其左邻域内梯度呈现从一个极值持续衰减的趋势，且满足以下条件时停止遍历：

确定下边界右候选点roc，即为此时i的值；

作为本发明的更进一步优化，所述的步骤4具体包括：

1)确定腕骨区域合理的最小长度MINH＝300，最大长度MAXH＝700；

2)根据候选点loc，roc分别计算候选区域长度，判断是否大于MINH且小于MAXH，同时剔除不规范的候选点；

3)根据候选点剩余情况，分类确定最终的坐标点bot；

4)根据上边界纵坐标top值，下边界纵坐标bot值计算腕骨区域左上角坐标tl、右下角坐标br，公式为：

5)根据两点坐标tl、br，绘制矩形并裁剪原始图像。

目前，许多图像切割方法，对图像质量要求较严格，而且不同人手腕形状各异，手掌关节摆放不正等因素，会降低切割准确率。本发明在保证腕骨区域分割准确的同时，还保证其具有较强的普适性。采用双向迭代扫描算法确定上边界，能够提升扫描速率，在进行大批量图像裁剪时，极大地减少了时间成本；在确定下边界时，考虑到手部关节位置变化会引起手掌两侧边缘曲率的变化，综合计算出左右两个边界点的位置，并剔除明显不符合情况的位置，减弱了手掌关节位置摆放不同对定位准确率的影响；最后利用得到的位置点信息，进行优化矫正，进一步保证切割区域的精确性。

本发明具有如下有益效果：

(1)准确分割出手腕图像中的腕骨区域；

(2)在确定上边界时，采用双向迭代扫描算法，极大提升了扫描速率；在确定下边界时，同时考虑了左右两个边界点的情况，并进行优化选择，减弱了手掌关节位置摆放不同对定位准确率的影响。

附图说明

图1是本发明步骤3中边缘集梯度变化示意图。

图2是本发明步骤4中判定流程图。

图3是本发明分割出的腕骨结果示意图。

图4是本发明分割方法总体流程图。

具体实施方式

下面结合附图进一步说明本发明的技术方案。

一种基于形状信息的腕骨区域分割方法，其特征在于，包括以下步骤：

步骤1：采用星型中值滤波，灰度图像二值化方法，依次对X光图像进行预处理；具体包括：

11)将原始图像转换为灰度图，公式为：

grey＝(299*R+587*G+114*B+500)/1000 (1-1)

其中，R、G、B代表三通道原始图像具体位置的三原色像素值，grey表示由此转换成的灰度值；

12)采用星型中值滤波器，对灰度图像进行降噪处理，公式为：

Gray(i,j)＝median[x(i,j),x(i,j+1),x(i,j-1),x(i-1,j),x(i+1,j)] (1-2)

其中，x代表经过上一步处理后的灰度图像，x(i,j)代表图像中坐标(i,j)位置处的灰度值，median表示对列表内的数值进行中值运算操作，Gray(i,j)表示图像坐标(i,j)位置经过滤波后的灰度值；

13)对图像进行二值化处理，将手腕和背景区分出来，手腕部分使用白色像素点填充，背景部分使用黑色像素点填充，得到二值化后的图像，二值化的公式为：

其中，0代表黑色像素填充，255代表白色像素点填充，threshold表示像素区分的阈值，其取值方法为：

首先计算灰度图像Gray的平均灰度值ths，其中mean表示对列表内的数值进行均值运算操作：

ths＝mean[Gray(i,j)] (1-4)

分别挑选出灰度值大于等于和小于ths值的所有像素点c1，c2；

c1＝{Gray(i,j)|Gray(i,j)＞＝ths} (1-5)

c2＝{Gray(i,j)|Gray(i,j)＜ths} (1-6)

然后计算阈值threshold，公式为：

判断此时fabs(threshold-ths)值是否小于0.001，若是则结束运算，否则运算ths＝threshold，重新挑选c1、c2，并根据(1-7)计算新的threshold值，直到满足以上条件；

14)对二值化图像进行8连通域的标记，使得同一连通区域内的所有像素相同，不同连通区域像素不同，并统计所有连通区域内白色像素点S_n，n代表标记为n的连通域。寻找白色像素点最多的一个连通区域m，即为手掌区域，满足S_m≥S_n，

N为连通区域的数量。其余的连通区域当作背景区域，使用0像素值填充：

其中，X表示经过预处理后的图像。

2.步骤2：根据腕骨区域分布特点，使用双向迭代扫描算法，确定腕骨区域上边界的纵坐标位置；具体包括：

21)求出连通区域S_m的最小外接矩阵，高和宽分别为BH，BW；预处理后的图像高和宽分别为XH，XW；

22)定义一个函数F(R_t)，接收图像中某一行t的所有像素值集合R_t作为参数，R_t＝{X(t,j)|j＝1,2,…,XW}，t∈{1,2,…,XH}，用于判断第t行是否满足以下条件：

p_i＝len(G(R_i)) (2-2)

其中，D值表示判定取样间隙，初始值为5，函数G(R_i)表示找出R_i集合中所有像素大于0的连续子序列段，p_i表示G(R_i)中序列段的数量；

23)从下往上，每隔50像素点对图像X进行扫描，起始位置为XH-BH/3，结束位置为XH-BH，当某行t，满足F(R_t)为TRUE时，XH-BH/3≤t≤XH-BH，跳转到24)，否则继续扫描；

24)从上往下，每隔5像素点对图像X进行扫描，起始位置为t，结束位置为t+50，当某行u，满足F(R_u)为TRUE时，t≤u≤t+50，结束扫描，否则继续扫描；

25)确定腕骨区域上边界的纵坐标位置top，公式为：

3.步骤3：根据手腕关节形状信息，计算手腕两侧边缘梯度变化，确定下边界左右候选点纵坐标位置；具体包括：

31)求出手腕图像下半部分左边缘集LC，右边缘集RC，公式为：

32)遍历左边缘集LC，并以step＝50为步长计算梯度grad_i，当某个点i，0<i<len(LC)-step，其左邻域内梯度呈现从一个极值持续衰减的趋势，且满足以下条件时停止遍历：

确定下边界左候选点loc，即为此时i的值；

33)遍历右边缘集RC，并以step＝50为步长计算梯度grad_i，当某个点i，0<i<len(RC)-step，其左邻域内梯度呈现从一个极值持续衰减的趋势，且满足以下条件时停止遍历：

确定下边界右候选点roc，即为此时i的值。

4.步骤4：综合上边界、左右候选点和腕骨区域信息，对坐标位置进行优化矫正，并根据得到的左上角、右下角坐标进行腕骨区域裁剪。具体包括：

41)确定腕骨区域合理的最小长度MINH＝300，最大长度MAXH＝700；

42)根据候选点loc，roc分别计算候选区域长度，判断是否大于MINH且小于MAXH，同时剔除不规范的候选点；

43)根据候选点剩余情况，分类确定最终的坐标点bot；

44)根据上边界纵坐标top值，下边界纵坐标bot值计算腕骨区域左上角坐标tl、右下角坐标br，公式为：

45)根据两点坐标tl、br，绘制矩形并裁剪原始图像。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (4)

1.一种基于形状信息的腕骨区域分割方法，其特征在于，包括以下步骤：

步骤1：采用星型中值滤波，灰度图像二值化方法，依次对X光图像进行预处理；对二值化图像进行8连通域的标记，使得同一连通区域内的所有像素相同，不同连通区域像素不同，并统计所有连通区域内白色像素点S_n，n代表标记为n的连通域；寻找白色像素点最多的一个连通区域m，即为手掌区域，满足S_m≥S_n，

N为连通区域的数量；其余的连通区域当作背景区域，使用0像素值填充：

其中，X表示经过预处理后的图像；

步骤2：根据腕骨区域分布特点，使用双向迭代扫描算法，确定腕骨区域上边界的纵坐标位置；具体包括：

21)求出连通区域S_m的最小外接矩阵，高和宽分别为BH，BW；预处理后的图像高和宽分别为XH，XW；

22)定义一个函数F(R_t)，接收图像中某一行t的所有像素值集合R_t作为参数，R_t＝{X(t,j)|j＝1,2,…,XW}，t∈{1,2,…,XH}，用于判断第t行是否满足以下条件：

D＝5

p_i＝len(G(R_i)) (2-2)

X(i,a-1)＝0,X(i,b+1)＝0}

其中，D值表示判定取样间隙，初始值为5，函数G(R_i)表示找出R_i集合中所有像素大于0的连续子序列段，p_i表示G(R_i)中序列段的数量；

23)从下往上，每隔50像素点对图像X进行扫描，起始位置为XH-BH/3，结束位置为XH-BH，当某行t，满足F(R_t)为TRUE时，XH-BH/3≤t≤XH-BH，跳转到24)，否则继续扫描；

24)从上往下，每隔5像素点对图像X进行扫描，起始位置为t，结束位置为t+50，当某行u，满足F(R_u)为TRUE时，t≤u≤t+50，结束扫描，否则继续扫描；

25)确定腕骨区域上边界的纵坐标位置top，公式为：

步骤3：根据手腕关节形状信息，计算手腕两侧边缘梯度变化，确定下边界左右候选点纵坐标位置；

步骤4：综合上边界、左右候选点和腕骨区域信息，对坐标位置进行优化矫正，并根据得到的左上角、右下角坐标进行腕骨区域裁剪。

2.根据权利要求1所述的一种基于形状信息的腕骨区域分割方法，其特征在于，所述的步骤1具体包括：

11)将原始图像转换为灰度图，公式为：

grey＝(299*R+587*G+114*B+500)/1000 (1-1)

其中，R、G、B代表三通道原始图像具体位置的三原色像素值，grey表示由此转换成的灰度值；

12)采用星型中值滤波器，对灰度图像进行降噪处理，公式为：

Gray(i,j)＝median[x(i,j),x(i,j+1),x(i,j-1),x(i-1,j),x(i+1,j)] (1-2)

其中，x代表经过上一步处理后的灰度图像，x(i,j)代表图像中坐标(i,j)位置处的灰度值，median表示对列表内的数值进行中值运算操作，Gray(i,j)表示图像坐标(i,j)位置经过滤波后的灰度值；

13)对图像进行二值化处理，将手腕和背景区分出来，手腕部分使用白色像素点填充，背景部分使用黑色像素点填充，得到二值化后的图像，二值化的公式为：

其中，0代表黑色像素填充，255代表白色像素点填充，threshold表示像素区分的阈值，其取值方法为：

首先计算灰度图像Gray的平均灰度值ths，其中mean表示对列表内的数值进行均值运算操作：

ths＝mean[Gray(i,j)] (1-4)

分别挑选出灰度值大于等于和小于ths值的所有像素点c1，c2；

c1＝{Gray(i,j)|Gray(i,j)＞＝ths} (1-5)

c2＝{Gray(i,j)|Gray(i,j)＜ths} (1-6)

然后计算阈值threshold，公式为：

判断此时fabs(threshold-ths)值是否小于0.001，若是则结束运算，否则运算ths＝threshold，重新挑选c1、c2，并根据(1-7)计算新的threshold值，直到满足以上条件；

14)对二值化图像进行8连通域的标记，使得同一连通区域内的所有像素相同，不同连通区域像素不同，并统计所有连通区域内白色像素点S_n，n代表标记为n的连通域；寻找白色像素点最多的一个连通区域m，即为手掌区域，满足S_m≥S_n，

N为连通区域的数量；其余的连通区域当作背景区域，使用0像素值填充：

其中，X表示经过预处理后的图像。

3.根据权利要求2所述的一种基于形状信息的腕骨区域分割方法，其特征在于，所述的步骤3具体包括：

31)求出手腕图像下半部分左边缘集LC，右边缘集RC，公式为：

32)遍历左边缘集LC，并以step＝50为步长计算梯度grad_P，当某个点P，0<P<len(LC)-step，其左邻域内梯度呈现从一个极值持续衰减的趋势，且满足以下条件时停止遍历：

确定下边界左候选点loc，即为此时P的值；

33)遍历右边缘集RC，并以step＝50为步长计算梯度grad_Q，当某个点Q，0<Q<len(RC)-step，其左邻域内梯度呈现从一个极值持续衰减的趋势，且满足以下条件时停止遍历：

确定下边界右候选点roc，即为此时Q的值。

4.根据权利要求3所述的一种基于形状信息的腕骨区域分割方法，其特征在于，所述的步骤4具体包括：

41)确定腕骨区域合理的最小长度MINH＝300，最大长度MAXH＝700；

42)根据候选点loc，roc分别计算候选区域长度，判断是否大于MINH且小于MAXH，同时剔除不规范的候选点；

43)根据候选点剩余情况，分类确定最终的下边界纵坐标bot；

44)根据上边界纵坐标top值，下边界纵坐标bot值计算腕骨区域左上角坐标tl、右下角坐标br，公式为：

tl＝(LC_bot,top)

45)根据两点坐标tl、br，绘制矩形并裁剪原始图像。

Images