Speed up collision detection / Bounding Box Tree for Point Cloud

1

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.

2

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
3

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.

4

It is simply because the function was added to the python interface after the 2019.1.0 release: 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).

5

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.

6

Hi plugged,

Thank you for providing the C++ code. I am using it but is seems somehow not working after I update mesh:

For example, I have mesh before deformation and the mesh after deformation using;
predeform = CollisionCounter.compute_collisions(mesh)
ALE.move(mesh, displacement)
postdeform=CollisionCounter.compute_collisions(mesh)
or I can do:
postdeform=CollisionCounter.compute_collisions(mesh)
along without the counting the predeform contacts.

The two postdeforms I get are somehow different. Could you please let me know why? Or it would just be helpful to have an example on how to use the code.

Thank you so much！

7

Hi,

here is an annotated MWE (you can replace the mesh with your mesh and the displacement with your displacement)

from fenics import *
import numpy as np
import collections

# define the mesh
mesh = ...

# count number of times a vertex is used in an element
cells = mesh.cells()
flat_cells = cells.flatten().tolist()
cell_counter = collections.Counter(flat_cells)
occurrences = np.array([cell_counter[i] for i in range(mesh.num_vertices())])
# occurrences[i] = j   means that vertex i is part of j elements

# create displacement and move mesh
VCG = VectorFunctionSpace(mesh, "CG", 1)
displacement = Function(VCG)
displacement.vector()[:] = ...
ALE.move(mesh, displacement)

# CollisionCounter.compute_collisions computes the amount of elements each mesh vertex
# collides with, i.e., collisions[i] = k  means that vertex i is part of k elements
self_intersections = False
collisions = CollisionCounter.compute_collisions(mesh)
if not (collisions == occurrences).all():
    self_intersections = True
    
# therefore, if collisions[i] == occurrences[i] for all i, there are no intersections,
# otherwise there are.

I use this code here cashocs/geometry.py at 3d0c2ebcf80c9d9bbe5ea2ddb50e47d486653d04 · sblauth/cashocs · GitHub

8

Hi plugged,

Thank you for the reply. But the problem is that I am using the code to run a loop in which the mesh is deformed in every step so I need to update the nodes in contact using like:

while displacement>1e-5:
…
ALE.move(mesh, displacement)
postdeform=CollisionCounter.compute_collisions(mesh)

in which the command CollisionCounter.compute_collisions(mesh) is used more than once. How can I get the correct result without creating a new mesh every time I deform the mesh？

Many thanks!

9

The code should work fine even when you deform your mesh. Are you “forgetting” to recompute / update the bounding box tree of the mesh? You can do so as follows:

mesh = ...
bbtree = mesh.bounding_box_tree()
...
while True:
    ...
    ALE.move(mesh, displacement)
    bbtree.build(mesh)
    postdeform = CollisionCounter.compute_collisions(mesh)

If this does not work, please provide a minimal working example, so that I can help with debugging your problem.

10

Thank you, it is working now :grinning: