基于明胶实验的弹创空腔序列图像处理技术*

沈 杰，王 涛，廖振强

(南京理工大学机械工程学院，南京 210094)

0 引言

空腔形成理论指出，高速飞行的杀伤元侵彻人体组织时，在组织中形成的瞬时空腔会经多次膨胀与压缩，使空腔周围组织在极短的时间内经历多次挤压与牵拉[1]。在瞬时空腔膨胀波及的范围内，人体组织受到猛烈扯动，发生明显的位移，甚至会撕裂和切断，该区域范围是组织器官受到严重创伤的区域[2]，因此研究瞬时空腔膨胀所能达到的最大空间范围具有重要意义。为了检测出弹创空腔最大膨胀波及范围的边缘轮廓，针对单帧图像表示的瞬时空腔轮廓精确提取困难，且不能有效反映最大膨胀波及范围等问题，文中提出序列瞬时空腔轮廓“或”运算叠加法，算法重点在最终的空腔轮廓叠加结果，对单幅图片处理的轮廓精度要求可以降低，而且在程序设计上更容易实现。

1 瞬时空腔最大膨胀波及范围的定义

如图1所示，设侵彻时刻从t1到t2时，弹创空腔轮廓由实线膨胀到虚线处，则在侵彻距离l1处，对应t1和t2时刻的空腔直径分别为Dl1t1和Dl1t2，在空腔膨胀周期内，此距离处的空腔直径不断变化，故在此侵彻距离上必然存在一个空腔直径的最大值Dl1max，同样在侵彻距离为l2处也必然存在一个空腔直径的最大值Dl2max，依次类推，即可得到一组不同侵彻距离处空腔直径的最大值，最后根据得出的这一组数值可以绘制出一个空腔轮廓，该轮廓即为瞬时空腔最大膨胀波及范围。

图1 瞬时空腔最大膨胀波及范围的定义

2 基于“或”运算的瞬时空腔最大膨胀波及范围轮廓提取

文献[3]提出了一种通过拍摄瞬时空腔在两个正交方向的投影图像，将空腔截面进行椭圆近似对空腔三维重建，为定量研究武器致伤效应提供了依据。但是，由于瞬时空腔最大膨胀波及范围是瞬时空腔在整个膨胀过程中所形成的，并不是某个时刻的瞬时空腔。传统的图像处理方法需对序列图像中的每帧图像进行精确的轮廓边缘检测，比较各个侵彻距离处的轮廓大小，得到一组最大直径值。因为要保留所有轮廓数据作比较，对每幅图片处理的轮廓精度要求较高。由于受实验环境的限制较大，得到的序列图片大多不够清晰，噪声较多，对每帧空腔图像进行精确的轮廓边缘检测尚无很好的通用算法。文献[4]探讨了创伤弹道空腔X射线投影图像的边缘检测技术，常规的边缘检测算子Sobel算子、Log算子、Canny算子不能有效的检测出封闭边缘。最后提出的基于简化M-S模型的改进水平集方法，也不适用于文中的高速摄影图像。

文中提出的算法先对单帧图像进行Canny算子边缘检测得到空腔轮廓，再对序列空腔轮廓图进行“或”运算叠加，在叠加的过程中进行适当的除噪处理，去除多余的伪边缘，尽可能的提高轮廓边缘检测的精度，最后通过滤波去噪和连接算法得到所需的轮廓，为空腔容积计算提供数据基础。

3 单帧图像瞬时空腔轮廓提取

3.1 Canny算子边缘检测

明胶经脱色处理后灰度均匀，瞬时空腔膨胀时轮廓边缘与明胶有较大的灰度差异。文中用信噪比大、检测精度高的Canny[5]算子检测边缘。由于Canny算子首先对图像进行高斯平滑，对噪声有一定的抑制能力。同时使用高阈值判断边缘点，用低阈值连接边缘点，避免了大量虚假边缘的产生，可以比较精确的检测出边缘。

图2 高速摄像拍摄的原始图像

图3 Canny算子边缘检测

杀伤元在明胶中飞行时，所形成的空腔内部也产生灰度变化，因此使用Canny算子检测出来的边缘除了空腔的轮廓边缘还含有很多非感兴趣边缘，如在瞬时空腔内部检测到的边缘，如图3。而瞬时空腔的膨胀范围只与其外轮廓有关。为了减少伪边缘带来的影响，文中在用Canny算子检测后的边缘图像基础上提出空腔外轮廓提取算法。

3.2 空腔外轮廓提取

将一幅经Canny算子检测后的边缘图像表示成大小为m×n的矩阵A:

1)产生一幅与边缘图像大小相等的临时图像Btemp，并将所有元素值置为0。

2)对第j列元素，从第1行向第m行检测，记录检测到的第一个值为1的元素，记下该元素的位置(iup，jup)。从第m行向第1行检测，记录检测到的第一个值为1的元素，记下该元素的位置(idown，jdown)。

3)在临时矩阵Btemp中将位置为(iup，jup)，(idown，jdown)的元素置为1。

4)遍历每列元素，重复2)、3)过程，最终得到的矩阵Btemp即为空腔轮廓外边缘图像。

图4 单帧图像轮廓提取

经外轮廓提取后，得到的单帧瞬时空腔外轮廓如图4所示。

4 空腔轮廓图像“或”运算

经轮廓提取后，得到了单帧图像的瞬时空腔轮廓图。将处理得到的序列空腔轮廓图像进行“或”运算叠加，多帧图像的瞬时空腔轮廓由于存在重叠、相邻，表现为轮廓变粗。

在杀伤元侵彻明胶时，由于明胶抖动或其他外界条件的影响，有时检测出来的单帧瞬时空腔轮廓并不精确，如图5(c)空腔下边缘轮廓并不能很好的检测出来，同时不可避免的还存在一些伪边缘。但是由于处理的不是单帧图像，而是连续多帧图像的“或”运算叠加，少数不精确的空腔轮廓可以由多数相对精确的轮廓弥补，随着序列图像帧数的增多，空腔轮廓继续变粗，表现得更突出。伪边缘由于随机性，在对瞬时空腔轮廓图进行“或”运算时，伪边缘重叠、相邻的概率小，呈现为一些细小的边缘，为了防止伪边缘被错误的加强，可设置一定的除噪周期，选取合适的结构元素对“或”运算叠加后的图像进行形态学开运算，消除不重要的图像细节，同时保证图像的全局几何形状不失真。

图6 连续30帧图像叠加后进行轮廓提取

5 瞬时空腔最大膨胀波及范围轮廓曲线滤波去噪

“或”运算叠加图经过开运算后伪边缘得到了抑制，再对其进行空腔轮廓提取，由空腔外轮廓提取算法可知，所得的轮廓曲线即为瞬时空腔最大膨胀波及范围轮廓。由于“开”操作以及Canny算子本身的局限，最终得到的空腔轮廓不连续且不可避免存在噪声点如图6(c)。提取得到的空腔轮廓由上、下两条轮廓边缘组成，每条轮廓上都存在相应的噪声点，噪声点表现为偏离空腔轮廓，传统的对曲线滤波方法有中值滤波、均值滤波、频域滤波，针对各种滤波方法所呈现出来的问题，文中提出一种基于斜率滤点的方法，取轮廓曲线上含有噪声点的一部分数据进行对比分析，结果表明基于斜率滤点的方法最有效。

5.1 中值滤波与均值滤波

通过一个采样窗口从原始数据中提取奇数个数据 y(i－ k)，y(i－ k+1)，…，y(i)，y(i+1)，…，y(i+k)，中值滤波是对这些数据按大小进行排序，用排序后的中间值作为输出值Y(i)，均值滤波是将提取数据的平均值作为输出值Y(i)。将采样窗口遍历所有数据。为了使新的数据与原来的数据个数相同，先对原来数据两端扩展k个数据，其值为端部的数据。

图7 中值滤波和均值滤波

中值滤波法可以去除一些噪声点，特别是当采样窗口中的噪声点个数小于采样窗口长度一半时，滤波效果显著，但对于宽度大的噪声曲线段，需加大窗口，但是，若盲目的加大窗口，对于无噪声段曲线，由于窗口的加大，会使曲线失真。对于此类无法预知宽度的噪声，中值滤波并不能完全适用。同样均值滤波有一定的滤波能力，使曲线趋于平滑，但是同样不能去除宽度大的噪声曲线段。

5.2 频域滤波

在信号处理中，噪声能量一般集中于高频，可设计低通滤波器，将高频段的噪声滤掉，使低频分量顺利通过而有效地阻止高频分量，通过频域处理去除噪声点。

1)空腔轮廓点是一系列的离散点，对离散点进行离散傅里叶变换，得到F(u):

2)建立低通滤波器H(u); (3)

3)G(u)=F(u)H(u); (4)

4)对G(u)进行离散傅里叶反变换;

图8 频域滤波

图8分别选用巴特沃斯低通滤波器和高斯低通滤波器滤波。由于噪声并不只是一些阶跃的点，有时也表现为偏移正常曲线的光滑曲线段，光用低通滤波的方法，并不能有效地去除噪声。

5.3 基于斜率滤波

文中提出基于斜率的滤波方法，分别对上、下两条轮廓曲线从曲线上某起始点向左右方向滤噪声点，可达到预期效果。

具体实现算法步骤如下:

1)对每条轮廓线，从连接良好，不存在噪声处的轮廓线上选取两相邻边缘点 Pi、Pj，j=i+1。

2)设置斜率阈值Δk。

3)计算Pi与左边相邻点Ph的斜率绝对值Ki，Pj与右边相邻点Pk的斜率绝对值Kj;

4)判断Ki是否大于 Δk，若是，则 Ph为噪声，去除。否则Ph为轮廓边缘点，保留，同时i=i－1;判断Kj是否大于 Δk，若是，则 Pk为噪声，去除。否则 Pk为轮廓边缘点，保留，同时j=j+1。

5)重复3)、4)过程，直到Pj为空腔轮廓最右边的边缘点，Pi为空腔轮廓最左边的边缘点。

6)通过插值的方法将原来噪声点的位置补上插值点实现轮廓边缘连接。

图9 基于斜率滤波

该算法不但能有效的去除偏离曲线的单个噪声点，而且对一些偏离曲线的大宽度噪声曲线段也能较好的去除，同时对于准确的边缘点不做处理，最大限度的保留轮廓曲线的真实形状。

6 实验结果与分析

对瞬时空腔剧烈膨胀成球形的序列图像处理如图10，由于拍摄条件好，空腔轮廓与明胶有较大的灰度差别，用σ=4的高斯滤波器先对图像平滑，在前5帧叠加时，没有带入伪边缘，叠加了10帧后，产生了细小的伪边缘，设定开运算周期为10，即每叠加10帧就对叠加后的边缘图进行开运算，随着帧数的增加，空腔轮廓加强越明显，伪边缘也不断地去除，对90帧空腔轮廓图像或运算结果如图10(c)，最终提取的空腔轮廓如图10(d)，该实验证明用文中提出的“或”运算叠加法能有效获得瞬时空腔的最大膨胀范围轮廓。

图10 瞬时空腔膨胀呈球形

图11 瞬时空腔膨胀呈锥形

瞬时空腔膨胀呈锥形的序列图像处理如图6。该实验中各帧图像包含的噪声较多，因此在图6(a)所示的30帧空腔轮廓“或”运算叠加结果中可以清楚地看到若干伪边缘存在，进行周期除噪后还有一些噪声点残留，使得图6(c)所示的最终提取的空腔轮廓并不完整，据此计算的空腔容积将产生较大的误差。用文中提出的基于斜率的曲线滤波算法得到的最终空腔轮廓如图11所示，对比可知，最终空腔轮廓得到了很好的改善。

7 结束语

文中针对杀伤元造成组织严重损伤区域的形成特点，提出了序列瞬时空腔轮廓“或”运算叠加处理方法，使瞬时空腔最大膨胀波及范围轮廓的提取重点落在最终叠加结果图的处理上，避免了需要精确检测单帧图像轮廓，通过基于斜率的滤波方法可以得到最终较为精确的轮廓，为空腔容积计算提供了数据基础。若对实验进行正交拍摄，采用此方法也可对弹创瞬时空腔最大膨胀波及范围进行三维重建。文中提出的方法可为枪弹的致伤效果提供一定的评价参考，也为各种杀伤武器的研究设计提供指导。

[1]安波，蒋浩征，李杨.高速投射物对生物致伤时瞬时空腔的理论研究[J].北京理工大学学报:英文版，2001，10(3):272－277.

[2]Fernando Spencer Netto，Dylan Pannell，Homer C Tien.Hollow-point ammunition and handguns:The potential for large temporary cavities[J]. Injury Extra，2008，39(2):50－52.

[3]林勇，王经瑾，宋征，等.明胶弹创空腔闪光 X射线投影图像三维重建[J].清华大学学报:自然科学版，2002，42(12):1576－1578.

[4]贺成，王涛，廖振强，等.创伤弹道空腔图像边缘检测技术研究[J].计算机工程与设计，2011，32(1):248 －250.

[5]Rafael C Gonzalez，Richand E Woods，Steven L Eddins.数字图像处理(MATLAB版)[M].阮秋琦，译.北京:电子工业出版社，2005.