How can I obtain the intersection line from a isosurface which is cut by other surface?

Hello,

I would like to obtain the representation of the intersection of an isosurface with other surface (in this example code, it is just a plane). The isosurface is extrated from a volume image via vtkMarchingCubes, vtkContourFilter or vtkFlyingEdges3D. And, in this simple example I use vtkCutter to try to obtain the intersection. However, after connect them (isosurface filter and vtkCutter) nothing is represented with the vtkPolyDataMapper + vtkActor, I only can see the background.

I’m pretty confused. I checked there should be an output from any of the isosurface filters (vtkMarchingCubes, vtkContourFilter or vtkFlyingEdges3D), because I can display the isosurface (with a mapper and an actor) before and after aplying other filters (like vtkStripper, vtkSmoothPolyDataFilter + vtkPolyDataNormals).But I have nothing when I connect them to the vtkCutter filter.

I attach the snippet code in order to show how I tried to connect the elements.
Please, could you help me to connect the isosurface to the cutting filter to obtain the intersecting lines?

Thank you very much in advance!

#include<vtkImageImport.h>
#include<vtkContourFilter.h>
#include<vtkMarchingCubes.h>
#include<vtkFlyingEdges3D.h>
#include<vtkPlane.h>
#include<vtkCutter.h>
#include<vtkPolyDataMapper.h>
#include<vtkNamedColors.h>
#include<vtkProperty.h>
#include<vtkActor.h>
#include<vtkRenderer.h>
#include<vtkRenderWindow.h>
#include<vtkRenderWindowInteractor.h>

int main(int argc, char *argv[])
{
    // Some missing code to load the volume image
    // [...]

    vtkSmartPointer<vtkImageImport> img_import_v1 = s_v_pixmap_synth_1->getImageImport();
    img_import_v1->Update();

    // I comment the alternative options, all of them have the same result
    //vtkSmartPointer<vtkContourFilter> isosurface = vtkSmartPointer<vtkContourFilter>::New();
    vtkSmartPointer<vtkMarchingCubes> isosurface = vtkSmartPointer<vtkMarchingCubes>::New();
    //vtkSmartPointer<vtkFlyingEdges3D> isosurface = vtkSmartPointer<vtkFlyingEdges3D>::New();
    isosurface->SetInputConnection(img_import_v1->GetOutputPort());
    //isosurface->ComputeGradientsOn();  // vtkFlyingEdges3D
    //isosurface->ComputeScalarsOff();   // vtkFlyingEdges3D
    isosurface->SetValue(0, 1150);

    vtkSmartPointer<vtkPolyDataMapper> isoMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    isoMapper->SetInputConnection(isosurface->GetOutputPort()); // iso
    isoMapper->ScalarVisibilityOff();

    // Create a plane to cut,here it cuts in the XZ direction (xz
    // normal=(1,0,0);XY =(0,0,1),YZ =(0,1,0)
    vtkSmartPointer<vtkPlane> plane = vtkSmartPointer<vtkPlane>::New();
    plane->SetOrigin(0, 0, 10);
    plane->SetNormal(0, 0, 1);

    // Create cutter plane
    vtkSmartPointer<vtkCutter> cutter = vtkSmartPointer<vtkCutter>::New();
    cutter->SetCutFunction(plane);
    cutter->SetInputConnection(isosurface->GetOutputPort()); // iso / iso Mapper (?)
    cutter->Update();

    vtkNew<vtkPolyDataMapper> cutterMapper;
    cutterMapper->SetInputConnection(cutter->GetOutputPort());
    cutterMapper->SetResolveCoincidentTopologyToPolygonOffset();
    cutterMapper->Update();

    // Create plane actor that should contain the intersction lines
    vtkSmartPointer<vtkNamedColors> colors = vtkSmartPointer<vtkNamedColors>::New();
    vtkSmartPointer<vtkActor> planeActor = vtkSmartPointer<vtkActor>::New();
    planeActor->GetProperty()->SetColor(colors->GetColor3d("Yellow").GetData());
    planeActor->GetProperty()->SetLineWidth(2);
    planeActor->GetProperty()->SetAmbient(1.0);
    planeActor->GetProperty()->SetDiffuse(0.0);
    planeActor->SetMapper(cutterMapper);
    
    // Create the RenderWindow and Renderer
    vtkSmartPointer<vtkRenderer> ren1 = vtkSmartPointer<vtkRenderer>::New();
    vtkSmartPointer<vtkRenderWindow> renWin = vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(ren1);
    vtkSmartPointer<vtkRenderWindowInteractor> iren = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    iren->SetRenderWindow(renWin);

    // This actor is only useful to display the original isosurface
    vtkSmartPointer<vtkActor> isoActor = vtkSmartPointer<vtkActor>::New();
    isoActor->SetMapper(isoMapper);  //isoMapper
    isoActor->GetProperty()->SetColor(colors->GetColor3d("Ivory").GetData());

    // Add the actors to the renderer, set the background and size
    //ren1->AddActor(isoActor);  // uncomment this to see the whole isosurface
    ren1->AddActor(planeActor);
    ren1->SetBackground(colors->GetColor3d("SlateGray").GetData());
    ren1->ResetCamera();

    renWin->SetSize(640, 480);
    renWin->SetWindowName("HeadBone");

    renWin->Render();
    iren->Start();
    return 0;
}

Just to clarify, I not looking for an orthogonal plane cut, it is just a simple example to illustrate and show it does not work as I built the chain. I try to do something more general, so the initial image volume (and it is associated isosurface) could be located and orientated in ay possible way, an the implicit cutting function too.
I tried to do this task in more complicated ways, but I was not able to make it work even in this simpler scenario.
Thank you in advance for any clue !

I reformulate muy question: How can I use vtkCutter on the output of an isosurface filter?
I guessed the isosurface filters (like vtkMarchingCubes, vtkContourFilter or vtkFlyingEdges3D) provide an output in vtkPolydata format, which could feed vtkCutter.
Do I need a triangulate filter between isosurface filter and vtkCutter? (I tried vtkStripper, but I got nothing)
What am I doing wrong?

Thanks for any clue / example.

A quick glance suggests you are close… Are you sure the plane is actually cutting through the isosurface? To debug it, render a plane at the same location as the implicit function vtkPlane (you could do this a variety of ways including using a vtkPlaneSource or vtkImplicitPlaneWidget2/vtkImplicitPlaneRepresentation or use ParaView).

1 Like

You are absolutely right, Will !
I had the plane in the wrong location to see any intersection output.
As soon as I followed your advise I could visualize the isosurface was under the cutting plane (well, the vtkPlaneSource clone).
Thank you so much for your help !!

I love happy endings!

1 Like