Creating a new segment from an open surface

Support
1

Hi,

I am trying to get a new segment from a part of a (closed) surface model (obtained from a segment in the Segment Editor) via Python scripting. I am working with a large closed surface, onto which I identify a part of the points and cells (mesh triangles) that fulfill a given set of conditions (see first image).

sense nom
sense nom752×1294 65.7 KB

I can generate a new surface isolating these points from the rest and creating a new vtkPolyData object, but I am left with open ends of the surface (image 2).

sense nom 2
sense nom 2700×500 22.8 KB

My goal is to generate a new segment from this open surface, but I understand that I need a closed surface to do so. If that is the case, what is the best way to close these open ends? I have tried to use vtkClipClosedSurface, but I haven’t found an automatic way to detect the clipping planes.

Thank you!

2

Maybe you could try using this:

edges = vtk.vtkFeatureEdges()
edges.BoundaryEdgesOn()

And the get all disconnected regions with vtkConnectivityFilter. Maybe you can use this data to initialize your algorithms.

This maybe a little bit related

1 Like
3

Thank you @mau_igna_06 for your reply. Following your advise, I have tried to get the boundary edges using the following code:

edges = vtk.vtkFeatureEdges()
edges.BoundaryEdgesOn()
edges.FeatureEdgesOff()
edges.ManifoldEdgesOff()
edges.SetInputData(circularSegment)
edges.Update()
output = edges.GetOutput()

In this example, the object circular segment is a the vtkPolyData containing the input surface from the original post. From here, I would expect that output would be a vtkPolyData containing two cells (vtkPolyLines?) corresponding to the two boundaries of the surface. Instead, I get all other edges of the original surface. In the image below you can see one of the open ends of the surface. On the left, the blue wireframe corresponds to the output object, while on the right I superposed the wireframe representation of the original surface in white.

sense nom
sense nom1824×664 259 KB

Is this expected behavior? Do you think there is something wrong with my original vtkPolyData? Also, I noticed that if I also introduce this:

edges.NonManifoldEdgesOff()

I get an empty object in return, which seems incorrect. Any ideas?

4

Do you calculated the normals of the model before processing it? if not try that

Try to know how many non-manifold edges your mesh has by using the corresponding option on and all the other off in the edges filter. If the mesh has lot of non manifold edges that can be a problem.

Try cleaning your mesh with vtkCleanPolyData before processing, that may help

Hope these ideas can help

Mauro

1 Like
5

I now introduced the normals computation, which didn’t have any noticeable effect in the edges computation. Then, I was experimenting a bit to see how many cells were computed for each of the 4 options (boundary, feature, manifold and non-manifold) and something strange is happening, at least if what I understand is correct. I have 0 boundary and feature edges, 53 manifold edges (these are the ones in the boundary!) and the rest (29537) are non-manifold… In the image I superposed the manifold edges form one of the open ends in red to the non-manifold in blue that we saw before.

sense nom
sense nom994×754 132 KB

At least I can continue to implement some filter to close those open ends, but it would be nice to have the expected behavior instead of this. If any ideas come up feel free to reply!

Thank you Mauro!

6

You are welcome.

I think to deep further into this you should ask here

Mauro

1 Like
7

Hi again,

In the end, I repeated the process computing the normals and passing the vtkCleanPolyData filter and everything works as expected now, I am not sure of what I was doing wrong before. In the end, the code looks something like this (I paste here from the vtkPolyData definition onwards):

# Define vtkPolyData from selected points (vtkPoints object) and polys (vtkCellArray object)
circularSegment = vtk.vtkPolyData()
circularSegment.SetPoints(pointsCircularSegmentModel)
circularSegment.SetPolys(cellArrayCircularSegmentModel)

# We can to compute the normals for all mesh triangles
normals = vtk.vtkPolyDataNormals()
normals.SetInputData(circularSegment)
normals.SetFeatureAngle(80)
normals.AutoOrientNormalsOn()
normals.UpdateInformation()
normals.Update()
circularSegment = normals.GetOutput()

# Clean the vtkPolyData. Without this, normal edges are seen as non-manifold and boundary edges as manifold. Normals have to be computed prior to this
cleanPolyData = vtk.vtkCleanPolyData()
cleanPolyData.SetInputData(circularSegment)
cleanPolyData.Update()
circularSegment = cleanPolyData.GetOutput()

# Now we can extract the boundary edges with the vtkFeatureEdges filter
edgesFilter = vtk.vtkFeatureEdges()
edgesFilter.BoundaryEdgesOn()
edgesFilter.FeatureEdgesOff()
edgesFilter.ManifoldEdgesOff()
edgesFilter.NonManifoldEdgesOff()
edgesFilter.SetInputData(circularSegment)
edgesFilter.Update()
boundaryEdges = edgesFilter.GetOutput()

# The vtkConnectivityFilter allows for separation of both boundaries (in my case I only have two groups, so I use first and final cell seeds for separation)
connectivityFilter = vtk.vtkConnectivityFilter()
connectivityFilter.SetExtractionModeToCellSeededRegions()
connectivityFilter.AddInputData(boundaryEdges)
connectivityFilter.AddSeed(0)
connectivityFilter.Update()
proximalBorder = connectivityFilter.GetOutput()

connectivityFilter2 = vtk.vtkConnectivityFilter()
connectivityFilter2.SetExtractionModeToCellSeededRegions()
connectivityFilter2.AddInputData(boundaryEdges)
connectivityFilter2.AddSeed(boundaryEdges.GetCell(boundaryEdges.GetNumberOfPoints() - 1).GetPointId(0))
connectivityFilter2.Update()
distalBorder = connectivityFilter2.GetOutput()

Now I continued the code introducing the vtkClipClosedSurface filter to close the surface at the edges (criteria for plane definition are a bit too adapted at my example, just want it to work altogether):

# We have to define planes for the vtkClipClosedSurface filter. First pool all points from both groups
proximalBorderPoints = np.ndarray([proximalBorder.GetNumberOfPoints(), 3])
distalBorderPoints = np.ndarray([distalBorder.GetNumberOfPoints(), 3])

for idx in range(proximalBorder.GetNumberOfPoints()):
    proximalBorderPoints[idx] = proximalBorder.GetPoints().GetPoint(idx)  
for idx in range(distalBorder.GetNumberOfPoints()):
    distalBorderPoints[idx] = distalBorder.GetPoints().GetPoint(idx)
    
# Get point with maximal z coordinate in proximal boundary and minimum z coordinate in distal boundary
proximalBorderZMax = np.amax(proximalBorderPoints[:, 2])
distalBorderZMin = np.amin(distalBorderPoints[:, 2])

# Cutting planes definition
proximalPlane = vtk.vtkPlane()
proximalPlane.SetNormal(0, 0, 1)
proximalPlane.SetOrigin(0, 0, proximalBorderZMax)
distalPlane = vtk.vtkPlane()
distalPlane.SetNormal(0, 0, -1)
distalPlane.SetOrigin(0, 0, distalBorderZMin)
clippingPlanes = vtk.vtkPlaneCollection()
clippingPlanes.AddItem(proximalPlane)
clippingPlanes.AddItem(distalPlane)

# Application of the vtkClipClosedSurface filter
clip = vtk.vtkClipClosedSurface()
clip.AddInputData(circularSegment)
clip.SetClippingPlanes(clippingPlanes)
clip.Update()
clippedCircularSegment = clip.GetOutput()

Now the object looks all right:

sense nom
sense nom994×754 49.8 KB

I am doing all of this in a Slicer context, and now I am trying to store it as a new segment in a vtkMRMLSegmentationNode. I am using the following line:

segmentationNode.AddSegmentFromClosedSurfaceRepresentation(clippedCircularSegment)

But this is not working (it returns an empty string). I have been able to define a new vtkMRMLModelNode and visualize the object in Slicer:

circularSegmentModelNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelNode")
circularSegmentModelNode.SetAndObservePolyData(clippedCircularSegment)

But I cannot add it as a new segment in the vtkMRMLSegmentationNode, which is what I would like to do, and I don’t really understand why… Any help with this?