VTK使用轮廓提取技术生成等值面

2022年2月24日 7753点热度 6人点赞 0条评论

参考之前记录的标量颜色映射，很多时候为了搞清楚标量在空间或者平面上的分布趋势，通常会通过提取等值面或者等值线的方式将这些数据可视化，这里会记录怎么使用VTK实现这样的功能。

工欲善其事必先利其器，实现这个功能的核心是vtkContourFIlter类，它提供了等值线（等值面）提取算法，只需要准备好数据，然后设置好轮廓参数，就可以实现了。关于过滤器的使用可以参考之前关于VTK渲染流水线的博客。

vtkContourFilter

这个过滤器的功能是从众多的标量值中找出标量值相等的位置，把这些位置连成线或者生成面，当然这个过程中会使用已有的标量值进行插值。它可以将任何数据集作为输入并输出等值面或者等值线（这里可以参考在标量颜色映射中对vtkDataSet类的介绍），输出的数据的形式却决于输入数据的维度，如果是3D单元数据那么输出的就是等值面，2D单元输出等值线，1D单元则输出一些点。

要使用这个过滤器必须为过滤器设置轮廓值（等值面的值），可以使用SetValue()来指定每一个轮廓值，或者使用GenerateValues()来生成一些均匀分布的轮廓值。

生成轮廓值需要耗费大量的时间，可以使用vtkScalarTree来加速这个操作，标量树用于快速定位包好等值表面的单元，尤其是对需要提取多个轮廓的时候效果极佳，如果要使用标量树需要调用UseScalarTreeOn()开启。

等值线提取

将之前标量颜色映射的示例拿过来，增加一个vtkContourFIlter。

//为三角形准备三个顶点
    vtkNew<vtkPoints> points;
    points->InsertNextPoint(1, 0, 0);
    points->InsertNextPoint(0, 0, 0);
    points->InsertNextPoint(0, 1, 0);

    //生成点的标量数据
    vtkNew<vtkIntArray> scalar;
    scalar->InsertNextTuple1(30);
    scalar->InsertNextTuple1(100);
    scalar->InsertNextTuple1(50);

    //使用顶点构建一个三角形 点使用ID传入
    vtkNew<vtkTriangle> triangle;
    triangle->GetPointIds()->SetId(0, 0);
    triangle->GetPointIds()->SetId(1, 1);
    triangle->GetPointIds()->SetId(2, 2);


    //将三角形传入单元集
    vtkNew<vtkCellArray> cellArray;
    cellArray->InsertNextCell(triangle);

    //为多边形数据添加单元集和点集
    vtkNew<vtkPolyData> polydata;
    polydata->SetPolys(cellArray);
    polydata->SetPoints(points);
    //为点数据添加标量
    polydata->GetPointData()->SetScalars(scalar);


    //添加一个过滤器
    vtkNew<vtkContourFilter> contourFilter;
    contourFilter->SetInputData(polydata);
    //设置等值线的值
    contourFilter->SetValue(0, 50);
    contourFilter->SetComputeScalars(true);
    contourFilter->Update();


    vtkNew<vtkPolyDataMapper> mapper;
    mapper->SetInputData(polydata);
    mapper->SetScalarRange(scalar->GetRange());
    mapper->Update();

    vtkNew<vtkPolyDataMapper> contourMapper;
    contourMapper->SetInputData(contourFilter->GetOutput());
    contourMapper->ScalarVisibilityOff();

    //为映射器添加一个颜色映射表
    //不添加颜色映射表会得到一个默认的颜色表
    //如果需要自定义颜色表需要调用SetTableValue()设置
 	vtkNew<vtkLookupTable> lookupTable;
    lookupTable->SetTableRange(scalar->GetRange());
    lookupTable->Build();
 	mapper->SetLookupTable(lookupTable);
    contourMapper->SetLookupTable(lookupTable);

    vtkNew<vtkScalarBarActor> barActor;
    barActor->SetLookupTable(lookupTable);
    barActor->SetNumberOfLabels(6);

    vtkNew<vtkActor> actor;
    actor->SetMapper(mapper);
    actor->GetProperty()->SetColor(255, 0, 0);

    vtkNew<vtkActor> contourActor;
    contourActor->SetMapper(contourMapper);
    contourActor->GetProperty()->SetColor(0,0,0);

    vtkNew<vtkRenderer> renderer;
    vtkNew<vtkRenderWindow> window;
    window->AddRenderer(renderer);

    vtkNew<vtkRenderWindowInteractor> interactor;
    interactor->SetRenderWindow(window);

    renderer->AddActor(actor);
    renderer->AddActor(contourActor);
    renderer->AddActor(barActor);
    window->Render();
    interactor->Start();

等值面轮廓提取

提取等值面与提取等值线的过程如出一辙，首先得提供一个三维点数据与标量数据。

(x^2 + (9/4) y^2 + z^2 - 1)^3 - x^2*z^3 - (9/80) y^2*z^3 =F, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}

我在网上找到了这么一个函数，据说这个F=0的时候，xyz的解会形成一个小心心。于是我生成了一个40*40*40的三维网格数据，然后每个网格的标量值使用这个函数进行赋值，最后再求F=0时的等值面。

vtkNew<vtkStructuredGrid> structuredGrid;
    vtkNew<vtkPoints> points;
    vtkNew<vtkFloatArray> scalars;
    unsigned int gridSize = 40;
    //构造一个点云，使用函数为点云构造标量值
    for (unsigned int i = 0; i < gridSize; i++)
    {
        for (unsigned int j = 0; j < gridSize; j++)
        {
            for (unsigned int k = 0; k < gridSize; k++)
            {
                points->InsertNextPoint(i, j, k);
                //实现函数 (x^2 + (9/4) y^2 + z^2 - 1)^3 - x^2*z^3 - (9/80) y^2*z^3 =F, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}
                double x = (i - (gridSize / 2.0)) / 10.0;
                double y = (j - (gridSize / 2.0)) / 10.0;
                double z = (k - (gridSize / 2.0)) / 10.0;
                scalars->InsertNextTuple1( pow(x*x + (9.0 / 4.0)*y*y + z*z- 1.0,3) - x*x*pow(z,3) - (9.0 / 80.0)*y*y * pow(z,3) );
            }
        }
    }



    //设置网格大小
    structuredGrid->SetDimensions(gridSize, gridSize, gridSize);
    //附加数据
    structuredGrid->SetPoints(points);
    structuredGrid->GetPointData()->SetScalars(scalars);


    //设置映射器
    vtkNew<vtkDataSetMapper> gridMapper;
    gridMapper->SetInputData(structuredGrid.Get());
    gridMapper->ScalarVisibilityOff();


    vtkNew<vtkActor> gridActor;
    gridActor->SetMapper(gridMapper.Get());
    gridActor->GetProperty()->SetOpacity(0.3);


    //设置等值线过滤器
    vtkNew<vtkContourFilter> contourFilter;
    contourFilter->SetInputData(structuredGrid.Get());
    //为标量为0时取一个等值面
    contourFilter->SetValue(0, 0);
    contourFilter->Update();

    //创建颜色表 这里只输入颜色值的数量，然后自动生成表
    vtkNew<vtkLookupTable> lut;
    lut->SetTableRange(scalars->GetRange());
    lut->Build();

    //为生成的等值面添加法向
    vtkNew<vtkPolyDataNormals> normals;
    normals->SetInputConnection(contourFilter->GetOutputPort());
    normals->SetFeatureAngle(45);
    normals->Update();

    vtkNew<vtkDataSetMapper> contourMaper;
    contourMaper->SetInputData(normals->GetOutput());
    contourMaper->SetLookupTable(lut);
    contourMaper->SetScalarRange(scalars->GetRange());
    vtkNew<vtkActor> contourActor;
    contourActor->SetMapper(contourMaper);


    vtkNew<vtkRenderer> renderer;
    vtkNew<vtkRenderWindow> renderWindow;
    renderWindow->AddRenderer(renderer);


    vtkNew<vtkScalarBarActor> barActor;
    barActor->SetLookupTable(lut);


    vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
    renderWindowInteractor->SetRenderWindow(renderWindow);


    
    renderer->AddActor(contourActor);
    renderer->AddActor(barActor);
    renderer->AddActor(gridActor);

    renderWindow->Render();

    renderWindowInteractor->Start();