VTK 从DICOM数据中重建三维

魔法师LQ

VTK读入DICOM数据

图像来自www.dicomlibrarymedDream

DICOM等医学影像本来将图像存储为多维的数组,因此从中恢复成3维的模型,而不是按照层次切片来查看是很自然的需求。

参考

[1] https://vtk.org/Wiki/VTK/Examples/Cxx/IO/ReadDICOMSeries

[2] https://lorensen.github.io/VTKExamples/site/Cxx/Medical/MedicalDemo4/

附录

https://lorensen.github.io/VTKExamples/site/Java/Medical/MedicalDemo4/

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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
import vtk.vtkNamedColors;
import vtk.vtkNativeLibrary;
import vtk.vtkRenderWindow;
import vtk.vtkRenderWindowInteractor;
import vtk.vtkRenderer;
import vtk.vtkMetaImageReader;
import vtk.vtkCamera;
import vtk.vtkColorTransferFunction;
import vtk.vtkPiecewiseFunction;
import vtk.vtkFixedPointVolumeRayCastMapper;
import vtk.vtkVolumeProperty;
import vtk.vtkVolume;

public class MedicalDemo4
{
// -----------------------------------------------------------------
// Load VTK library and print which library was not properly loaded
static
{
if (!vtkNativeLibrary.LoadAllNativeLibraries())
{
for (vtkNativeLibrary lib : vtkNativeLibrary.values())
{
if (!lib.IsLoaded())
{
System.out.println(lib.GetLibraryName() + " not loaded");
}
}
}
vtkNativeLibrary.DisableOutputWindow(null);
}
// -----------------------------------------------------------------


public static void main(String args[])
{

//parse command line arguments
if (args.length != 1)
{
System.err.println("Usage: java -classpath ... Filename(.mhd) e.g FullHead.mhd");
return;
}
String inputFilename = args[0];

vtkNamedColors colors = new vtkNamedColors();

double Bgcolor[] = new double[4];

colors.GetColor("SteelBlue", Bgcolor);

// Create the renderer, render window and interactor.
vtkRenderer ren = new vtkRenderer();
vtkRenderWindow renWin = new vtkRenderWindow();
renWin.AddRenderer(ren);
vtkRenderWindowInteractor iren = new vtkRenderWindowInteractor();
iren.SetRenderWindow(renWin);

// The following reader is used to read a series of 2D slices (images)
// that compose the volume. The slice dimensions are set, and the
// pixel spacing. The data Endianness must also be specified. The reader
// uses the FilePrefix in combination with the slice number to construct
// filenames using the format FilePrefix.%d. (In this case the FilePrefix
// is the root name of the file: quarter.)

vtkMetaImageReader reader = new vtkMetaImageReader();
reader.SetFileName(inputFilename);

// The volume will be displayed by ray-cast alpha compositing.
// A ray-cast mapper is needed to do the ray-casting.
vtkFixedPointVolumeRayCastMapper volumeMapper = new vtkFixedPointVolumeRayCastMapper();
volumeMapper.SetInputConnection(reader.GetOutputPort());

// The color transfer function maps voxel intensities to colors.
// It is modality-specific, and often anatomy-specific as well.
// The goal is to one color for flesh (between 500 and 1000)
// and another color for bone (1150 and over).
vtkColorTransferFunction volumeColor = new vtkColorTransferFunction();
volumeColor.AddRGBPoint(0, 0.0, 0.0, 0.0);
volumeColor.AddRGBPoint(500, 1.0, 0.5, 0.3);
volumeColor.AddRGBPoint(1000, 1.0, 0.5, 0.3);
volumeColor.AddRGBPoint(1150, 1.0, 1.0, 0.9);

// The opacity transfer function is used to control the opacity
// of different tissue types.
vtkPiecewiseFunction volumeScalarOpacity = new vtkPiecewiseFunction();
volumeScalarOpacity.AddPoint(0, 0.00);
volumeScalarOpacity.AddPoint(500, 0.15);
volumeScalarOpacity.AddPoint(1000, 0.15);
volumeScalarOpacity.AddPoint(1150, 0.85);

// The gradient opacity function is used to decrease the opacity
// in the "flat" regions of the volume while maintaining the opacity
// at the boundaries between tissue types. The gradient is measured
// as the amount by which the intensity changes over unit distance.
// For most medical data, the unit distance is 1mm.
vtkPiecewiseFunction volumeGradientOpacity = new vtkPiecewiseFunction();
volumeGradientOpacity.AddPoint(0, 0.0);
volumeGradientOpacity.AddPoint(90, 0.5);
volumeGradientOpacity.AddPoint(100, 1.0);

// The VolumeProperty attaches the color and opacity functions to the
// volume, and sets other volume properties. The interpolation should
// be set to linear to do a high-quality rendering. The ShadeOn option
// turns on directional lighting, which will usually enhance the
// appearance of the volume and make it look more "3D". However,
// the quality of the shading depends on how accurately the gradient
// of the volume can be calculated, and for noisy data the gradient
// estimation will be very poor. The impact of the shading can be
// decreased by increasing the Ambient coefficient while decreasing
// the Diffuse and Specular coefficient. To increase the impact
// of shading, decrease the Ambient and increase the Diffuse and Specular.
vtkVolumeProperty volumeProperty = new vtkVolumeProperty();
volumeProperty.SetColor(volumeColor);
volumeProperty.SetScalarOpacity(volumeScalarOpacity);
volumeProperty.SetGradientOpacity(volumeGradientOpacity);
volumeProperty.SetInterpolationTypeToLinear();
volumeProperty.ShadeOn();
volumeProperty.SetAmbient(0.4);
volumeProperty.SetDiffuse(0.6);
volumeProperty.SetSpecular(0.2);

// The vtkVolume is a vtkProp3D (like a vtkActor) and controls the position
// and orientation of the volume in world coordinates.
vtkVolume volume = new vtkVolume();
volume.SetMapper(volumeMapper);
volume.SetProperty(volumeProperty);
double c[] = new double[3];
c=volume.GetCenter();

ren.AddViewProp(volume);

// Set up an initial view of the volume. The focal point will be the
// center of the volume, and the camera position will be 400mm to the
// patient's left (which is our right).

vtkCamera camera = new vtkCamera();
camera.SetViewUp (0, 0, -1);
camera.SetPosition (c[0], c[1] - 400, c[2]);
camera.SetFocalPoint (c[0], c[1], c[2]);
camera.Azimuth(30.0);
camera.Elevation(30.0);
camera.Dolly(0.75);

// Actors are added to the renderer. An initial camera view is created.
// The Dolly() method moves the camera towards the FocalPoint, thereby enlarging the image.

ren.SetActiveCamera(camera);

// Set a background color for the renderer and set the size of the
// render window (expressed in pixels).
ren.SetBackground(Bgcolor);

// Note that when camera movement occurs (as it does in the Dolly()
// method), the clipping planes often need adjusting. Clipping planes
// consist of two planes: near and far along the view direction. The
// near plane clips out objects in front of the plane; the far plane
// clips out objects behind the plane. This way only what is drawn
// between the planes is actually rendered.
ren.ResetCameraClippingRange ();

renWin.SetSize(300, 300);
renWin.Render();

iren.Initialize();
iren.Start();
}
}