断层数据三维重建就是基于一些列二维断层图像还原出被检物体的三维结构。其基本思想就是由一些列二维断层图像构成的数据集,再由此数据集形成三维空间采用数据集,进而采用OpenGL对三维空间数据集进行渲染,以还原出三维实物。采用的算法是MC 算法,它是在三维数据场中提取等值面,由等值面来反映物体原貌。这种方法基础,但是有待优化。以下在Visual Studio2010下给出一个实例,该实例代表了一种三维重建的思想,具有这一类问题的代表性。
第一步:新建一个单文档的MFC工程,命名为SurfaceReconstruction;
第二步:设置工程属性:
1.在SurfaceReconstructionView.h的头文件中天际以下两个头文件:
#include "GL/gl.h"
#include "GL/glu.h"
2.打开属性设置对话框:
1).设置程序为静态运行方式,配置属性->常规->MFC使用->在静态库中使用MFC,然后单击确定;
2).连接OpenGL库文件路径,配置属性->链接器->常规->附加库路径,输入OpenGL的库文件所在的路径,单击确定;
3).连接OpenGL库文件,再选择配置属性->链接器->输入->附加依赖库,在弹出的对话框中输入OpenGL的库文件OpenGL32.lib和glu32.lib
4).编译运行,检查配置是否正确,若有错误,则返回检查,重新配置。
第三步:添加视图类成员和成员函数
1.在SurfaceReconstructionView.h的开头定义两个结构体,结构体如下:
typedef struct{
double x, y, z;
} POINT3D; //定义点的三维坐标
typedef struct{
POINT3D p[3]; //三角片三个顶点坐标
POINT3D n[3]; //三角片三个顶点法向量
} TRIANGLE;
2.在SurfaceReconstructionView.cpp文件的开始定义表示边索引表和三角片索引表的两个数组:
//边索引表
int EdgeTable[256] = {
0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 };
//三角片索引表
int TriTable[256][16] = {
{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
{3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
{9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
{10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
{2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
{11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
{11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
{9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
{6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
{6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
{3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
{10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
{10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
{0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
{3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
{10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
{7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
{0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
{7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
{7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
{10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
{7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
{6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
{8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
{1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
{10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
{10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
{9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
{7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
{6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
{6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
{9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
{1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
{0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
{5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
{11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
{2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
{1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
{9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
{5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
{9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
{11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
{1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
{4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
{0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
{1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
};
3.为SurfaceReconstructionView类添加变量:
CClientDC* pDrawDC; // 用于指向当前DC的指针,用作与RC建立关联
BOOL bRotate; // 控制图形是否旋转
float RotateAngle; // 控制选择的角度
int m_nTriangles; // 记录构成等值面的三角片数目
BOOL m_bMCSurface; // 值为TRUE时,表示基于MC方法重建
HANDLE hTriAngles; // 由三角片构成的等值面句柄
TRIANGLE* pTriangles; // 指向构成等值面的三角片集合的指针
<pre name="code" class="cpp"><pre name="code" class="cpp" style="color: rgb(85, 85, 85); font-size: 15px; line-height: 35px;">BOOL m_bValidSurface; <span style="font-family: 'microsoft yahei';">// 是否抽取出三维体表面数据</span>
4.为SurfaceReconstructionView类添加成员函数
BOOL PixelformatSetting(void); // 用于设置像素格式
void GLsetting(void); // 创建渲染描述表
// 基于MC方法提取等值面
int MarchingCubes(TRIANGLE* pTriangles, LPBYTE pDataField, int nSliceWidth, int nSliceHeight, int nSliceNum, double CubeLengthX, double CubeLengthY, double CubeLengthZ,BYTE Threshold);
// 计算等值点的坐标和法向量
POINT3D CalculateIntersection(int x, int y, int z, double* gx, double* gy, double* gz, POINT3D* pNorm, double CubeLengthX, double CubeLengthY, double CubeLengthZ, int nEdgeNo);
<pre name="code" class="cpp">HANDLE LoadBMP(LPCTSTR lpFileName); //载入位图文件 <pre name="code" class="cpp" style="color: rgb(85, 85, 85); font-size: 15px; line-height: 35px;">WORD BMPColorNum(LPBYTE lpbmInfo); <span style="font-family: 'microsoft yahei';">// 计算颜色表中的颜色项数</span>
<span style="font-family: 'microsoft yahei';"></span><pre name="code" class="cpp">void DrawSurface(void); //绘制三维物体表面
5.为SurfaceReconstructionView类添加5个消息响应函数:
afx_msg int OnCreate(LPCREATESTRUCT lpCreateStruct); //响应WM_CREATE消息,消息宏ON_WM_CREATE()
afx_msg void OnDestroy();//响应WM_DESTROY消息,消息宏ON_WM_DESTROY()
afx_msg void OnSize(UINT nType, int cx, int cy);//响应WM_SIZE消息,消息宏ON_WM_SIZE()
afx_msg void OnTimer(UINT_PTR nIDEvent);//响应WM_SIZE消息,消息宏ON_WM_SIZE()
afx_msg void OnLButtonDown(UINT nFlags, CPoint point);//响应WM_LBUTTONDOWN消息,消息宏ON_WM_LBUTTONDOWN()
第四步:添加工具栏按钮及其响应函数:
1.添加工具栏按钮,打开资源视图中的Toolbar文件夹,双击IDR_MAINFRAME_256资源,在弹出的工具栏编辑窗口中创建一个图标按钮,如下图所示:
2.为新建的工具按钮设置ID,右击新建的按钮打开属性对话框,在ID栏中为其设置ID值为ID_BUTTON_MC;
3.为按钮MC手动添加消息响应函数OnButtonMC(),方法如下:
1)在SurfaceReconstructionView.h的头文件中,在CSurfaceReconstructionView类中添加公有属性函数:
afx_msg void OnButtonMC(); //按钮MC的消息响应函数
2)在SurfaceReconstructionView.cpp文件中在消息宏映射(BEGIN_MESSAGE_MAP之后和END_MESSAGE_MAP之前的代码都是消息映射宏)中添加:
ON_COMMAND(ID_BUTTON_MC, &CSurfaceReconstructionView::OnButtonMC) //按钮MC的消息响应宏
3)在SurfaceReconstructionView.cpp文件中手动添加函数体:
//工具按钮MC的消息响应函数
void CSurfaceReconstructionView::OnButtonMC()
{
}
一般手动添加消息响应函数或其它类成员函数的方法类似。
第五步:完善对所有声明函数的定义:
1.HANDLE CSurfaceReconstructionView::LoadBMP(LPCTSTR lpFileName)
/********************************************************************************************
*LoadBMP(LPCTSTR lpFileName)
*功能:从文件装载位图数据到内存
*参数:LPCTSTR lpFileName,指定位图文件的文件名(路径名)
*返回值:HANDLE,位图数据缓冲区句柄,包括位图信息块和像素阵列
*********************************************************************************************/
HANDLE CSurfaceReconstructionView::LoadBMP(LPCTSTR lpFileName)
{
HANDLE hBMP; //载入数据所在缓冲区的句柄
HANDLE hFile;
BITMAPFILEHEADER bmfHeader; //文件头
UINT nNumColors; //位图颜色表的颜色数
HANDLE hBMPtmp;
LPBITMAPINFOHEADER lpbmInfo; //指向信息头的指针
DWORD dwOffBits; //像素阵列偏移量
DWORD dwRead;
if((hFile=CreateFile((LPCWSTR)lpFileName,GENERIC_READ,FILE_SHARE_READ,NULL,OPEN_EXISTING,
FILE_ATTRIBUTE_NORMAL|FILE_FLAG_SEQUENTIAL_SCAN,NULL))==INVALID_HANDLE_VALUE)
{
return NULL; //若打开文件失败,则返回NULL
}
//为位图信息头和颜色表分配初始内存,后续可根据实际需要扩大内存
hBMP=GlobalAlloc(GMEM_MOVEABLE,(DWORD)(sizeof(BITMAPINFOHEADER)+256*sizeof(RGBQUAD)));
if(!hBMP)
{
CloseHandle(hFile);
return NULL;
}
lpbmInfo=(LPBITMAPINFOHEADER)GlobalLock(hBMP); //指向信息头的指针
if(!lpbmInfo)
{
goto ErrorExit_NoUnlock;
}
//从文件读入位图文件头BITMAPFILEHEADER
if(!ReadFile(hFile,(LPBYTE)&bmfHeader,sizeof(BITMAPFILEHEADER),&dwRead,NULL))
{
goto ErrorExit;
}
if(sizeof(BITMAPFILEHEADER) !=dwRead)
{
goto ErrorExit;
}
//读入位图信息头BITMAPINFOHEADER
if(!ReadFile(hFile,(LPBYTE)lpbmInfo,sizeof(BITMAPINFOHEADER),&dwRead,NULL))
{
goto ErrorExit;
}
if(sizeof(BITMAPINFOHEADER)!=dwRead)
{
goto ErrorExit;
}
//确定颜色大小
if(!(nNumColors=(UINT)lpbmInfo->biClrUsed))
{
if(lpbmInfo->biBitCount !=24) //如果非真彩色,则根据biBitCount计算颜色项数
{
nNumColors=1<<lpbmInfo->biBitCount;
}
}
if(lpbmInfo->biClrUsed == 0)
{
lpbmInfo->biClrUsed=nNumColors;
}
//计算像素阵列占用空间,字节对齐
if(lpbmInfo->biSizeImage == 0)
{
lpbmInfo->biSizeImage=((((lpbmInfo->biWidth*(DWORD)lpbmInfo->biBitCount)+31)&~31)>>3)*lpbmInfo->biHeight;
}
//重新根据图像实际大小分配内存,用于存放信息头、颜色表和像素阵列
GlobalUnlock(hBMP);
hBMPtmp=GlobalReAlloc(hBMP,lpbmInfo->biSize+nNumColors*sizeof(RGBQUAD)+lpbmInfo->biSizeImage,0);
if(!hBMPtmp)
{
goto ErrorExit_NoUnlock;
}
else
{
hBMP=hBMPtmp;
}
lpbmInfo=(LPBITMAPINFOHEADER)GlobalLock(hBMP);
//读入颜色表
ReadFile(hFile,(LPBYTE)(lpbmInfo)+lpbmInfo->biSize,nNumColors*sizeof(RGBQUAD),&dwRead,NULL);
//计算像素阵列偏移量
dwOffBits=lpbmInfo->biSize+nNumColors*sizeof(RGBQUAD);
if(bmfHeader.bfOffBits!=0L)
{
SetFilePointer(hFile,bmfHeader.bfOffBits,NULL,FILE_BEGIN);
}
//读入图像像素阵列数据
if(ReadFile(hFile,(LPBYTE)lpbmInfo+dwOffBits,lpbmInfo->biSizeImage,&dwRead,NULL))
{
goto SuccessExit;
}
ErrorExit:
GlobalUnlock(hBMP);
ErrorExit_NoUnlock:
GlobalFree(hBMP);
CloseHandle(hFile);
return NULL;
SuccessExit:
CloseHandle(hFile);
GlobalUnlock(hBMP);
return hBMP;
}
2.WORD CSurfaceReconstructionView::BMPColorNum(LPBYTE lpbmInfo)
/********************************************************************************************
*BMPColorNum(LPBYTE lpbmInfo)
*功能:计算颜色表的颜色项数
*参数:LPBYTE lpbmInfo,位图缓冲区地址,即指向位图信息头的指针
*返回值:DOWRD,颜色表中的颜色项数
*********************************************************************************************/
WORD CSurfaceReconstructionView::BMPColorNum(LPBYTE lpbmInfo)
{
WORD wBitCount;
//根据信息头的biBitCount计算颜色表项数
wBitCount=((LPBITMAPINFOHEADER)lpbmInfo)->biBitCount;
switch(wBitCount)
{
case 1:
return 2;
case 4:
return 16;
case 8:
return 256;
default:
return 0;
}
}
3.void CSurfaceReconstructionView::OnButtonMC()
void CSurfaceReconstructionView::OnButtonMC()
{
if( !m_bValidSurface ) //尚未提取三维表面数据
{
HANDLE h2DImg; //读入二维断层图像的句柄
DWORD dwBMPSize, //BMP文件的信息头、颜色表、像素阵列大小
dwLineBytes; //每一行的字节数
char FileName[100]; //BMP文件名
LPBITMAPINFOHEADER lpbmInfo; //指向位图数据的指针,即指向位图信息头
int ImgWidth; //二维断层图像宽度
int ImgHeight; //二维断层图像的高度
int SliceSize; //二维断层图像的大小
HANDLE hDataField; //三维数据场句柄
LPBYTE pDataField; //指向三维数据场的指针
int nSliceNum=145; //切片数目,本例使用145张CT切片图像
int i;
int x,y;
//获得当前项目所在的目录
CString sPath,sPathFileName;
GetModuleFileName(NULL,sPath.GetBufferSetLength(MAX_PATH+1),MAX_PATH);
sPath.ReleaseBuffer();
int pos=sPath.ReverseFind('\\');
sPath=sPath.Left(pos);
//读入断层图像序列构成离散数据场
for(i=0;i<nSliceNum;i++)
{
//载入当前项目目录文件夹T4中的二维断层图像
sprintf_s(FileName,"\\T4\\%d.bmp",i);
sPathFileName=sPath+(CString)FileName;
h2DImg=LoadBMP(sPathFileName);
//获取指向位图信息头的指针
lpbmInfo=(LPBITMAPINFOHEADER)GlobalLock(h2DImg);
//获取图像宽、高,每行字节数等信息,从而为三维数据场分配存放空间
if(i==0)
{
ImgHeight=lpbmInfo->biHeight;
ImgWidth=lpbmInfo->biWidth;
SliceSize=ImgHeight * ImgWidth; //断层图像的大小,即像素值
dwLineBytes=(DWORD)WIDTHBYTES((lpbmInfo->biWidth) * (lpbmInfo->biBitCount)); //保证每行的字节数是4的整数倍
dwBMPSize = (DWORD)GlobalSize( h2DImg );
//为三维数据场分配存放空间
hDataField = GlobalAlloc( GHND, ImgHeight * ImgWidth * nSliceNum * sizeof(BYTE) );
pDataField = (LPBYTE)GlobalLock( hDataField ); //指向三维数据场的指针
}
//将图像数据放入pDataField所指向的数据场
for (y=0;y<ImgHeight;y++)
{
for(x=0;x<ImgWidth;x++)
{
*(pDataField + i*SliceSize + ImgWidth*y + x) =
*((LPBYTE)lpbmInfo + dwBMPSize - dwLineBytes - y*dwLineBytes + x);
}
}
GlobalUnlock( h2DImg );
GlobalFree( h2DImg );
}
//为将要抽取的等值面分配空间
hTriAngles = GlobalAlloc( GHND, ImgWidth * ImgHeight * nSliceNum * sizeof(TRIANGLE) );
pTriangles = (TRIANGLE*)GlobalLock( hTriAngles );//指向等值面三角片集合的指针
//从数据场抽取等值面,等值面阈值设为60,并获取构成等值面的三角片数目
m_nTriangles = MarchingCubes( pTriangles, pDataField, ImgWidth, ImgHeight, nSliceNum, 1, 1, 1, 60 );
GlobalUnlock( hDataField );
GlobalFree( hDataField ); //释放数据场占用的内存空间
if( m_nTriangles )
{
m_bMCSurface = TRUE;
m_bValidSurface = TRUE;
Invalidate(); //该函数将使用户界面重绘,这将调用函数OnDraw()
}
}
}
3.int CSurfaceReconstructionView::MarchingCubes(TRIANGLE* pTriangles, LPBYTE pDataField, int nSliceWidth, int nSliceHeight, int nSliceNum, double CubeLengthX, double CubeLengthY, double CubeLengthZ,BYTE Threshold)
/********************************************************************************************
*MarchingCubes(TRIANGLE* pTriangles, LPBYTE pDataField, int nSliceWidth, int nSliceHeight,
int nSliceNum, double CubeLengthX, double CubeLengthY, double CubeLengthZ,BYTE Threshold)
*功能:基于MC方法抽取等值面
*参数:TRIANGLE* pTriangle,指向存放三角片(顶点及其法向量)内存的指针
LPBYTE pDataField,指向由二维切片图像构成的离散三维数据场的指针
int nSliceWidth, 二维切片图像的宽度(以像素为单位)
int nSliceHeight, 二维切片图像的高度(以像素为单位)
int nSliceNum,切片数目
double CubeLengthX, X方向立方体长
double CubeLengthY, Y方向立方体长
double CubeLengthZ,Z方向立方体长
BYTE Threshold, 指定的等值面阈值
*返回值:int,三角片数目
*********************************************************************************************/
int CSurfaceReconstructionView::MarchingCubes(TRIANGLE* pTriangles, LPBYTE pDataField, int nSliceWidth, int nSliceHeight, int nSliceNum, double CubeLengthX, double CubeLengthY, double CubeLengthZ,BYTE Threshold)
{
int GridIntensity[8]; //一个立方体的8个顶点灰度值
POINT3D Vertices[12]; //一个立方体的12条边上的等值点坐标
POINT3D Normals[12]; //一个立方体的12条边上的等值点法向量
double gx[8], gy[8], gz[8]; //一个立方体的8个顶点在X、Y、Z方向的梯度
int nTriangles = 0; //三角片总数目,初始化为0
int nCubeNumX = nSliceWidth - 1; // X方向立方体数目
int nCubeNumY = nSliceHeight - 1; // Y方向立方体数目
int nCubeNumZ = nSliceNum - 1; // Z方向立方体数目
int TableIndex; //索引表的索引值
int nPointsPerSlice = nSliceWidth * nSliceHeight; //一层切片图像的采样点数
int x,y,z;
for(z=0;z<nCubeNumZ;z++)
{
for(y=0;y<nCubeNumY;y++)
{
for(x=0;x<nCubeNumX;x++)
{
//获取当前立方体的8个顶点灰度值
GridIntensity[0] = pDataField[z*nPointsPerSlice + y*nSliceWidth + x];
GridIntensity[1] = pDataField[z*nPointsPerSlice + y*nSliceWidth + (x+1)];
GridIntensity[2] = pDataField[(z+1)*nPointsPerSlice + y*nSliceWidth + (x+1)];
GridIntensity[3] = pDataField[(z+1)*nPointsPerSlice + y*nSliceWidth + x];
GridIntensity[4] = pDataField[z*nPointsPerSlice + (y+1)*nSliceWidth + x];
GridIntensity[5] = pDataField[z*nPointsPerSlice + (y+1)*nSliceWidth + (x+1)];
GridIntensity[6] = pDataField[(z+1)*nPointsPerSlice + (y+1)*nSliceWidth + (x+1)];
GridIntensity[7] = pDataField[(z+1)*nPointsPerSlice + (y+1)*nSliceWidth + x];
TableIndex = 0;
//据立方体顶点灰度值和给定等值面阈值的比较确定索引表的索引值
if (GridIntensity[0] < Threshold)
TableIndex |= 1;
if (GridIntensity[1] < Threshold)
TableIndex |= 2;
if (GridIntensity[2] < Threshold)
TableIndex |= 4;
if (GridIntensity[3] < Threshold)
TableIndex |= 8;
if (GridIntensity[4] < Threshold)
TableIndex |= 16;
if (GridIntensity[5] < Threshold)
TableIndex |= 32;
if (GridIntensity[6] < Threshold)
TableIndex |= 64;
if (GridIntensity[7] < Threshold)
TableIndex |= 128;
//若该立方体和等值面不相交,即该立方体内无三角片,则跳过该立方体继续处理后续立方体
if (EdgeTable[TableIndex] == 0)
{
continue;
}
//确定立方体内三角片的组成
//计算当前立方体的8个顶点的梯度
int xl, xr, yl, yr, zl, zr; // 在X,Y,Z方向,当前顶点(x,y,z)的前、后坐标值
xl = x - 1; xr = x + 2;
yl = y - 1; yr = y + 2;
zl = z - 1; zr = z + 2;
//防止越界
if( xl < 0 ) xl = x;
if( xr >= nSliceWidth ) xr = nSliceWidth - 1;
if( yl < 0 ) yl = y;
if( yr >= nSliceHeight ) yr = nSliceHeight - 1;
if( zl < 0 ) zl = z;
if( zr >= nSliceNum ) zr = nSliceNum - 1;
//通过灰度差分计算立方体的8个顶点在X,Y,Z方向的梯度
gx[0] = (double)(GridIntensity[1] - pDataField[z*nPointsPerSlice + y*nSliceWidth + xl]) / 2.0;
gx[1] = (double)(pDataField[z*nPointsPerSlice + y*nSliceWidth + xr] - GridIntensity[0]) / 2.0;
gx[2] = (double)(pDataField[(z+1)*nPointsPerSlice + y*nSliceWidth + xr] - GridIntensity[3]) / 2.0;
gx[3] = (double)(GridIntensity[2] - pDataField[(z+1)*nPointsPerSlice + y*nSliceWidth + xl]) / 2.0;
gx[4] = (double)(GridIntensity[5] - pDataField[z*nPointsPerSlice + (y+1)*nSliceWidth + xl]) / 2.0;
gx[5] = (double)(pDataField[z*nPointsPerSlice + (y+1)*nSliceWidth + xr] - GridIntensity[4]) / 2.0;
gx[6] = (double)(pDataField[(z+1)*nPointsPerSlice + (y+1)*nSliceWidth + xr] - GridIntensity[7]) / 2.0;
gx[7] = (double)(GridIntensity[6] - pDataField[(z+1)*nPointsPerSlice + (y+1)*nSliceWidth + xl]) / 2.0;
gy[0] = (double)(GridIntensity[4] - pDataField[z*nPointsPerSlice + yl*nSliceWidth + x]) / 2.0;
gy[1] = (double)(GridIntensity[5] - pDataField[z*nPointsPerSlice + yl*nSliceWidth + (x+1)]) / 2.0;
gy[2] = (double)(GridIntensity[6] - pDataField[(z+1)*nPointsPerSlice + yl*nSliceWidth + (x+1)]) / 2.0;
gy[3] = (double)(GridIntensity[7] - pDataField[(z+1)*nPointsPerSlice + yl*nSliceWidth + x]) / 2.0;
gy[4] = (double)(pDataField[z*nPointsPerSlice + yr*nSliceWidth + x] - GridIntensity[0]) / 2.0;
gy[5] = (double)(pDataField[z*nPointsPerSlice + yr*nSliceWidth + (x+1)] - GridIntensity[1]) / 2.0;
gy[6] = (double)(pDataField[(z+1)*nPointsPerSlice + yr*nSliceWidth + x+1] - GridIntensity[2]) / 2.0;
gy[7] = (double)(pDataField[(z+1)*nPointsPerSlice + yr*nSliceWidth + x] - GridIntensity[3]) / 2.0;
gz[1] = (double)(GridIntensity[2] - pDataField[zl*nPointsPerSlice + y*nSliceWidth + (x+1)]) / 2.0;
gz[2] = (double)(pDataField[zr*nPointsPerSlice + y*nSliceWidth + x+1] - GridIntensity[1]) / 2.0;
gz[3] = (double)(pDataField[zr*nPointsPerSlice + y*nSliceWidth + x] - GridIntensity[0]) / 2.0;
gz[4] = (double)(GridIntensity[7] - pDataField[zl*nPointsPerSlice + (y+1)*nSliceWidth + x]) / 2.0;
gz[5] = (double)(GridIntensity[6] - pDataField[zl*nPointsPerSlice + (y+1)*nSliceWidth + (x+1)]) / 2.0;
gz[6] = (double)(pDataField[zr*nPointsPerSlice + (y+1)*nSliceWidth + (x+1)] - GridIntensity[5]) / 2.0;
gz[7] = (double)(pDataField[zr*nPointsPerSlice + (y+1)*nSliceWidth + x] - GridIntensity[4]) / 2.0;
//由索引值查边索引表确定和等值面相交的边,
//并调用函数CalculateIntersection()计算等值点坐标和法向量
//将等值点坐标存放在Vertices[]中,法向量存放在Normals[]中
if (EdgeTable[TableIndex] & 1)
Vertices[0] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 0 );
if (EdgeTable[TableIndex] & 2)
Vertices[1] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 1 );
if (EdgeTable[TableIndex] & 4)
Vertices[2] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 2 );
if (EdgeTable[TableIndex] & 8)
Vertices[3] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 3 );
if (EdgeTable[TableIndex] & 16)
Vertices[4] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 4 );
if (EdgeTable[TableIndex] & 32)
Vertices[5] = CalculateIntersection(x, y, z,gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 5 );
if (EdgeTable[TableIndex] & 64)
Vertices[6] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 6 );
if (EdgeTable[TableIndex] & 128)
Vertices[7] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 7 );
if (EdgeTable[TableIndex] & 256)
Vertices[8] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 8 );
if (EdgeTable[TableIndex] & 512)
Vertices[9] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 9 );
if (EdgeTable[TableIndex] & 1024)
Vertices[10] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 10 );
if (EdgeTable[TableIndex] & 2048)
Vertices[11] = CalculateIntersection(x, y, z, gx, gy, gz, Normals, CubeLengthX, CubeLengthY, CubeLengthZ, 11 );
//根据索引值查找三角片索引表,确定三角片组成,
//并将三角片的顶点坐标和法向量存放在pTriangles所指向的三角片存放空间
for ( int i = 0; TriTable[TableIndex][i] != -1; i += 3)
{
pTriangles[nTriangles].p[0] = Vertices[TriTable[TableIndex][i]];
pTriangles[nTriangles].p[1] = Vertices[TriTable[TableIndex][i+1]];
pTriangles[nTriangles].p[2] = Vertices[TriTable[TableIndex][i+2]];
pTriangles[nTriangles].n[0] = Normals[TriTable[TableIndex][i]];
pTriangles[nTriangles].n[1] = Normals[TriTable[TableIndex][i+1]];
pTriangles[nTriangles].n[2] = Normals[TriTable[TableIndex][i+2]];
nTriangles++;
}
}
}
}
return nTriangles;
}
4.POINT3D CSurfaceReconstructionView::CalculateIntersection(int x, int y, int z, double* gx, double* gy, double* gz, POINT3D* pNorm, double CubeLengthX, double CubeLengthY, double CubeLengthZ, int nEdgeNo)
/********************************************************************************************
*CalculateIntersection(int x, int y, int z, double* gx, double* gy, double* gz,
POINT3D* pNorm, double CubeLengthX, double CubeLengthY, double CubeLengthZ, int nEdgeNo)
*功能:计算等值点的坐标和法向量
*参数:int x, 立方体原点的X坐标
int y, 立方体原点的Y坐标
int z, 立方体原点的Z坐标
double* gx, 指向立方体顶点X方向梯度值的指针
double* gy, 指向立方体顶点Y方向梯度值的指针
double* gz, 指向立方体顶点Z方向梯度值的指针
POINT3D* pNorm, 指向等值点法向量的指针
double CubeLengthX, X方向立方体边长
double CubeLengthY, Y方向立方体边长
double CubeLengthZ, Z方向立方体边长
int nEdgeNo, 立方体中当前边的边号
*返回值:POINT3D,等值点的三维坐标
*********************************************************************************************/
POINT3D CSurfaceReconstructionView::CalculateIntersection(int x, int y, int z, double* gx, double* gy, double* gz, POINT3D* pNorm, double CubeLengthX, double CubeLengthY, double CubeLengthZ, int nEdgeNo)
{
double x1, y1, z1, x2, y2, z2; //构成当前边的两个顶点的坐标(x1, y1, z1)和(x2, y2, z2)
int index1, index2; //构成当前边的两个顶点的索引(编号)
int v1x = x, v1y = y, v1z = z; //一条边的顶点位置(v1x,v1y,v1z),初始化为(x,y,z)
int v2x = x, v2y = y, v2z = z; //一条边的顶点位置(v2x,v2y,v2z),初始化为(x,y,z)
double nx1, nx2, ny1, ny2, nz1, nz2; //构成当前边的两顶点的法向量
POINT3D intersection; //当前边和等值面的交点,即等值点
POINT3D Norm; //等值点的法向量
//据边号确定该边的两顶点位置和索引
switch (nEdgeNo)
{
case 0:
v2x += 1;
index1 = 0;
index2 = 1;
break;
case 1:
v1x += 1;
v2x += 1;
v2z += 1;
index1 = 1;
index2 = 2;
break;
case 2:
v1x += 1;
v1z += 1;
v2z += 1;
index1 = 2;
index2 = 3;
break;
case 3:
v2z += 1;
index1 = 0;
index2 = 3;
break;
case 4:
v1y += 1;
v2x += 1;
v2y += 1;
index1 = 4;
index2 = 5;
break;
case 5:
v1x += 1;
v1y += 1;
v2x += 1;
v2y += 1;
v2z += 1;
index1 = 5;
index2 = 6;
break;
case 6:
v1x += 1;
v1y += 1;
v1z += 1;
v2y += 1;
v2z += 1;
index1 = 6;
index2 = 7;
break;
case 7:
v1y += 1;
v2y += 1;
v2z += 1;
index1 = 4;
index2 = 7;
break;
case 8:
v2y += 1;
index1 = 0;
index2 = 4;
break;
case 9:
v1x += 1;
v2x += 1;
v2y += 1;
index1 = 1;
index2 = 5;
break;
case 10:
v1x += 1;
v1z += 1;
v2x += 1;
v2y += 1;
v2z += 1;
index1 = 2;
index2 = 6;
break;
case 11:
v1z += 1;
v2y += 1;
v2z += 1;
index1 = 3;
index2 = 7;
break;
}
//计算当前边两顶点的坐标
x1 = v1x * CubeLengthX;
y1 = v1y * CubeLengthY;
z1 = v1z * CubeLengthZ;
x2 = v2x * CubeLengthX;
y2 = v2y * CubeLengthY;
z2 = v2z * CubeLengthZ;
//采用中点选择法,由当前边两顶点的坐标计算等值点坐标
intersection.x = (x2 + x1) / 2.0;
intersection.y = (y1 + y2) / 2.0;
intersection.z = (z1 + z2) / 2.0;
//计算当前边两顶点的法向量
nx1 = gx[index1] / CubeLengthX;
ny1 = gy[index1] / CubeLengthY;
nz1 = gz[index1] / CubeLengthZ;
nx2 = gx[index2] / CubeLengthX;
ny2 = gy[index2] / CubeLengthY;
nz2 = gz[index2] / CubeLengthZ;
//采用中点选择法,由当前边两顶点的法向量计算等值点的法向量
Norm.x = (nx1 + nx2) / 2.0;
Norm.y = (ny1 + ny2) / 2.0;
Norm.z = (nz1 + nz2) / 2.0;
//法向量归一化处理
double squaroot = sqrt( Norm.x * Norm.x + Norm.y * Norm.y + Norm.z * Norm.z );
if( squaroot <= 1e-6 )
squaroot = 1.0;
Norm.x /= squaroot;
Norm.y /= squaroot;
Norm.z /= squaroot;
*(pNorm + nEdgeNo) = Norm; //存放等值点法向量
return intersection; //返回等值点坐标
}
5.void CSurfaceReconstructionView::GLsetting(void)
void CSurfaceReconstructionView::GLsetting(void)
{
HGLRC hRC; //渲染描述表句柄
pDrawDC = new CClientDC(this); //使pDrawDC指向当前DC
ASSERT(pDrawDC !=NULL); //断言pDrawDC不为空
if(!PixelformatSetting()) //设置像素格式
return;
//创建和当前DC兼容的RC,并和当前DC关联
hRC = wglCreateContext(pDrawDC->GetSafeHdc());
wglMakeCurrent(pDrawDC->GetSafeHdc(),hRC);
//定义和设置环境光、漫反射光、光源位置
GLfloat LightAmbient[] = {0.6f, 0.6f, 0.6f, 1.0f };
GLfloat LightDiffuse[] = {1.0f, 1.0f, 1.0f, 1.0f };
GLfloat LightPosition[] = {0.0f, 0.0f, -1.0f, 0.0f };
glLightfv( GL_LIGHT1, GL_AMBIENT, LightAmbient );
glLightfv( GL_LIGHT1, GL_DIFFUSE, LightDiffuse );
glLightfv( GL_LIGHT1, GL_POSITION, LightPosition );
}
6.void CSurfaceReconstructionView::DrawSurface(void)
void CSurfaceReconstructionView::DrawSurface(void)
{
double px, py, pz; //顶点坐标
double nx, ny, nz; //顶点法向量
int i, j;
glEnable( GL_LINE_SMOOTH ); //接受线段平滑处理
glColor3f( 1.0f, 1.0f, 1.0f ); //指定绘制线条颜色为白色
glPolygonMode( GL_FRONT_AND_BACK, GL_LINE ); //前后面显示为边线框架方式
//绘制基于MC方法抽取的等值面
if( m_bMCSurface )
{
glEnable( GL_DEPTH_TEST ); //启动深度测试
glEnable( GL_LIGHTING); //启动光照
glEnable( GL_LIGHT1 ); //使用GL_LIGHT1光照
glTranslatef( -50.0f, -60.0f, -750.0f ); //移动物体到显示区
glRotated( RotateAngle, 1.0, 1.0, 1.0 ); //绕从原点指向(1.0, 1.0, 1.0)的射线旋转RotateAngle度
glBegin( GL_TRIANGLES ); //绘制三角形
for( i = 0; i < m_nTriangles; i++ )
{
for( j = 0; j < 3; j++ )
{
px = pTriangles[i].p[j].x;
py = pTriangles[i].p[j].y;
pz = pTriangles[i].p[j].z;
nx = pTriangles[i].n[j].x;
ny = pTriangles[i].n[j].y;
nz = pTriangles[i].n[j].z;
glVertex3f( (GLfloat)px, (GLfloat)py, (GLfloat)pz );
glNormal3f( (GLfloat)nx, (GLfloat)ny, (GLfloat)nz );
}
}
glEnd();
glDisable( GL_LIGHT1 ); //取消GL_LIGHT1
glDisable( GL_LIGHTING); //关闭光照
glDisable( GL_DEPTH_TEST ); //关闭深度测试
}
glDisable( GL_LINE_SMOOTH ); //取消线段平滑处理
}
7.void CSurfaceReconstructionView::OnDraw(CDC* /*pDC*/)
void CSurfaceReconstructionView::OnDraw(CDC* /*pDC*/)
{
CSurfaceReconstructionDoc* pDoc = GetDocument();
ASSERT_VALID(pDoc);
if (!pDoc)
return;
// TODO: 在此处为本机数据添加绘制代码
static BOOL bBusy = FALSE; //定义开关变量
if( bBusy ) return; //开关,使绘制完成后才可更新缓存
bBusy = TRUE;
glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT ); //清除深度缓存和颜色缓存
glMatrixMode( GL_MODELVIEW); //启动模型矩阵
glLoadIdentity(); //初始化为单位矩阵
//若已获得三维体表面数据,则调用DrawSurface()完成表面绘制
if( m_bValidSurface )
{
DrawSurface();
}
glFinish(); //绘制完成
SwapBuffers( wglGetCurrentDC() ); //更新缓存
bBusy = FALSE;
}
8.CSurfaceReconstructionView::~CSurfaceReconstructionView()
//析构函数释放已分配的空间
CSurfaceReconstructionView::~CSurfaceReconstructionView()
{
if( m_bValidSurface )
{
GlobalUnlock( hTriAngles );
GlobalFree( hTriAngles );
m_bValidSurface = FALSE;
}
}
9.int CSurfaceReconstructionView::OnCreate(LPCREATESTRUCT lpCreateStruct)
int CSurfaceReconstructionView::OnCreate(LPCREATESTRUCT lpCreateStruct)
{
if (CView::OnCreate(lpCreateStruct) == -1)
return -1;
// TODO: Add your specialized creation code here
GLsetting();
return 0;
}
10.void CSurfaceReconstructionView::OnDestroy()
void CSurfaceReconstructionView::OnDestroy()
{
CView::OnDestroy();
// TODO: Add your message handler code here
HGLRC hRC;
hRC = ::wglGetCurrentContext(); //获取当前RC句柄
::wglMakeCurrent( NULL, NULL ); //解除当前RC和DC的关联,并把当前RC非当前化
if( hRC )
{
::wglDeleteContext( hRC ); //删除RC
}
if( pDrawDC )
{
delete pDrawDC;
}
}
11.void CSurfaceReconstructionView::OnSize(UINT nType, int cx, int cy)
void CSurfaceReconstructionView::OnSize(UINT nType, int cx, int cy)
{
CView::OnSize(nType, cx, cy);
// TODO: Add your message handler code here
if( cy > 0 )
{
glMatrixMode( GL_PROJECTION ); // 启动投影矩阵
glLoadIdentity(); // 初始化为单位阵
// 设置透视投影变换,指定视场
glFrustum( -1.0, 1.0, -1.0*cy/cx, 1.0*cy/cx, 3.0, 1000.0 );
// 设置视口,即定义显示范围
glViewport( 0, 0, cx, cy );
}
RedrawWindow(); //显示更新
}
12.void CSurfaceReconstructionView::OnLButtonDown(UINT nFlags, CPoint point)
void CSurfaceReconstructionView::OnLButtonDown(UINT nFlags, CPoint point)
{
// TODO: Add your message handler code here and/or call default
if( m_bValidSurface )
{
if( bRotate )
{
SetTimer(1, 100, NULL ); //设置定时器1,100毫秒触发一次
}
else
{
KillTimer(1); //移除定时器1
}
bRotate = !bRotate;
}
CView::OnLButtonDown(nFlags, point);
}
13.void CSurfaceReconstructionView::OnTimer(UINT_PTR nIDEvent)
void CSurfaceReconstructionView::OnTimer(UINT_PTR nIDEvent)
{
// TODO: Add your message handler code here and/or call default
//若定时器1计时到达预定时刻,则旋转角度增加10度
if( nIDEvent == 1 )
{
RotateAngle += 10.0f;
Invalidate(FALSE); //使重新绘制
}
CView::OnTimer(nIDEvent);
}
第六步:将包含二维断层图像的T4文件夹放在工程Debug(可执行文件.exe所在的)路径下,编译运行程序,在弹出的MFC窗口中单击M按钮即可出现三维重建之后的脊椎骨图形。如下图:
总结,整个绘图的思路是:先建立OpenGL的绘图环境,然后读入二维图像创建三维数据场将其存放在hTriangles所指向的结构体指针(实际是数组)中,通过OpenGL在DrawSurface函数中进行体绘制。
参考文献:张俊华.<医学图像三维重建和可视化——VC++实现>.科学出版社