2

I have two convex hulls. Let's assume they are given as scipy.spatial.ConvexHulls

import numpy as np


points1 = np.random.rand((10, 3))
points2 = np.random.rand((10, 3))

hull1 = ConvexHull(points1)
hull2 = ConvexHull(points2)

I would like the convex hull that is the intersection of these two convex hulls, but could not find a built in method to do this.

I assume this can be done manually somehow by using scipy.spatial.HalfspaceIntersection by using half spaces defined by hull1 to cut off hull2, but still having trouble doing it, and can't believe this is not already implemented somewhere.

Note that I don't mind if scipy is not used.

martineau
  • 119,623
  • 25
  • 170
  • 301
Gulzar
  • 23,452
  • 27
  • 113
  • 201
  • Convex hulls are convex polyhedrons. Look for a utility to intersect two convex polyhedrons, or just two polyhedrons. –  Dec 22 '21 at 17:40

2 Answers2

4

I would try pycddlib, which implements the double description of polyhedra. The double description of a polyhedron is:

  • V-description: description by vertices
  • H-description: description by system of linear inequalities ("H" for "hyperplanes")

You probably have the vertices of your two convex polyhedra. Convert to the H-descriptions, then combine the two systems of linear inequalities, and then convert to the V-representation.

Here is an example.

import numpy as np
import pyvista as pv
import cdd as pcdd
from scipy.spatial import ConvexHull

# take one cube
cube1 = pv.Cube()
# take the same cube but translate it 
cube2 = pv.Cube() 
cube2.translate((0.5, 0.5, 0.5))

# plot 
pltr = pv.Plotter(window_size=[512,512])
pltr.add_mesh(cube1)
pltr.add_mesh(cube2)
pltr.show()

enter image description here

# I don't know why, but there are duplicates in the PyVista cubes;
# here are the vertices of each cube, without duplicates
pts1 = cube1.points[0:8, :]
pts2 = cube2.points[0:8, :]

# make the V-representation of the first cube; you have to prepend
# with a column of ones
v1 = np.column_stack((np.ones(8), pts1))
mat = pcdd.Matrix(v1, number_type='fraction') # use fractions if possible
mat.rep_type = pcdd.RepType.GENERATOR
poly1 = pcdd.Polyhedron(mat)

# make the V-representation of the second cube; you have to prepend
# with a column of ones
v2 = np.column_stack((np.ones(8), pts2))
mat = pcdd.Matrix(v2, number_type='fraction')
mat.rep_type = pcdd.RepType.GENERATOR
poly2 = pcdd.Polyhedron(mat)

# H-representation of the first cube
h1 = poly1.get_inequalities()

# H-representation of the second cube
h2 = poly2.get_inequalities()

# join the two sets of linear inequalities; this will give the intersection
hintersection = np.vstack((h1, h2))

# make the V-representation of the intersection
mat = pcdd.Matrix(hintersection, number_type='fraction')
mat.rep_type = pcdd.RepType.INEQUALITY
polyintersection = pcdd.Polyhedron(mat)

# get the vertices; they are given in a matrix prepended by a column of ones
vintersection = polyintersection.get_generators()

# get rid of the column of ones
ptsintersection = np.array([
    vintersection[i][1:4] for i in range(8)    
])

# these are the vertices of the intersection; it remains to take
# the convex hull
ConvexHull(ptsintersection)
Stéphane Laurent
  • 75,186
  • 15
  • 119
  • 225
0

I have searched for a suitable solution. Tried pycddlib variant proposed here, but it failed because of numerical inconsistency. Also I have tried pyvista library for mesh intersections, but this solution was not exact enough and also crushed in some cases. I have succeeded with Pymeshlab library:

https://pymeshlab.readthedocs.io/en/latest/index.html

Here you can easily create a mesh from a list of points and simplices, then there is a filter meshing_re_orient_faces_coherentely which helps to orient the faces properly. Then you should take some face and translate its center of mass a little bit in the direction of normal (I took 1e-4 of normal length) and check if this point lies into the hull (https://stackoverflow.com/a/55484399/22359849 this answer helps to compute the signed distance, negative is inside). If the point is out of the hull you should use meshing_invert_face_orientation filter. This filter should automatically set the proper orientation, but for me it did not work right. Now, when all the faces are oriented properly your mesh is ready. You need to add two properly oriented meshes to a single MeshSet and apply generate_boolean_intersection filter. The new mesh could be extracted, it would be a answer.

Arsenii Raiko
  • 1
  • 2