Clipping data with manually constructed closed surface polyData.

Hi everyone,

I have the following task I want to accomplish. I have results I would like to clip to a specified shape. The shape in my case is constructed via number of 2D polygon surfaces, that combined define an enclosed volume. The results are larger compared to the defined volume.

My initial try was to use implicit data set, in which I fed in the appended 2D polygon surfaces (see attached code example). The issue is that the clipping picks up only the slice of the data, which coincides with the boundaries I am trying to impose. Cutting with infinite planes is not an option, because the region can be concave (like L-shaped).

I suspect that the reason why clipping does not show the volume is that there are no reasonable
function values that would allow to differentiate inside/outside. Consequently, I am currently thinking about coding-up a function that would be able to provide values both inside/outside of this region. This seems not a straightforward task though.

Therefore I wanted to ask:
*) Does anyone now an existing took within VTK that would allow to accomplish what I am trying to do?

Best,
Ugis

CMakeLists.txt (629 Bytes)
minEx_clipBasedOnPolyData.cxx (8.0 KB)

If you can combine the parallel contours into a closed surface then you can use Boolean mesh operations implemented in vtkbool to cut the other surface. VTK has built-in Boolean mesh operations, too, but they often fail - even if the inputs are valid, well-behaving meshes.

1 Like

@Charles_Gueunet @Tiffany_Chhim

Hi @lassoan,

thanks a lot for the reference to vtkbool, it could be helpful. But for that to work, I would actually need to convert my boundary “shell” into volume mesh first, so I can use “intersection” operation. I have been looking in that direction, but have not found a way yet.

So sub-question would be - do you (or anyone else) know how to define volume based on encompassing 2D polygon surfaces in VTK and mesh it (including concave volumes)? I tried to pass all the points of the surface to vtkDelaunay3D filter, but that does not work (as it is explicitly stated in the documentation) for concave volumes. In other words, it produces convex encompassing volume. In meshing tools (such as Gmsh) one would define a surface loop, but looking through VTK classes, it seems for me that such concept does not exist in VTK.

Thanks a lot for the input!
/Ugis

Hello @ugis ,

The VTK boolean operation expect two watertight sufraces, so you do not need to compute a volume mesh.
If you can deal with GPLv3 code, in my experience the best boolean operation are thoses of the vespa project, based on CGAL. You may also find tools to remesh, smooth and deform your surfaces (as long as they are triangulated first).

Ok, I see. I guess my intuition was too much “CAD-oriented”. I will give the boolean operations a try and see how far I get.

The 2D planar contours do not fully specify a closed 3D surface. When you slice concave objects then they show up as several contours and there is no way of knowing which one to connect to which (you can make some guesses based on spatial proximity, but there is no definitive answer), the intersection may have holes in it, there may be no continuation of a contour on the next slice, etc.

For this reason, I would recommend to avoid representing concave structures by planar contour sets and represent them using an image or closed surface instead. If you have absolutely no other choice (e.g., you already have lots of data set with planar contour set that you need to use) then you can find VTK-based solutions for specific application domains that work reasonably well.

For example, for medical images we developed a solution for reconstructing smooth 3D surfaces from planar contour sets in DICOM RT Structure Sets. The method can handle holes, branching, generate smooth end-capping, etc.) - it is available in SlicerRT extension of 3D Slicer and it is automatically applied when you import DICOM RT studies. It may work for other domains as well.

Right, I understand. I think that luckily in my case the problem is simpler, because I have a closed surface without holes (perhaps I was not perfectly clear about that before). So then based on this

my hope was that boolean operations would give me what I need.

I have adjusted my example:

  1. I create L-shaped closed volume by appending square polygon data to each other;
  2. I create an example data set, which is mostly larger than the closed surface;
  3. I slice the data set and refine the result, which creates planar poly data;
  4. I try to apply boolean operations on it, to limit the data extent to within the closed volume.
    CMakeLists.txt (696 Bytes)
    minEx_clipBasedOnPolyData.cxx (11.5 KB)

As mentioned in other responses before, the VTK native filter does not behave too well, while vtkbool produced much nicer results. However, no boolean operation mode gave me what I hoped for. I hoped for planar output within the defined volume. The closest I got was another closed volume basically clipped by the plane I wanted to clip.

Am I doing something wrong with the boolean operation or am I misunderstanding the way it is supposed to work? Perhaps boolean operation does not work if one of the surfaces is only planar?

1 Like

Unfortunately, even if your closed surface does not have holes, many of your slice intersections may still have holes in them:

So, in general, reconstructing a non-convex object from a set of parallel contours is a complex task, because even for very simple shapes you in some planes you usually get multiple contours, some of them with holes, and the contours can branch and merge between the planes.

Hi Andras,

Thanks a lot for your efforts trying to clarify the reconstruction problem for me, I really appreciate it and also your work towards finding solution for this problem. And apologies if I am being unable to convey my problem appropriately clearly.

I have been thinking quite a bit and I do not see where in my case the reconstruction would come in. My interest lies in two scenarios (both, I think, is only clipping):

  1. Clipping planar results (according to the second example I posted), in which as inputs I have one continuous plane (no holes) and one closed surface. I am trying to clip the plane to be within the closed surface;
  2. Clipping iso-surface results, in which instead of plane, I have an isosurface (this one could have holes, but I would not want to reconstruct it, but clip it). And I am trying to clip the isosurface to be only within the closed surface.

I have really hard time to understand, where the surface reconstruction would come in, because I construct the closed clipping surface explicitly. Is there really some reconstruction going on, when I am appending my planar polyData to each other?

I am a bit confused about why boolean operations did not work out, because it seemed theoretically applicable, especially according to @Charles_Gueunet answer, that I do not need volume meshes, but only water-tight surfaces. Perhaps appending polyData to each other is not enough to construct a closed surface?

Thanks again for all the input!

Surface reconstruction is only needed if you want to reconstruct a closed surface from a set of planar contours, but it seems that you don’t need this.

I don’t know exactly what you mean by volume meshes. If you have a structured grid of small cubes a.k.a. 3D image a.k.a. image volume then it simplifies things, as you may perform Boolean operations in the image domain. If you have a volumetric mesh of various volumetric cells (e.g., tetrahedra, wedges, …) then it is probably better to extract the outer contour using VTK geometry filter because I am not aware of any Boolean operation implementations for VTK volumetric meshes.

Appending polydata does not create a closed surface. You need to at least use VTK clean polydata filter to merge coincident points and make sure that you did not have any faces inside the outer closed surface.

If you provide a few screenshots or sketches of the problems then it makes easier to understand what you are trying to achieve.

For me, C++ code samples are not that helpful, because I would need to create a project and build it and it may turn out that there are build errors, etc. so I would need to spend even more time with it. Therefore, I rather don’t even try to do it. If you provide standalone Python code snippet that I can just copy-paste into a Python console in a few seconds then I can try it without risking investing tens of minutes into something that may not end up being useful.

Hi again @lassoan

Thanks a lot about willingness to try my example. Sorry about the delay, I was kept busy with some other things, and it also took some time to get a VTK-Python environment with vtkBool up and running. But here comes hopefully explanatory sketch and also the python code.

So, let us look at the sketch/screenshots I made:


What I have first is the closed surface, constructed by appending vtkPolyData (see red L-shaped object). Then I have two data types I want to clip, planar data (upper row) and isosurface data (lower row). What I would like to do, is to use boolean operation to clip out the part of the data that lies within the surface I have defined (see “expected output”, two images to the right - wireframe of “bounding box” given only to further clarify that the data lies within the defined surface). However, the vtkBool (I have tried all operation modes) does not give what I want. So, the question(s) would be i) if I am misunderstanding the way vtkBool should be used and ii) is there another already implemented way in VTK to achieve this.

To generate the “expected output”, I did two passes of clipping by planes and then appending final result. This is something that would definitely work, but in general I can have much more complex shapes than the one I use in this example and the bookkeeping would become rather tedious.

Here is the python code (one for each data type) I used to generate the attempt captured in the illustration:
python_minEx_cliBasedOnPolyData_isosurface.py (12.3 KB)
python_minEx_cliBasedOnPolyData_planar.py (12.4 KB)
To switch over different views, just press any character key on the keyboard.