体绘制算法

三维绘制包括面绘制和体绘制。体绘制可以提供最接近人眼视觉的绘制效果。体绘制的方法包括：

RayCasting
Splatting
Shear Warp
Texture based volume rendering

其中RayCasting是最常用的方法之一。包括四个步骤：Ray casting, Sampling, Shading and Composition。具体的描述如下：



在RayCasting算法中，主要涉及到的是transfer function的概念。Transfer Function通常被用来对不同密度或者不同位置的体素单元加上颜色和透明度属性。Transfer Function的实例如下：



在VTK中，实现volume Rendering的程序示例如下：

#include <vtkSmartPointer.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkVolume16Reader.h>
#include <vtkVolume.h>
#include <vtkVolumeRayCastMapper.h>
#include <vtkGPUVolumeRayCastMapper.h>
#include <vtkVolumeRayCastCompositeFunction.h>
#include <vtkVolumeProperty.h>
#include <vtkColorTransferFunction.h>
#include <vtkPiecewiseFunction.h>
#include <vtkCamera.h>
#include <vtkDICOMImageReader.h>
#include <vtkImageShiftScale.h>

int main (int argc, char *argv[])
{

vtkSmartPointer< vtkDICOMImageReader > reader = vtkSmartPointer< vtkDICOMImageReader >::New();
reader->SetDirectoryName("E:\\Coding\\ReadDICOMSeries\\data\\DicomSeries");
reader->SetDataScalarTypeToUnsignedShort();
reader->SetDataSpacing(1, 1, 1);
reader->SetDataByteOrderToLittleEndian();

//图像数据预处理。常见两种操作。一种是类型转换，通过 vtkImageShiftScale 将不同类型的数据集转换为 VTK 可以处理的数据；
//另一种是剔除冗余数据，通过 vtkStripper 放置无效的旧单元的存在，提高绘制速度。
vtkImageShiftScale *shiftScale = vtkImageShiftScale::New();
shiftScale->SetInputConnection(reader->GetOutputPort());
shiftScale->SetOutputScalarTypeToUnsignedChar();

// The volume will be displayed by ray-cast alpha compositing.
// A ray-cast mapper is needed to do the ray-casting, and a
// compositing function is needed to do the compositing along the ray.
vtkSmartPointer<vtkVolumeRayCastCompositeFunction> rayCastFunction =
vtkSmartPointer<vtkVolumeRayCastCompositeFunction>::New();

//vtkSmartPointer<vtkVolumeRayCastMapper> volumeMapper = vtkSmartPointer<vtkVolumeRayCastMapper>::New();
vtkSmartPointer<vtkGPUVolumeRayCastMapper> volumeMapper = vtkSmartPointer<vtkGPUVolumeRayCastMapper>::New();
volumeMapper->SetInputConnection(shiftScale->GetOutputPort());
//volumeMapper->SetVolumeRayCastFunction(rayCastFunction);// 如果不用GPU绘制，就使用这句话
volumeMapper->SetBlendModeToComposite();
volumeMapper->SetSampleDistance(0.1);

// 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).
vtkSmartPointer<vtkColorTransferFunction>volumeColor =
vtkSmartPointer<vtkColorTransferFunction>::New();
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.
vtkSmartPointer<vtkPiecewiseFunction> volumeScalarOpacity =
vtkSmartPointer<vtkPiecewiseFunction>::New();
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.
vtkSmartPointer<vtkPiecewiseFunction> volumeGradientOpacity =
vtkSmartPointer<vtkPiecewiseFunction>::New();
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.
vtkSmartPointer<vtkVolumeProperty> volumeProperty =
vtkSmartPointer<vtkVolumeProperty>::New();
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.
vtkSmartPointer<vtkVolume> volume =
vtkSmartPointer<vtkVolume>::New();
volume->SetMapper(volumeMapper);
volume->SetProperty(volumeProperty);

// Create the renderer, the render window, and the interactor. The renderer
// draws into the render window, the interactor enables mouse- and
// keyboard-based interaction with the scene.
vtkSmartPointer<vtkRenderer> ren = vtkSmartPointer<vtkRenderer>::New();
vtkSmartPointer<vtkRenderWindow> renWin = vtkSmartPointer<vtkRenderWindow>::New();
if (NULL == renWin)
{
std::cout << "error" << std::endl;
//return 0;
}
renWin->AddRenderer(ren);
vtkSmartPointer<vtkRenderWindowInteractor> iren =
vtkSmartPointer<vtkRenderWindowInteractor>::New();
iren->SetRenderWindow(renWin);

// Finally, add the volume to the renderer
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 = ren->GetActiveCamera();
double *c = volume->GetCenter();
camera->SetFocalPoint(c[0], c[1], c[2]);
camera->SetPosition(c[0] + 400, c[1], c[2]);
camera->SetViewUp(0, 0, -1);

// Increase the size of the render window
renWin->SetSize(640, 480);

// Interact with the data.
iren->Initialize();
iren->Start();

return EXIT_SUCCESS;
}

结果如下：



参考文档：

http://www.mccauslandcenter.sc.edu/CRNL/tools/surface-rendering-with-matlab
https://graphics.ethz.ch/teaching/former/scivis_07/Notes/Handouts/03-raycasting.pdf
http://en.wikipedia.org/wiki/Volume_ray_casting
http://www-pequan.lip6.fr/~tierny/stuff/teaching/tierny_intro_vol_rend09.pdf
http://web.cs.wpi.edu/~matt/courses/cs563/talks/powwie/p1/ray-cast.htm
http://en.wikipedia.org/wiki/Volume_rendering#Volume_ray_casting
http://www.kitware.com/source/home/post/160