VTK笔记——医学图像等值面提取（Marching Cubes）

什么是等值面

引用维基百科的解释：

等值曲面是一种曲面。在空间里，假若，每一点都有一个设定的值。这值可能是压力、温度、速度、密度。那么，一个等值曲面所包含的每一个点，其设定值是一样的。换句话说，以三维空间为定义域的连续函数，其每一个水平集都是一个等值曲面。

应用计算机图形学，我们可以简易地显示出等值曲面的线框图或明暗图。在计算流体力学里，数据视觉化方法时常会用等值曲面来表示流体（液体或气体）流过物体时的瞬时状态。这是工程师研究发展新科技的一个利器。他们可以观察一个系统在任何时间的状态，从而发现其中奥秘。例如，等值曲面可以代表超音速飞行的单独震波。或者，我们可以制造几个等值曲面来代表，当空气流过飞机翅膀时，随着时间演变的一系列压力值。

面对着一大堆三维空间的数据，一个明智又受欢迎的选择，就是采用等值曲面为数据视觉化的主要形式。简单的多边形造型渲染的等值曲面，不需要用到很多的中央处理单元的资源，就能够迅速的计算出所要显示的图形。

在医学影像里，三维的电脑断层扫描用等值曲面来代表一个密度值区的部位。这样，我们可以将内部器官、骨头、等等，这些结构视觉化。

简单的来说，我们可以使用等值面技术可视化感兴趣的区域，比如皮肤，骨骼等。

等值面生成

等值面可视化的三种技术，体绘制，移动立方体（Marching Cubes）和部分立方体（Dividing Cubes）.
这篇笔记主要记录的是移动立方体（Marching Cubes）.

Marching Cubes的原理

Marching Cubes（移动立方体）方法是由W.E.Lorenson和H.E.Cline在1987年提出来的。由于这一方法原理简单，易于实现，目前已经得到了较为广泛的应用，成为三维数据等值面生成的经典算法，Marching Cubes算法又简称为MC算法。

在Marching Cubes方法中，假定原始数据是离散的三维空间规则数据，一个体元定义为由相邻层上的8个顶点组成的一个长方体。为了在三维数据中构造等值面，应先给定所求等值面的值，该方法的基本原理是逐个处理所有的体元，将体元各顶点处的值与给定的阈值进行比较，首先找出与等值面相交的体元，然后通过插值求等值面与体元棱边的交点，并将各交点连成三角形来构成等值面片，所有体元中的三角形集合就构成了等值面。由于这一方法是逐个处理所有的体元，因此被称为Marching Cubes方法。

关于Marching Cubes原理的详细解释。有兴趣可以参考Lorensen, W. E. and Cline, H. E., “Marching Cubes: A High Resolution 3D Surface Construction Algorithm,” Computer Graphics, vol. 21, no. 3, pp. 163-169, July 1987.

vtk实现类

在vtk中实现Marching Cubes的算法类有很多，如vtkMarchingCubes,vtkSliceCubes,vtkMarchingSquares,vtkImageMarchingCubes等,它们是针对特定数据集类型而定义的。
vtkMarchingCubes
vtkMarchingCubes用于提取体数据。

generate isosurface(s) from volume

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

vtkContourFilter
vtkContourFilter根据特定的数据集类型自动创建一个最快的子类方法。

generate isosurfaces/isolines from scalar values

vtkContourFilter is a filter that takes as input any dataset and generates on output isosurfaces and/or isolines. The exact form of the output depends upon the dimensionality of the input data. Data consisting of 3D cells will generate isosurfaces, data consisting of 2D cells will generate isolines, and data with 1D or 0D cells will generate isopoints. Combinations of output type are possible if the input dimension is mixed.

To use this filter you must specify one or more contour values. You can either use the method SetValue() to specify each contour value, or use GenerateValues() to generate a series of evenly spaced contours. It is also possible to accelerate the operation of this filter (at the cost of extra memory) by using a vtkScalarTree. A scalar tree is used to quickly locate cells that contain a contour surface. This is especially effective if multiple contours are being extracted. If you want to use a scalar tree, invoke the method UseScalarTreeOn().

下面示例一个生成头部数据皮肤和骨骼的方法。

The example uses FullHead.mhd which references FullHead.raw.gz.

MarchingCubes.cxx

#include "vtkSmartPointer.h"
#include "vtkMetaImageReader.h"
#include "vtkMarchingCubes.h"
#include "vtkContourFilter.h"
#include "vtkStripper.h"
#include "vtkPolyDataNormals.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkProperty.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkOutlineFilter.h"
#include "vtkNamedColors.h"

#include <array>

int main(int argc, char* argv[])
{
	if (argc < 2)
	{
		std::cout << "Usage: " << argv[0] << " FullHead.mhd" << std::endl;
		return EXIT_FAILURE;
	}

	auto colors = vtkSmartPointer<vtkNamedColors>::New();
	std::array<unsigned char, 4> skinColor{ {255, 125, 64} };
	colors->SetColor("SkinColor", skinColor.data());

	auto reader = vtkSmartPointer<vtkMetaImageReader>::New();
	reader->SetFileName(argv[1]);

#if 0
	// An isosurface, or contour value of 500 is known to correspond to the
	// skin of the patient.
	auto skinExtractor = vtkSmartPointer<vtkMarchingCubes>::New();
	skinExtractor->SetInputConnection(reader->GetOutputPort());
	skinExtractor->SetValue(0, 500);

	// An isosurface, or contour value of 1150 is known to correspond to the
	// bone of the patient.
	auto boneExtractor = vtkSmartPointer<vtkMarchingCubes>::New();
	boneExtractor->SetInputConnection(reader->GetOutputPort());
	boneExtractor->SetValue(0, 1150);

	// The triangle stripper is used to create triangle strips from the
	// isosurface; these render much faster on may systems.
	auto skinStripper = vtkSmartPointer<vtkStripper>::New();
	skinStripper->SetInputConnection(skinExtractor->GetOutputPort());

	auto boneStripper = vtkSmartPointer<vtkStripper>::New();
	boneStripper->SetInputConnection(boneExtractor->GetOutputPort());
#else
	// An isosurface, or contour value of 500 is known to correspond to the
	// skin of the patient.
	auto skinExtractor = vtkSmartPointer<vtkContourFilter>::New();
	skinExtractor->SetInputConnection(reader->GetOutputPort());
	skinExtractor->SetValue(0, 500);

	// An isosurface, or contour value of 1150 is known to correspond to the
	// bone of the patient.
	auto boneExtractor = vtkSmartPointer<vtkContourFilter>::New();
	boneExtractor->SetInputConnection(reader->GetOutputPort());
	boneExtractor->SetValue(0, 1150);

	auto skinNormals = vtkSmartPointer<vtkPolyDataNormals>::New();
	skinNormals->SetInputConnection(skinExtractor->GetOutputPort());
	skinNormals->SetFeatureAngle(60.0);

	auto boneNormals = vtkSmartPointer<vtkPolyDataNormals>::New();
	boneNormals->SetInputConnection(boneExtractor->GetOutputPort());
	boneNormals->SetFeatureAngle(60.0);

	// The triangle stripper is used to create triangle strips from the
	// isosurface; these render much faster on may systems.
	auto skinStripper = vtkSmartPointer<vtkStripper>::New();
	skinStripper->SetInputConnection(skinNormals->GetOutputPort());

	auto boneStripper = vtkSmartPointer<vtkStripper>::New();
	boneStripper->SetInputConnection(boneNormals->GetOutputPort());
#endif

	auto skinMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
	skinMapper->SetInputConnection(skinStripper->GetOutputPort());
	skinMapper->ScalarVisibilityOff();
	
	auto boneMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
	boneMapper->SetInputConnection(boneStripper->GetOutputPort());
	boneMapper->ScalarVisibilityOff();

	auto skin = vtkSmartPointer<vtkActor>::New();
	skin->SetMapper(skinMapper);
	skin->GetProperty()->SetDiffuseColor(colors->GetColor3d("SkinColor").GetData());
	skin->GetProperty()->SetSpecular(.3);
	skin->GetProperty()->SetSpecularPower(20);
	skin->GetProperty()->SetOpacity(.5);

	auto bone = vtkSmartPointer<vtkActor>::New();
	bone->SetMapper(boneMapper);
	bone->GetProperty()->SetDiffuseColor(colors->GetColor3d("Ivory").GetData());

	// An outline provides context around the data.
	auto outlineData = vtkSmartPointer<vtkOutlineFilter>::New();
	outlineData->SetInputConnection(reader->GetOutputPort());

	auto outlineMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
	outlineMapper->SetInputConnection(outlineData->GetOutputPort());

	auto outline = vtkSmartPointer<vtkActor>::New();
	outline->SetMapper(outlineMapper);
	outline->GetProperty()->SetColor(1, 0, 0);

	auto renderer = vtkSmartPointer<vtkRenderer>::New();
	renderer->AddActor(skin);
	renderer->AddActor(bone);
	renderer->AddActor(outline);

	auto renWindow = vtkSmartPointer<vtkRenderWindow>::New();
	renWindow->SetSize(640, 480);
	renWindow->AddRenderer(renderer);

	auto interactor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
	interactor->SetRenderWindow(renWindow);

	renWindow->Render();
	interactor->Initialize();
	interactor->Start();

	return EXIT_SUCCESS;
}

Ref

等值曲面
等值面
vtkMarchingCubes Class Reference
vtkContourFilter Class Reference
VTKExamples/Cxx/Medical/MedicalDemo2
An Implementation of the Marching Cubes[1] Algorithm