Marching Cubes算法再回顾

本文主要是介绍Marching Cubes算法再回顾，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

1,确定包含等值面的体元

首先介绍一下 体元的概念，体元是三维图像中由相邻的八个体素点组成的正方体方格，英语也叫 Cube，体元中角点函数值分为两种情况，一种是大于等于给定等值面的值 C0 ,则将角点设为 1 称该角点在等值面内部，否则设为0，在等值面之外，

一般来说，会出现一个角点在内，一个角点在外，则角点之间的连线(也就是体元的边)必然与等值面相交，根据这个原理就能判断等值面与哪些体元相交。

——————————————————————————————————————

三维空间中，平行且相邻的两个二维图像（每个图像中的正方形四个像素顶点组成一个基本的像素图像单元）组成一个基本的三维图像单元。，下图中由6个这样的基本三维图像单元：

vtkImageData结构由尺寸、间距和原点来定义。尺寸标注是沿着每个主轴的体素或像素的数量。原点是数据的第一个切片的左下角的世界坐标位置。间距是沿三个主要轴的像素之间的距离。

原点是数据集左下角的世界坐标位置。

尺寸是沿着三个主要轴的体素或像素的数量。

间距是体素的高度、长度和宽度，或相邻像素之间的距离，这取决于是将数据视为相同的方框还是连续函数中的样本点。
——————————————————————————————————————————

Marching Cubes算法根据一个立方体的8个顶点，判断这8个顶点的每个顶点在等值面的内部还是外部（每个顶点只有“在等值面内”和“在等值面外”这两种状态，设为0和1）从而根据这8个顶点的状态，建立一个包含共256种状态的查找表（根据平面对称性、中心对称性，256种最终降到15种）。

顶点值高于等值在表面的内部，等于等值在表面上，低于等值在表面外。

体元的每个顶点有两种状态，总共有256种，可以制作一个查找表(look up table)

但由于反转状态不变，所以可以减少一半，为128种。

再根据旋转不变形，又可以减少到15种情况。

可以认为这15种情况类似于基，经过旋转，反转可以得到256种状态对应的结果Triangulated Cubes：

根据每个顶点的状态，我们可以为每类制作一个8位索引Cube Numbering：

索引指向边表，给出了边的交叉情况。

相交边的编码——通过编码记录对应的cube，相交边的编号

二进制：00000010

十进制：2

Table[256]表示哪些边有交点

Table[2]=0x103=0000 0010 0000 0011

表示0,1,9号边上有交点

为了避免每次转化成二进制进行解码，可以直接记录与哪些边有交点，之后直接查表即可。

Table2[256][16]

Table2[2]=(0, 1, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

因为每个正方体中最多有1个或4个三角形，所以线性插值就足够了。

2，确定等值面与体元边界的交点

找到含有等值面的体元之后，接下来就是确定等值面与体元边界的交点，体元间的数值都是呈线性变化，求交点时一般采用的是线性插值，如 Case0 中等值面的两个端点一个在外为( 标记0) ，一个在内 ( 标记为1 ) 则交点为0.5；

3，求等值面的法向量

以上步骤 1，2，3 为实现 MC 算法步骤流程，但利用 VTK ，不需要这么繁琐，主要算法步骤都已经封装到 vtkMarchingCube 类中，使用 vtkMarchingCube 时，需要设置三个参数：

SetValue(int i,double value) 设置第i 个等值面的值b，(提醒一下，医学图像中的灰度值范围不是 0-256 而是0-65326，但大部分取值范围都在0-1000）。
SetNumberofContours(int number)，设置等值面的个数
ComputerNormalsOn() 设置计算等值面的法向量，提高渲染质量；

上面这张图显示的就是 vtk 呈像的基本流程，下面是仿照官网写的用面绘制来对图像重建的代码部分：

#include<vtkRenderWindow.h>
#include<vtkRenderWindowInteractor.h>
#include<vtkDICOMImageReader.h>
#include<vtkMarchingCubes.h>
#include<vtkPolyDataMapper.h>
#include<vtkStripper.h>
#include<vtkActor.h>
#include<vtkProperty.h>
#include<vtkCamera.h>
#include<vtkOutlineFilter.h>
#include<vtkOBJExporter.h>
#include<vtkRenderer.h>
#include<vtkMetaImageReader.h>
#include<vtkInteractorStyleTrackballCamera.h>#include<iostream>
#include<string.h>
//需要进行初始化，否则会报错
#include <vtkAutoInit.h> 
#include<vtkRenderingVolumeOpenGL2ObjectFactory.h>
#include<vtkRenderingOpenGL2ObjectFactory.h>using namespace std;
int main()
{///Marching Cube; vtkObjectFactory::RegisterFactory(vtkRenderingOpenGL2ObjectFactory::New());vtkObjectFactory::RegisterFactory(vtkRenderingVolumeOpenGL2ObjectFactory::New());vtkSmartPointer<vtkRenderer> ren = vtkSmartPointer<vtkRenderer>::New();vtkSmartPointer<vtkRenderWindow> renWin = vtkSmartPointer<vtkRenderWindow>::New();//WINDOW；renWin->AddRenderer(ren);vtkSmartPointer<vtkRenderWindowInteractor> iren = vtkSmartPointer<vtkRenderWindowInteractor>::New();//wininteratcor;iren->SetRenderWindow(renWin);vtkSmartPointer<vtkDICOMImageReader> reader = vtkSmartPointer<vtkDICOMImageReader>::New();reader->SetDirectoryName("E:/DIcom_Data/DICOM");reader->SetDataByteOrderToLittleEndian();reader->Update();/*vtkDICOMImageReader *reader = vtkDICOMImageReader::New();reader->SetDirectoryName("E:/Coding Pra/VTK/VTK_Examples_StandardFormats_Input_DicomTestImages/DICOM");reader->SetDataByteOrderToLittleEndian();reader->Update();*/cout << "读取数据完毕" << endl;cout << "The width is" << reader->GetWidth() << endl;cout << "The height is" << reader->GetHeight() << endl;cout << "The depth is" << reader->GetPixelSpacing() << endl;cout << "The Output port is" << reader->GetOutputPort() << endl;vtkSmartPointer<vtkMarchingCubes> marchingcube = vtkSmartPointer<vtkMarchingCubes>::New();marchingcube->SetInputConnection(reader->GetOutputPort());//获得读取的数据的点集；marchingcube->SetValue(0, 200);//Setting the threshold;marchingcube->ComputeNormalsOn();//计算表面法向量;vtkSmartPointer<vtkStripper> Stripper = vtkSmartPointer<vtkStripper>::New();Stripper->SetInputConnection(marchingcube->GetOutputPort());//获取三角片vtkSmartPointer<vtkPolyDataMapper> Mapper = vtkSmartPointer<vtkPolyDataMapper>::New();//将三角片映射为几何数据；Mapper->SetInputConnection(Stripper->GetOutputPort());Mapper->ScalarVisibilityOff();//vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();//Created a actor;actor->SetMapper(Mapper);//获得皮肤几何数据actor->GetProperty()->SetDiffuseColor(1, .49, .25);//设置皮肤颜色；actor->GetProperty()->SetSpecular(0.3);//反射率；actor->GetProperty()->SetOpacity(1.0);//透明度；actor->GetProperty()->SetSpecularPower(20);//反射光强度；actor->GetProperty()->SetColor(1, 0, 0);//设置角的颜色；actor->GetProperty()->SetRepresentationToWireframe();//线框；//vtkSmartPointer<vtkCamera> camera = vtkSmartPointer<vtkCamera>::New();//Setting the Camera;//camera->SetViewUp(0, 0, -1);//设置相机向上方向；//camera->SetPosition(0, 1, 0);//位置：世界坐标系，相机位置；//camera->SetFocalPoint(0, 0, 0);//焦点，世界坐标系，控制相机方向；//camera->ComputeViewPlaneNormal();//重置视平面方向，基于当前的位置和焦点；vtkSmartPointer<vtkOutlineFilter> outfilterline = vtkSmartPointer<vtkOutlineFilter>::New();outfilterline->SetInputConnection(reader->GetOutputPort());vtkSmartPointer<vtkPolyDataMapper> outmapper = vtkSmartPointer<vtkPolyDataMapper>::New();outmapper->SetInputConnection(outfilterline->GetOutputPort());vtkSmartPointer<vtkActor> OutlineActor = vtkSmartPointer<vtkActor>::New();OutlineActor->SetMapper(outmapper);OutlineActor->GetProperty()->SetColor(0, 0, 0);//线框颜色ren->AddActor(actor);ren->AddActor(OutlineActor);//ren->SetActiveCamera(camera);//设置渲染器的相机；ren->ResetCamera();ren->ResetCameraClippingRange();//camera->Dolly(1.5);//使用Dolly()方法延伸着视平面法向移动相机；ren->SetBackground(1, 1, 1);//设置背景颜色；renWin->SetSize(1000, 600);vtkInteractorStyleTrackballCamera *style = vtkInteractorStyleTrackballCamera::New();iren->SetInteractorStyle(style);renWin->Render();iren->Initialize();iren->Start();vtkSmartPointer<vtkOBJExporter> porter = vtkSmartPointer<vtkOBJExporter>::New();porter->SetFilePrefix("E:/ceshi/aaa/regist_after/polywrite.obj");//重建图像输出porter->SetInput(renWin);porter->Write();return EXIT_SUCCESS;
}

上面就是 VTK 基于 Marching Cube算法实现的重建效果：

体绘制重建

体绘制时分为两部分：

1，定义 vtkVoluemRayCastMapper 对象

体绘制中最常用的方法；vtkVolumeRayCastMapper() 光线投影，体绘制时，首先定义一个Mapper 然后接受两个输入：

SetInput(vtkImageDate *) 用于设置输入图像数据；
SetVolumeRayCastFunction(vtkVolumeRayCastFunction *) 用于设置光线投影函数类型；

2，利用 vtkVolumeProperty 定义体绘制属性；

SetScalarOpacity() 设置灰度不透明函数；
SetColor() 颜色传输函数；

3，定义 vtkVolume 对象接收 Mapper对象和 Property 对象

SetMapper()接受 Mapper 对象；
SetProperty() 接受 Property 对象；

vtk 中体绘制核心就是改变 Mapper 和 vtkVolumeRayCastFunction() ，上面中vtkColumeRayCastMapper 只是 VolumeMapper 其中的一种，且投影函数类 vtkVolumeRayCastFunction 一共有三个子类：

vtkVolumeRayCastCompositeFunction
vtkVolumeRayCasMIPFunction、
vtkVolumeRayCastIsosurfaceFunction，
因此，其细分的话vtk中的体绘制也不止一种

而下面这个是最常用到的（`vtkVolumeRayCastMapper + vtkVolumeRayCastCompositeFunction）

//体绘制#include<vtkRenderWindowInteractor.h>
#include<vtkDICOMImageReader.h>
#include<vtkCamera.h>
#include<vtkActor.h>
#include<vtkRenderer.h>
#include<vtkVolumeProperty.h>
#include<vtkProperty.h>
#include<vtkPolyDataNormals.h>
#include<vtkImageShiftScale.h>
#include "vtkVolumeRayCastMapper.h"
#include<vtkPiecewiseFunction.h>
#include<vtkColorTransferFunction.h>
#include<vtkVolumeRayCastCompositeFunction.h>
#include<vtkRenderWindow.h>
#include<vtkImageCast.h>
#include<vtkVolumeRayCastCompositeFunction.h>
#include<vtkOBJExporter.h>
#include<vtkOutlineFilter.h>
#include<vtkPolyDataMapper.h>#include<vtkInteractorStyleTrackballCamera.h>
#include<vtkRenderingVolumeOpenGL2ObjectFactory.h>
#include<vtkRenderingOpenGL2ObjectFactory.h>
#include<vtkMetaImageReader.h>#include<vtkLODProp3D.h>//体绘制加速//Gpu光照映射
#include<vtkGPUVolumeRayCastMapper.h>#include<iostream>int main()
{vtkObjectFactory::RegisterFactory(vtkRenderingOpenGL2ObjectFactory::New());vtkObjectFactory::RegisterFactory(vtkRenderingVolumeOpenGL2ObjectFactory::New());//定义绘制器；vtkRenderer *aRenderer = vtkRenderer::New();//指向指针；vtkSmartPointer<vtkRenderWindow> renWin = vtkSmartPointer<vtkRenderWindow>::New();renWin->AddRenderer(aRenderer);vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();iren->SetRenderWindow(renWin);//读取数据；/*vtkDICOMImageReader *reader = vtkDICOMImageReader::New();reader->SetDirectoryName("E:/Coding Pra/VTK/VTK_Examples_StandardFormats_Input_DicomTestImages/DICOM");reader->SetDataByteOrderToLittleEndian();*/vtkSmartPointer<vtkDICOMImageReader> reader = vtkSmartPointer<vtkDICOMImageReader>::New();reader->SetDirectoryName("E:/DIcom_Data/DICOM");reader->SetDataByteOrderToLittleEndian();//图像数据预处理，类型转换:通过 vtkimageCast 将不同类型数据集转化为 vtk 可以处理的数据集；vtkImageCast *cast_file = vtkImageCast::New();cast_file->SetInputConnection(reader->GetOutputPort());cast_file->SetOutputScalarTypeToUnsignedShort();cast_file->Update();//透明度映射函数定义；vtkPiecewiseFunction *opacityTransform = vtkPiecewiseFunction::New();opacityTransform->AddPoint(0, 0.0);opacityTransform->AddPoint(20, 0.0);opacityTransform->AddPoint(200, 1.0);opacityTransform->AddPoint(300, 1.0);//颜色映射函数定义,梯度上升的vtkColorTransferFunction *colorTransformFunction = vtkColorTransferFunction::New();colorTransformFunction->AddRGBPoint(0.0, 0.0, 0.0, 0.0);colorTransformFunction->AddRGBPoint(64.0, 0.0, 0.0, 0.0);colorTransformFunction->AddRGBPoint(128.0, 1.0, 0.0, 0.0);colorTransformFunction->AddRGBPoint(192.0, 1.0, 0.0, 0.0);colorTransformFunction->AddRGBPoint(255.0, 1.0, 0.0, 0.0);vtkPiecewiseFunction *gradientTransform = vtkPiecewiseFunction::New();gradientTransform->AddPoint(0, 0.0);gradientTransform->AddPoint(20, 2.0);gradientTransform->AddPoint(200, 0.1);gradientTransform->AddPoint(300, 0.1);//体数据属性；vtkVolumeProperty *volumeProperty = vtkVolumeProperty::New();volumeProperty->SetColor(colorTransformFunction);volumeProperty->SetScalarOpacity(opacityTransform);volumeProperty->SetGradientOpacity(gradientTransform);volumeProperty->ShadeOn();//应用volumeProperty->SetInterpolationTypeToLinear();//直线间样条插值；volumeProperty->SetAmbient(0.4);//环境光系数；volumeProperty->SetDiffuse(0.6);//漫反射；volumeProperty->SetSpecular(0.2);volumeProperty->SetSpecularPower(10);//高光强度；计算光照效应；利用 vtkBolumeRayCaseMapper进行计算；//vtkVolumeRayCastMapper *volunemapper = vtkVolumeRayCastMapper::New();//vtkVolumeRayCastCompositeFunction *compositeFunction = vtkVolumeRayCastCompositeFunction::New();//光纤映射类型定义：vtkSmartPointer<vtkVolumeRayCastCompositeFunction> compositecast =vtkSmartPointer<vtkVolumeRayCastCompositeFunction>::New();//Mapper定义,vtkSmartPointer<vtkVolumeRayCastMapper> hiresMapper = vtkSmartPointer<vtkVolumeRayCastMapper>::New();hiresMapper->SetInputData(cast_file->GetOutput());hiresMapper->SetVolumeRayCastFunction(compositecast);vtkSmartPointer<vtkLODProp3D> prop = vtkSmartPointer<vtkLODProp3D>::New();prop->AddLOD(hiresMapper,volumeProperty,0.0);////volunemapper->SetVolumeRayCastFunction(compositeFunction);//载入体绘制方法；//volunemapper->SetInputConnection(cast_file->GetOutputPort());//vtkFixedPointVolumeRayCastMapper *fixedPointVolumeMapper = vtkFixedPointVolumeRayCastMapper::New()//fixedPointVolumeMapper->SetInput()vtkVolume *volume = vtkVolume::New();volume->SetMapper(hiresMapper);volume->SetProperty(volumeProperty);//设置体属性；double volumeView[4] = { 0,0,0.5,1 };vtkOutlineFilter *outlineData = vtkOutlineFilter::New();//线框；outlineData->SetInputConnection(reader->GetOutputPort());vtkPolyDataMapper *mapOutline = vtkPolyDataMapper::New();mapOutline->SetInputConnection(outlineData->GetOutputPort());vtkActor *outline = vtkActor::New();outline->SetMapper(mapOutline);outline->GetProperty()->SetColor(0, 0, 0);//背景纯黑色；aRenderer->AddVolume(volume);aRenderer->AddActor(outline);aRenderer->SetBackground(1, 1, 1);aRenderer->ResetCamera();//重设相机的剪切范围；aRenderer->ResetCameraClippingRange();renWin->SetSize(800, 800);renWin->SetWindowName("测试");vtkRenderWindowInteractor *iren2 = vtkRenderWindowInteractor::New();iren2->SetRenderWindow(renWin);//设置相机跟踪模式vtkInteractorStyleTrackballCamera *style = vtkInteractorStyleTrackballCamera::New();iren2->SetInteractorStyle(style);renWin->Render();iren2->Initialize();iren2->Start();vtkOBJExporter *porter = vtkOBJExporter::New();porter->SetFilePrefix("E:/ceshi/aaa/regist_after/esho.obj");porter->SetInput(renWin);porter->Write();porter->Update();return EXIT_SUCCESS;}