专利名称：基于超声波测量的残肢骨骼与皮肤特征的提取方法
技术领域：
本发明涉及一种残肢骨骼与皮肤特征的提取方法，尤其涉及一种基于超声波测量的残肢骨骼与皮肤特征的提取方法。
背景技术：
在医疗康复技术上，对于残缺肢体，主要是通过安装人工假体，部分恢复原有生理功能。假肢主要由接受腔、胫管、假体组成，接受腔是人工肢体与患者残体的结合部分，其结构合理性直接影响患者使用假体的舒适度。制造接受腔首先要对残肢取型，根据残肢数据制造接受腔。传统的残肢取型采用石膏复型法，这种方法存在以下缺点①精度低，石膏复型劳动强度大；②石膏复型只对皮肤外轮廓取型，无法对内部骨骼取型，骨骼和软组织的生物力学特性被忽略。
现在的残肢轮廓测量技术还有接触式(包括身体测量法、浸水法和圆周测量、口型圈测量法、阴型扫描法、阴型数字测量法)和非接触式(包括电子扫描式、激光扫描式、测面影像测量等)等数字化采集技术，尽管这些数字化采集技术可以得到数字化的残肢外形，但也无法对内部骨骼进行测量。
利用现代医学影像技术如X射线投影技术、断层扫描CT技术以及磁共振MRI等可以对内部组织进行测量。常用的扫描技术是CT和MRI，它在不破坏物体的前提下，根据投影，重建物体特定断面上的无重叠二维图像，但这种方法价格昂贵，还会对人体产生电离辐射。
发明内容
本发明的目的是提供一种利用超声波来测量残肢骨骼及皮肤形成二维图像进行特征提取的方法，解决了石膏复型法精度低、CT扫描等方法对人体有辐射的不足之处，可以同时对皮肤的残肢外形和内部骨骼进行测量。
本发明所采用的技术方案是采用超声波测量装置、B超诊断仪、超声影像工作站和图形处理系统组成超声波测量系统，超声波测量装置上的探头与B超诊断仪连接，B超诊断仪通过数据线与超声影像工作站和图形处理系统连接，其特征在于该方法包括以下步骤a.图像采集探头围绕残体的旋转中心以每0°＜α°≤90°的间隔做圆周运动，同时由超声影像工作站对超声测量图进行采集、存储，得到k＝360°/α帧B超图像，然后探头垂直下移，残体的旋转中心也随之下移，探头又围绕新的旋转中心以每0°＜α°≤90°的间隔做旋转测量，获得k＝360°/α帧B超图像，以此类推，即得到被测残体各水平面的多帧二维截面断层图像；b.图像二维重建第一步先确定残体的旋转中心，通过测量标准圆柱件确定其圆心位置，如果圆心与探头聚焦中心重合，则该圆心即是旋转中心；第二步将k帧图像按各自测量时的旋转角α绕旋转中心旋转到新坐标系，根据坐标旋转原理，其变换表达式为 上式中，(x，y)是原图像坐标系的横纵坐标点，(u，v)是旋转后的图像坐标系的横纵坐标点，再把上式写成矩阵形式为 其中
这里，Ra是旋转变换矩阵，(x，y)是单帧图像的坐标系，(u，v)是复合图像的坐标系；第三步对旋转到新位置的k帧图像进行叠加，设单帧测量图像为Ak，复合图像为Bk+1，Bk+1=Σl=1k+1Alc=Σl=1k+1aij(l)c]]>其中aij(l)表示m×n矩阵Al，m＝768，n＝576，1≤i≤m，1≤j≤n，1≤l≤k+1；第四步对像素值取整，表示为Ck+1＝{Bk+1}，其中{x}表示超过x的最小整数；第五步去噪，阈值设为128， Ck+1表示对复合后的图像像素取整和去噪后的图像像素矩阵，cij表示经过复合，去噪后的的图像像素在第i行，第j列上的灰度值；c.对图像进行处理，提取骨骼和皮肤的边缘特征用函数f(x，y)表示图像中任意一点的灰度，则该点的灰度变化梯度为▿f(x,y)=∂f∂x∂f∂yT=GxGyT]]>是梯度算子，T表示转置，Gx和Gy表示在x和y方向上的灰度值的变化程度，梯度是一个矢量，其幅值为|▿f|=[Gx2+Gy2]1/2]]>表示某点图像的边缘强度，满足下面三个条件既可找到图像的边缘点1.以该点为中心的3×3邻域中的边缘强度极大值小于某阀值，2.该点的边缘强度大于沿该点梯度方向的两个相邻像素点的边缘强度，
3.这一点的相位值φ(x,y)=arctan(GyGx)]]>值小于45°，其中，图像矩阵每个点f(x，y)，离散化成f(i，j)，记为fi，j，梯度Gx，Gy离散化为Gx＝[fi-1，j-1+2fi，j-1+fi+1，j-1]-[fi-1，j+1+2fi，j+1+fi-1，j+1]Gy＝[fi+1，j-1+2fi+1，j+fi+1，j+1]-[fi-1，j-1+2fi-1，j+fi-1，j+1]算子的模板为Gx=-101-202-101,Gy=121000-1-2-1]]>经过以上步骤就完成对图像的骨骼和皮肤的特征提取。
本发明的有益效果是本发明利用超声波成像原理，对残肢的骨骼及皮肤进行超声波测量，可以同时对皮肤外轮廓和内部骨骼进行测量，其中的图像采集，既可获得不同层面的图像，而且在同一层面上可得到多角度的图像，提高了测量的精确度，与同类技术相比，极大地降低了成本、减少了辐射，图像的重建及骨骼和皮肤的特征提取为残肢的三维CAD的建立奠定了基础。
图1是超声测量系统的结构示意图；图2是图像复合原理图；图3是内轮廓精度检验，其中图3a是内圆柱复合图像，图3b是轮廓精度检验；图4是外轮廓精度检验，其中图4a是外轮廓复合图像，图4b是外轮廓精度检验；图5是重建的猪腿中间部位的截面图像，其中5a是一截面的图像，5b是一截面的图像，5c是一截面的图像；图6是图像特征提取得到的骨骼与皮肤轮廓。
图中，1.超声波测量装置，2.探头，3.储液桶，4.B超诊断仪，5.超声影像工作站，6.图形处理系统。
具体实施例方式
下面结合附图对本发明作详细说明。
参见图1，采用超声波测量装置1、B超诊断仪4、超声影像工作站5和图形处理系统6组成超声波测量系统，超声波测量装置1上的探头2与B超诊断仪4连接，B超诊断仪4通过数据线与超声影像工作站5和图形处理系统6连接，其特征在于该方法包括以下步骤a.图像采集探头2围绕残体的旋转中心以每0°＜α°≤90°的间隔做圆周运动，同时由超声影像工作站5对超声测量图进行采集、存储，得到k＝360°/α帧B超图像，然后探头2垂直下移，残体的旋转中心也随之下移，探头2又围绕新的旋转中心以每0°＜α°≤90°的间隔做旋转测量，获得k＝360°/α帧B超图像，以此类推，即得到被测残体各水平面的多帧二维截面断层图像，参见图2；b.图像二维重建第一步先确定残体的旋转中心，通过测量标准圆柱件确定其圆心位置，如果圆心与探头聚焦中心重合，则该圆心即是旋转中心；第二步将k帧图像按各自测量时的旋转角α绕旋转中心旋转到新坐标系，根据坐标旋转原理，其变换表达式为
上式中，(x，y)是原图像坐标系的横纵坐标点，(u，v)是旋转后的图像坐标系的横纵坐标点，再把上式写成矩阵形式为 其中 这里，Ra是旋转变换矩阵，(x，y)是单帧图像的坐标系，(u，v)是复合图像的坐标系；第三步对旋转到新位置的k帧图像进行叠加，设单帧测量图像为Ak，复合图像为Bk+1，Bk+1=Σl=1k+1Alc=Σl=1k+1aij(l)c]]>其中aij(l)表示m×n矩阵Al，m＝768，n＝576，1≤i≤m，1≤j≤n，1≤l≤k+1；第四步对像素值取整，表示为Ck+1＝{Bk+1}，其中{x}表示超过x的最小整数；第五步去噪，阈值设为128， Ck+1表示对复合后的图像像素取整和去噪后的图像像素矩阵，cij表示经过复合，去噪后的的图像像素在第i行，第j列上的灰度值；c.对图像进行处理，提取骨骼和皮肤的边缘特征用函数f(x，y)表示图像中任意一点的灰度，则该点的灰度变化梯度为▿f(x,y)=∂f∂x∂f∂yT=GxGyT]]>是梯度算子，T表示转置，Gx和Gy表示在x和y方向上的灰度值的变化程度，梯度是一个矢量，其幅值为|▿f|=[Gx2+Gy2]1/2]]>表示某点图像的边缘强度，满足下面三个条件既可找到图像的边缘点1.以该点为中心的3×3邻域中的边缘强度极大值小于某阀值，
2.该点的边缘强度大于沿该点梯度方向的两个相邻像素点的边缘强度，3.这一点的相位值φ(x,y)=arctan(GyGx)]]>值小于45°，其中，图像矩阵每个点f(x，y)，离散化成f(i，j)，记为fi，j，梯度Gx，Gy离散化为Gx＝[fi-1，j-1+2fi，j-1+fi+1，j-1]-[fi-1，j+1+2fi，j+1+fi-1，j+1]Gy＝[fi+1，j-1+2fi+1，j+fi+1，j+1]-[fi-1，j-1+2fi-1，j+fi-1，j+1]算子的模板为Gx=-101-202-101,Gy=121000-1-2-1]]>经过以上步骤就完成对图像的骨骼和皮肤的特征提取。
通过测量标准圆柱件来检验超声波测量的精度和图像特征提取的精度1)内轮廓精度分析用猪肉包埋一个直径为30mm的氟乙烯圆柱棒，将其置于储液桶3中，设定旋转角度为9°对其做超声波测量，采用上述图像二维重建的方法对得到的k＝40帧图像进行重建，图3a所示是重建的内圆柱横截面图像，提取内圆柱轮廓，测得直径大小为30.2mm。取旋转中心为圆心，直径为30作一标准圆，如图3b所示，虚线是标准圆，实线是提取的内圆柱轮廓线，轮廓线与标准圆十分吻合。两者的最大偏差在底部，最大偏差与尺寸误差见表1。
2)外轮廓精度分析测量一个直径为90mm聚胺脂圆柱件，采用上述图像二维重建的方法对其图像进行重建，结果如图4a所示，对图像进行轮廓提取，测量其轮廓直径大小为89mm，取旋转中心为圆心，直径为90作一标准圆，如图4b所示，虚线是标准圆，实线是提取的外圆柱轮廓线，轮廓线与标准圆十分吻合。两者的最大偏差处在左下侧，最大偏差与尺寸误差见表1。
表1 内外轮廓精度分析结果
综合1)、2)可得由超声波测量及特征提取得到的图像轮廓与实物轮廓的最大偏差距离δ＝±1.01mm，尺寸相对误差ε＝1.1％，如表1所示。上述的误差结果对残肢的CAD三维重建并不影响，超声波测量产生δ和ε是可接受的。
实施例1选用从活体剥离的新鲜猪前腿(平均半径约240mm，高度约160mm)作为测量对象，方法如下a.图像采集先将猪前腿浸没在储满水的储液桶3中，下端置于托盘上并用胶泥固定，标定旋转中心，测量旋转角设定为9°，探头下移距离为5mm，在每一水平方向上，探头以每9°的间隔逆时针绕着猪腿一周作旋转测量，得到k＝360°/9°＝40帧图像，像素大小是768×576，垂直方向上，探头每次下移5mm，一共测24个水平面，这样共可测得24×40帧图像，参见图5；
b.图像二维重建第一步先确定残体的旋转中心，通过测量标准圆柱件确定其圆心位置，如果圆心与探头聚焦中心重合，则该圆心即是旋转中心；第二步将k帧图像按各自测量时的旋转角9°绕旋转中心旋转到新坐标系，根据坐标旋转原理，其变换表达式为 上式中，(x，y)是原图像坐标系的横纵坐标点，(u，v)是旋转后的图像坐标系的横纵坐标点，再把上式写成矩阵形式为 其中 这里，Ra是旋转变换矩阵，(x，y)是单帧图像的坐标系，(u，v)是复合图像的坐标系；第三步对旋转到新位置的40帧图像进行叠加，设单帧测量图像为Ak，复合图像为Bk+1，Bk+1=Σl=1k+1Alc=Σl=1k+1aij(l)c]]>其中aij(l)表示m×n矩阵Al，m＝768，n＝576，1≤i≤m，1≤j≤n，1≤l≤k+1；第四步对像素值取整，表示为Ck+1＝{Bk+1}，其中{x}表示超过x的最小整数；第五步去噪，阈值设为128， Ck+1表示对复合后的图像像素取整和去噪后的图像像素矩阵，cij表示经过复合，去噪后的的图像像素在第i行，第j列上的灰度值；c.对图像进行处理，提取骨骼和皮肤的边缘特征用函数f(x，y)表示图像中任意一点的灰度，则该点的灰度变化梯度为▿f(x,y)=∂f∂x∂f∂yT=GxGyT]]>是梯度算子，T表示转置，Gx和Gy表示在x和y方向上的灰度值的变化程度，梯度是一个矢量，其幅值为|▿f|=[Gx2+Gy2]1/2]]>表示某点图像的边缘强度，满足下面三个条件既可找到图像的边缘点1.以该点为中心的3×3邻域中的边缘强度极大值小于某阀值，2.该点的边缘强度大于沿该点梯度方向的两个相邻像素点的边缘强度，3.这一点的相位值φ(x,y)=arctan(GyGx)]]>值小于45°，其中，图像矩阵每个点f(x，y)，离散化成f(i，j)，记为fi，j，梯度Gx，Gy离散化为Gx＝[fi-1，j-1+2fi，j-1+fi+1，j-1]-[fi-1，j+1+2fi，j+1+fi-1，j+1]Gy＝[fi+1，j-1+2fi+1，j+fi+1，j+1]-[fi-1，j-1+2fi-1，j+fi-1，j+1]算子的模板为Gx=-101-202-101,Gy=121000-1-2-1]]>经过以上步骤就完成对图像的骨骼和皮肤的特征提取。
实施例2被测物体同实施例1，方法如下a.图像采集，先将猪前腿浸没在储满水的储液桶3中，下端置于托盘上并用胶泥固定，标定旋转中心，测量旋转角设定为90°，探头下移距离为3mm，在每一水平方向上，探头以每90°的间隔逆时针绕着猪腿一周作旋转测量，得到K＝360°/90°＝4帧图像，像素大小是768×576，垂直方向上，探头每次下移3mm，一共测45个水平面，这样共可测得45×4帧图像。
以下步骤同实施例1。
实施例3被测物体同实施例1，方法如下a.图像采集，先将猪前腿浸没在储满水的储液桶3中，下端置于托盘上并用胶泥固定，标定旋转中心，测量旋转角设定为45°，探头下移距离为4mm，在每一水平方向上，探头以每45°的间隔逆时针绕着猪腿一周作旋转测量，得到K＝360°/45°＝8帧图像，像素大小是768×576，垂直方向上，探头每次下移4mm，一共测35个水平面，这样共可测得35×8帧图像。
以下步骤同实施例1。
实施例4被测物体同实施例1，方法如下a.图像采集，先将猪前腿浸没在储满水的储液桶3中，下端置于托盘上并用胶泥固定，标定旋转中心，测量旋转角设定为30°，探头下移距离为5mm，在每一水平方向上，探头以每30°的间隔逆时针绕着猪腿一周作旋转测量，得到K＝360°/30°＝12帧图像，像素大小是768×576，垂直方向上，探头每次下移5mm，一共测25个水平面，这样共可测得25×12帧图像。
以下步骤同实施例1。
实施例5被测物体同实施例1，方法如下a.图像采集，先将猪前腿浸没在储满水的储液桶3中，下端置于托盘上并用胶泥固定，标定旋转中心，测量旋转角设定为60°，探头下移距离为6mm，在每一水平方向上，探头以每60°的间隔逆时针绕着猪腿一周作旋转测量，得到K＝360°/60°＝6帧图像，像素大小是768×576，垂直方向上，探头每次下移6mm，一共测8个水平面，这样共可测得8×6帧图像。
以下步骤同实施例1。
权利要求
1.一种基于超声波测量的残肢骨骼与皮肤特征的提取方法，采用超声波测量装置(1)、B超诊断仪(4)、超声影像工作站(5)和图形处理系统(6)组成超声波测量系统，超声波测量装置(1)上的探头(2)与B超诊断仪(4)连接，B超诊断仪(4)通过数据线与超声影像工作站(5)和图形处理系统(6)连接，其特征在于该方法包括以下步骤a.图像采集探头(2)围绕残体的旋转中心以每0°＜α°≤90℃的间隔做圆周运动，同时由超声影像工作站(5)对超声测量图进行采集、存储，得到k＝360°/α帧B超图像，然后探头(2)垂直下移，残体的旋转中心也随之下移，探头(2)又围绕新的旋转中心以每0°＜α°≤90°的间隔做旋转测量，获得k＝360°/α帧B超图像，以此类推，即得到被测残体各水平面的多帧二维截面断层图像；b.图像二维重建第一步先确定残体的旋转中心，通过测量标准圆柱件确定其圆心位置，如果圆心与探头聚焦中心重合，则该圆心即是旋转中心；第二步将k帧图像按各自测量时的旋转角α绕旋转中心旋转到新坐标系，根据坐标旋转原理，其变换表达式为 上式中，x，y是原图像坐标系的横纵坐标点，u，v是旋转后的图像坐标系的横纵坐标点，再把上式写成矩阵形式为uv=Rα·xy]]>其中Rα=cosαsinα-sinαcosα]]>这里，Ra是旋转变换矩阵，x，y是单帧图像的坐标系，u，v是复合图像的坐标系；第三步对旋转到新位置的k帧图像进行叠加，设单帧测量图像为Ak，复合图像为Bk+1，Bk+1=Σl=1k+1Alc=Σl=1k+1aij(l)c]]>其中αij(l)表示m×n矩阵Al，m＝768，n＝576，1≤i≤m，1≤j≤n，1≤l≤k+1；第四步对像素值取整，表示为Ck+1＝{Bk+1}，其中{x}表示超过x的最小整数；第五步去噪，阈值设为128， Ck+1表示对复合后的图像像素取整和去噪后的图像像素矩阵，cij表示经过复合，去噪后的的图像像素在第i行，第j列上的灰度值；c.对图像进行处理，提取骨骼和皮肤的边缘特征用函数fx，y表示图像中任意一点的灰度，则该点的灰度变化梯度为▿f(x,y)=∂f∂x∂f∂yT=CxCyT]]>是梯度算子，T表示转置，Gx和Gy表示在x和y方向上的灰度值的变化程度，梯度是一个矢量，其幅值为|▿f|=[Gx2+Gy2]1/2]]>表示某点图像的边缘强度，满足下面三个条件既可找到图像的边缘点1.以该点为中心的3×3邻域中的边缘强度极大值小于某阀值，2.该点的边缘强度大于沿该点梯度方向的两个相邻像素点的边缘强度，3.这一点的相位值φ(x,y)=arctan(GyGx)]]>值小于45°，其中，图像矩阵每个点fx，y，离散化成fi，j，记为fi，j，梯度Gx，Gy离散化为Gx＝[fi-1，j-1+2fi，j-1+fi+1，j-1]-[fi-1，j+1+2fi，j+1+fi-1，j+1]Gy＝[fi+1，j-1+2fi+1，j+fi+1，j+1]-[fi-1，j-1+2fi-1，j+fi-1，j+1]算子的模板为Gx=-101-202-101]]>Gy=121000-1-2-1]]>经过以上步骤就完成对图像的骨骼和皮肤的特征提取。
全文摘要
一种基于超声波测量的残肢骨骼与皮肤特征的提取方法，采用超声波测量装置、B超诊断仪、超声影像工作站和图形处理系统组成超声波测量系统，超声波测量装置上的探头与B超诊断仪连接，B超诊断仪通过数据线与超声影像工作站和图形处理系统连接，该方法为首先在每一水平面对残肢进行旋转测量，获得多帧超声测量图像，然后对多帧测量图像进行复合，重建其二维截面图像，进而利用边缘检测法对骨骼与皮肤进行特征提取。本发明利用超声波成像原理，对残肢进行超声波测量，解决了目前无法对残肢骨骼和皮肤同时取型的问题，为残肢的三维CAD建立奠定了基础，而且极大地降低了成本、减少了辐射，为假肢康复工程和医学图像的重建提供了一种新的手段。
文档编号A61B8/00GK1608589SQ20041007316
公开日2005年4月27日 申请日期2004年10月11日 优先权日2004年10月11日
发明者赵万华, 卢秉恒, 郑淑娴, 范晓逶 申请人:西安交通大学