New Isocontouring filter for discrete, labled volumetric data: vtkSurfaceNets3D

Development
1

There is a new isocontouring algorithm available in VTK master: vtkSurfaceNets3D. Unlike classic algorithms such as marching cubes or flying edges that operate on continuous scalar fields, vtkSurfaceNets3D operates on discrete, labeled volumetric data (e.g., segmentation maps, label maps). It simultaneously extracts labeled objects with fully shared boundaries. It is fast because it has been threaded, using concepts from flying edges such as edge-by-edge processing.

Note that besides being much faster, this algorithm produces better results than vtkDiscreteMarchingCubes and vtkDiscreteFlyingEdges which produce output which is not fully shared across object boundaries. For example, while vtkDiscreteMarchingCubes shares output points, it may produce coincident triangles. Also, surface nets has a built-in (optional) constrained smoothing step that is quite powerful and produces great results.

The sequential Surface Nets algorithm was originally proposed by Sarah Frisken now at Brigham/Harvard Medical School. (Her recent paper “SurfaceNets for Multi-Label Segmentations with Preservation of Sharp Boundaries”, J. Computer Graphics Techniques, 2022 provides details.) Originally patented, the algorithm is now unencumbered and extremely useful for creating complex geometry such as that found in anatomical atlases (see https://www.openanatomy.org/ for example). While just introduced in VTK, the plan is to make this fast surface nets implementation available in 3D Slicer and other applications.

Here is an example: segmented volume of 317^2 x 835 unsigned char voxel values with 105 labeled objects. On my 20-thread laptop, execution times are around 0.1 second to extract the objects and smooth them (~2.8 million triangles produced). The three images are: left) extracted voxelized surface net; center) after smoothing; and 3) using alternative windowed sinc smoothing.

image
image1788Ă—885 118 KB

I would encourage anyone with interest to try it out and provide feedback. Make sure you build with a non-sequential SMP backend for best performance. I prefer TBB because it does a better job of load balancing.

9 Likes
2

Will and all,

this looks absolutely fabulous! Dr. Frisken’s publication caught my eye, and I was very glad to discover that you have not only adopted SurfaceNets but also spent the time parallelizing it. Never having used VTK before, I was impressed that I was able to build VTK from source and even get it to run, using the python wrappers, in the course of a few hours - this speaks to your dedication to making it accessible much more than to my skills. We are currently working on a number (>200) of decently sized synchrotron µCT reconstructions of mouse mandibles (typical size ~1500x900x3500 voxels, 14-20 classes/labels/contour values). I managed to load one of our label matrices using vtkTIFFReader, and looking at https://gitlab.kitware.com/will.schroeder/surface-nets-testing/-/blob/main/RunSN.cxx I cobbled a Jupyter script together that actually worked, and was blindingly fast, too. Color me impressed.

But I seem to have exhausted my streak of getting things to work. Here are a few things that I would like to do, hoping that you might be able to point me in the right direction.

  1. Extract the BoundaryLabels Array from the output, and parse it, for example to determine the surface areas of the various interfaces. I have a rough sense on how I might be able to extract the vertices and faces (connectivity) from the vtkPolyData I can pull out via the .GetPolys() method. But I can’t figure out how to get the BoundaryLabels, and do not understand how these labels map onto the polys (in my case, triangles). Because I am less familiar with vtk, my instinct is to try and export to numpy arrays (and ultimately Matlab, where most of our current code lives). Any suggestions on how to do this are very welcome.

  2. One thing that we’d ultimately like to do is 3D print some of the models using a multi-material printer. For this, I was hoping to find a method more elegant and efficient than using .SetOutputStyleToSelected(), .InitializeSelectedLabelsList(), and AddSelectedLabel(label) and iterating over the label to save separate STLs. Any suggestions on how to transform the output into a format such as 3MF or a similar industry-standard format that supports multi-material printing would be greatly appreciated.

Thank you very much,

Derk

PS: And one more question. Using a very simple rendering pipeline on my Mac (OS X 12.6.3), I can generate beautiful windows and interact with the render. However, I cannot figure out how to close the window (the red button results in a spinning pizza) short of restarting the Python kernel. The keyboard commands I read about either do nothing or also result in Python becoming unresponsive. How do I get the renderer/interactor to stop? Here is my Jupyter script:

# set up tiff reader and import
reader = vtkTIFFReader()
reader.SetFileName('/path/to/SurfaceNetsTest.tiff')
reader.Update()
LabelMatrix = reader.GetOutput()

# vtkSurfaceNets3D class methods can be found here: https://vtk.org/doc/nightly/html/classvtkSurfaceNets3D.html#a51c2dd330853912b207b944831cc8d92
SN = vtkSurfaceNets3D()
# create instance of SurfaceNets3D
SN.SetInputData(LabelMatrix) 
# this set the number of classes ander their values for (eventual) coloring, I think
SN.GenerateValues(14,0,13)
# the next two are parameters that we need to explore
SN.SetNumberOfIterations(25)
SN.SetConstraintScale(2)
SN.Modified()
SN.Update()

# create color lookup table
lut = vtkLookupTable()
lut.SetNumberOfColors(14);
lut.SetHueRange(0.667, 0.0);
lut.Build();

# create PolyDataMapper
pdm = vtkPolyDataMapper()
pdm.SetInputData(SN.GetOutput())
pdm.SetLookupTable(lut)
pdm.SetScalarRange(0,13)
pdm.SetColorModeToMapScalars()

# Create Actor
actor = vtkActor();
actor.SetMapper(pdm);

# Create Renderer
ren = vtkRenderer();
ren.AddActor(actor);

# Create Window
renWin = vtkRenderWindow();
renWin.SetSize(500,500);
renWin.AddRenderer(ren);

# Create Interactor
iren= vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)

# Render w/interactive controls
renWin.Render();
iren.Start();
3

Derk I’m very happy that you tried this out and have had some success. We are in the middle of tuning this class etc. and all feedback is welcome. LIkely there will be some blood on the floor before we get it right, but we are pretty excited with the possibilities and will keep at it. It’s interesting, we are now working with several groups acquiring data in similar fashion using a variety of imaging modalities, and this seems to be an increasingly important visualization algorithm. Anyway, let me try to answer some of your questions in the following.

  1. The datasets in VTK consist of a geometric component (e.g., points GetPoints()) and a topological component polygons (e.g., GetPolys()). Then associated with points and polygons are point data (e.g., polydata.GetPointData()) and cell data (polydata.GetCellData()). which is basically an array of data arrays associated with the points or cells. I believe you want (polydata.GetCellData().GetArray(“BoundaryLabels”)). There is a 1-to-1 association between the nth triangle and the nth 2-tuple cell data. If you need help getting at the data in the boundary labels data array let us know and we’ll help.

  2. I agree that this is an inefficient interface. Just today I was conferring with some colleagues on how to efficiently extract labeled “objects” from the output. We have some ideas and are going to try and implement a “vtkAtlasModeller” filter to do so shortly. I’m behind the times with STL / 3MF, maybe someone else on this list can help. However, you probably noticed the STL writer etc so there is a way, albeit an inefficient one. We’re coming to realize that there are some auxiliary classes to write. We’re already using vtkPackLabels to manipulate the label maps, and it’s clear there are more to add. Anyway, if you are willing to bleed on the edge, when we have a prototype vtkAtlasModeller to try out, I’ll point you to a branch and you can provide feedback.

As far as the VTK building process, I am just amazed how well the community keeps this going. I just try to stay out of their way and throw an algorithm over the fence now and then :slight_smile: I’m curious, are you running with TBB enabled? or using a different SMPTools backend.

Did you try hitting the “e” key in the window? By default this works for me.

Finally, we have a SurfaceNets publication under review that we have selectively shared with researchers. I can’t say more at this point on a public forum until the review process is complete, but if you want to take a look and promise to keep it within your organization I can share it. Contact me at my email will.schroeder@kitware.com if you are interested.

4

Hello Will and all

Following Derk with some very enthusiastic feedback. We have been experimenting with your new filter here and it is indeed amazingly fast ! The fact that it includes a smoothing algorithm saves us a lot of time. And also the fact that it can operate on all the labels at once !
Before that, we were doing something (probably stupid) like

for all labels in a label map :
1. Extract one label and set it to 255;
2. Perform a slow gaussian Blur
3. Genereate a mesh with the Flying Edges at 128
4. Add the mesh to the renderWindow

Now it can be
LabelMap > Surface Nets > Add the mesh to the renderWindow
Easier to read and to maintain :slight_smile:

However we are also facing the same difficulties as Derk : How to extract properly and efficiently a mesh corresponding to a specific label ?

Please find attached a full example that demonstrates 2 different solutions to generate the meshes and render them, but we have the feeling here that there are more efficient solutions :sweat_smile:
Maybe someone can guide us ? Or maybe that’s something that should be integrated within the Surface Nets implementation itself, the possibility to output several meshes instead of one…

Anyway the program is meant to be used with the attached label-map.vti
Once build, just run it with one of the following line :

SimpleSurfaceNets.exe C:/Path/to/label-map.vti one-by-one 2 3
SimpleSurfaceNets.exe C:/Path/to/label-map.vti all-at-once 1
SimpleSurfaceNets.exe C:/Path/to/label-map.vti all-at-once 1 2 3 4 5

Thank you so much for your work

Simon

PS: for Derk, I think the program can help you if you are not too worried about performance, the extraction of all the meshes takes more time than the execution of the Surface Nets :wink:

label-map.vti (800.8 KB)
CMakeLists.txt (524 Bytes)
main.cpp (10.0 KB)

5

I have been waiting for a multi-label meshing algorithm for years. This sounds exactly like what I want. We wrote a very crude mesher in our “DREAM.3D” application years ago but it is rudimentary at best. An example of the data that we want to mesh is the reconstructed microstructure from a 3D set of slices. As an example of the output, take a look at http://www.dream3d.io. The image at the top of the page, are you saying that this algorithm would give me a mesh where each of the colored regions could be extracted (through the BoundaryCells) property? Currently we do the meshing, then run a laplacian smoothing algorithm to make it look better but the results are sub par at best most of the time.

What branch of VTK do I need to be on to try out this new algorithm?

And does there happen to be a stand-alone version of the algorithm by chance?

Thanks for this contribution
Mike Jackson

6

Thanks for the questions and feedback. I see that the filter does not support some of your workflows well. I will extend the API in the current filter next week. In the future, we are also considering creating a vtkAxisModeller which would essentially be a geometry/modelling kernel, with lots of advanced “atlas” functionality, but that won’t happen for a while. Here’s some additional elaboration if you are interested.

-------- More information -----

  • As Derk indicated, the current API is “.SetOutputStyleToSelected(), .InitializeSelectedLabelsList(), and AddSelectedLabel(label)” to extract labels/regions one-by-one (or to extract a subset of labels). I think another output style, maybe something like .SetOutputStyleToPartitioned(), to produce a grouped dataset (i.e., vtkPartitionedDataSet) would work. Each label would be placed into a separate vtkPolyData contained within an output vtkPartitionedDataSet. Of course to provide surface closure, shared points and triangles would be duplicated (does this make sense? comments?) Then it’d be quick and easy to extract vtkPolyData one by one and process them. I think Derk’s idea of writing a 3MF writer is of interest to me, but since I am not yet familiar with this, and to save some time, does anyone have any code that we can reuse?

  • Credit for the original algorithm goes to Sarah Frisken "SurfaceNets for Multi-Label Segmentations with Preservation of Sharp Boundaries” (see this paper for algorithm concepts and details). At one time the algorithm was patented, but now that it’s off patent we wanted to include it into VTK and make it really fast. We adapted some algorithmic concepts from Flying Edges, and mixed in vtkSMPTools for threading; the resulting performance is very good IMHO (best if you use VTK_SMP_IMPLEMENTATION_TYPE=TBB; typical speed ups range from 1-2 orders of magnitude depending on the data, number of labels, etc.). Note that vtkDiscreteFlyingEdges will be faster than vtkSurfaceNets3D for very small numbers of labels because FE assumes manifold isosurfaces, which gives it a performance boost over SN which assumes its dealing with non-manifold, joined surfaces. But for large numbers of labels SN is much faster, since FE makes a separate, complete pass for each label over the input annotation volume.

  • As far as future work goes, I occasionally see “bumps” in the surface near the intersection of multiple materials. I suspect that the smoothing stencils could be improved in some cases (it’s on my todo list). Also, the smoothing parameters are likely not tuned well so this will need to be experimented with. I also mentioned the vtkAxisModeller - it leverages the best kept secret of SNs, the BoundaryLabels 2-tuple cell data (indicating what regions are on either side of a triangle). You can probably imagine how this adjacency information could be used to navigate anatomical atlases. especially when combined with anatomical / ontology information.

  • If you are a glutton for punishment, the repository here has some standalone little programs for testing SNs, and comparing it against other algorithms. Dr. Frisken’s initial MultiMaterial SN code is also included (MM*) if you want to try that out. The mini-program RunSN.cxx lets you run and view results. Five datasets, some small, some quite large 2048^3 are also included. You will have to build (using the primitive CMakeLists.txt file) against a recent VTK version, either in VTK master or a version later than 9.2.20230420.

1 Like
7

Hi Will,

I think your suggestion regarding output makes perfect sense. Regarding 3MF, note that there is a library: https://lib3mf.readthedocs.io/en/master/ that may greatly accelerate things for those comfortable in C++.

It would be very convenient if vtkPolyData was able to export/import homogenous meshes as an array of vertices and an array of faces, i.e. a connectivity matrix. I was able to extract this information from an instance of vtkSurfaceNets3D (SN) to numpy arrays, but it seems both slow and clumsy:

Off  = SN.GetOutput().GetPolys().GetOffsetsArray()
Conn = SN.GetOutput().GetPolys().GetConnectivityArray()
Faces = np.empty([Off.GetNumberOfTuples()-1,3],dtype='int32')
for row in range(Off.GetNumberOfTuples()-1):
    off = int(Off.GetTuple(row)[0])
    for col in range(3):
        Faces[row,col] = int(Conn.GetTuple(off+col)[0])
        
Verts = numpy_support.vtk_to_numpy(SN.GetOutput().GetPoints().GetData())
BoundaryLabels = np.array(SN.GetOutput().GetCellData().GetArray("BoundaryLabels"))

Iterating over tuples to extract the connectivity data seems to be the slow step. Any advice on how to accelerate is welcome. The numpy convenience functions included with the python wrappers do not seem to work with the vtkTypeInt64Array that one gets with .GetConnectivityArray().

I did have another question regarding vtkSurfaceNets. As you pointed out, the algorithm passes through just once, no matter how many labels are selected for output. So I thought that cycling through labels as show here:

SN.SetOutputStyleToSelected()
SN.InitializeSelectedLabelsList()

for i in range(1,numclasses,1):
    print('Extracting SurfaceNet for '+keys.Label[I])
    if i>1:
        SN.DeleteSelectedLabel(i-1)
        
    SN.AddSelectedLabel(i)
    SN.Modified()
    SN.Update()
    stlWriter.SetInputData(SN.GetOutput())
    stlWriter.SetFileTypeToBinary()
    stlWriter.SetFileName('/path/to/file'+keys.Label[i]+'.stl')
    stlWriter.Write()

would not require re-calculating everything. But that is what seems to be happening:

2023-04-26 07:26:34.132 (1599.634s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2334  INFO| Executing Surface Nets 3D
2023-04-26 07:26:45.679 (1611.180s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2408  INFO| Extracted: 29666048 points, 29749004 quads
2023-04-26 07:26:45.686 (1611.188s) [         13BF1E9]   vtkSurfaceNets3D.cxx:1828  INFO| Smoothing output
2023-04-26 07:26:45.754 (1611.255s) [         13BF1E9]vtkConstrainedSmoothing:460   INFO| Executing constrained smoothing filter
2023-04-26 07:26:54.562 (1620.063s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2044  INFO| Transforming output mesh type to: 1
2023-04-26 07:26:55.650 (1621.151s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2459  INFO| Triangulated to produce: 59498008 triangles
2023-04-26 07:26:56.407 (1621.908s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2472  INFO| Selected: 5609800 cells
2023-04-26 07:26:59.466 (1624.967s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2334  INFO| Executing Surface Nets 3D
2023-04-26 07:27:10.887 (1636.389s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2408  INFO| Extracted: 29666048 points, 29749004 quads
2023-04-26 07:27:10.894 (1636.395s) [         13BF1E9]   vtkSurfaceNets3D.cxx:1828  INFO| Smoothing output
2023-04-26 07:27:10.976 (1636.477s) [         13BF1E9]vtkConstrainedSmoothing:460   INFO| Executing constrained smoothing filter
2023-04-26 07:27:19.734 (1645.235s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2044  INFO| Transforming output mesh type to: 1
2023-04-26 07:27:20.769 (1646.270s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2459  INFO| Triangulated to produce: 59498008 triangles
2023-04-26 07:27:21.188 (1646.689s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2472  INFO| Selected: 8748100 cells
2023-04-26 07:27:23.611 (1649.112s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2334  INFO| Executing Surface Nets 3D
2023-04-26 07:27:33.439 (1658.940s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2408  INFO| Extracted: 29666048 points, 29749004 quads
2023-04-26 07:27:33.443 (1658.944s) [         13BF1E9]   vtkSurfaceNets3D.cxx:1828  INFO| Smoothing output
2023-04-26 07:27:33.476 (1658.977s) [         13BF1E9]vtkConstrainedSmoothing:460   INFO| Executing constrained smoothing filter
2023-04-26 07:27:41.872 (1667.373s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2044  INFO| Transforming output mesh type to: 1
2023-04-26 07:27:42.829 (1668.330s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2459  INFO| Triangulated to produce: 59498008 triangles
2023-04-26 07:27:43.146 (1668.647s) [         13BF1E9]   vtkSurfaceNets3D.cxx:2472  INFO| Selected: 2958346 cells

I am still new to this entire pipeline thinking in VTK, so any pointers as to why this may be happening is welcome.

8

Simon,

thank you for including your code. I am afraid I am not that comfortable in C++ , so I’ll just have to work it out in python. So far, it’s been slow, as you say, but getting there…

9

I like that idea, that would save us some time and be very convenient for our workflow. Some of the cells will indeed be duplicated so that’s a small drawback. But on the other hand, sometime we would like to “hide” a specific part of the mesh linked to one label because the user is editing another label and want clear feedback on the 3D view. In that case, having only one mesh with shared triangles makes it difficult to hide just the specific part of the mesh properly. Because we have to choose one label or the other one for the shared triangles, sometimes we get big holes on the part of the mesh that is still visible (don’t know if that’s understandable, I can send some snapshots …). With multiple closed mesh it will be really easy to control the opacity/color/visibility for each part independently.

Once again, thanks for the work and dedication, looking forward to the next changes and glad to test that if needed :slight_smile:

10

Hi Will,

playing around some more, I realized that decimating the output mesh(es) might be useful in some circumstances. It would probably be a good idea to decimate using the default output style in which points at interior boundaries are shared. However, doing that, the color information encoded in BoundaryLabels seems to get lost. Would it make sense to offer decimation as an option in vtkSurfaceNets3D or is the VTK philosophy to keep these filter steps separate? If so, is there a way to make sure that vtkDecimatePro also acts on the BoundaryLabels produced by vtkSurfaceNets3D?

11

Decimation or remeshing is definitely needed as part of a region extraction workflow. In VTK we try to break these up into separate filters / classes so they can be more easily reused and combined into other workflows.

What I could see happening is a very specialized decimation algorithm that operates on “patches” (patches being the area of contact between two regions). The patches are bounded by a complex network of edges along which the patches join (sort of like the curves at which soap bubbles touch each other). Within each patch the adjacency 2-tuple remains unchanged so cell data would be preserved. To preserve compatibility, the decimation would have to reduce the patch edges and carry these changes into the connected patches etc. It would be fun to do, but tricky :slight_smile:

And if you are a glutton for punishment, I pushed a draft MR with a method to extract a region given its label. The method is:

void ExtractRegion(vtkIdType labelId, vtkPolyData *regionData,
bool boundaryFaces=false, bool cleanPoints=false);

(Note that the cleanPoints boolean has no effect yet, it is a placeholder for this ongoing work.) There is also a test TestSurfaceNets3D4.py which show how it works.

Running this extraction, I learned a lot. I am clearly seeing the effects of different regions in contact with another - I’m not happy with what I am seeing in some cases. (This shows up after smoothing is enabled, the contact patch edges shows up clearly.) I am pretty sure that the smoothing can be improved and I’ll work on it. In the mean time, it’s possible to extract the unsmoothed SurfaceNet regions (set smoothing iterations=0, or smoothing off) and then use a separate filter like vtkWindowedSincPolyDataFilter or vtkConstrainedSmoothingFilter to produce results.

Here’s the difference. When you smooth the SurfaceNet all together, the regions in contact with one another exert smoothing forces that cause slight “ridges” where the region edges join. The good news is that points and triangles are shared across boundaries. On the other hand, when you “extract” a region, and then smooth it, the region is decoupled from the SurfaceNet and you don’t see the effects of adjacent regions on the regions surface. The bad news is that the objects may not longer touch (after smoothing) and points and triangles are no longer shared. I’m reworking the API so that both workflows are possible. Also threading this stuff so the extraction process is as fast as the basic SurfaceNets algorithm.

Hopefully this makes some sense :slight_smile:

2 Likes
12

FYI- I just updated the MR which supports extraction of regions into separate vtkPolyData and/or vtkPartitionedDataSet (which is a collection of vtkPolyData, one per labeled region). It’s threaded and quite fast.

There are two ways to do this: by setting the output style and optionally using AddSelectedLabel(); or by invoking ExtractRegion() or vtkExtractRegions() (these latter two methods do not cause filter reexecution). Note that there are options to extract just cells on the boundary of their region (i.e., touching the background label/region), and/or to clean the region’s points (i.e., points and cells are renumbered so that only points used by the regions cells are produced on output).

Look at the test TestSurfaceNets3D4.py so see how this works.

13

Will,

I came across this post when I was looking for a python library that went from VTK to 3MF. Still not there and since I got my management lobotomy, I’m not up to writing something myself.

But good news. You can save as VRML (*.wrl) and Stratasys GrabCAD can read that format. We are a Stratasys reseller and also an Ansys reseller, so my tool converts simulation results through VTK to VRML and we 3D Print color objects that way.

I just use the following in Python:

myplt = pyvista.Plotter()
#add what I want to plot
myplt.export_vrml(“myfile.wrl”)

You can save as VTK in your program then test in Paraview as well.

Hope this helps.

Eric
Principal & Co-Owner
PADT, Inc

14

Hello,
where can I find the input data set you mention?
Best, Andreas

15

Check out this link to a preprint paper. It has some information in the appendix about obtaining the data and reproducing results:
https://www.researchgate.net/publication/372871904_A_High-Performance_SurfaceNets_Discrete_Isocontouring_Algorithm

Or if you want to go directly to the data and build instructions:
https://drive.google.com/drive/folders/10Rf7TjP9ZXAGRPbwJa6CKoCHv4DygicU

16

Thanks @will.schroeder! This is really nice. The only thing I would propose to change is that the BoundaryLabels does not sort the labels. This makes it impossible to directly extract the triangles with correct orientation. Instead if you simply comment out this line the orientation of triangles/quads is preserved and no post-processing to generate consistent normals is necessary.

I created an issue here: https://gitlab.kitware.com/vtk/vtk/-/issues/19156

17

Just to close the loop on discourse: I’ve commented on the issue. I very much appreciate the insightful feedback, thank you.

18

Hi Will. Is there a public license for the code, and for the data sets?

19

I’m not sure I know what you are asking for, but here goes. vtkSurfaceNets3D uses the standard 3-clause BSD license and is part of the VTK distribution. I believe that the paper describes the sources for the data, you’d have to go to these sources to track down the licenses. If I am not answering your question, please follow up and I’ll do what I can to help.