VTK随笔十一：VTK图形处理（表面重建、点云配准、纹理映射）

秀木易风

于 2024-08-31 07:30:00 发布

阅读量1.6k

点赞数 21

分类专栏： VTK随笔文章标签： c++ qt

本文链接：https://blog.csdn.net/u012959478/article/details/141643268

版权

VTK随笔专栏收录该内容

13 篇文章 3 订阅

订阅专栏

一、表面重建

通过三维扫描仪所获取的实际物体的空间点云数据仅仅表示物体的几何形状，而无法表达其内部的拓扑结构。拓扑结构对于实际图形处理以及可视化具有更重要的意义。因此这就需要利用表面重建技术将点云数据转换为面模型，通常为三角网格模型。

1、三角剖分

三角剖分将一些散乱的点云数据剖分为一系列的三角形网格。最常用的三角剖分技术是Delaunay三角剖分。Delaunay三角剖分具有许多优良的性质，如最大化最小角特性，即在所有可能的三角剖分中，其所生成的三角形的最小角的角度最大。所以，Delaunay 三角剖分无论从哪个区域开始构建，最终生成的三角网格是唯一的。
VTK的 vtkDelaunay2D 类实现了二维三角剖分。该类的输入数据为一个 vtkPointSet 或其子类表示的三维空间点集，其输出为一个三角网格 vtkPolyData数据。虽然输入的是三维数据，但是算法仅使用 XY平面数据进行平面三角剖分，而忽略Z方向数据。当然，也可以为vtkDelaunay2D 设置一个投影变换从而在新的投影平面上进行三角剖分。需要注意的是，在不加任何限制的条件下，该类生成的平面三角网格为一个凸包。

示例代码：

#include <QApplication>
#include <vtkSmartPointer.h>
#include <vtkRenderer.h>
#include <vtkGenericOpenGLRenderWindow.h>
#include <QVTKOpenGLNativeWidget.h>
#include <vtkDelaunay2D.h>
#include <vtkVertexGlyphFilter.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>

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

    unsigned int gridSize = 10;
    vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
    for(unsigned int x = 0; x < gridSize; x++)
    {
        for(unsigned int y = 0; y < gridSize; y++)
        {
            points->InsertNextPoint(x, y, vtkMath::Random(0.0, 3.0));
        }
    }

    vtkSmartPointer<vtkPolyData> polydata = vtkSmartPointer<vtkPolyData>::New();
    polydata->SetPoints(points);

    vtkSmartPointer<vtkDelaunay2D> delaunay = vtkSmartPointer<vtkDelaunay2D>::New();
    delaunay->SetInputData(polydata);
    delaunay->Update();

    vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter = vtkSmartPointer<vtkVertexGlyphFilter>::New();
    glyphFilter->SetInputData(polydata);
    glyphFilter->Update();

    vtkSmartPointer<vtkPolyDataMapper> pointsMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    pointsMapper->SetInputData(glyphFilter->GetOutput());

    vtkSmartPointer<vtkActor> pointsActor = vtkSmartPointer<vtkActor>::New();
    pointsActor->SetMapper(pointsMapper);
    pointsActor->GetProperty()->SetPointSize(3);
    pointsActor->GetProperty()->SetColor(1,0,0);

    vtkSmartPointer<vtkPolyDataMapper> triangulatedMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    triangulatedMapper->SetInputData(delaunay->GetOutput());

    vtkSmartPointer<vtkActor> triangulatedActor = vtkSmartPointer<vtkActor>::New();
    triangulatedActor->SetMapper(triangulatedMapper);

    vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
    renderer->AddActor(pointsActor);
    renderer->AddActor(triangulatedActor);
    renderer->SetBackground(1.0, 1.0, 1.0);

    vtkSmartPointer<vtkGenericOpenGLRenderWindow> renderWindow = vtkSmartPointer<vtkGenericOpenGLRenderWindow>::New();
    renderWindow->AddRenderer(renderer);

    QVTKOpenGLNativeWidget w;
    w.resize(640, 320);
    w.setRenderWindow(renderWindow);
    w.setWindowTitle("PolyDataDelaunay2D");
    w.show();

    return a.exec();
}

运行效果：

vtkDelaunay2D 还支持加入边界限制。用户需要设置另外一个 vtkPolyData 数据，其内部的线段、闭合或者非闭合的线段集合将作为边界条件控制三角剖分的过程。其中组成这些边界的点的索引必须与原始点集数据一致。加入边界条件后，最后的剖分结果可能不再满足Delaunay准则。

另外，VTK的 vtkDelaunay3D 类可实现三维三角剖分。该类的使用方法与vtkDelaunay2D基本一致，不同的是，三维三角剖分得到的结果并非三角网格，而是四面体网格。因此其输出数据的类型为 vtkUnstructuredGrid，在未加入边界条件下的三维三角剖分通常也为一个凸包。

2、等值面提取

等值面(线)提取是一种常用的可视化技术，常应用于医学、地质、气象等领域，例如在医学图像处理中，由于CT、MRI等图像分辨率越来越高，虽然体绘制技术可以清晰地对数据内部结构进行可视化，但是其计算量和效率却制约了其使用。此时可通过等值面提取技术，仅提取感兴趣的一个或者几个组织轮廓，并生成网格模型以供后续的处理和显示。

示例代码：

#include <QApplication>
#include <vtkSmartPointer.h>
#include <vtkRenderer.h>
#include <vtkGenericOpenGLRenderWindow.h>
#include <QVTKOpenGLNativeWidget.h>
#include <vtkMetaImageReader.h>
#include <vtkMarchingCubes.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>

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

    vtkSmartPointer<vtkMetaImageReader> reader = vtkSmartPointer<vtkMetaImageReader>::New();
    reader->SetFileName("D:/data/HeadMRVolume.mhd");
    reader->Update();

    double isoValue = 200;

    vtkSmartPointer<vtkMarchingCubes> surface = vtkSmartPointer<vtkMarchingCubes>::New();
    surface->SetInputData((vtkDataObject*)reader->GetOutput());
    surface->ComputeNormalsOn();
    surface->SetValue(0, isoValue);
    //surface->GenerateValues(5, 150,200);

    // vtkSmartPointer<vtkContourFilter> surface =
    // vtkSmartPointer<vtkContourFilter>::New();
    // surface->SetInput(reader->GetOutput());
    // surface->ComputeNormalsOn();
    // surface->SetValue(0, isoValue);

    vtkSmartPointer<vtkPolyDataMapper> surfMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    surfMapper->SetInputConnection(surface->GetOutputPort());
    surfMapper->ScalarVisibilityOff();

    vtkSmartPointer<vtkActor> surfActor = vtkSmartPointer<vtkActor>::New();
    surfActor->SetMapper(surfMapper);
    surfActor->GetProperty()->SetColor(1.0, 0.0, 0.0);

    vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
    renderer->AddActor(surfActor);
    renderer->SetBackground(1.0, 1.0, 1.0);

    vtkSmartPointer<vtkGenericOpenGLRenderWindow> renderWindow = vtkSmartPointer<vtkGenericOpenGLRenderWindow>::New();
    renderWindow->AddRenderer(renderer);

    QVTKOpenGLNativeWidget w;
    w.resize(640, 480);
    w.setRenderWindow(renderWindow);
    w.setWindowTitle("PolyDataMarchingCubes");
    w.show();

    return a.exec();
}

运行效果：

3、点云重建

虽然 Delaunay 三角剖分算法可以实现网格曲面重建，但是其应用主要在二维剖分，在三维空间网格生成中遇到了问题。因为在三维点云曲面重建中，Delaunay条件不再满足，不仅基于最大最小角判断的对角线交换准则不再成立，而且基于外接圆判据的Delaunay 三角化也不能保证网格质量。
vtkSurfaceReconstructionFilter 则实现了一种隐式曲面重建的方法，即将曲面看作一个符号距离函数的等值面，曲面内外的距离值的符号相反，而零等值面即为所求的曲面。该方法需要对点云数据进行网格划分，然后估算每个点的切平面和方向，并以每个点与最近的切平面距离来近似表面距离。这样即可得到一个符号距离的体数据，使用vkContourFilter 来提取零等值面即可得到相应的网格。

示例代码：

#include <QApplication>
#include <vtkSmartPointer.h>
#include <vtkRenderer.h>
#include <vtkGenericOpenGLRenderWindow.h>
#include <QVTKOpenGLNativeWidget.h>
#include <vtkPolyDataReader.h>
#include <vtkSurfaceReconstructionFilter.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkContourFilter.h>
#include <vtkVertexGlyphFilter.h>
#include <vtkCamera.h>

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

    vtkSmartPointer<vtkPolyDataReader> reader = vtkSmartPointer<vtkPolyDataReader>::New();
    reader->SetFileName("D:/data/fran_cut.vtk");
    reader->Update();

    vtkSmartPointer<vtkPolyData> points = vtkSmartPointer<vtkPolyData>::New();
    points->SetPoints(reader->GetOutput()->GetPoints());

    vtkSmartPointer<vtkSurfaceReconstructionFilter> surf = vtkSmartPointer<vtkSurfaceReconstructionFilter>::New();
    surf->SetInputData(points);
    surf->SetNeighborhoodSize(20);
    surf->SetSampleSpacing(0.005);
    surf->Update();

    vtkSmartPointer<vtkContourFilter> contour = vtkSmartPointer<vtkContourFilter>::New();
    contour->SetInputConnection(surf->GetOutputPort());
    contour->SetValue(0, 0.0);
    contour->Update();

    double leftViewport[4] = {0.0, 0.0, 0.5, 1.0};
    double rightViewport[4] = {0.5, 0.0, 1.0, 1.0};

    vtkSmartPointer<vtkVertexGlyphFilter> vertexGlyphFilter = vtkSmartPointer<vtkVertexGlyphFilter>::New();
    vertexGlyphFilter->AddInputData(points);
    vertexGlyphFilter->Update();

    vtkSmartPointer<vtkPolyDataMapper> vertexMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    vertexMapper->SetInputData(vertexGlyphFilter->GetOutput());
    vertexMapper->ScalarVisibilityOff();

    vtkSmartPointer<vtkActor> vertexActor = vtkSmartPointer<vtkActor>::New();
    vertexActor->SetMapper(vertexMapper);
    vertexActor->GetProperty()->SetColor(1.0, 0.0, 0.0);

    vtkSmartPointer<vtkRenderer> vertexRenderer = vtkSmartPointer<vtkRenderer>::New();
    vertexRenderer->AddActor(vertexActor);
    vertexRenderer->SetViewport(leftViewport);
    vertexRenderer->SetBackground(1.0, 1.0, 1.0);

    vtkSmartPointer<vtkPolyDataMapper> surfMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    surfMapper->SetInputData(contour->GetOutput());
    surfMapper->ScalarVisibilityOff();

    vtkSmartPointer<vtkActor> surfActor = vtkSmartPointer<vtkActor>::New();
    surfActor->SetMapper(surfMapper);
    surfActor->GetProperty()->SetColor(1.0, 0.0, 0.0);

    vtkSmartPointer<vtkRenderer> surfRenderer = vtkSmartPointer<vtkRenderer>::New();
    surfRenderer->AddActor(surfActor);
    surfRenderer->SetViewport(rightViewport);
    surfRenderer->SetBackground(1.0, 1.0, 1.0);

    vtkSmartPointer<vtkGenericOpenGLRenderWindow> renderWindow = vtkSmartPointer<vtkGenericOpenGLRenderWindow>::New();
    renderWindow->AddRenderer(vertexRenderer);
    renderWindow->AddRenderer(surfRenderer);

    QVTKOpenGLNativeWidget w;
    w.resize(640, 320);
    w.setRenderWindow(renderWindow);
    w.setWindowTitle("PolyDataSurfaceReconstruction");
    w.show();

    surfRenderer->GetActiveCamera()->SetPosition(0, -1, 0);
    surfRenderer->GetActiveCamera()->SetFocalPoint(0, 0, 0);
    surfRenderer->GetActiveCamera()->SetViewUp(0, 0, 1);
    surfRenderer->GetActiveCamera()->Azimuth(30);
    surfRenderer->GetActiveCamera()->Elevation(30);
    surfRenderer->ResetCamera();
    vertexRenderer->SetActiveCamera(surfRenderer->GetActiveCamera());

    return a.exec();
}

运行效果：

二、点云配准

在计算机逆向工程中，通过三维扫描等实物数字化技术可以获取各种点云数据。但是受测量环境和设备的影响，在一次测量的情况下，难以获取实物整体的点云数据，因此需要多次从不同角度进行测量。但不同的测量数据之间可能会存在平移错误和旋转错位等问题。这就需要使用点云配准技术来对测量点云数据进行局部配准和整合以得到完整的模型数据。另外，在外科手术导航技术中，图像标记点数据与人体表面标记点数据的配准是一个关键步骤对于手术定位的精度有着重要的影响。这通常需要使用基于标记点的配准技术(也属于点云配准的一种)。因此，点云配准即是对一组源点云数据应用一个空间变换，使得变换后的数据与目标点云数据能够一一映射，使两组数据之间的平均距离误差最小。

vtkLandmarkTransform 实现了标记点配准算法，使得两个点集在配准后的平均距离最小。通过 SetSourceLandmarks()和 SetTargetLandmarks()函数分别设置源标记点集和目标标记点集需要注意的是，源标记点集和目标标记点集序号要对应，如1号源点要映射到1号目标点。

点云数据配准最经典的方法是迭代最近点算法(IterativeClosest Points，ICP)。ICP 算法是一个迭代的过程，每次迭代中对于源数据点P找到目标点集O中的最近点，然后基于最小二乘原理求解当前的变换T。通过不断迭代直至收敛，即完成了点集的配准。

示例代码：

#include <QApplication>
#include <vtkSmartPointer.h>
#include <vtkRenderer.h>
#include <vtkGenericOpenGLRenderWindow.h>
#include <QVTKOpenGLNativeWidget.h>
#include <vtkPolyDataReader.h>
#include <vtkTransform.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkTransformPolyDataFilter.h>
#include <vtkVertexGlyphFilter.h>
#include <vtkIterativeClosestPointTransform.h>
#include <vtkLandmarkTransform.h>

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

    vtkSmartPointer<vtkPolyDataReader> reader =  vtkSmartPointer<vtkPolyDataReader>::New();
    reader->SetFileName("D:/data/fran_cut.vtk");
    reader->Update();
    vtkSmartPointer<vtkPolyData> original = reader->GetOutput();

    vtkSmartPointer<vtkTransform> translation = vtkSmartPointer<vtkTransform>::New();
    translation->Translate(0.2, 0.0, 0.0);
    translation->RotateX(30);

    vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter1 = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
    transformFilter1->SetInputData(reader->GetOutput());
    transformFilter1->SetTransform(translation);
    transformFilter1->Update();

    vtkSmartPointer<vtkPolyData> source = vtkSmartPointer<vtkPolyData>::New();
    source->SetPoints(original->GetPoints());

    vtkSmartPointer<vtkPolyData> target = vtkSmartPointer<vtkPolyData>::New();
    target->SetPoints(transformFilter1->GetOutput()->GetPoints());

    vtkSmartPointer<vtkVertexGlyphFilter> sourceGlyphFilter = vtkSmartPointer<vtkVertexGlyphFilter>::New();
    sourceGlyphFilter->SetInputData(source);
    sourceGlyphFilter->Update();

    vtkSmartPointer<vtkVertexGlyphFilter> targetGlyphFilter = vtkSmartPointer<vtkVertexGlyphFilter>::New();
    targetGlyphFilter->SetInputData(target);
    targetGlyphFilter->Update();

    vtkSmartPointer<vtkIterativeClosestPointTransform> icpTransform = vtkSmartPointer<vtkIterativeClosestPointTransform>::New();
    icpTransform->SetSource(sourceGlyphFilter->GetOutput());
    icpTransform->SetTarget(targetGlyphFilter->GetOutput());
    icpTransform->GetLandmarkTransform()->SetModeToRigidBody();
    icpTransform->SetMaximumNumberOfIterations(20);
    icpTransform->StartByMatchingCentroidsOn();
    icpTransform->Modified();
    icpTransform->Update();

    vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter2 = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
    transformFilter2->SetInputData(sourceGlyphFilter->GetOutput());
    transformFilter2->SetTransform(icpTransform);
    transformFilter2->Update();

    vtkSmartPointer<vtkPolyDataMapper> sourceMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    sourceMapper->SetInputConnection(sourceGlyphFilter->GetOutputPort());

    vtkSmartPointer<vtkActor> sourceActor = vtkSmartPointer<vtkActor>::New();
    sourceActor->SetMapper(sourceMapper);
    sourceActor->GetProperty()->SetColor(0,1,0);
    sourceActor->GetProperty()->SetPointSize(3);

    vtkSmartPointer<vtkPolyDataMapper> targetMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    targetMapper->SetInputConnection(targetGlyphFilter->GetOutputPort());

    vtkSmartPointer<vtkActor> targetActor = vtkSmartPointer<vtkActor>::New();
    targetActor->SetMapper(targetMapper);
    targetActor->GetProperty()->SetColor(1,0,0);
    targetActor->GetProperty()->SetPointSize(3);

    vtkSmartPointer<vtkPolyDataMapper> solutionMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    solutionMapper->SetInputConnection(transformFilter2->GetOutputPort());

    vtkSmartPointer<vtkActor> solutionActor = vtkSmartPointer<vtkActor>::New();
    solutionActor->SetMapper(solutionMapper);
    solutionActor->GetProperty()->SetColor(0,0,1);
    solutionActor->GetProperty()->SetPointSize(3);

    vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
    renderer->SetBackground(1.0, 1.0, 1.0);
    renderer->AddActor(sourceActor);
    renderer->AddActor(targetActor);
    renderer->AddActor(solutionActor);

    vtkSmartPointer<vtkGenericOpenGLRenderWindow> renderWindow = vtkSmartPointer<vtkGenericOpenGLRenderWindow>::New();
    renderWindow->AddRenderer(renderer);

    QVTKOpenGLNativeWidget w;
    w.resize(640, 480);
    w.setRenderWindow(renderWindow);
    w.setWindowTitle("PolyDataICP");
    w.show();

    return a.exec();
}

运行效果：

三、纹理映射

纹理映射(Texture Mapping)是将纹理空间中的纹理像素映射到屏幕空间中的像素的过程。纹理生成过程实质上是将所定义的纹理映射为某种三维物体表面的属性，并参与后续的光照计算。在三维图形中，纹理映射运用得非常广泛，尤其是描述具有真实感的物体。比如绘制一面砖墙，就可以使用一幅具有真实感的图像或者照片作为纹理贴到一个矩形上，这样一面逼真的砖墙就绘制好了。

实现纹理映射主要是建立纹理空间与模型空间、模型空间与屏幕空间之间的映射关系。

而对于无参数化曲面的纹理映射技术，通常需要将纹理空间到模型空间的映射分解为两个简单映射。这里需要引入一个包围景物的中介三维曲面作为中介映射媒介，主要实现步骤如下。

先将二维纹理空间映射为一个简单的三维物体表面，例如球面、圆柱面等;然后将上述中介物体表面的纹理映射到模型表面，例如以模型表面法线与中介模型的交点作为映射点。这样即可实现由纹理空间到模型空间的映射。

VTK中定义了多个类实现纹理空间到模型空间的映射，例如vtkTextureMapToPlane通过一个平面建立纹理空间到模型空间的映射关系；vtkTextureMapToCylinder 通过圆柱面建立映射关系；vtkTextureMapToSphere 通过球面建立映射关系。vtkTexture 则实现加载纹理。另外，vtkTransformTextureCoords也是一个非常有用的类，可以实现纹理坐标的平移和缩放，例如，如果要实现重复纹理，只需通过vtkTransformTextureCoords::SetScale()将纹理坐标每个方向进行放大，如由[0,1]变换到[0,10]即可。

示例代码：

#include <QApplication>
#include <vtkSmartPointer.h>
#include <vtkRenderer.h>
#include <vtkGenericOpenGLRenderWindow.h>
#include <QVTKOpenGLNativeWidget.h>
#include <vtkBMPReader.h>
#include <vtkTexture.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkTextureMapToCylinder.h>

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

    vtkSmartPointer<vtkBMPReader> texReader = vtkSmartPointer<vtkBMPReader>::New();
    texReader->SetFileName("D:/data/masonry.bmp");

    vtkSmartPointer<vtkTexture> texture = vtkSmartPointer<vtkTexture>::New();
    texture->SetInputConnection(texReader->GetOutputPort());

    vtkSmartPointer<vtkXMLPolyDataReader> modelReader = vtkSmartPointer<vtkXMLPolyDataReader>::New();
    modelReader->SetFileName("D:/data/cow.vtp");

    vtkSmartPointer<vtkTextureMapToCylinder> texturemap = vtkSmartPointer<vtkTextureMapToCylinder>::New();
    texturemap->SetInputConnection(modelReader->GetOutputPort());

    vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    mapper->SetInputConnection(texturemap->GetOutputPort());

    vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper( mapper );
    actor->SetTexture( texture );

    vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
    renderer->AddActor(actor);
    renderer->SetBackground(1.0, 1.0, 1.0);

    vtkSmartPointer<vtkGenericOpenGLRenderWindow> renderWindow = vtkSmartPointer<vtkGenericOpenGLRenderWindow>::New();
    renderWindow->AddRenderer(renderer);

    QVTKOpenGLNativeWidget w;
    w.resize(640, 480);
    w.setRenderWindow(renderWindow);
    w.setWindowTitle("TextureMap");
    w.show();

    return a.exec();
}

运行效果：