CN102609980A - 混凝土ct图像三维重构方法 - Google Patents

Abstract

本发明公开了一种混凝土CT图像三维重建的方法，首先用CT机对混凝土试样进行断层扫描；其次光线投影算法对混凝土CT图像进行三维重建，得到骨料、砂浆、孔洞的真实图像；最后利用可视化工具VTK对混凝土CT图像进行三维立体重构。本发明的积极效果是：实时性强、运算效率高、能清晰可靠的分辨骨料、砂浆及孔洞。三维重建结果得到了混凝土真实的细观结构模型，可以与混凝土细观数值模拟结果进行试验比较，并校正数值分析模型。

Description

混凝土CT图像三维重构方法

技术领域

本发明涉及混凝土CT图像处理技术，尤其涉及一种根据对象的二维断层图像数据重构其三维图像的方法。

背景技术

混凝土是一种工程应用相当广泛的不均匀材料。由粗骨料、细骨料、水泥水化物、未水化水泥颗粒、孔隙及裂纹等组成的非均质复介材料。研究各种空间尺度缺陷对混凝土材料的变形、破坏和稳定性所起的作用，或者作用的大小，是目前混凝土材料研究所面临的最深刻而又最艰难的课题之一。近些年X射线混凝土CT试验成为混凝土、岩石材料细观破裂过程的热点研究课题，X射线CT观察混凝土裂纹演化过程最大优点在于它的无损探测性能和较高分辨率。对于实验后的得到混凝土的CT图像数据分析，采用直接图像分析法。直接图像分析是指直接对混凝土受力后的结构变化通过图像分析。然而混凝土CT图像中的许多萌生细微裂纹都无法直接从图像中直接准确的判断是否存在；另外由于CT机的扫描原理，所得到混凝土CT图像只是横切面图片，对于混凝土体柱受力后图像分析，无法直接进行，要想通过二维断层图像认识真实物体，就必须完成从二维图像到三维图像的恢复工作，需要对混凝土CT图像进行三维重构。三维重建结果得到的混凝土真实的细观结构模型，可以与混凝土细观数值模拟结果进行试验比较，并校正数值分析模型。

体绘制算法用于完成三维结果图像的生成和显示。光线投影算法是最常用的以图像空间为序的典型直接体绘制算法。

利用光线投射法绘制出的图像质量很高、效果很好，但是此方法存在很大的问题，那就是处理的数据量比较大，对内存的要求比较高，由于要对每条光线进行重采样计算，所以绘制速度比较慢。如何进行算法的改进以加快绘制速度已经成为目前光线投射算法研究中最需要解决的问题，也是体绘制技术中最主要的研究方向。目前，人们最常用的加速技术主要从两个方面考虑：一方面从图像空间的相关性着手使发射出的光线数目尽量减少，另一方面从对象空间的相关性着手，省略灰度、梯度方向相同或相近体素的采样，从而使获得的采样点数目尽量减少。因此，提供一种在不影响绘制图像质量的情况下，从这两方面考虑对传统的光线投射算法进行改进，使绘制速度得到很大提高的方法成目前亟待解决的技术问题。

发明内容

针对现有光线投射算法中存在的问题和不足，本发明的目的在于提供一种混凝土图像重构方法，利用对象的二维断层图像数据重构其三维图像。运用三维重构可视化技术，对混凝土CT图像进行立体重构，真实反映出混凝土柱应力后结构的变化。为分析混凝土的裂纹扩展及破坏形态，揭示出混凝土材料的损伤破裂机理提供了一种重要的辅助分析手段和方法。

本发明的技术方案是：一种混凝土CT图像三维重构的方法，包括步骤如下：

1)采用医用CT或工业CT扫描混凝土试件断面，获得混凝土静力压缩CT二维断层图像数据；

2)运用三维重构可视化技术对混凝土静力压缩CT二维断层图像进行立体重构，生成混凝土CT三维显示图像；

3)对混凝土CT三维显示图像进行三维优化，最终完成混凝土CT图像三维重构。

本发明进一步的特征在于：

所述步骤2)中，运用三维重构可视化技术对混凝土CT二维断层图像进行立体重构，采用光线投影算法按照设定视线的方向，从混凝土静力压缩CT二维屏幕上的每个像素点出发，发出一条射线，沿着该射线按照一定的步长进行等距采样；对距离某个采样点最近的八个体素的属性值做三线性插值运算从而计算出，得到该采样点的属性值(颜色值和不透明度值)，对得到的每个采样点的属性值进行逐点合成，从而得到所对应的屏幕上该像素点的颜色值和不透明度值，最终生成混凝土CT三维显示图像。

上述对得到的每个采样点的属性值按照由后向前的图像合成和由前向后的图像合成的图像合成算法进行逐点合成。

所述步骤3)中，图像进行三维重建化过程包括下述步骤：

①数据预处理优化：使用包围盒技术对三维空间数据场进行裁剪，在平行光线的投影下，利用规则网格数据场特点，使用Graham方法实现用一个能把三维体数据场中的有用信息都包围起来的最小长方体包围盒，把目标体数据限定在其中；把包围盒投影到成像屏幕上，形成一个投影多边形；对所得到的最优投影多边形进行一次扫描转换实现投射光线的生成；

②投射光线体素化及重采样优化：通过图像空间坐标计算转换为物体空间坐标，接着在物体空间坐标系OXYZ中利用Bresenham方法推广到三维，对投射光线进行体素化和求交处理。

所述步骤①中，用Graham求解最小包围盒凸壳的算法过程如下：

a)先计算平面凸集中所有点的y坐标，记录y值最小点，之后把它与凸集中所有剩余点进行连接形成线段，接着计算所得到的线段与水平线之间的夹角；然后对这些点先按照夹角大小进行排序，如果夹角相同，就按照与y值最小点的距离排序，依次连接得到的序列点集p1，p2，…，pn，从而形成一个多边形；

b)通过判断某两个点是不是在某一条指定线段的同侧，把p3到pn-1中不是凸壳顶点的点删除；

c)最后把经过删除所得到的凸壳顶点进行输出。

所述步骤②中，投射光线体素化及重采样优化，坐标转换是通过下述迭代关系式实现的：

p2′.x＝p2.x+r₀；p2′.y＝p2.y+r₃；p2′.z＝p2.z+r₆

其中，(p2.x，p2.y，p2.z)是点p1在物体空间坐标系中对应的点p1′的坐标；r₀、r₃、r₆分别表示像素在x，y，z轴上的物体空间位置。

与现有技术相比，本发明的积极效果是：本发明所生成的混凝土CT三维图像提供立体的混凝土骨料、砂浆、孔洞和裂纹区域结构，为混凝土细观破裂研究提供分析研究观测平台。本发明方法实时性强、运算效率高，能清晰的分辨骨料、砂浆及孔洞。

本发明具有下述特点：

1)混凝土图像重构方法是在光线投射算法基础上提出了一种加速的光线投射算法，极大的提高了传统光线投射法图像绘制的速度，对图像的绘制质量几乎没有影响。并且，由于大部分的加速工作都是在预处理阶段完成的，从而使图像的交互速度也有了很大的提高。

2)混凝土三维重建结果得到了混凝土真实的细观结构，可以与混凝土细观数值模拟结果进行试验比较，校正数值分析模型。

3)三维重建技术为混凝土CT图像立体细观研究提供了一定的观测技术平台。

4)混凝土真实细观结构的三维重建结果可以作为导入到大型工程软件(如ANSYS)的前期技术基础，从而可以建立起能更加准确表征混凝土非均质特性的三维有限元模型，为土木工程学者进一步研究混凝土材料的力学特性和细观破坏机理提供更好的支持。

附图说明

下面结合附图对本发明作进一步的详细说明。

图1为本发明方法步骤流程图。

图2(a)-图2(e)为本发明Graham方法求解投影多边形凸壳过程；

图3为本发明坐标转换成像平面相邻点；

图4为本发明排序两交点；

图5为本发明光栅化直线方法Bresenham方法；

图6为本发明体素化投射光线。

具体实施方式

本发明的混凝土图像重构方法，如图1所示，该方法包括：

第一步，采用医用CT或工业CT扫描混凝土试件断面，获得混凝土静力压缩CT图像二维断层图像数据；

第二步，运用三维重构可视化技术对混凝土静力压缩CT二维断层图像进行立体重构，生成三维显示图像；具体方法如下：

运用三维重构可视化技术对混凝土CT二维断层图像进行立体重构，采用光线投影算法按照设定视线的方向，从混凝土静力压缩CT二维屏幕上的每个像素点出发，发出一条射线，沿着该射线按照一定的步长进行等距采样；对距离某个采样点最近的八个体素的属性值做三线性插值运算从而计算出，得到该采样点的属性值(颜色值和不透明度值)，对得到的每个采样点的属性值按照从后到前或从前到后的顺序进行逐点合成，从而得到所对应的屏幕上该像素点的颜色值和不透明度值，最终混凝土CT生成三维显示图像。

第三步，对混凝土CT三维显示图像进行三维优化，从而完成混凝土CT图像三维重构。

①数据预处理优化。

A.使用包围盒技术对三维空间数据场进行必要的裁剪，这样可以避免很多不必要的数据计算，从而提供绘制速度。

运用包围盒技术对其进行裁剪的过程为：在平行光线的投影下，利用规则网格数据场的一些特点，用一个能把三维体数据场中的有用信息都包围起来的最小长方体包围盒，把目标体数据限定在其中。这样就可以利用包围盒对数据场进行有效的裁剪，从而极大的缩小了数据场的规模。我们假设OUVW是图像空间坐标系，OXYZ是物体空间坐标系。假如用Δx，Δy，Δz表示三维数据场中X，Y，Z三个坐标轴上每个体素的步长间隔，对体素进行遍历就可以确定对应的数据场的最小包围盒的长、宽、高，它们分别是：

(max_x-min_x)*ΔX，(max_y-min_y)*Δy，(max_z-min_z)*Δz，

如果体素在包围盒内部，就可以对其进行采样；如果它在包围盒外部就忽略，即认为该体素为空体元。

要对体数据进行显示必须完成体视见变换，即把体数据从物体空间OXYZ转换到图像空间OUVW。但是有些体素的体视见变换是不必要的，为了提高算法速度，必须减少这些不必要的变换。把包围盒投影到成像屏幕上，形成一个投影多边形，只从投影多边形中的像素发出光线进行投射，就可以保证所发射的每条光线都能与三维空间的体数据相交。这种方法不仅能减少了一些不必要的体视见变换，同时也通过减少所包含的像素点使发出的投射光线数目与求交运算也大幅度的下降，并且还能有效的保证数据场中非空体元的遗漏。这一过程使用Graham方法实现。

如图2(a)-图2(e)所示，用Graham求解最小包围盒凸壳的算法过程如下：

(1)先计算平面凸集中所有点的y坐标，把y值最小的点记作p1，之后把它与凸集中所有剩余点进行连接形成线段，接着计算所得到的线段与水平线之间的夹角。然后对这些点进行词典式排序，即先按照夹角大小进行排序，如果夹角相同就按照与p1的距离排序，依次连接得到的序列点集p1，p2，…，pn，从而形成一个多边形。因为p1是凸壳边界的起点，且p1p2和p1pn与水平线的夹角分别最小和最大，所以p2和pn也肯定是凸壳的顶点。

(2)通过判断某两个点是不是在某一条指定线段的同侧，把p3到pn-1中不是凸壳顶点的点进行删除。

(3)最后把经过删除所得到的凸壳顶点进行输出。

B.对所得到的最优投影多边形进行一次扫描转换实现投射光线的生成。进行扫描转换的目的就是用一些常用的转换算法按照显示器要求的模式进行图像的精确组织，从而可以用视频技术来显示所得到的光栅图像，采用活化边表和排序Y桶的有序边表算法来实现投影多边形的扫描转换，从而记录在图像空间OUVW中扫描到的有效像素点位置。其具体算法如下：

(1)数据准备：

采用位于y+1/2的中心扫描线，对于多边形的每条边求出与其相交的最高扫描线。把该多边形的边存入与这一扫描线相对应的y桶。把x交点的初值，多边形边穿过的扫描线条数Δy以及相邻扫描线之间的x增量Δx存入链表。

(2)数据的扫描转换：

对于每条扫描线，检查相应的y桶是否有新的边，把新的边移入活化边表。在活化边表中对应交点安x递增顺序排序，即当x₁≤x₂时，x₁位于x₂之前。从已x排序的表中成对地取出交点。在扫描线y上激活x整数值以满足

的像素。活化边表中每条边的Δy减1。若Δy＜0，从活化边表中删除这条边。计算新的x截距交点x新＝x旧+Δx。对下一条扫描线重复上述过程。

②投射光线体素化及重采样优化。

这一阶段的优化是与上一节中投影多边形的扫描转换相结合一起完成的。通过投影多边形的扫描转换我们记录了在图像空间OUVW中扫描得到的有效像素点的位置，首先需要把这些坐标进行一下转换运算，即通过图像空间坐标计算转换为物体空间坐标，接着在物体空间坐标系OXYZ中利用Bresenham方法推广到三维，对投射光线进行体素化和求交处理。

A.坐标转换。用R表示旋转矩阵，T表示平移矩阵，则可以假设从图像空

间坐标系到物体空间坐标系的转换矩阵是M＝R*T。

$R=\left(\begin{array}{cccc}{r}_0& r_1& r_2& 0\\ r_3& r_4& r_4& 0\\ r_6& r_7& r_8& 0\\ 0& 0& 0& 1\end{array}\right)$ $T=\left(\begin{array}{cccc}1& 0& 0& t_1\\ 0& 1& 0& t_2\\ 0& 0& 1& t_3\\ 0& 0& 0& 1\end{array}\right)$

在坐标转换过程中可以利用坐标之间的相关性进行优化，由它们之间的相关性得到迭代关系式，使坐标转换之间的运算量减少，从而极大的提高了图像的绘制速度。

在图3中可知，p1的坐标为(u1，v1，w1)，p2的坐标为(u2，v2，w2)，它们是相邻的两个像素点，都位于图像空间坐标系中的同一扫描线上，由于p1和p2之间存在相关性，所以p2的坐标也可以表示为(p1.u+1，p1.v，p1.w)。这样就可以运用下面的迭代关系式快速的把p2的坐标转换为物体空间坐标系上的坐标：

p2′.x＝p2.x+r₀；p2′.y＝p2.y+r₃；p2′.z＝p2.z+r₆

其中(p2.x，p2.y，p2.z)是点p1在物体空间坐标系中对应的点p1′的坐标。r₀、r₃、r₆分别表示像素在x，y，z轴上的物体空间位置。

经过这样的迭代处理，使进行坐标转换处理时的运算量比之前普通的矩阵运算有了大幅度的减少，极大的提高了绘制速度。对应处在同一扫描线上的所有像素点全都采用这样的迭代关系式来进行坐标转换的计算。

按照视线方向从点P2发射一条光线，利用投射光线与平面族求交代替直接与体元求交算法，同时结合上一阶段找到的体数据场的包围盒就可以非常快速的求出此投射光线与三维体数据场的起始和终止交点，有效的剪裁了投射光线，这样就可以避免后面阶段的盲目采样，从而减少了无用采样点的数目。在求交的过程中，我们只需要考虑光线与长方体包围盒六个表面M1，M2，……，M6的相交情况。

因为入射光线方向是任意的，它与长方体包围盒相交的入点和出点有以下三种情况：1.交点落在长方体包表面的内点处；2.交点落在长方体棱边的内点处；3.交点落在长方体顶点处。如果求的的两个交点重合，那么该光线对应的屏幕上的像素点颜色值设为背景颜色值；否则需要利用投射光线的方向对这两个交点进行排序，从而确定出入点。在此我们利用判定两向量同向的方法对这两个交点进行排序。

在图4中，我们假设M(x_M，y_M，z_M)、N(x_N，y_N，z_N)是当前入射光线与长方体包围盒Box的两个交点，其中

表示投射光线的方向，并且它不是零向量。现在对M，N进行排序，从而确定出入点。我们先不妨假设M是入点，N是出点，从入点到出点的向量

与投射光线的向量

的关系如下：即要断向量 $\stackrel{\→}{\mathrm{MN}}=(x_B-x_M,y_N-y_M,z_N-z_M)$ 是否与向量同向。本文采用符合函数sign(x)进行判断：

$\mathrm{sign}\left(x\right)=\left\{\begin{array}{ccc}1& \mathrm{if}& x>0\\ -1& \mathrm{if}& x<0\\ 0& \mathrm{if}& x=0\end{array}\right.$

从向量

的三个分量中找出一个非零向量，我们不妨设x_N-x_M≠0，如果sign(x_N-x_m)＝sign(x_L)，则可以得到M是入点，N是出点；否则M是出点，N是入点。

确定了投射光线与长方体包围盒的起止交点后，对入点到出点之间的线段做体素化处理，即把物体的几何表示形式转换为对应的体素表示形式。如图5中，本发明利用最有效的光栅化直线方法Bresenham方法，将它推广至三维进行投射光线的体素化。

Bresenham算法循环的步长沿着向量

的三个分量中最大分量绝对值的方向有入点M向出点N前进，设

max(|x_N-x_M|，|y_N-y_M|，|z_N-z_M|)＝|x_N-x_M|

即沿着x分量的方向前进，其中前进的步长为sign(x_N-x_M)×Δx；同样可以利用类似于二维Bresenham光栅化直线的方法判断出其他两个分量方向的前进情况，y、z方向上的前进步长分别是sign(y_N-y_M)×Δy、sign(z_N-z_M)×Δz，其中Δx、Δy、Δz分别是相邻两个体素在三个坐标轴上的距离。如图6所示。

投射光线体素化过程完成之后，就可以对该光线上的体素进行重采样。体绘制是通过对体数据的不断采样和合成来绘制图像的，因此在重采样阶段进行一些优化处理，可以有效的提高图像的绘制速度。在重采样阶段有很多优化方法，本发明采用自适应采样方法。两个相邻体素的灰度或梯度方向往往相同或相近，利用这种相邻体素间的相关性，可以在物体空间中避免对具有相同性质的体素进行不必要的重复采样，而利用体数据本身的性质进行自适应采样。采样步长的大小可以依据采样点周围体数据属性值的均匀程度来确定。如果周围体数据比较均匀，采样步长我们就取大一点；采样点周围的体数据不够均匀，那么采样步长就取小一点，从而有效的避免了对属性值相同的体素的重复采样。

采样步长Δt的大小是根据采样点周围数据场是否均匀自适应确定的。而数据场是否均匀主要是用体素的特征值表示的，本发明主要考虑体素的颜色值和不透明度值。Δt的大小可以使用投射光线方向上体素个数加以度量。由于视线的方向是任意的，则穿过体素的投射光线的方向也是任意的，所以，用一种保守的度量方式即正方体表示法来反应采用点周围数据场的均匀性范围。当前采样点位于该正方形的中心，该采用点与下一个采用点之间的距离也就是采样步长为Δt，正方体的边长是Δt+1，则可以用下式来确定正方体区域的均匀性：

F(x+i，y+j，z+k)＝F(x，y，z)， $\&ForAll;i,j,k\&Element;[-\mathrm{\&Delta;t},\mathrm{\&Delta;t}]$

其中F表示的是体素的属性值，即特征值。经过计算，满足上式的Δt就是我们所求的该点的采样步长。通过对每个体元进行周围数据场均匀性的判断，从而确定其采样步长Δt。一般，周围体数据比较均匀时，得到的采样步长就比较大；采样点周围的体数据不够均匀时，得到的采样步长比较小，这样就避免了对属性值相同的体素的重复采样，最后对求的此采用点进行三线性插值运算计算出该采样点的属性值(颜色值和不透明度值)，从而有效的提高了图像的绘制速度。

第三步，图像合成优化。

图像合成的目的是按照设定视线的方向，从屏幕上的每个像素点出发，发出一条射线，沿着这条射线按照一定的步长进行等距采样，并计算每个采样点的属性值，运用一定的规则对各个采样点的属性值进行合成，从而得到所对应的屏幕上该像素点的颜色值和不透明度值，当屏幕上每个像素点的颜色值和不透明度值都得到之后，最终生成显示图像。在像素体绘制算法中，有两种不同的图像合成算法：由后向前的图像合成和由前向后的图像合成。本发明采用由前向后合成法进行图像的合成，利用由前向后合成法进行图像合成时，不透明度值O肯定是慢慢增大的，当它增大到接近于1时，光线后面的体素对该像素点的图像将不会有任何贡献，所以就可以结束计算。由前向后合成法可以剔除掉很多无用的计算，有效的提高了绘制速度。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (6)

1.混凝土CT图像三维重构方法，其特征在于，该方法包括下述步骤：

第一步，采用医用CT或工业CT扫描混凝土试件断面，获得混凝土静力压缩CT二维断层图像数据；

第二步，运用三维重构可视化技术对混凝土静力压缩CT二维断层图像进行立体重构，生成混凝土CT三维显示图像；

第三步，对混凝土CT三维显示图像进行三维优化，最终完成混凝土CT图像三维重构。

2.根据权利要求1所述的混凝土CT图像三维重构方法，其特征在于，所述第二步骤中，运用三维重构可视化技术对混凝土CT二维断层图像进行立体重构，采用光线投影算法按照设定视线的方向，从混凝土静力压缩CT二维屏幕上的每个像素点出发，发出一条射线，沿着该射线按照一定的步长进行等距采样；对距离某个采样点最近的八个体素的属性值做三线性插值运算，得到该采样点的属性值，对得到的每个采样点的属性值进行逐点合成，从而得到所对应的屏幕上该像素点的颜色值和不透明度值，最终生成混凝土CT三维显示图像。

3.根据权利要求2所述的一种混凝土CT图像三维重构方法，其特征在于，所述对得到的每个采样点的属性值按照由后向前的图像合成和由前向后的图像合成的图像合成算法进行逐点合成。

4.根据权利要求1所述的混凝土CT图像三维重构方法，其特征在于，所述第三步骤中，对混凝土CT三维显示图像进行三维优化过程包括下述步骤：

①数据预处理优化：使用包围盒技术对三维空间数据场进行裁剪，在平行光线的投影下，利用规则网格数据场特点，使用Graham方法实现用一个能把三维体数据场中的有用信息都包围起来的最小长方体包围盒，把目标体数据限定在其中；把包围盒投影到成像屏幕上，形成一个投影多边形；对所得到的最优投影多边形进行一次扫描转换实现投射光线的生成；

②投射光线体素化及重采样优化：通过图像空间坐标计算转换为物体空间坐标，接着在物体空间坐标系OXYZ中利用Bresenham方法推广到三维，对投射光线进行体素化和求交处理。

5.根据权利要求4所述的一种混凝土CT图像三维重构方法，其特征在于，所述步骤①中，用Graham求解最小包围盒凸壳的算法过程如下：

a)先计算平面凸集中所有点的y坐标，记录y值最小点，之后把它与凸集中所有剩余点进行连接形成线段，接着计算所得到的线段与水平线之间的夹角；然后对这些点先按照夹角大小进行排序，如果夹角相同，就按照与y值最小点的距离排序，依次连接得到的序列点集p1，p2，…，pn，从而形成一个多边形；

b)通过判断某两个点是不是在某一条指定线段的同侧，把p3到pn-1中不是凸壳顶点的点删除；

c)最后把经过删除所得到的凸壳顶点进行输出。

6.根据权利要求4所述的一种混凝土CT图像三维重构方法，其特征在于，所述步骤②中，投射光线体素化及重采样优化，坐标转换是通过下述迭代关系式实现的：

p2′.x＝p2.x+r₀；p2′.y＝p2.y+r₃；p2′.z＝p2.z+r₆

其中，(p2.x，p2.y，p2.z)是点p1在物体空间坐标系中对应的点p1′的坐标；r₀、r₃、r₆分别表示像素在x，y，z轴上的物体空间坐标位置。

Images