Article Info

1.东南大学生物科学与医学工程学院，江苏南京 210018；2.南京中医药大学附属医院，江苏南京 210029；3.南京市口腔医院南京大学医学院附属口腔医院，江苏南京 210001

[摘要]为了更全面便捷地衡量动物临床试验骨质疏松情况，开发了一套基于MircoCT图像的用于测量骨小梁参数的软件系统。在保留骨结构完好的前提下，使用12组成年猪额骨及颌骨样本进行MicroCT扫描成像，并根据二维断层图像，经过导入数据、三维可视化、圈定感兴趣区、二值化处理等步骤后进行骨密度、骨小梁厚度等参数测量。测量结果与BoneJ进行对比，结果验证了本系统具有较高的测量准确性（各项参数相对误差均＜5%），可视化效果直观，在动物骨小梁参数测量方面具有较高的应用价值。

[关键词]骨小梁；距离变换；三维可视化；骨密度；骨小梁厚度

骨小梁是松质骨的主要结构，其形态结构参数是研究松质骨病变，判断骨质疏松情况的重要参考依据。研究动物骨小梁在试验前后的变化情况对于临床医学具有重要意义。MircoCT（又称微型CT或显微CT），是指可测定体素空间分辨率小于100的X线三维成像系统，具有空间分辨率高、焦点小、辐射剂量少等特点，是研究动物骨小梁结构的有力辅助工具[1]。而对于动物骨小梁的结构参数测量，除了比利时SkyScan有配合测量其MircoCT图像的对应骨分析软件以外，目前国内对动物骨小梁结构多数还在使用显微切片或扫描电镜观察等二维图像分析方法，缺少对动物骨小梁完整三维结构测量和观察的辅助工具[2]。因此，针对MircoCT扫描骨小梁图像，开发了一套基于C++的骨小梁参数测量系统，能够测得骨密度、骨小梁厚度、骨小梁间隙、骨小梁数量等多个参数，可以对动物骨病变情况有更全面的衡量。

1 测量原理

骨密度，全称骨骼矿物质密度，是反映骨质疏松程度，衡量骨骼健康情况最为重要的一个参数。使用定量CT（QCT）所测密度为该区域的体密度，单位为g/cm3。

CT成像中当X线束进入物体时，将受到物体对X线的吸收和散射，并且以吸收为主。X线束透过物体后剩余能量将被检测器所接受作为投影数据，投影数据经过CT重建后获得最终CT扫描图像，因此扫描图像的CT值可以直接反映物体对X线束的吸收情况。而物体对X线的吸收程度与物体密度、原子系数以及X线的吸收能量等参数密切相关[3]。在同一CT相同电压电流的拍摄条件下，可以保证X线的吸收能量相同。使用专业制造与骨成分近似的模体对比拍摄，可以保证原子系数相同。在控制其他变量的情况下可视为物体对X线的吸收程度与物体的密度成正比。即得到以下公式：

其中，CTLevel表示当前测量区域平均CT值，BMDh表示模体高密度部分密度，BMD1表示模体低密度部分密度，CTh表示模体高密度部分测量CT值，CT1表示模体低密度部分测量CT值。计算结果的单位通常用mg/cm3表示。其中，对于比较好的骨密度模体，有不止两个高低密度值的情况，也可以采用最小二乘法来求得目标区域骨密度。

骨小梁厚度指的是骨小梁的平均厚度，是衡量骨小梁微观形态最主要的参数之一，它与骨小梁间隙和骨小梁数量这两个参数在计算上是紧密联系的。结合骨小梁间隙（Th.Sp）与骨小梁数量（Th.b）这两个参数可以判断骨质疏松程度。当发生骨质疏松时，骨小梁厚度值会降低，而骨小梁间隙增大，骨小梁数量减小。

计算骨小梁厚度可采用Hildebrand和Ruegsegger[4]（1997）所提出的算法。这种算法基于MicroCT拍摄的3D图像分析，不同于过去借助形状建模的算法，这种独立于模型的算法更加精确。根据该算法，对物体上一个点的局部厚度可以定义为满足以下两个条件的最大球体的直径：① 球体包含该点（该点不一定是球心）；② 球体边界在物体内。

基于这个定义，一段骨小梁结构的厚度即骨小梁内骨架上所有点局部厚度的平均值。而图像骨架可以使用距离变换法（Distance Transform）求出[5]。距离变换法指的是对二值图像，计算每个前景像素最近的背景像素距离，将其变换为灰度图像，即距离图像的一种算法。

通过距离变换求得灰度图像的峰值即原二值图像的骨架，对骨架上的点遍历计算局部最大球体直径，求和平均即得到骨小梁厚度。计算骨小梁间隙的过程与骨小梁厚度相同，只是需要先将图像取反再进行运算。而骨小梁平均数量则是由以下公式得出[6]：

其中，Tb.Th表示骨小梁厚度，Tb.Sp表示骨小梁间隙，因此骨小梁数量常用单位为 1/mm。

结构模型指数（Structure Model Index，SMI）表达的是骨小梁的三维形态相对而言接近圆盘状或圆柱状的程度。当出现骨质疏松时，骨小梁的形态结构将会从盘状向杆状转变，而这个参数就能用于衡量骨小梁的形变程度。一个理想的圆盘其SMI值为0，一个理想的圆柱其SMI值为3，一个理想的球体其SMI值为4。相应地，一个理想的圆柱状孔洞SMI值为-3，一个理想的圆柱状孔洞SMI值为-4。因为SMI涉及到表面凸面曲率的测量，如果目标是凹面，则其测量值为负数。

计算SMI是基于三维体素模型的膨胀求得的。首先将目标松质骨进行二值化，计算其体积与表面积，然后对二值化后物体进行三维膨胀，再次计算体积与表面积（这个步骤与计算形状指数Tb.Pf的过程相同），然后SMI可用如下公式求出[7]：

其中，

表示膨胀前后表面积之差，V表示膨胀前骨小梁体积，S表示膨胀前骨小梁表面积。

需要说明的是对于封闭孔洞而言，其表面是凹陷的，因此求得的SMI也是负数，由于对这样的封闭空间膨胀后物体的表面积变小了，所求得的表面积差S'也就是负数。因此对于一个包含封闭孔洞超过50%以上空间的物体（如骨小梁），最终SMI参数将是一个负数。也可以说，参数SMI与体积百分比关系紧密。由于人为圈定目标区域产生的边和角改变了测量物体的体积，也会增大测量的SMI值。

2 测量方法

实验仪器：本次实验使用的MircoCT为东南大学生物科学与医学工程学院利用国家自然科学基金科学仪器研究专项资金研发的高分辨显微CT（Hiscan M1000），最大分辨率为20 μm，最大扫描范围可达到60 mm×200 mm（多次扫描）。图像重建及后处理软件为苏州海斯菲德信息技术有限公司提供的MCT_Recon_Sever（版本号1.0.10）。

材料制备：从正常猪体内使用锯子分离获得颌骨12组，然后采用浓度4%的福尔马林浸泡，保持其骨结构不变，骨矿含量不流失（由于骨组织与非骨组织在图像灰度值上有很大差异，因此在扫描前不需要特意剔除骨表面的其他组织）。用锯子将颌骨样本切割为长度不超过15 cm的骨块，使其能够完全放入CT的聚丙烯样本槽。

扫描过程：将12组猪颌骨样本依次置于聚丙烯样本槽中，将扫描X 射线的能量设置为400 kV，电流设置为200 μA，扫描旋转360°，采集720张投影图，每张投影图的曝光时间为50 ms，每个样本平均扫描时间3min。使用MCT_Recon_ Sever对每个样本重建，平均获得800张大小为1874×1874像素的断层重建数据。重建图像分辨率为25 μm，Z方向间隔25 μm。

（1）首先导入CT断层扫描图像数据，导入数据将以三视图方式显示。在三视图中选择目标区域，可以对该区域数据进行三维体绘制显示。如果数据是已经分割完成的感兴趣区（Region of Interest，ROI）数据，则可以直接进入二值化处理环节。在三维可视化环节，通过旋转位移调整物体，可以直观地观察骨小梁形态，选择最合适的角度获得ROI的二维断层图像。然后在二维断层图像上手动圈定ROI区域。圈定时程序会根据手动圈定结果对中间所有层进行插值运算获得所有层的圈定结果。圈定ROI完成后，将灰度图像进行二值化并保存数据（灰度图像数据保留在内存中不会删除），然后按照上文所述算法进行不同的参数测量。

导入数据：由于图像分辨率较高，单张DCM图像大小可达到5～6 MB，全部数据导入将占用1 GB以上内存导致程序卡顿，因此在导入数据过程中并非将图像数据全部读入内存，而是使用了“预加载”技术，只读入DCM图像中的tag信息，在实际需要该层图像数据时再读取真正的图像数据，极大提升了加载速度。同时考虑到内存空间限制，对于过大的数据会进行自动降采样减小内存占用。

（2）三维可视化的操作界面，见图2。为了直观地进行ROI圈定和观察目标松质骨的骨质疏松情况，需要先对图像进行三维可视化。使用OpenGL库进行体绘制来实现三维可视化，其具体方法是导入数据经采样和插值后，将每个体素构造为独立的理想化物理模型，考虑其介质属性，按照phone光照模型调节传递函数对每个体素分配光强和不透明度，然后沿视线观察方向积分，最后形成半透明的三维投影图像。

（3）圈定ROI的操作界面，见图3。计算骨密度、骨小梁厚度等参数，需要保证计算的图像范围只包含目标松质骨而不含有其他结构或空腔，否则会导致计算结果偏小。因此需要对ROI进行准确地圈定。圈定ROI时，首先将三维显示的松质骨结构按照垂直屏幕方向获得二维断层图像，浏览断层图像选择ROI起始层和终止层，在起始层和终止层上再次手动圈定区域，将会在中间每一层按照渐变插值计算出该层的圈定区域，若中间层的圈定效果不满意，可以在该层手动圈定，然后按照新的圈定结果重新计算每一层的圈定区域。圈定完成后，程序将ROI区域以外的数据清零，保证进行计算时不受到其他结构影响，并自动将当前图像保存为新的DCM序列便于再次测量和查看结果。

（4）二值化处理的操作界面，见图4。在计算骨结构参数前，需要将骨组织和非骨组织区分，故先将图像进行二值化处理。具体过程为通过手动对断层图像调窗，将非骨部分数据置为0，将骨部分数据置为1，同时将调窗结果同步至三维可视化窗口，直观显示调窗效果，便于判断是否将骨组织和非骨组织区分开。图像二值化数据可以保存至硬盘便于再次查看和测量。

（5）参数测量的操作界面，见图5。参数测量分为“骨密度测量”和“骨小梁参数测量”两个部分，其中“骨密度测量”使用未二值化的原始图像数据，通过计算区域灰度平均值，再导入模体图像数据，计算模体图像灰度值并输入模体密度，按照第一节所述原理进行密度计算。“骨小梁参数测量”使用原始图像二值化后数据，经过去噪（消除图像中由于CT扫描造成的非骨部分的噪点），按照第一节所述原理对所有骨小梁参数进行测量。最后输出测量结果并生成测量报告。

3 测量结果

猪颌骨三维效果示意图，见图6。在右侧三视图中框选出松质骨的大致区域后，可以看到完整松质骨的三维可视化效果，左下角标记了当前物体的旋转情况，确定好旋转角度后将按照垂直屏幕方向生成新的断层图像。相比在CT断层图像上选择ROI的方法，使用三维可视化进行ROI选择效果更加清楚直观，操作更加方便准确。

ImageJ是一个基于Java的公共的图像处理软件[8]，它是由美国国立卫生研究院（National Institutes of Health，NIH）的Wayne Rasband等人开发，可运行于Microsoft Windows、Mac OS、Mac OS X、Linux和Sharp Zaurus等多种平台，在世界范围的科学研究领域中广泛应用。其中BoneJ属于ImageJ中专门测量骨结构参数的插件[9]，由 Doube M、K?osowski MM、Arganda-Carreras I等人研发。其测量结果具有一定的权威性和参考价值，本文采用BoneJ测量来验证本系统的测量准确度。

测量结果见表1，与BoneJ测量结果对比，12组标本骨密度的平均相对误差为0.0%，骨小梁厚度的平均相对误差为4.275%，骨小梁间隙的平均相对误差为1.127%，骨小梁数量的平均相对误差为1.952%，骨小梁结构指数的平均相对误差为3.954%。

分析误差原因，主要有如下几个方面：① 计算骨表面积方法不同，BoneJ中使用三角面片法计算骨表面积，而本系统使用体素统计计算表面积。表面积计算的误差直接影响到参数计算结果；② 圈定ROI时有部分骨组织被截断，形成边际效应影响测量结果。

4 结论

本测量系统能够对MicroCT断层扫描获得的动物骨小梁图像进行骨密度、骨小梁厚度、骨结构指数等参数进行全面准确的测量，并提供骨小梁结构的三维可视化效果，为临床研究动物骨质疏松程度、骨健康情况提供了一个更全面的参考工具。

表1 猪颌骨各项参数测量结果

[1]高鹏,戎军艳,廖琪梅,等.MicroCT系统软件平台的设计与实现[J].医疗卫生装备,2014,35(12):22-24.

[2]付鑫,马剑雄,董宝康.骨微结构检测方法研究进展[J].中国骨与关节外科,2009,2(6):509-515.

[3]Ma XQ,Overton TR.Automated Image Analysis for Bone Density Measurement Using Computed Tomography[J].IEEE Trans Med Imaging,1991,10(4):611-616.

[4]Hildebrand T,Rüegsegger P.A new method for the modelindependent assessment of thickness in three-dimensional images[J].J Microsc,1997,185(1):67-75.

[5]Borgefors G.Distance transformation in arbitrary dimensions[J].Comput Vis Graph Image Process,1984,27(3):321-345.

[6]Hahn M,Vogel M,Pompesius-Kempa M,et al.Trabecular bone pattern factor-a new parameter for simple quantification of bone microarchitecture[J].Bone,1992,13(4):327-330.

[7]Hildebrand T,Rüegsegger P.Quantifcation of Bone Microarchirecture with the Structure Model Index[J].Comput Methods Biomech Biomed Engin,1997,1(1):15-23.

[8]Schneider CA,Rasband WS,Eliceiri KW.NIH Image to ImageJ:25 years of image analysis[J].Nat Methods,2012,9(7):671-675.

[9]Doube M,K?osowski MM,Arganda-Carreras I,et al.BoneJ:Free and extensible bone image analysis in ImageJ[J].Bone,2010,47 (6):1076-1079.

Analysis of the Application Effects of Trabecular Bone Measurement System Development Based on MircoCT

Abstract：In order to measure osteoporosis cases in animal clinical trials,the research developed a trabecular bone parameters measurement software system based on MircoCT images.On the premise of keeping skeletal structure intact,the paper carried on MicroCT scan on 12 samples of adult pig frontal bones and maxillary bones.Based on the two-dimensional sectional images,measurement of parameters such as bone density and trabecular thickness was conducted through the process of data import,threedimensional visualization,and delineation of region of interests,and binary image processing.The results of the measurement was compared with BoneJ,which confrmed the accuracy of measurement system (relative error ＜5% for each parameter).The system also allowed intuitive visualization and had a high application value in terms of animal trabecular parameter measurement.

Key words：trabecular;distance transform;three-dimensional visualization;bone density;trabecular thickness

WU Pei-ze1,LUO Shou-hua1,CHEN Gong2,LIU Hui-fen3,CHEN Bin3
1.Department of Biological Science and Medical Engineering,Southeast University,Nanjing Jiangsu 210018,China;2.The Affiliated Hospital of Nanjing University of Traditional Chinese Medicine,Nanjing Jiangsu 210029,China;3.Nanjing Stomatological Hospital,Medical School of Nanjing University,Nanjing Jiangsu 210001,China

基金项目：国家自然科学基金（61127002，61179035）；中央高校基本业务费重点项目培育计划（021414310017）。

基于MicroCT的骨小梁参数测量系统的应用效果分析

1 测量原理

2 测量方法

3 测量结果

4 结论