Speed up collision detection / Bounding Box Tree for Point Cloud

Hello,

I’ve got a quick question regarding the dolfin implementation of BoundingBoxTree. In the C++ interface (see Bitbucket) there seems to be a possibility to generate a bounding box tree for point clouds… However this interface does not seem to be available for python, whenever I try using

from dolfin import BoundingBoxTree, Point
bbtree = BoundingBoxTree()
bbtree.build(Point(0.,0.), 2)

or

bbtree.build([Point(0.,0.)], 2)

I get the error, that this is not supported.
Is there any way that this can be achieved?

Background: I am trying to speed up a collision / intersection test, where I currently use the following code

bbtree = BoundingBoxTree()
bbtree.build(mesh)
cells = mesh.cells()
cells_flat = cells.flatten().tolist()
counter = Counter(cells_flat)
for j in range(coordinates.shape[0]):
        x = Point(coordinates[j])
        intersections = len(bbtree.compute_entity_collisions(x))
        occurences = counter[j]

        if intersections > occurences:
            self_intersections = True
            break

So if anyone has suggestions how to speed up this test, this would also be very helpful.

Thanks a lot in advance.

The following runs for me ( '2019.2.0.dev0'):

from dolfin import BoundingBoxTree, Point 
bbtree = BoundingBoxTree() 
bbtree.build([Point(0.,0.)], 2)  
Computed bounding box tree with 1 nodes for 1 points.

As your code is not complete (not a minimal example illustrating what you want to speed up), I can’t really help you to speed up the code any further.

1 Like

Okay, yes my example was a bit incomplete, let me be more thorough now.

I am working on shape optimization problems. Due to the deformations of my mesh, there is the possibility that some mesh elements may intersect with others, which I want to check for so that I can avoid it.

The method I thought would be efficient is the following. I go over all nodes of my mesh and compute the entity collisions with cells of the very same mesh. This will return all the cells which the node is part of the vertex set as well as the cells where it lies in the interior / on the boundary (and is not a vertex of the cell). As I can precompute the number of cells which have node i as a vertex, I can then just compare these quantities - and if the collision test returns more cells than the precomputed quantity, there is an intersection, which is what I want to detect. A MWE for this is the following:

from dolfin import *
from collections import Counter

# initialization
mesh = UnitSquareMesh(10, 10)
bbtree = mesh.bounding_box_tree()
coordinates = mesh.coordinates()

# precompute number of cells, for which a given node is part of the vertices
cells = mesh.cells()
cells_flat = cells.flatten().tolist()
counter = Counter(cells_flat)

self_intersections = False
for i in range(coordinates.shape[0]):
    # Compute the number of collisions with cells for node i
    x = Point(coordinates[i])
    collision_cells = bbtree.compute_entity_collisions(x)
    no_collisions = len(collision_cells)

    # get the precomputed value
    no_vertex_cells = counter[i]

    if no_collisions > no_vertex_cells:
        self_intersections = True
        break
    # sanity check
    elif no_collisions < no_vertex_cells:
        raise Exception("Error with the collision detection")

Obviously, as I am using UnitSquareMesh and have no deformations here, no intersections will be found. But in the case where the mesh would be obtained by a deformation, there could be intersections and the test would show this.

I have noticed that my implementation of this check can still take a rather long time, which is why I wanted to speed this up. And one idea for this was to also use a bounding box tree for the point cloud, which I could not get to work. I am using fenics 2019.1.0 via conda-forge, and get the following behavior when I add these lines to my MWE

mytree = BoundingBoxTree()
mytree.build([Point(0.0, 0.0)], 2)
Traceback (most recent call last):
  File "xxx/miniconda3/envs/fenics19/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3441, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-3-3cad0bfbbc2d>", line 2, in <module>
    mytree.build([Point(0.0, 0.0)], 2)
TypeError: build(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfin.cpp.geometry.BoundingBoxTree, arg0: dolfin.cpp.mesh.Mesh) -> None
    2. (self: dolfin.cpp.geometry.BoundingBoxTree, arg0: dolfin.cpp.mesh.Mesh, arg1: int) -> None
Invoked with: <dolfin.cpp.geometry.BoundingBoxTree object at 0x7f4677d5f7f0>, [<dolfin.cpp.geometry.Point object at 0x7f4678168930>], 2

So my question is: Is there any way to speed up this detection of intersections of meshes. And if the answer involves a bonding box tree for a point cloud: Why does this not work with version 2019.1.0?

Thanks a lot for your time.

It is simply because the function was added to the python interface after the 2019.1.0 release: fenics-project / DOLFIN / commit / 8f32523665e5 — Bitbucket

I agree that using a boundingboxtree of pointclouds is the easiest way to do such a computation.

There are plenty of smoothers (Riesz representations and similar methods) to avoid such behavior, see for instance: https://arxiv.org/pdf/2001.10058.pdf (Section 5.1.1-> onwards plus the references).

Okay, thanks a lot.
I’ve implemented an interface for the point cloud bounding box tree, but apparently it gives wrong results when using compute_entity_collisions.

I have resolved my problem by writing the test in C++. For reference, my implementation can be found here: cashocs/geometry.py at 3d0c2ebcf80c9d9bbe5ea2ddb50e47d486653d04 · sblauth/cashocs · GitHub

Thanks a lot for your help, though.