VTK – Cut Image Data By Plane

The post shows how to cut image by a plane smoothly.

Header And Macros

#include <vtkVersion.h>
#include <vtkImageData.h>
#include <vtkImageMapper3D.h>
#include <vtkImageStencil.h>
#include <vtkImageStencilData.h>
#include <vtkImageToImageStencil.h>
#include <vtkJPEGReader.h>
#include <vtkPointData.h>
#include <vtkSmartPointer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkInteractorStyleImage.h>
#include <vtkRenderer.h>
#include <vtkImageActor.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkPolyData.h>
#include <vtkSTLReader.h>
#include <vtkProperty.h>
#include <vtkImageProperty.h>
#include <vtkSphereSource.h>
#include <vtkPolyDataMapper.h>
#include <vtkImageDataGeometryFilter.h>
#include <vtkImageMarchingCubes.h>
#include <vtkPolyDataNormals.h>
#include <vtkFillHolesFilter.h>
#include <vtkCleanPolyData.h>
#include <vtkTriangleFilter.h>
#include <vtkAppendPolyData.h>
#include <vtkImageMarchingCubes.h>
#include <vtkPlane.h>
#include <vtkImageWeightedSum.h>
#include <vtkConeSource.h>
#include <vtkPolyDataConnectivityFilter.h>
#include <vtkXMLImageDataWriter.h>
#include <vtkXMLImageDataReader.h>

#include "../vtkLearn/tool.h"

#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();

unsigned char inval = 255;
unsigned char outval = 0;
double spacingVal = 0.1;

Read And Write Image Data

    vtkSPtrNew( writer, vtkXMLImageDataWriter );
    writer->SetInputData( _sphereImage );
    writer->SetFileName( "image.vti" );
    writer->Write();

    vtkSPtrNew( reader, vtkXMLImageDataReader );
    reader->SetFileName( "image.vti" );
    reader->Update();

    vtkSPtr<vtkImageData> sphereImage = reader->GetOutput();

Convert Image To PolyData

vtkSPtr<vtkPolyData> imageToPolyData( vtkSPtr<vtkImageData> image )
{
    vtkSPtrNew( marchingCubes, vtkImageMarchingCubes );
    marchingCubes->SetInputData( image );
    marchingCubes->SetValue( 0, (inval + outval)/2 );
    marchingCubes->SetNumberOfContours( 1 );
    marchingCubes->Update();

    auto fillHoles = vtkSmartPointer<vtkFillHolesFilter>::New();
    fillHoles->SetInputConnection( marchingCubes->GetOutputPort() );
    fillHoles->SetHoleSize(1000.0);

    auto normals = vtkSmartPointer<vtkPolyDataNormals>::New();
    normals->SetInputConnection(fillHoles->GetOutputPort());
    // Make the triangle winding order consistent
    normals->ConsistencyOn();
    normals->SplittingOff();
    normals->GetOutput()->GetPointData()->SetNormals(
        marchingCubes->GetOutput()->GetPointData()->GetNormals());
    normals->Update();

    return normals->GetOutput();
}

Convert PolyData To Image

vtkSPtr<vtkImageData> ConvertPolydataToImage(vtkPolyData *polydata)
{
    double bounds[6];
    polydata->GetBounds(bounds);

    int dim[3];
    for (int i = 0; i < 3; i++)
    {
        dim[i] = static_cast<int>( ceil((bounds[2 * i + 1] - bounds[2 * i]) / spacingVal) );
    }
    double origin[3];
    origin[0] = bounds[0] + spacingVal / 2;
    origin[1] = bounds[2] + spacingVal / 2;
    origin[2] = bounds[4] + spacingVal / 2;

    vtkSPtrNew(imageData, vtkImageData)
    imageData->SetSpacing(spacingVal, spacingVal, spacingVal);
    imageData->SetDimensions(dim);
    imageData->SetOrigin(origin);
    imageData->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
    imageData->SetExtent( 0, dim[0] - 1, 0, dim[1] - 1, 0, dim[2] - 1 );

    // fill the imageData with foreground voxels
    vtkIdType count = imageData->GetNumberOfPoints();
    for (vtkIdType i = 0; i < count; ++i)
    {
        imageData->GetPointData()->GetScalars()->SetTuple1(i, inval);
    }

    // polygonal data --> imageData stencil:
    vtkSPtrNew(pd2stenc, vtkPolyDataToImageStencil);
    pd2stenc->SetInputData(polydata);
    pd2stenc->SetOutputOrigin(origin);
    pd2stenc->SetOutputSpacing(spacingVal, spacingVal, spacingVal);
    pd2stenc->SetOutputWholeExtent(imageData->GetExtent());
    pd2stenc->Update();

    // cut the corresponding white imageData and set the background:
    vtkSPtrNew(imgstenc, vtkImageStencil);
    imgstenc->SetInputData(imageData);
    imgstenc->SetStencilConnection(pd2stenc->GetOutputPort());
    imgstenc->ReverseStencilOff();
    imgstenc->SetBackgroundValue( outval );
    imgstenc->Update();

    imageData->DeepCopy(imgstenc->GetOutput());
    return imageData;
}

Cut Image

void CutImage(vtkImageData *image, vtkPlane *cutPlane )
{
    if( nullptr == image || nullptr == cutPlane )
    {
        return ;
    }
    int negtiveCount = 0;
    PointStruct spaceVals( image->GetSpacing() );
    double maxDistance = sqrt(spaceVals[0] * spaceVals[0] + spaceVals[1] * spaceVals[1] + spaceVals[2] * spaceVals[2]);
    double threshold = 0.01;
    cout << "maxDistance: " << maxDistance << endl;
    auto isoSurfaceValue = (inval+outval)*0.5;
    for (int i = 0; i < image->GetNumberOfPoints(); i++)
    {
        double pointPos[3];
        image->GetPoint(i, pointPos);
        double scaleVal = image->GetPointData()->GetScalars()->GetTuple1(i);
        if (scaleVal < isoSurfaceValue ) //(outval + 0.1))  // speed up
            continue;

        double val = cutPlane->EvaluateFunction(pointPos);
        if( abs( val ) < maxDistance ) // smooth cut plane
        {
            auto value = (1 + val / maxDistance) * isoSurfaceValue;
            if( value > inval ){ value = inval; }
            if( value < outval ){ value = outval; }

            image->GetPointData()->GetScalars()->SetTuple1(i, value);
        }
        else if( val < 0 ) // clip part
        {
            image->GetPointData()->GetScalars()->SetTuple1(i, outval);
            negtiveCount++;
        }
    }
}

int main(int, char *[])
{
    setbuf( stdout, nullptr );

    vtkSPtrNew( sphereSource, vtkSphereSource );
    sphereSource->SetRadius(1);
    sphereSource->SetPhiResolution(10);
    sphereSource->SetThetaResolution(10);
    sphereSource->Update();

    vtkSPtr<vtkPolyData> spherePd = sphereSource->GetOutput();

    vtkSPtr<vtkImageData> sphereImage = ConvertPolydataToImage( spherePd );

    auto originImage = imageToPolyData( sphereImage );

    vtkSPtrNew( cutPlane, vtkPlane );
    cutPlane->SetOrigin( 0, 0, 0 );
    cutPlane->SetNormal( 0, -1, 0 );
    CutImage( sphereImage, cutPlane );

    auto newImage = imageToPolyData( sphereImage );

    vtkSPtrNew( leftMapper, vtkPolyDataMapper );
    leftMapper->SetInputData( originImage );
    leftMapper->Update();
    leftMapper->SetScalarVisibility( false );

    vtkSPtrNew( leftActor, vtkActor );
    leftActor->SetMapper( leftMapper );
    leftActor->GetProperty()->SetColor( 1, 1, 1 );
    // ============= Finish: show image data ===================


    vtkSPtrNew( rightMapper, vtkPolyDataMapper );
    rightMapper->SetInputData( newImage);
    rightMapper->SetScalarVisibility( false );
    rightMapper->Update();

    vtkSPtrNew( rightActor, vtkActor );
    rightActor->SetMapper( rightMapper );

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

    // Setup renderers
    vtkSPtrNew( leftRenderer, vtkRenderer );
    leftRenderer->SetViewport( leftViewport );
    leftRenderer->AddActor( leftActor );
    leftRenderer->SetBackground(.6, .5, .4);
    leftRenderer->ResetCamera();

    vtkSPtrNew( rightRenderer, vtkRenderer );
    rightRenderer->SetViewport(rightViewport);
    rightRenderer->AddActor( rightActor );
    rightRenderer->SetBackground(.4, .5, .6);
    rightRenderer->SetActiveCamera( leftRenderer->GetActiveCamera() );

    // Setup render window
    vtkSPtrNew( renderWindow, vtkRenderWindow );
    renderWindow->AddRenderer( leftRenderer );
    renderWindow->AddRenderer( rightRenderer );
    renderWindow->SetSize( 900, 300 );

    // Setup render window interactor
    vtkSPtrNew( renderWindowInteractor, vtkRenderWindowInteractor );

    // Render and start interaction
    renderWindowInteractor->SetRenderWindow(renderWindow);
    renderWindowInteractor->Start();

    return EXIT_SUCCESS;
}