added minkowski sum and intersection

pycapacity/objects.py

@@ -118,6 +118,29 @@ def find_from_point_cloud(self, points):
         self.H, self.d = vertex_to_hspace(points)
         self.vertices, self.face_indices = hspace_to_vertex(self.H,self.d)
 
+    # minkowski sum of two polytopes
+    def __add__(self, p):
+        if self.vertices is None:
+            self.find_vertices()
+        if p.vertices is None:
+            p.find_vertices()   
+        vertices_sum = []
+        for v1 in self.vertices.T:
+            for v2 in p.vertices.T:
+                vertices_sum.append(v1+v2)
+        P_sum = Polytope()
+        P_sum.find_from_point_cloud(points=np.array(vertices_sum).T)
+        return P_sum
+
+    # intersecting two polytopes
+    def __sub__(self, p):
+        if self.H is None or self.d is None:
+            self.find_halfplanes()
+        if p.H is None or p.d is None:
+            p.find_halfplanes()
+        H_int = np.vstack((self.H,p.H))
+        d_int = np.vstack((self.d,p.d))
+        return Polytope(H=H_int,d=d_int)
 
 
 

tests/testing_objects.py

@@ -52,4 +52,23 @@ def test_polytope_from_vertex_to_face():
     p = Polytope(vertices=points)
     p.find_faces()
     # check if all faces in the face representation are in the vertex representation
-    assert np.all(np.isin(p.faces,p.vertices))
+    assert np.all(np.isin(p.faces,p.vertices))
+
+# test minowski sum of two polytopes
+# testing for a cube
+def test_polytope_minkowski_sum():
+    p1 = Polytope(H = np.vstack((np.eye(3),-np.eye(3))), d = np.vstack((np.ones((3,1)),np.ones((3,1)))))
+    p2 = Polytope(H = np.vstack((np.eye(3),-np.eye(3))), d = np.vstack((2*np.ones((3,1)),2*np.ones((3,1)))))
+    p_sum = p1 + p2
+    assert np.all(p_sum.H@p1.vertices - p_sum.d[:,None]  < 1e-5) and  np.all(p_sum.H@p2.vertices - p_sum.d[:,None]  < 1e-5)
+
+# test intersection of two polytopes
+# testing for a cube
+def test_polytope_intersection():
+    p1 = Polytope(H = np.vstack((np.eye(3),-np.eye(3))), d = np.vstack((np.ones((3,1)),np.ones((3,1)))))
+    p2 = Polytope(H = np.vstack((np.eye(3),-np.eye(3))), d = np.vstack((2*np.ones((3,1)),2*np.ones((3,1)))))
+    p_int = p1 - p2
+    p1.find_halfplanes()
+    p2.find_halfplanes()
+    p_int.find_vertices()
+    assert np.all(p1.H@p_int.vertices - p1.d[:,None]  < 1e-5) and  np.all(p2.H@p_int.vertices - p2.d[:,None]  < 1e-5)