Marching Cubes算法再回顾

2024-01-08 10:37:25
1,确定包含等值面的体元

首先介绍一下?体元的概念,体元是三维图像中由相邻的八个体素点组成的正方体方格,英语也叫?Cube,体元中角点函数值分为两种情况,一种是大于等于给定等值面的值?C0?,则将角点设为?1?称该角点在等值面内部,否则设为0,在等值面之外,

一般来说,会出现一个角点在内,一个角点在外,则角点之间的连线(也就是体元的边)必然与等值面相交,根据这个原理就能判断等值面与哪些体元相交。

——————————————————————————————————————

三维空间中,平行且相邻的两个二维图像(每个图像中的正方形四个像素顶点组成一个基本的像素图像单元)组成一个基本的三维图像单元。,下图中由6个这样的基本三维图像单元:

vtkImageData结构由尺寸、间距和原点来定义。尺寸标注是沿着每个主轴的体素或像素的数量。原点是数据的第一个切片的左下角的世界坐标位置。间距是沿三个主要轴的像素之间的距离。

原点是数据集左下角的世界坐标位置。

尺寸是沿着三个主要轴的体素或像素的数量。

间距是体素的高度、长度和宽度,或相邻像素之间的距离,这取决于是将数据视为相同的方框还是连续函数中的样本点。
——————————————————————————————————————————

Marching Cubes算法根据一个立方体的8个顶点,判断这8个顶点的每个顶点在等值面的内部还是外部(每个顶点只有“在等值面内”和“在等值面外”这两种状态,设为0和1)从而根据这8个顶点的状态,建立一个包含共256种状态的查找表(根据平面对称性、中心对称性,256种最终降到15种)。

顶点值高于等值在表面的内部,等于等值在表面上,低于等值在表面外。

体元的每个顶点有两种状态,总共有256种,可以制作一个查找表(look up table)

但由于反转状态不变,所以可以减少一半,为128种。

再根据旋转不变形,又可以减少到14种情况。

可以认为这14中类似于基,经过旋转,反转可以得到256种状态对应的结果Triangulated Cubes:

根据每个顶点的状态,我们可以为每类制作一个8位索引Cube Numbering:

索引指向边表,给出了边的交叉情况。

相交边的编码——通过编码记录对应的cube,相交边的编号

二进制:00000010

十进制:2

Table[256]表示哪些边有交点

Table[2]=0x103=0000 0010 0000 0011

表示0,1,9号边上有交点

为了避免每次转化成二进制进行解码,可以直接记录与哪些边有交点,之后直接查表即可。

Table2[256][16]

Table2[2]=(0, 1, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

因为每个正方体中最多有1个或4个三角形,所以线性插值就足够了。

2,确定等值面与体元边界的交点

找到含有等值面的体元之后,接下来就是确定等值面与体元边界的交点,体元间的数值都是呈线性变化,求交点时一般采用的是线性插值,如?Case0?中等值面的两个端点 一个在外为( 标记0) ,一个在内 ( 标记为1?) 则交点为0.5;

3,求等值面的法向量

以上步骤?1,2,3?为实现?MC?算法步骤流程,但利用 VTK ,不需要这么繁琐,主要算法步骤都已经封装到?vtkMarchingCube?类中,使用?vtkMarchingCube?时,需要设置三个参数:

  • SetValue(int i,double value)?设置第i 个等值面的值b,(提醒一下,医学图像中的灰度值范围不是?0-256?而是0-65326,但大部分取值范围都在0-1000)。
  • SetNumberofContours(int number),设置等值面的个数
  • ComputerNormalsOn()?设置计算等值面的法向量,提高渲染质量;

dab9fdfd6e307b681d5c9054c536065e.png

上面这张图显示的就是 vtk 呈像的基本流程,下面是仿照官网写的用面绘制来对图像重建的代码部分:

#include<vtkRenderWindow.h>
#include<vtkRenderWindowInteractor.h>
#include<vtkDICOMImageReader.h>
#include<vtkMarchingCubes.h>
#include<vtkPolyDataMapper.h>
#include<vtkStripper.h>
#include<vtkActor.h>
#include<vtkProperty.h>
#include<vtkCamera.h>
#include<vtkOutlineFilter.h>
#include<vtkOBJExporter.h>
#include<vtkRenderer.h>
#include<vtkMetaImageReader.h>
#include<vtkInteractorStyleTrackballCamera.h>
 
 
#include<iostream>
#include<string.h>
//需要进行初始化,否则会报错
#include <vtkAutoInit.h> 
#include<vtkRenderingVolumeOpenGL2ObjectFactory.h>
#include<vtkRenderingOpenGL2ObjectFactory.h>
 
using namespace std;
int main()
{
    ///Marching Cube; 
 
    vtkObjectFactory::RegisterFactory(vtkRenderingOpenGL2ObjectFactory::New());
    vtkObjectFactory::RegisterFactory(vtkRenderingVolumeOpenGL2ObjectFactory::New());
 
    vtkSmartPointer<vtkRenderer> ren = vtkSmartPointer<vtkRenderer>::New();
    vtkSmartPointer<vtkRenderWindow> renWin = vtkSmartPointer<vtkRenderWindow>::New();//WINDOW;
 
    renWin->AddRenderer(ren);
 
    vtkSmartPointer<vtkRenderWindowInteractor> iren = vtkSmartPointer<vtkRenderWindowInteractor>::New();//wininteratcor;
    iren->SetRenderWindow(renWin);
 
    vtkSmartPointer<vtkDICOMImageReader> reader = vtkSmartPointer<vtkDICOMImageReader>::New();
    reader->SetDirectoryName("E:/DIcom_Data/DICOM");
    reader->SetDataByteOrderToLittleEndian();
    reader->Update();
 
 
    /*vtkDICOMImageReader *reader = vtkDICOMImageReader::New();
    reader->SetDirectoryName("E:/Coding Pra/VTK/VTK_Examples_StandardFormats_Input_DicomTestImages/DICOM");
    reader->SetDataByteOrderToLittleEndian();
    reader->Update();*/
 
    cout << "读取数据完毕" << endl;
    cout << "The width is" << reader->GetWidth() << endl;
    cout << "The height is" << reader->GetHeight() << endl;
    cout << "The depth is" << reader->GetPixelSpacing() << endl;
    cout << "The Output port is" << reader->GetOutputPort() << endl;
 
    
    vtkSmartPointer<vtkMarchingCubes> marchingcube = vtkSmartPointer<vtkMarchingCubes>::New();
    marchingcube->SetInputConnection(reader->GetOutputPort());//获得读取的数据的点集;
    marchingcube->SetValue(0, 200);//Setting the threshold;
    marchingcube->ComputeNormalsOn();//计算表面法向量;
 
    vtkSmartPointer<vtkStripper> Stripper = vtkSmartPointer<vtkStripper>::New();
    Stripper->SetInputConnection(marchingcube->GetOutputPort());//获取三角片
 
    vtkSmartPointer<vtkPolyDataMapper> Mapper = vtkSmartPointer<vtkPolyDataMapper>::New();//将三角片映射为几何数据;
    Mapper->SetInputConnection(Stripper->GetOutputPort());
    Mapper->ScalarVisibilityOff();//
 
 
    vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();//Created a actor;
    actor->SetMapper(Mapper);//获得皮肤几何数据
    actor->GetProperty()->SetDiffuseColor(1, .49, .25);//设置皮肤颜色;
    actor->GetProperty()->SetSpecular(0.3);//反射率;
    actor->GetProperty()->SetOpacity(1.0);//透明度;
    actor->GetProperty()->SetSpecularPower(20);//反射光强度;
    actor->GetProperty()->SetColor(1, 0, 0);//设置角的颜色;
    actor->GetProperty()->SetRepresentationToWireframe();//线框;
 
    //vtkSmartPointer<vtkCamera> camera = vtkSmartPointer<vtkCamera>::New();//Setting the Camera;
    //camera->SetViewUp(0, 0, -1);//设置相机向上方向;
    //camera->SetPosition(0, 1, 0);//位置:世界坐标系,相机位置;
    //camera->SetFocalPoint(0, 0, 0);//焦点,世界坐标系,控制相机方向;
    //camera->ComputeViewPlaneNormal();//重置视平面方向,基于当前的位置和焦点;
 
    vtkSmartPointer<vtkOutlineFilter> outfilterline = vtkSmartPointer<vtkOutlineFilter>::New();
    outfilterline->SetInputConnection(reader->GetOutputPort());
    vtkSmartPointer<vtkPolyDataMapper> outmapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    outmapper->SetInputConnection(outfilterline->GetOutputPort());
    vtkSmartPointer<vtkActor> OutlineActor = vtkSmartPointer<vtkActor>::New();
    OutlineActor->SetMapper(outmapper);
    OutlineActor->GetProperty()->SetColor(0, 0, 0);//线框颜色
 
    ren->AddActor(actor);
    ren->AddActor(OutlineActor);
    //ren->SetActiveCamera(camera);//设置渲染器的相机;
    ren->ResetCamera();
    ren->ResetCameraClippingRange();
 
    //camera->Dolly(1.5);//使用Dolly()方法延伸着视平面法向移动相机;
    ren->SetBackground(1, 1, 1);//设置背景颜色;
    renWin->SetSize(1000, 600);
 
 
    vtkInteractorStyleTrackballCamera *style = vtkInteractorStyleTrackballCamera::New();
    iren->SetInteractorStyle(style);
 
    renWin->Render();
    iren->Initialize();
    iren->Start();
 
    vtkSmartPointer<vtkOBJExporter> porter = vtkSmartPointer<vtkOBJExporter>::New();
    porter->SetFilePrefix("E:/ceshi/aaa/regist_after/polywrite.obj");//重建图像输出
    porter->SetInput(renWin);
    porter->Write();
 
 
    return EXIT_SUCCESS;
}

33a9534cf4e61f81f4d5a1ff644aa7f5.png

上面就是 VTK 基于?Marching Cube算法实现的重建效果:

体绘制重建

体绘制时分为两部分:

1,定义 vtkVoluemRayCastMapper 对象

体绘制中最常用的方法 ;vtkVolumeRayCastMapper()?光线投影,体绘制时,首先定义一个Mapper?然后接受两个输入:

  • SetInput(vtkImageDate *)?用于设置输入图像数据;
  • SetVolumeRayCastFunction(vtkVolumeRayCastFunction *)?用于设置光线投影函数类型;
2,利用 vtkVolumeProperty 定义体绘制属性;
  • SetScalarOpacity()?设置灰度不透明函数;
  • SetColor()?颜色传输函数;
3, 定义 vtkVolume 对象接收 Mapper对象和 Property 对象
  • SetMapper()接受 Mapper 对象;
  • SetProperty()?接受 Property 对象;

vtk 中体绘制 核心就是改变?Mapper?和?vtkVolumeRayCastFunction()?,上面中vtkColumeRayCastMapper?只是?VolumeMapper?其中的一种,且投影函数类?vtkVolumeRayCastFunction?一共有三个子类:

  • vtkVolumeRayCastCompositeFunction
  • vtkVolumeRayCasMIPFunction、
  • vtkVolumeRayCastIsosurfaceFunction
  • 因此,其细分的话vtk中的体绘制也不止一种

而下面这个是最常用到的(`vtkVolumeRayCastMapper?+?vtkVolumeRayCastCompositeFunction

//体绘制
 
#include<vtkRenderWindowInteractor.h>
#include<vtkDICOMImageReader.h>
#include<vtkCamera.h>
#include<vtkActor.h>
#include<vtkRenderer.h>
#include<vtkVolumeProperty.h>
#include<vtkProperty.h>
#include<vtkPolyDataNormals.h>
#include<vtkImageShiftScale.h>
#include "vtkVolumeRayCastMapper.h"
#include<vtkPiecewiseFunction.h>
#include<vtkColorTransferFunction.h>
#include<vtkVolumeRayCastCompositeFunction.h>
#include<vtkRenderWindow.h>
#include<vtkImageCast.h>
#include<vtkVolumeRayCastCompositeFunction.h>
#include<vtkOBJExporter.h>
#include<vtkOutlineFilter.h>
#include<vtkPolyDataMapper.h>
 
 
 
#include<vtkInteractorStyleTrackballCamera.h>
#include<vtkRenderingVolumeOpenGL2ObjectFactory.h>
#include<vtkRenderingOpenGL2ObjectFactory.h>
#include<vtkMetaImageReader.h>
 
#include<vtkLODProp3D.h>
 
 
//体绘制加速
 
//Gpu光照映射
#include<vtkGPUVolumeRayCastMapper.h>
 
#include<iostream>
 
int main()
{
 
    vtkObjectFactory::RegisterFactory(vtkRenderingOpenGL2ObjectFactory::New());
    vtkObjectFactory::RegisterFactory(vtkRenderingVolumeOpenGL2ObjectFactory::New());
 
 
    //定义绘制器;
    vtkRenderer *aRenderer = vtkRenderer::New();//指向指针;
    vtkSmartPointer<vtkRenderWindow> renWin = vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(aRenderer);
 
    vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
    iren->SetRenderWindow(renWin);
 
    //读取数据;
    /*vtkDICOMImageReader *reader = vtkDICOMImageReader::New();
    reader->SetDirectoryName("E:/Coding Pra/VTK/VTK_Examples_StandardFormats_Input_DicomTestImages/DICOM");
    reader->SetDataByteOrderToLittleEndian();*/
 
 
    vtkSmartPointer<vtkDICOMImageReader> reader = vtkSmartPointer<vtkDICOMImageReader>::New();
    reader->SetDirectoryName("E:/DIcom_Data/DICOM");
    reader->SetDataByteOrderToLittleEndian();
 
 
 
 
    //图像数据预处理,类型转换:通过 vtkimageCast 将不同类型数据集转化为 vtk 可以处理的数据集;
    vtkImageCast *cast_file = vtkImageCast::New();
    cast_file->SetInputConnection(reader->GetOutputPort());
    cast_file->SetOutputScalarTypeToUnsignedShort();
    cast_file->Update();
 
 
    //透明度映射函数定义;
    vtkPiecewiseFunction *opacityTransform = vtkPiecewiseFunction::New();
    opacityTransform->AddPoint(0, 0.0);
    opacityTransform->AddPoint(20, 0.0);
    opacityTransform->AddPoint(200, 1.0);
    opacityTransform->AddPoint(300, 1.0);
 
    //颜色映射函数定义,梯度上升的
    vtkColorTransferFunction *colorTransformFunction = vtkColorTransferFunction::New();
    colorTransformFunction->AddRGBPoint(0.0, 0.0, 0.0, 0.0);
    colorTransformFunction->AddRGBPoint(64.0, 0.0, 0.0, 0.0);
    colorTransformFunction->AddRGBPoint(128.0, 1.0, 0.0, 0.0);
    colorTransformFunction->AddRGBPoint(192.0, 1.0, 0.0, 0.0);
    colorTransformFunction->AddRGBPoint(255.0, 1.0, 0.0, 0.0);
 
    vtkPiecewiseFunction *gradientTransform = vtkPiecewiseFunction::New();
    gradientTransform->AddPoint(0, 0.0);
 
    gradientTransform->AddPoint(20, 2.0);
    gradientTransform->AddPoint(200, 0.1);
    gradientTransform->AddPoint(300, 0.1);
 
 
 
    //体数据属性;
    vtkVolumeProperty *volumeProperty = vtkVolumeProperty::New();
    volumeProperty->SetColor(colorTransformFunction);
    volumeProperty->SetScalarOpacity(opacityTransform);
    volumeProperty->SetGradientOpacity(gradientTransform);
    volumeProperty->ShadeOn();//应用
    volumeProperty->SetInterpolationTypeToLinear();//直线间样条插值;
    volumeProperty->SetAmbient(0.4);//环境光系数;
    volumeProperty->SetDiffuse(0.6);//漫反射;
    volumeProperty->SetSpecular(0.2);
    volumeProperty->SetSpecularPower(10);//高光强度;
 
 
    计算光照效应;利用 vtkBolumeRayCaseMapper进行计算;
    //vtkVolumeRayCastMapper *volunemapper = vtkVolumeRayCastMapper::New();
    //vtkVolumeRayCastCompositeFunction *compositeFunction = vtkVolumeRayCastCompositeFunction::New();
 
 
    //光纤映射类型定义:
    vtkSmartPointer<vtkVolumeRayCastCompositeFunction> compositecast =
        vtkSmartPointer<vtkVolumeRayCastCompositeFunction>::New();
 
    //Mapper定义,
    vtkSmartPointer<vtkVolumeRayCastMapper> hiresMapper = 
        vtkSmartPointer<vtkVolumeRayCastMapper>::New();
    hiresMapper->SetInputData(cast_file->GetOutput());
    hiresMapper->SetVolumeRayCastFunction(compositecast);
 
 
    vtkSmartPointer<vtkLODProp3D> prop = vtkSmartPointer<vtkLODProp3D>::New();
    prop->AddLOD(hiresMapper,volumeProperty,0.0);
 
    //
    //volunemapper->SetVolumeRayCastFunction(compositeFunction);//载入体绘制方法;
    //volunemapper->SetInputConnection(cast_file->GetOutputPort());
 
    //vtkFixedPointVolumeRayCastMapper *fixedPointVolumeMapper = vtkFixedPointVolumeRayCastMapper::New()
    //fixedPointVolumeMapper->SetInput()
 
 
 
    vtkVolume *volume = vtkVolume::New();
    volume->SetMapper(hiresMapper);
    volume->SetProperty(volumeProperty);//设置体属性;
 
    double volumeView[4] = { 0,0,0.5,1 };
 
    vtkOutlineFilter *outlineData = vtkOutlineFilter::New();//线框;
    outlineData->SetInputConnection(reader->GetOutputPort());
    vtkPolyDataMapper *mapOutline = vtkPolyDataMapper::New();
    mapOutline->SetInputConnection(outlineData->GetOutputPort());
    vtkActor *outline = vtkActor::New();
    outline->SetMapper(mapOutline);
    outline->GetProperty()->SetColor(0, 0, 0);//背景纯黑色;
 
    aRenderer->AddVolume(volume);
    aRenderer->AddActor(outline);
    aRenderer->SetBackground(1, 1, 1);
    aRenderer->ResetCamera();
 
 
    //重设相机的剪切范围;
    aRenderer->ResetCameraClippingRange();
    renWin->SetSize(800, 800);
    renWin->SetWindowName("测试");
 
    vtkRenderWindowInteractor *iren2 = vtkRenderWindowInteractor::New();
    iren2->SetRenderWindow(renWin);
 
    //设置相机跟踪模式
    vtkInteractorStyleTrackballCamera *style = vtkInteractorStyleTrackballCamera::New();
    iren2->SetInteractorStyle(style);
 
    renWin->Render();
    iren2->Initialize();
 
    iren2->Start();
 
    vtkOBJExporter *porter = vtkOBJExporter::New();
    porter->SetFilePrefix("E:/ceshi/aaa/regist_after/esho.obj");
    porter->SetInput(renWin);
    porter->Write();
    porter->Update();
 
 
    return EXIT_SUCCESS;
 
}

0c8404965dbab9c8abcab0cd04762657.png

上面是体绘制的结果,相对来说体绘制需要计算资源更大些, vtk 在这方面有所考虑,提供了vtKGPUVolumeRayCastMapper?GUP 加速的光线投射算法。

文章来源:https://blog.csdn.net/weixin_55035144/article/details/135449958
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。