如何使用vtk进行三维形态学腐蚀膨胀

Question

想将.txt格式的点云数据读入，然后作三维形态学开闭操作，目标去掉图中植物根系外周的小土屑
原来的数据查资料后读成了vtkPolyData,但OpenClose3D先前报错说需要vtkImageData类型
用vtkPolyData绘制出来的原数据如下（3D）：

查阅各种资料写出的数据类型转换加形态学操作代码如下：

```c++
#include <vtkAutoInit.h>
VTK_MODULE_INIT(vtkRenderingOpenGL2);
VTK_MODULE_INIT(vtkInteractionStyle);
#include <vtkActor.h>
#include <vtkImageActor.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPoints.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkVertexGlyphFilter.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkImageOpenClose3D.h>
#include <vtkInteractorStyleImage.h>
#include <vtkSmartPointer.h>
#include <vtkImageAccumulate.h>
#include <vtkImageCast.h>
#include <vtkImageData.h>
#include <vtkImageMapper3D.h>
#include <vtkPointData.h>
#include <vtkImageThreshold.h>

#include <vtkImageStencil.h>




#include <sstream>

int main(int argc, char* argv[])
{
    vtkNew<vtkNamedColors> colors;

    // Verify input arguments
    if (argc != 2)
    {
        std::cout << "Usage: " << argv[0] << " Filename(.txt) e.g. TeapotPoints.txt"
            << std::endl;
        return EXIT_FAILURE;
    }
    // Get all data from the file
    std::string filename = argv[1];
    std::ifstream filestream(filename.c_str());

    std::string line;
    vtkNew<vtkPoints> points;
    char* p;
    while (std::getline(filestream, line))
    {
        double x, y, z;
        std::stringstream linestream;
        linestream << line;
        p = strtok((char*)line.data(), ",");
        x = atof(p);
        p = strtok(NULL, ",");
        y = atof(p);
        p = strtok(NULL, ",");
        z = atof(p);

        points->InsertNextPoint(x, y, z);
    }

    filestream.close();

    vtkNew<vtkPolyData> polyData;

    polyData->SetPoints(points);

    // 将vtkPolyData转换成vtkImageData，先前OpenClose3D报错说需要vtkImageData类型
    vtkSmartPointer<vtkImageData> whiteImage =
        vtkSmartPointer<vtkImageData>::New();
    double bounds[6];
    polyData->GetBounds(bounds);
    double spacing[3]; // desired volume spacing
    spacing[0] = 0.5;
    spacing[1] = 0.5;
    spacing[2] = 0.5;
    whiteImage->SetSpacing(spacing);

    // compute dimensions
    int dim[3];
    for (int i = 0; i < 3; i++)
    {
        dim[i] = static_cast<int>(ceil((bounds[i * 2 + 1] - bounds[i * 2]) / spacing[i]));
    }
    whiteImage->SetDimensions(dim);
    whiteImage->SetExtent(0, dim[0] - 1, 0, dim[1] - 1, 0, dim[2] - 1);

    double origin[3];
    origin[0] = bounds[0] + spacing[0] / 2;
    origin[1] = bounds[2] + spacing[1] / 2;
    origin[2] = bounds[4] + spacing[2] / 2;
    whiteImage->SetOrigin(origin);
    whiteImage->AllocateScalars(VTK_UNSIGNED_CHAR, 1);

    // fill the image with foreground voxels:
    unsigned char inval = 255;
    unsigned char outval = 0;
    vtkIdType count = whiteImage->GetNumberOfPoints();
    for (vtkIdType i = 0; i < count; ++i)
    {
        whiteImage->GetPointData()->GetScalars()->SetTuple1(i, inval);
    }

    // polygonal data --> image stencil:
    vtkSmartPointer<vtkPolyDataToImageStencil> pol2stenc =
        vtkSmartPointer<vtkPolyDataToImageStencil>::New();
    pol2stenc->SetInputData(polyData);
    pol2stenc->SetOutputOrigin(origin);
    pol2stenc->SetOutputSpacing(spacing);
    pol2stenc->SetOutputWholeExtent(whiteImage->GetExtent());
    pol2stenc->Update();

    // cut the corresponding white image and set the background:
    vtkSmartPointer<vtkImageStencil> imgstenc =
        vtkSmartPointer<vtkImageStencil>::New();
    imgstenc->SetInputData(whiteImage);
    imgstenc->SetStencilConnection(pol2stenc->GetOutputPort());
    imgstenc->ReverseStencilOff();
    imgstenc->SetBackgroundValue(outval);
    imgstenc->Update();



    vtkNew<vtkImageOpenClose3D> openClose;
    openClose->SetInputData(whiteImage);

    openClose->SetOpenValue(0);
    openClose->SetCloseValue(255);
    openClose->SetKernelSize(5, 5, 3);
    openClose->ReleaseDataFlagOff();
    openClose->GetOutput();
    openClose->GetCloseValue();
    openClose->GetOpenValue();

 

    vtkNew<vtkImageActor> openCloseActor;
    openCloseActor->GetMapper()->SetInputConnection(openClose->GetOutputPort());

    // Define viewport ranges
    // (xmin, ymin, xmax, ymax)
    double originalViewport[4] = { 0.0, 0.0, 0.5, 1.0 };
    double openCloseViewport[4] = { 0.5, 0.0, 1.0, 1.0 };

    // Setup renderers
  

    vtkNew<vtkRenderer> openCloseRenderer;
    openCloseRenderer->SetViewport(openCloseViewport);
    openCloseRenderer->AddActor(openCloseActor);
    openCloseRenderer->ResetCamera();
    openCloseRenderer->SetBackground(.4, .5, .7);
    openCloseRenderer->SetBackground(colors->GetColor3d("RoyalBlue").GetData());

    vtkNew<vtkRenderWindow> renderWindow;
    renderWindow->SetSize(600, 300);
    renderWindow->AddRenderer(openCloseRenderer);
    renderWindow->SetWindowName("ImageOpenClose3D");

    vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
    vtkNew<vtkInteractorStyleImage> style;

    renderWindowInteractor->SetInteractorStyle(style);

    renderWindowInteractor->SetRenderWindow(renderWindow);
    renderWindow->Render();
    renderWindowInteractor->Initialize();

    renderWindowInteractor->Start();

    return EXIT_SUCCESS;
}

运行结果：

好多代码似懂非懂不知问题在哪里，求解答

Answer 1

你这个需要将点云生成一个闭合的图形，然后转换为vtkImageData。
点云直接转vtkImageData可能就不会生成正确的vtkImageData；