vtk 提取等值面并显示

marchingcube是提取等值面比较通用的算法，本文利用vtk 的marching cube接口提取等值面，并通过其绘制管线把等值面绘制出来。

其原理请参考下文：

1.等值面的定义及其三角面片近似

等值面是空间中的一张曲面，在该曲面上函数F(x，y，z)的值等于某一给定值。准确地讲，是指在某一网格空间中，假若每一结点保存着三变量函数F(x，y，z)，而且网格单元在x，y，z方向上的连续采样值为F(x，y，z)，则对于某一给定值Fi，等值面是由所有满足

S={(x，y，z) | F(x，y，z)=Fi}

的点组成的一张曲面。

按照此严格定义下得到的等值面表达式如下：

F(x，y，z)=a0+a1x+a2y+a3Z+a4xy+a5yz+a6xz+a7xyz

可以看出等值面是三次代数曲面，提取过程复杂而且不利于显示。为了简化等值面的提取，W.ELorenson和H.E.CIine在1987年提取了一种等值面的简化提取方法，该方法首先找到等值面经过的六面体网格，求出该六面体与等值面的交点，将这些交点按照一定的拓扑连接关系连接起来，作为等值面在该六面体网格中的近似表示。如图所示，为了叙述的方便以及程序的编写，对所有六面体的顶点采用如下的编码规定：

假设某一个六面体中顶点3的值比要提取的等值面值小，而其他的7个顶点比等值面值大，则显然等值面一定经过该六面体，而且该六面体内的等值面可以采用如图(4.1.1.2)所示的三角面片进行近似。

该方法虽然在处理时是针对数据场中的某一个六面体单元，似乎等值面的构造是相互独立的。但是处理完所有的六面体单元后，这些等值面之间确是统一，相互关联和连续的。由于它是以六面体为单位，一个接一个的处理，因此被称为MarehingCubes方法。

2.六面体单元网格与等值面的位置关系

MarhcingCubeS方法是采用三角面片来近似六面体内的实际等值面，假设我们得到了一个六面体与等值面的所有交点，这些交点之间该进行怎么样的拓扑连接呢?怎么样的拓扑连接才会保证不同六面体之间拓扑连接的一致性呢?为此，我们将根据六面体8个顶点的值与等值面的位置关系，制定不同的拓扑连接关系。

就一个六面体而言，假设某一个顶点的值大于(或者等于)给定的等值面值C，我们就将该点标记位置1，表示该顶点位于等值面之内(或者之上)。而如果某一个顶点的值小于给定的等值面值C，我们就将该点标记位置0，表示该顶点位于等值面之外。如果六面体中某一条边的一个顶点在等值面之内，另一个顶点在等值面之外，那么该边一定与等值面相交。根据这一方法，就可以确定等值面是否与当前处理的六面体相交。

每一个六面体有8个顶点,而每一个顶点都有0和l两种可能的状态,所以每一个六面体单元根据8个顶点0和1的分布,共有2/=256种不同的状态。如果去定义256种可能的拓扑连接,是十分繁琐的,而且容易出现错误。为此,我们利用两种不同的对称性去简化这256种可能的状态。

第一种对称性是利用等值面与8个交点的相对位置关系,如果将等值面的值和8个顶点的物理值的大小关系颠倒过来,六面体内等值面与六面体8个顶点的拓扑关系不会发生改变。也就是说,如果将一个六面体中其物理值大于给定值C的顶点和小于给定值C的顶点所标记的0,1值互换,生成的等值面是相同的(如图(4.1.1.3)。图(a)中顶点3的值大于等值面的值,而其他的7个顶点的值都小于等值面的值,图b()中顶点3的值小于等值面的值,而其他的7个顶点的值都小于等值面的值,在六面体内等值面连接都是一样的,所以可以归纳为一起。

第二种对称性是利用六面体8个顶点的旋转对称性,进一步对可能的状态进行组合”如图(4.1.1.4)所示,在这两种六面体内三角面的拓扑连接相对与图中的黑点来说,是一样的,因此,我们将它们组合成一类。经过上述两种对称性处理之后,我们得到了如图（4.1.1.5)所示的15种类型。

3.等值面与六面体单元边的求交

给定一个六面体单元,如何确定该立方体的类型以及等值面与该六面体的哪条边有交点,交点坐标的插值计算？

为了解决上面提出的问题,我们设计了一个六面体状态表,如图(4.1.1.6)所示,该状态表中的每一位可以表示出该单元体中一个顶点的0和l两种状态。根据这一状态表,我们就可以得到当前单元体的一个索引号。

得到当前单元体的索引号之后,我们就可以在插值边标记数组中知道该单元体的哪些边与等值面有交点。

Marchingcubes算法中进行了如下假设:即六面体网格单元在x,少,:方向是函数值是呈线性变化的。基于该假设,等值面与当前单元体边的交点就可以通过该边两个端点函数值的线性插值得到,插值公式如下:

P=P1+(isovalue一V1)(P2一P1)/(V2一V1)

其中P代表等值点坐标,P1、P2代表两个端点的坐标,V1、V2代表两个端点的物理量值,isovalue代表当前提取的数值。求出等值面与当前单元体所有交点之后,就可以将这些交点按照事先定义号的拓扑连接关系连接成三角面片,作为等值面在该单元中的近似表示。

4.近似等值面三角面片顶点的法向计算

为了利用图形硬件显示提取的等值面图像,必须给出形成等值面的三角面片的法向矢量方向。对于等值面上的每一点,其沿面的切线方向的剃度分量应该为零,因此,该点的梯度矢量的方向就代表了等值面在该点的法向方向。对于规则数据场,经常采用中心差分计算出数据体各个顶点处的梯度,然后通过当前单元体边的两个顶点的梯度进行线性插值,得到三角面片个顶点的梯度,即三角面片个顶点的法向方向。得到各个顶点的法向量之后,即可用光照模型计算出各数据点出的漫反射分量,更加突出提取的等值面边界。

假设数据体中数据点的数值以f(xi,yj,zk)表示,采用中心差分方法就可以求出该数据点处的梯度值,即

5.MacrhingCubes算法提取等值面的流程

(1)将原始数据经过预处理之后,读入特定的数组中;

(2)从网格数据体中提取一个单元体,成为当前单元体”同时获取该单元体的所有信息,例如8个顶点的值,坐标位置等;

(3)将当前单元体8个顶点的函数值与给定等值面值C进行比较,得到该单元体的状态表;

(4)根据当前单元体的状态表索引,找出与等值面相交的单元体棱边,并采用线性插值的方法,计算出各个交点的位置坐标;

(5)利用中心差分法,求出当前单元体8个顶点的法向量,在采用线性插值的方法,得到三角面片各个顶点的法向

(6)根据各个三角面片顶点的坐标,顶点法向量进行等值面图象的绘制.

提取及显示方法：

#include "stdafx.h"

#include <vtkSmartPointer.h>

#include <vtkVolume16Reader.h>

#include <vtkMarchingCubes.h>

#include <vtkVectorNorm.h>

#include <vtkPolyDataMapper.h>

#include "vtkActor.h"

#include "vtkCamera.h"

#include "vtkPolyData.h"

#include "vtkPolyDataMapper.h"

#include "vtkRenderWindow.h"

#include "vtkRenderWindowInteractor.h"

#include "vtkRenderer.h"

#include <vtkMergePoints.h>

1. vtkImageData *image = vtkImageData::New();  
2.     image->SetDimensions(30, 30, 30);  
3.     image->SetSpacing(5, 5, 5);  
4.     image->SetOrigin(0, 0, 0);  
5.     image->SetScalarTypeToInt();  
6.     image->SetNumberOfScalarComponents(1);  
7.     image->AllocateScalars();  
8.     int *ptr = (int *)image->GetScalarPointer();  
9.   
10.     for(int i=0; i<30*30*30; i++)  
11.     {  
12.         if(i%9 == 0)  
13.             *ptr++ = 0;  
14.         else  
15.             *ptr++ = 2;  
16.     }  
17.   
18.     vtkMarchingCubes *iso = vtkMarchingCubes::New();  
19.     iso->SetInput(image);  
20.     iso->SetNumberOfContours(1);  
21.     iso->SetValue(0,1);  
22.     iso->ComputeGradientsOn();  
23.     iso->ComputeNormalsOff();  
24.     iso->ComputeScalarsOff();  
25.   
26.     vtkVectorNorm *gradiant = vtkVectorNorm::New();  
27.     gradiant->SetInputConnection(iso->GetOutputPort());  
28.   
29.     vtkDataSetMapper *cubeMapper = vtkDataSetMapper::New();  
30.     cubeMapper->SetInputConnection(gradiant->GetOutputPort());  
31.     vtkActor *cubeActor = vtkActor::New();  
32.     cubeActor->SetMapper(cubeMapper);  
33.   
34.     vtkCamera *camera = vtkCamera::New();  
35.     camera->SetPosition(1,1,1);  
36.     camera->SetFocalPoint(0,0,0);  
37.   
38.     vtkRenderer *renderer = vtkRenderer::New();  
39.     vtkRenderWindow *renWin = vtkRenderWindow::New();  
40.     renWin->AddRenderer(renderer);  
41.   
42.     vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();  
43.     iren->SetRenderWindow(renWin);  
44.   
45.     renderer->AddActor(cubeActor);  
46.     renderer->SetActiveCamera(camera);  
47.     renderer->ResetCamera();  
48.     renderer->SetBackground(1,1,1);  
49.   
50.     renWin->SetSize(512,512);  
51.     renWin->Render();  
52.     iren->Start();  
53.   
54.     image->Delete();  
55.     cubeMapper->Delete();  
56.     cubeActor->Delete();  
57.     camera->Delete();  
58.     renderer->Delete();  
59.     renWin->Delete();  
60.     iren->Delete();