VTK学习笔记之使用vtkMarchingCubes

前段时间学习了一下VTK的基础使用方法，由于关于这方面的中文教程比较少，再加上英文水平有不是很高，所以学习过程相当缓慢。中途又由于期中考试耽搁了。

最近接到一个大规模散乱数据场三维重现的一个小程序，所以又开始研究起VTK。 以下分两个部分：

一，Marching Cubes方法

 在三维规则数据场中构造等值面是计算机视觉也是VTK当中的一个重要的问题，Marching Cubes方法就是解决这个问题的一个成熟的被广泛使用的方法，简称MC方法。VTK当中的vtkMarchingCubes就是对这个方法的代码实现。

下面结合唐泽圣的《三维数据场可视化》和vtk的《The Design and Implementation of an Object-OrientedToolkit for 3D Graphics and Visualization》两本书的内容介绍一下MC方法。

在MC方法中, 假定原始数据是离散的三维空间规则数据场。用于医疗诊断的断层扫描仪( CT ) 及核磁共振仪( MRI ) 等产生的图象均属于这一类型。为了在这一数据场中构造等值面, 用户应先给出所求等值面的值, 设为C0。MC 方法首先找出该等值面经过的体元的位置, 求出该体元内的等值面并计算出相关参数, 以便由常用的图形软件包或图形硬件提供的面绘制功能绘制出等值面。由于这一方法是逐个体元依次处理的, 因此被称为Marching
Cubes方法。

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

我们假定某体元一个角点的函数值大于( 或等于) 给定的等值面的值C0 , 则将该角点赋值为1, 并称该角点位于等值面之内( 或之上) 。如果该角点的函数值小于等值面的值C0 , 则将该角点赋值为零, 并称该角点位于等值面之外。显然, 如果某体元中一条边的一个角点在等值面之内, 而另一个角点在等值面之外, 那么, 该边必然与所求等值面相交。根据这一原理就可以判断所求等值面将与哪些体元相交, 或者说将穿过哪些体元。

由于每个顶点都有0和1这两个可能的值，所以一个体元中可能的等值面情况是256（2的8次方）种。但是可以用两种对称性将256中情况降到15中情况。第一，若将体元中值等于1的顶点的值和值等于0的顶点的值进行对换，等值面的情况不变，这样可能的情况变成原来的一半。第二，立方体是中心对称的，根据这个特性情况再次减少，降为15种。如下图所示。



2) 求等值面与体元边界的交点

当三维离散数据场的密度较高, 也就是说当每个体元很小时, 可以假定函数值沿体元边界呈线性变化。因此, 等值面与体元边界的交点可以通过该边两端点函数值的线性插值求出。求出了等值面与体元边界的交点以后, 就可以将这些交点连接成三角形或多边形, 形成等值面的一部分。

即在1)的15种情况下，根据体元每个顶点的实际值和等值面的值，确定等值面与体元的棱的交点所在的具体位置。确定的方法是简单的线性插值。例如某棱的两个顶点的值分别是0和1，等值面的值是0.5，则等值面与该棱的交点就在棱的中点位置。

3) 求等值面的法向

为了利用图形硬件显示等值面图象, 必须给出形成等值面的各三角面片的法向。这个一般直接使用vtkMarchingCubes编程可以先不关注。

二，使用vtkMarchingCubes

帮助文档中关于vtkMarchingCubes的详细描述为：

vtkMarchingCubes is a filter that takes as input a volume (e.g., 3D structured point set) and generateson output one or more isosurfaces. One or more contour values must be specifiedto generate the isosurfaces. Alternatively,
you can specify a min/max scalarrange and the number of contours to generate a series of evenly spaced contourvalues.

即：vtkMarchingCubes是一个过滤器，该过滤器接受一个体数据的输入（三维规则的点集），输出一个或者多个等值面。使用它，必须给等值面赋值。或者，你可以设定的一个scalar的范围和等值面的间隔来决定这个范围内的一系列的等值面的值。

几个重要的方法是：

void SetValue (int i, double value) 设置第i个等值面的值为value

void SetNumberOfContours (int number) 设置等值面的个数

第一种用法是从三维文件当中读取数据，用MC方法划分等值面，获得结果。第一个例子是VTK自带的例子headBone.tcl，使用TCL脚本语言编写的。

package require vtk
package require vtkinteraction
package require vtktesting

# Create the RenderWindow, Renderer and both Actors
#
vtkRenderer ren1
vtkRenderWindow renWin
renWin AddRenderer ren1
vtkRenderWindowInteractor iren
iren SetRenderWindow renWin
vtkLight lgt

# create pipeline
#
vtkMergePoints locator
locator SetDivisions 32 32 46
locator SetNumberOfPointsPerBucket 2
locator AutomaticOff

vtkVolume16Reader v16
v16 SetDataDimensions 64 64
[v16 GetOutput] SetOrigin 0.0 0.0 0.0
v16 SetDataByteOrderToLittleEndian
v16 SetFilePrefix "$VTK_DATA_ROOT/Data/headsq/quarter"
v16 SetImageRange 1 93
v16 SetDataSpacing 3.2 3.2 1.5

vtkMarchingCubes iso
iso SetInputConnection [v16 GetOutputPort]
iso SetValue 0 1150
iso ComputeGradientsOn
iso ComputeScalarsOff
iso SetLocator locator

vtkVectorNorm gradient
gradient SetInputConnection [iso GetOutputPort]

vtkDataSetMapper isoMapper
isoMapper SetInputConnection [gradient GetOutputPort]
isoMapper ScalarVisibilityOn
isoMapper SetScalarRange 0 1200
isoMapper ImmediateModeRenderingOn

vtkActor isoActor
isoActor SetMapper isoMapper
set isoProp [isoActor GetProperty]
eval $isoProp SetColor $antique_white

vtkOutlineFilter outline
outline SetInputConnection [v16 GetOutputPort]
vtkPolyDataMapper outlineMapper
outlineMapper SetInputConnection [outline GetOutputPort]
vtkActor outlineActor
outlineActor SetMapper outlineMapper
set outlineProp [outlineActor GetProperty]
#eval $outlineProp SetColor 0 0 0

# Add the actors to the renderer, set the background and size
#
ren1 AddActor outlineActor
ren1 AddActor isoActor
ren1 SetBackground 1 1 1
ren1 AddLight lgt
renWin SetSize 250 250
ren1 SetBackground 0.1 0.2 0.4

ren1 ResetCamera
set cam1 [ren1 GetActiveCamera]
$cam1 Elevation 90
$cam1 SetViewUp 0 0 -1
$cam1 Zoom 1.5
eval lgt SetPosition [$cam1 GetPosition]
eval lgt SetFocalPoint [$cam1 GetFocalPoint]

ren1 Res
4000
etCameraClippingRange

# render the image
#
iren AddObserver UserEvent {wm deiconify .vtkInteract}

renWin Render

# prevent the tk window from showing up then start the event loop
wm withdraw .

输出结果是：



第二种用法是自己构造规则的三维数据场作为输入，这种用法用在了将多轮廓线重构转换成等值线的划分方法中，

这里就不做介绍。VTK中可以表示三维规则数据场的一个重要的结构是vtkImageData。

vtkImageData *image = vtkImageData::New();
image->SetDimensions(30, 30, 30);
image->SetSpacing(5, 5, 5);
image->SetOrigin(0, 0, 0);
image->SetScalarTypeToInt();
image->SetNumberOfScalarComponents(1);
image->AllocateScalars();
int *ptr = (int *)image->GetScalarPointer();

for(int i=0; i<30*30*30; i++)
{
if(i%9 == 0)
*ptr++ = 0;
else
*ptr++ = 2;
}

vtkMarchingCubes *iso = vtkMarchingCubes::New();
iso->SetInput(image);
iso->SetNumberOfContours(1);
iso->SetValue(0,1);
iso->ComputeGradientsOn();
iso->ComputeNormalsOff();
iso->ComputeScalarsOff();

vtkVectorNorm *gradiant = vtkVectorNorm::New();
gradiant->SetInputConnection(iso->GetOutputPort());

vtkDataSetMapper *cubeMapper = vtkDataSetMapper::New();
cubeMapper->SetInputConnection(gradiant->GetOutputPort());
vtkActor *cubeActor = vtkActor::New();
cubeActor->SetMapper(cubeMapper);

vtkCamera *camera = vtkCamera::New();
camera->SetPosition(1,1,1);
camera->SetFocalPoint(0,0,0);

vtkRenderer *renderer = vtkRenderer::New();
vtkRenderWindow *renWin = vtkRenderWindow::New();
renWin->AddRenderer(renderer);

vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
iren->SetRenderWindow(renWin);

renderer->AddActor(cubeActor);
renderer->SetActiveCamera(camera);
renderer->ResetCamera();
renderer->SetBackground(1,1,1);

renWin->SetSize(300,300);
renWin->Render();
iren->Start();

image->Delete();
cubeMapper->Delete();
cubeActor->Delete();
camera->Delete();
renderer->Delete();
renWin->Delete();
iren->Delete();

输出结果为：