Here are some preliminary functions. If you just want to read the document, then you can ignore this first cell.

#preliminary functions

def hyperplane_arrangement_to_fan(hyperplane_arrangement):
    '''
    Input: Linear hyperplane arrangement
    Output: Maximal cones of the induced fan. Each maximal cone is the closure of a connected component
    of R^d after removal of the hyperplane arrangement
    '''
    n = len(hyperplane_arrangement)
    p = list(powerset(range(n)))
    d = hyperplane_arrangement[0].ambient_dim()
    cones = []
    for entry in p:
        signs = [1 if i in entry else -1 for i in range(n)]
        inequalities = []
        for i in range(n):
            hyperplane = hyperplane_arrangement[i]
            A = hyperplane.Hrepresentation()[0].A()
            ieq = [0]+[signs[i]*a for a in A]
            inequalities.append(ieq)
        C = Polyhedron(ieqs = inequalities)
        if C.dimension() == d:
            cones.append(C)
    return cones

def intersecting_edges(polytope, u):
    '''
    Input: A polytope and a normal vector u
    Output: Set of edges of polytope that has nonempty intersection with the hyperplane orthogonal to u
    '''
    u_perp = Polyhedron(eqns= [ [0]+list(u) ])
    intersecting_edges = []
    for edge in polytope.faces(1):
        if edge.as_polyhedron().intersection(u_perp).dimension()>-1:
            intersecting_edges.append(edge)
    return intersecting_edges

Computing the intersection body of the cube.

We consider the \(3\)-dimensional cube \([-1,1]^3\) and compute the radial function of its intersection body step by step.

P = polytopes.cube()

The combinatorial type of the intersection of \(P\) with a hyperplane through the origin is fixed for each region of the normal fan of the following zonotope Z. You can visualize this fan by looking at the dual polytope.

Z = sum([Polyhedron(vertices = [-vector(v),vector(v)]) for v in P.vertices()])
dualZ = Z.polar()
dualZ.show(frame=False, point = False)

We now construct the associated hyperplane arrangement \(H\).

H = list(set([Polyhedron(eqns = [[0]+list(vertex)]) for vertex in P.vertices()]))
Fan = hyperplane_arrangement_to_fan(H)

Let us focus now on one of the open chambers of this hyperplane arrangement, and compute the radial function in here.

We sample a vector in the interior of this region.

C = Fan[1]
sample_vector = sum([vector(ray) for ray in C.rays()])
sample_vector

(-2, 2, 2)

We now collect all facets of \(P\) that intersect the hyperplane perp orthogonal to sample_vector. Note that in this case \(Q = P\cap\)perp is a hexagon: indeed we are intersecting \(6\) facets of \(P\).

perp = Polyhedron(eqns= [ [0]+list(sample_vector) ])

intersecting_facets = []
if P.intersection(perp).dimension()==2:
    for facet in P.facets():
        if facet.as_polyhedron().intersection(perp).dimension() == 1:
            intersecting_facets.append(facet)
intersecting_facets

[A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices]

The vertices of such a hexagon arise as the intersection of perp with some edges of the above facets of \(P\). We now compute these edges.

intersecting_edges_of_facets = []
for facet in intersecting_facets:
    edges = intersecting_edges(facet.as_polyhedron(),sample_vector)
    intersecting_edges_of_facets.append(edges)
intersecting_edges_of_facets

[[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
  A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
 [A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
  A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
 [A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
  A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
 [A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
  A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
 [A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
  A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
 [A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
  A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices]]

Let \(x \in C\) arbitrary. By the theory we know that \(x^\perp\) intersects \(P\) in the set of edges in intersecting_edges_of_facets. Therefore, for each such edge \(\hbox{conv}(a,b)\) we write the corresponding vertex \(v\) of \(P\cap x^\perp\) as:

\[v = \frac{(b\cdot x)\,a - (a\cdot x)\,b}{(b-a)\cdot x} \,.\]

x = vector(var('x1,x2,x3'))
x1, x2, x3 = x
v = []
for edges_of_facet in intersecting_edges_of_facets:
    intersecting_points_of_facets = []
    for edge in edges_of_facet:
        [a,b] = list(edge.vertices())
        a = vector(a)
        b = vector(b)
        intersection_point = ((b.inner_product(x)*a - a.inner_product(x)*b )/ (b-a).inner_product(x))
        intersecting_points_of_facets.append(intersection_point)
    v.append(intersecting_points_of_facets)
v

[[((x2 - x3)/x1, -1, 1), (-1, -1, (x1 + x2)/x3)],
 [(-(x2 - x3)/x1, 1, -1), (-1, (x1 + x3)/x2, -1)],
 [(-1, (x1 + x3)/x2, -1), (-1, -1, (x1 + x2)/x3)],
 [(1, -(x1 + x3)/x2, 1), ((x2 - x3)/x1, -1, 1)],
 [(1, 1, -(x1 + x2)/x3), (-(x2 - x3)/x1, 1, -1)],
 [(1, -(x1 + x3)/x2, 1), (1, 1, -(x1 + x2)/x3)]]

We are now ready to compute the volume of the hexagon. This will be done by subdividing the polygon into simplices, and by computing the volume of each simplex as a determinant.

determinants = []
for simplex in v:
    M = matrix(
        simplex + [x]
    )
    det_M = M.determinant()
    if det:
        sgn = sign(det_M(x1 = sample_vector[0], x2 = sample_vector[1], x3 = sample_vector[2]))
        determinants.append([sgn,M.determinant()])
determinants

[[-1, -x1*((x1 + x2)/x3 - 1) - x2 - (x2 - x3)*((x1 + x2)*x2/x3 + x3)/x1 - x3],
 [1, x1*((x1 + x3)/x2 - 1) - (x2 + (x1 + x3)*x3/x2)*(x2 - x3)/x1 + x2 + x3],
 [1,
  x1*((x1 + x2)*(x1 + x3)/(x2*x3) - 1) + x2 + (x1 + x2)*x2/x3 + (x1 + x3)*x3/x2 + x3],
 [-1, -x1*((x1 + x3)/x2 - 1) + (x2 + (x1 + x3)*x3/x2)*(x2 - x3)/x1 - x2 - x3],
 [1, x1*((x1 + x2)/x3 - 1) + x2 + (x2 - x3)*((x1 + x2)*x2/x3 + x3)/x1 + x3],
 [1,
  x1*((x1 + x2)*(x1 + x3)/(x2*x3) - 1) + x2 + (x1 + x2)*x2/x3 + (x1 + x3)*x3/x2 + x3]]

Finally, by adding all the determinants with a correction factor, we obtain the radial function of \(IP\) restricted to the region \(C\).

correcting_factor = 1/(factorial(P.dimension()-1) * sum([x_i^2 for x_i in x]))
rho = correcting_factor*sum([sng*det for [sng,det] in determinants])
rho.factor()

(x1^2 + 2*x1*x2 + x2^2 + 2*x1*x3 - 2*x2*x3 + x3^2)/(x1*x2*x3)

From the radial function we can compute the algebraic boundary of \(IP\). In \(C\) this is the zero set of the following polynomial.

p, q = rho.numerator_denominator()
alg_boundary = q-p
alg_boundary

x1*x2*x3 - x1^2 - 2*x1*(x2 + x3) - x2^2 + 2*x2*x3 - x3^2