使用VTK调用Marching Cubes和Volume Rendering

发表于   |   分类于   |     |   阅读次数   |  

前一段时间在可视化的实验课上，我第一次接触了VTK (The Visualization Toolkit)这一开源项目。在本篇文章中，将给出如何通过VTK来调用 Marching Cubes 和 Volume Rendering的方法。要说明的一点是，以下内容是基于 Ubuntu 16.04 的平台上运行的，其中VTK的版本为最新的7.0.0。至于 Ubuntu 下 VTK的编译安装将不属于本文的讨论范围，有兴趣的同学可以参考网上的相关教程进行环境搭建。

调用 Marching Cubes

Marching Cubes.cxx 文件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
#include <vtkVersion.h>
#include <vtkSmartPointer.h>
#include <vtkMarchingCubes.h>
#include <vtkVoxelModeller.h>
#include <vtkSphereSource.h>
#include <vtkImageData.h>
#include <vtkDICOMImageReader.h>
#include <vtkActor.h>
#include <vtkPolyDataMapper.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
int main(int argc, char *argv[]) {
    vtkSmartPointer <vtkImageData> volume =
            vtkSmartPointer<vtkImageData>::New();
    double isoValue;
    if (argc < 3) {
        vtkSmartPointer <vtkSphereSource> sphereSource =
                vtkSmartPointer<vtkSphereSource>::New();
        sphereSource->SetPhiResolution(20);
        sphereSource->SetThetaResolution(20);
        sphereSource->Update();
        double bounds[6];
        sphereSource->GetOutput()->GetBounds(bounds);
        for (unsigned int i = 0; i < 6; i += 2) {
            double range = bounds[i + 1] - bounds[i];
            bounds[i] = bounds[i] - .1 * range;
            bounds[i + 1] = bounds[i + 1] + .1 * range;
        }
        vtkSmartPointer <vtkVoxelModeller> voxelModeller =
                vtkSmartPointer<vtkVoxelModeller>::New();
        voxelModeller->SetSampleDimensions(50, 50, 50);
        voxelModeller->SetModelBounds(bounds);
        voxelModeller->SetScalarTypeToFloat();
        voxelModeller->SetMaximumDistance(.1);
        voxelModeller->SetInputConnection(sphereSource->GetOutputPort());
        voxelModeller->Update();
        isoValue = 0.5;
        volume->DeepCopy(voxelModeller->GetOutput());
    } else {
        vtkSmartPointer <vtkDICOMImageReader> reader =
                vtkSmartPointer<vtkDICOMImageReader>::New();
        reader->SetDirectoryName(argv[1]);
        reader->Update();
        volume->DeepCopy(reader->GetOutput());
        isoValue = atof(argv[2]);
    }
    vtkSmartPointer <vtkMarchingCubes> surface =
            vtkSmartPointer<vtkMarchingCubes>::New();
#if VTK_MAJOR_VERSION <= 5
    surface->SetInput(volume);
#else
    surface->SetInputData(volume);
#endif
    surface->ComputeNormalsOn();
    surface->SetValue(0, isoValue);
    vtkSmartPointer <vtkRenderer> renderer =
            vtkSmartPointer<vtkRenderer>::New();
    renderer->SetBackground(.1, .2, .3);
    vtkSmartPointer <vtkRenderWindow> renderWindow =
            vtkSmartPointer<vtkRenderWindow>::New();
    renderWindow->AddRenderer(renderer);
    vtkSmartPointer <vtkRenderWindowInteractor> interactor =
            vtkSmartPointer<vtkRenderWindowInteractor>::New();
    interactor->SetRenderWindow(renderWindow);
    vtkSmartPointer <vtkPolyDataMapper> mapper =
            vtkSmartPointer<vtkPolyDataMapper>::New();
    mapper->SetInputConnection(surface->GetOutputPort());
    mapper->ScalarVisibilityOff();
    vtkSmartPointer <vtkActor> actor =
            vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(mapper);
    renderer->AddActor(actor);
    renderWindow->Render();
    interactor->Start();
    return EXIT_SUCCESS;
}

CMakeLists.txt文件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
cmake_minimum_required(VERSION 2.8)
PROJECT(MarchingCubes)
find_package(VTK REQUIRED)
include(${VTK_USE_FILE})
add_executable(MarchingCubes MACOSX_BUNDLE MarchingCubes)
if(VTK_LIBRARIES)
    target_link_libraries(MarchingCubes ${VTK_LIBRARIES})
else()
    target_link_libraries(MarchingCubes vtkHybrid vtkWidgets)
endif()

这是已经写好的完整项目 MarchingCubes。下载以后进行解压，并且：

1
2
$ cd MarchingCubes/build
$ cmake ..

生成项目并且运行

1
2
$ make
$ ./MarchingCubes

运行结果为：

如果使用外部 Marching Cubes Data Sets数据集中的 MarchingMan 数据文件，则将其解压至 build 文件夹中

1
$ ./MarchingCubes Marching\ Man 1000

此处是在运行时传入两个参数：一个是数据文件所在目录 Marching Man，另一个是 ISO 感光度的值1000。并且，需要注意的是，由于解压后的数据文件夹为 Marching Man ，中间含有空格，因此须用反斜杠 \ 进行转义。当然，你也可以将其重命名为 MarchingMan ，这样在运行传参的时候就不需要对空格进行转义了，即可写为：

1
$ ./MarchingCubes MarchingMan 1000

最后的运行结果为：

如果使用 Mummies 中的 Male 数据，则

1
$ ./MarchingCubes Mummies/Male 100

如果使用 Mummies 中的 Female 数据，则

1
$ ./MarchingCubes Mummies/Female 1000

如果使用 PlasticSkull 中的数据，则

1
$ ./MarchingCubes PlasticSkull 1000

调用 Volume Rendering

由于在 VTK 中 Volume Rendering 的调用方式有多种，这里我们选取其中的 SmartVolumeMapper 进行讲解。

MinIntensityRendering.cxx

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
#include <vtkSmartPointer.h>
#include <vtkCamera.h>
#include <vtkFiniteDifferenceGradientEstimator.h>
#include <vtkImageClip.h>
#include <vtkPiecewiseFunction.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkStructuredPoints.h>
#include <vtkStructuredPointsReader.h>
#include <vtkVolume.h>
#include <vtkVolumeProperty.h>
#include <vtkFixedPointVolumeRayCastMapper.h>
#include <vtkColorTransferFunction.h>
int main(int argc, char *argv[]) {
    if (argc < 2) {
        std::cerr << "Required arguments: vtkFile" << std::endl;
        return EXIT_FAILURE;
    }
    std::string filename = argv[1]; //  "/Data/ironProt.vtk";
    // Create the renderers, render window, and interactor
    vtkSmartPointer <vtkRenderWindow> renWin =
            vtkSmartPointer<vtkRenderWindow>::New();
    vtkSmartPointer <vtkRenderWindowInteractor> iren =
            vtkSmartPointer<vtkRenderWindowInteractor>::New();
    iren->SetRenderWindow(renWin);
    vtkSmartPointer <vtkRenderer> ren =
            vtkSmartPointer<vtkRenderer>::New();
    renWin->AddRenderer(ren);
    // Read the data from a vtk file
    vtkSmartPointer <vtkStructuredPointsReader> reader =
            vtkSmartPointer<vtkStructuredPointsReader>::New();
    reader->SetFileName(filename.c_str());
    reader->Update();
    // Create a transfer function mapping scalar value to opacity
    vtkSmartPointer <vtkPiecewiseFunction> oTFun =
            vtkSmartPointer<vtkPiecewiseFunction>::New();
    oTFun->AddSegment(0, 1.0, 256, 0.1);
    vtkSmartPointer <vtkColorTransferFunction> cTFun =
            vtkSmartPointer<vtkColorTransferFunction>::New();
    cTFun->AddRGBPoint(0, 1.0, 1.0, 1.0);
    cTFun->AddRGBPoint(255, 1.0, 1.0, 1.0);
    // Need to crop to actually see minimum intensity
    vtkSmartPointer <vtkImageClip> clip =
            vtkSmartPointer<vtkImageClip>::New();
    clip->SetInputConnection(reader->GetOutputPort());
    clip->SetOutputWholeExtent(0, 66, 0, 66, 30, 37);
    clip->ClipDataOn();
    vtkSmartPointer <vtkVolumeProperty> property =
            vtkSmartPointer<vtkVolumeProperty>::New();
    property->SetScalarOpacity(oTFun);
    property->SetColor(cTFun);
    property->SetInterpolationTypeToLinear();
    vtkSmartPointer <vtkFixedPointVolumeRayCastMapper> mapper =
            vtkSmartPointer<vtkFixedPointVolumeRayCastMapper>::New();
    mapper->SetBlendModeToMinimumIntensity();
    mapper->SetInputConnection(clip->GetOutputPort());
    vtkSmartPointer <vtkVolume> volume =
            vtkSmartPointer<vtkVolume>::New();
    volume->SetMapper(mapper);
    volume->SetProperty(property);
    ren->AddViewProp(volume);
    renWin->Render();
    iren->Start();
    return EXIT_SUCCESS;
}

CMakeLists.txt

1
2
3
4
5
6
7
8
9
10
11
12
13
14
cmake_minimum_required(VERSION 2.8)
 
PROJECT(MinIntensityRendering)
 
find_package(VTK REQUIRED)
include(${VTK_USE_FILE})
 
add_executable(MinIntensityRendering MACOSX_BUNDLE MinIntensityRendering)
 
if(VTK_LIBRARIES)
	target_link_libraries(MinIntensityRendering ${VTK_LIBRARIES})
else()
	target_link_libraries(MinIntensityRendering vtkHybrid vtkWidgets)
endif()

在这里，同样提供已经写好的 MinIntensityRendering 项目工程文件和运行所需要的 VTK 文件 ironProt.vtk 。下载解压，将 ironProt.vtk 放入 build 文件夹中，然后

1
2
3
4
$ cd MinIntensityRendering/build
$ cmake ..
$ make
$ ./MinIntensityRendering ironProt.vtk

运行结果：

坚持原创技术分享，您的支持将鼓励我继续创作！

热评文章