常见三维重建算法及其实现（含源码下载链接）

主要整理部分本科毕业设计

这里三维重建输入的是断层图像。

对图像首先需要进行一些常见操作，直方图均衡化、常见滤波、图像锐化、边缘提取、二值化等等，把常见操作集成在一个MATLAB GUI中。当初也是从网上下的demo改的，这里就不要积分了。

网址：http://download.csdn.net/detail/bizer_csdn/9644508

这里主要注意一下图像二值化操作（后面MC算法需要输入二值化图像）

请看这篇博客：http://blog.csdn.net/bizer_csdn/article/details/52717940

代码编写核心就是就是用一个比较大高斯卷积核与图像卷积后作为二值化的阈值。

下面开始讲讲三维重建把，我这里主要针对医学断层图像。

三维重建的方法大概有两种：直接将体素投影到计算机显示平面的方法被称为体绘制（VolumeRendering），也可以被称为直接绘制；而通过几何单元（一般近似选取三角面片）来拟合物体的表面，这种方法被称为面绘制（SurfaceFiting），也可以称为间接绘制。其中，体绘制是直接分析光线穿过三维体数据场时的变化情况，并且得到最终绘制结果；而面绘制方法主要利用计算机图形学技术，借助光照模型，把物体的表面绘制在计算机屏幕上。

MATLAB自带isosurface其实就是一种等值面绘制技术，具体操作请看：

https://cn.mathworks.com/help/matlab/visualize/techniques-for-visualizing-scalar-volume-data.html

下面我们仔细讲一讲：

首先，体绘制，大概有以下几种：

光线投射算法（Ray casting）是以图像空间为序的体绘制算法，并被广泛地运用。在光线投射过程中，该算法需要累加计算每一条穿过体数据的光线的不透明度和颜色的属性值，并且投射结束的条件设定为光线不透明或者光线已经穿过体素空间，最后合成当前的不透明度值和颜色值，并把其作为该像素的属性写入帧缓冲区。

光线投射算法大致有三步，首先，输入原始数据，并可以对其进行一些预处理操作；其次，是计算光强的过程，即根据一定规则对不同的数据赋予不同的不透明度值和颜色值，在VTK编程中，有专门的成员函数负责分段线性化工作，主要依靠vtkSmartPointer<vtkPiecewiseFunction>和vtkSmartPointer<vtkColorTransferFunction>来实现；最后，就是图像与合成显示工作。

纹理映射算法（Texture Mapping），也称为纹理贴图，由于纹理映射算法中的纹理生成步骤都是在专有的图形硬件中生成，所以它的绘制速度也比光线投射算法快，但质量缺没光线投射算法高。所谓的纹理映射，是指依靠计算机图形学技术描述真实物体表面细节（纹理可以理解为物体表面细节）。二维纹理映射和三维纹理映射是纹理映射算法的典型代表算法。

人们把物体纹理放入一个称为纹理空间的坐标系中（纹理空间可以理解为一个二维图像），纹理映射，就是建立从物体表面上的点到纹理空间中的点的映射关系，这个过程就像“贴墙纸”。最后，根据映射关系，纹理图案就可以贴到物体表面上。

三维重建的输入数据可以是二维图像序列，也可以是三维体数据。纹理映射首先要根据输入数据生成纹理，然后才是纹理绘制过程。如果用二维图像序列作为输入，按着不同图形硬件要求，对图像进行裁剪或者插值；然后通过转换函数，分别建立灰度对不透明度和颜色的对应关系，使原来的图像序列变为二维纹理切片序列，构成纹理空间。如果三维体数据作为输入，仍按着照图形硬件要求对体素进行裁剪或者插值，此后，再按某个正交方向对体数据切片，生成二维图像序列，下面步骤就完全与二维图像输入的情况一样，通过转换函数生成二维纹理序列。在VTK编程中，也是依靠vtkSmartPointer<vtkPiecewiseFunction>和vtkSmartPointer<vtkColorTransferFunction>来实现灰度值对不透明度和颜色的转换。

对于安装VTK的朋友，VTK的官网有demo，这里不贴出来了。

好了，下面开始重头戏，面绘制。

面绘制主要通过平面元来逼近物体表面，目前面绘制主要有切片级重建和体素级重建两种，划分依据就是绘制过程中处理元素的级别不同。由于，切片级重建重建方法有一些技术难点尚待解决，本文只介绍介绍体素级重建方法。

体素级重建方法是体数据内的体素（这里“体素”定义，是三维数据场中相邻两层中相邻的八个顶点组成的立方体，与前文体绘制中“体素”定义完全不同）为单位进行表面跟踪，小面片重建是在体素内完成的。最后，借助于计算机图形学，通过光照模型，明暗处理后，显示出小面片组成的物体。

表面重建最经典的就是Marching cubes（MC）算法了

这篇文章是对算法的编程进行了介绍（opengl），很经典,一定要看（ps:后面MC改进我是在VTK平台下实现,VTK实现与该文章大同小异）

http://paulbourke.net/geometry/polygonise/

下面接单介绍：

MC算法，通过等值面提取（IsosurfaceExtraction）重建物体表面，其本质就是从一个三维数据场中抽取一个等值面。等值面定义为三维空间中具有相同属性的点的集合，其数学描述为：

   

此外，还需要了解下三线性插值，这里就不介绍了。

标准MC算法的步骤如下：

（1）确定体素中等值面的剖分方式

首先需要把体素的八个顶点分为两类，如果顶点数据大于或等于等值面的标量值，记该点的状态为“1”；如果顶点数据小于等值面的标量值，记为该点的状态为“0”。 由于MC算法是基于沿着立方体的边数据场呈连续性变化这一假设上而提出的，也就是说立方体边数据场是线性变化的。所以，如果体素中一条边上的两个顶点分别大于、小于等值面的值（即该棱边上两个顶点状态分别为“1”和“0”），则在该边上必有也仅有一点是这条边与等值面的交点，然后可以通过线性插值方法计算该边上等值面的顶点和该顶点的方向量。

因为每个体素有八个顶点，而每个顶点又有“1”和“0”两种状态，所以一共有2的8次方种可能，即256种组合。MC算法把256种组合化简了最基本的15种情况。首先，把体素顶点状态置反，即“1”变“0”， “0”变“1”，而对体素内三角面片的剖分方式并没有影响，这一性质也被称为互补对称性，这样256种组合就可以化简为128种。然后，对于上述情况某种确定体素而言，如果这经过旋转后，这个体素与另一情况的体素是一致的，即表现为顶点位置及其状态相同，那么这两种体素情况就又可以合并为同一种，这一性质也被称为旋转对称性，于是128种情况化简为最终的15种情况。

下图所示为MC算法的15种基本情况，其中黑点表示标记为1的角点，标记为0的角点在图中立方体上为未作记号的顶点。

在VTK编程处理中，定义了一个构型索引（index），index中存储了体素中八个角点的所有状态（“1”变“0”），如图所示：

 

为了方便求出每个立方体的索引，实际中让索引从0开始，并根据立方体每个顶点的标量值的状态，决定索引与1,2,4,8,16,32,64,128中依次哪几个（这些数字称之为掩码）进行“或”运算得到。具体核心代码如下所示：

for ( ii=0, index = 0; ii < 8; ii++)
{
   if ( s[ii] >= value )
   {
     index |= CASE_MASK[ii];  // static int CASE_MASK[8] = {1,2,4,8,16,32,64,128};
   }
}

另外VTK通过构造“三角剖分”查找表中确定三角剖分形式，图3·13举例了查找表其中一种构型的构造过程，左为立方体顶点和其边的编号，右为基本构型中的其中一种。

对于上图右边这一基本构型，我们如下所示的表来表示其三角剖分形式：

其中数字3,11,2分别表示第3号边、第11号边和第2号边上有三角面片的顶点，而-1表示缺省。由于每个立方体最多可能产生5个三角面片，所以这里一共需要16位数字来表示。

（2）求取等值面与体素边界的交点

当体素内三角面片剖分模式确定后，三角面片顶点位置是通过线性插值方法计算得到的，这是因为MC算法假设立方体边数据场呈线性变化。对于体素某一棱边，设其两个端点（p1、p2）的标量数据为v1、v2，那么插值得到的交点（即三角面片的顶点）为p：

（3）计算等值面顶点处的法向量

在计算机图形学中，为了生成具有真实感三维图形，就需要选择合适的局部光照模型进行光照计算，所以，MC算法除了要给出三角面片各顶点的坐标外，还是给出该顶点处的法向量。

在计算出三角面片各顶点处坐标后，虽然可以通过叉积的方法计算该顶点处的法向量，但这种做法计算量较大。MC注意到，对于等值面上任意一点而言，该点处等值面的切线方向在等值面上无变化率可言，而梯度方向是方向导数变化率最大的方向，因此，该点的梯度方向就是该点在等值面上法向量的方向。前文已经讲述MC算法假设立方体边数据场呈线性变化，为了方便求出三角面片顶点处的法向量，可以先分别求出该三角面片顶点所在体素棱边上的两个角点的梯度，然后再次通过线性插值，计算该三角面片顶点的法向量。所谓的中心差分，就是计算当前角点在某个方向上前向和后向两个角点的标量的差分，体素角点的中心差分公式如下：

此时，法向量为：

设立方体顶点处的法向量为N1、N2（标量数据仍为V1、V2），再次利用线性插值得到的三角面片顶点处法向量为：

以上分析了MC算法的基本理论，下面逐步说明MC算法的编程步骤：

     1) 分层读入三维离散数据场的数据；

     2) 逐层扫描三维数据场，把相邻两层的八个顶点组成的立方体作为当前体素；

     3) 分别比较当前体素中八个角点的标量值与等值面的值，计算出当前体素的索引（index）；

     4) 排除index为0或者255的体素；

     5) 利用线性插值方法求出三角面片的各顶点坐标；

     6) 利用中心差分法，先计算体素各角点处的法向量，然后再次利用线性插值，

计算三角面片各顶点处的法向量；

     7) 将三角面片各顶点坐标及相应的法向量送入表面绘制环节。

MC算法的伪代码描述如下：

 for(k=1;k<Nz;k++)

{

       读入数据点值

        for(j=1;j<Ny,j++)

        {

           读入数据点值

          for(i=1;j<Nx;i++)

          {

             a．读入数据点值

             b. 由当前体素的八个角点(i,j,k)，(i+1,j,k)，(i,j+1,k)，(i+1,j+1,k)，(i,j,k+1)，(i+1,j,k+1)，(i,j+1,k+1)和(i+1,j+1,k+1),判定每个角点上标量值与等值面的值相对大小，并由此                  计算出当前体素的索引下标index，并且排除index为0或255的体素

             c. 通过当前index对应“三角剖分”插值表，确定当前三角剖分形式

             d. 由当前三角剖分形式，可知当前体素哪条棱边（以及确定该棱边对应的哪两个顶点）会与三角面片相交

             e. 由线性插值计算出体素棱边上三角面片顶点的坐标

             f. 由中心差分法计算体素各角点处的法向量

            g. 再次利用线性插值计算三角面片各顶点处的法向量

            h. 把三角顶点和法向量送入表面绘制步骤

           }

       }

}

实现就是这样，如果使用VTK的话，参考一下vtkMarchingCubes这个类

VTK实现MC源码是这几个文件：vtkMarchingCubes.h,vtkMarchingCubes.cxx,vtkMarchingCubesCases.h,vtkMarchingCubesCases.cxx

下面看一下结果，可以发现MC去噪声能力有待提高，并且MC算法还存在二义性问题。

下面对MC方法的改进：

算法基于文章 "何晖光， 田捷， 赵明昌， 杨骅． 基于分割的三维医学图像表面重建算法[J]． 软件学报， 2002， 02： 219-226"

对于图像二值化算法，采用的是前文自适应二值化方法，在此处就不再赘述。另外采用二值化数据重建表面，可以不用考虑等值面的标量值，这一点也是有优势的。

标准算法是遍历所有体素，实际中有些体素并不需要遍历，如上文MC的15种基本构型中的第0号构型，是可以直接忽略的。文献采用的算法把标准MC算法遍历所有体素改为一种基于区域增长的立方体检测方法，并不需要遍历所有体素，这能够加快重建速度。

除了如前所述，第0号构型不需要考虑外，对于连通区域，如果一个立方体中的三角剖分的形式确定后，那么与它相邻的立方体中的三角面片只可能按着按着三维空间内，前、后、右、左、上、下这六个方向其中的一个或几个方向去延伸。比如，对于图3.11中15种基本构型中的第1号构型，它延伸的方向可能会向前、向左和向下，其他方位相邻的立方体可以不用考虑；

基于以上分析，当立方体确定后，与之相应的是该立方体的延伸方向也随之确定。因此在标准MC算法的索引（index）和“三角剖分”查找表的基础上，在构建一个立方体邻接表，来确定当前立方体的邻接方向。于此同时，算法可以用一个标记Flag来表示立方体是否被处理过（1代表处理过，0代表未处理过），并在程序初始化时候，让立方体的标记Flag置零。

MC采用线性插值的方法计算三角面片的顶点坐标和当前顶点的法向量，文献采用的算法直接默认三角面片的顶点在立方体的边的中点，与之相应，法向量计算也改为对立方体边的法向量直接求平均，这样做可以进一步减少计算量。

所以，算法在标准MC算法基础上，在添加一些数据结构，其基本数据结构包括标准MC算法中的“三角剖分”查找表triCases、一个立方体邻接表neiborCases、一个表示立方体的结构体Cube、一个队列CubeQueue和一个由若干个该结构体组成的向量List。“三角剖分”查找表的构造过程已经在第3.1.1节中阐释，立方体邻接表与其相似，按着前、后、右、左、上、下这六个访问构造，其中1表示此处方向邻接，0表示此处方向不邻接，具体过程这里不再赘述。而对于立方体，本文用图3.13中第0号顶点表示当前立方体坐标，其结构体如下所示：

struct Cube{
	        int oi;
	        int oj;
	        int ok;
	        int index;
	        bool Flag;
        Cube(int i=0,int j=0,int k=0,int id=0,bool f=0):oi(i),oj(j),ok(k),index(id),Flag(f){}
};

其中oi,oj,ok分别表示前文中第0号顶点的坐标，index即是前文所示的立方体的构型索引，Flag是立方体的访问标记，并通过构造函数，将所有值初始化为零。

程序是在VTK的基础上实现的，队列和向量的实现可以通过C++中STL，很方便的实现。

算法开始时，在遍历体素过程中一旦发现构型索引不为0或255（即立方体index不为0或255）直接停止遍历，把当前立方体作为种子点压入队列CubeQueue中，并置当前立方体的标记Flag为1，随后通过当前立方体的构型索引(index)和“三角剖分”查找表确定当前立方体的三角剖分形式，并按前文所述方法计算各个三角面片的顶点以及各个顶点处的法向量；然后，通过邻接表，把该立方体相邻的立方体压入到队列中。随后循环重复上述步骤，直至队列为空时候停止。

算法伪代码描述如下：

输入：分割后的数据（用于表面追踪），原始三维数据（用于计算法向量）

输出：三维模型

辅助数据结构：struct Cube；std::queue<Cube> CubeQueue；

std::vector<Cube>List；

Begin：

（1）       初始化List；初始化CubeQueue，置其为空；

（2）       然后检测数据集，选取一个构型索引index不为0或255的立方体作为种子点，并将其压入队列CubeQueue中；

while(!CubeQueue.empty()){

（3）       从CubeQueue中取出队首元素，得到前文图中第0号顶点的坐标（i，j，k），并令List中该立方体的编号idx=i+jOffset+kOffset；

（4）       if(List[idx].Flag==1)  continue；

（5）       置List[idx].Flag=1；

（6）       根据当前立方体的八个顶点状态得到其构型索引index；

（7）       由构型索引index查找“三角剖分”查找表triCases得到当前立方体的三角剖分形式，同时采用本节前文所述的方法计算各三角面片的顶点坐标和顶点处的法向量；并将其输出备用；

（8）       通过构型索引，查找立方体邻接表neiborCases将邻接的立方体压入到队列中去；

}

End

再看一下结果，因为有二值图像输入，噪声少很多了：

代码我放在了github上，欢迎下载指正（VTK6.2）：

https://github.com/bizerfr/new-marching-cubes

另外，还有把一些常见三维重建算法与上文算法集成在一个简单VTK,QT程序中：

https://github.com/bizerfr/Simple_Qt_design_3D_reconstruction