Find one point per pyramidal region

Question

Consider 3 dimensional space which has been partitioned by at least two planes which go through the origin. If there are n planes then the number of distinct pyramidal regions this creates is exactly 2 - n + n^2 as long as they are in "general position", a term I make specific below.

The challenge is write code that will return one point per pyramidal region. That is to return 2 - n + n^2 points in 3d space, each of which is in a different pyramidal region.

Notice that a plane through the origin can be defined by a single vector starting at the origin which is orthogonal to the plane.

Your code should take a list of planes in as input, each plane being described by a vector from the origin which is orthogonal to it. The vector itself is just represented as a space separated list of three numbers. Your should output a list of 2 - n + n^2 points. You can assume the planes are in general position, that is that no three of the vectors defining the planes are coplanar.

I will verify that the points are really all in different pyramidal regions using the following code.

Let us assume the variable planes contains an array of vectors, each one of which defines the plane. For example, let n = 4, the planes could equal

[[ 0.44060338 -0.55711491 -0.70391167]
 [ 0.15640806  0.89972402  0.40747172]
 [-0.48916566  0.86280671  0.12759912]
 [ 0.25920378 -0.94828262  0.18323068]]

We can then choose a number of points and put them in a 2d array called testpoints. To determine if they are all in distinct pyramidal regions we can call the following Python function which just looks to see which side of each plane the points are on and checks that each point is on at least one different side from every other point.

def all_in_distinct_pyramidal_regions(testpoints, planes):
    signs = np.sign(np.inner(testpoints, hyperplanes))
    return (len(set(map(tuple,signs)))==len(signs))

You can create some random test planes by first sampling points on a sphere and then using those points to define the hyperplanes. Here is a Python function that does this sampling.

import numpy as np
def points_on_sphere(N, norm=np.random.normal):
    """
    http://en.wikipedia.org/wiki/N-sphere#Generating_random_points
    """
    normal_deviates = norm(size=(N, 3))
    radius = np.sqrt((normal_deviates ** 2).sum(axis=0))
    points = normal_deviates / radius
    return points

You can use any language you choose for which there is a free and easy way to run the code in Linux. You can use any standard libraries of that language. I will however want to test your code so please give full instructions for how to run it.

Code to produce 10 random points on a sphere can be run using http://ideone.com/cuaHy5 .

Here is a picture of the 4 planes described above along with 2 - 4 + 16 = 14 points, one in each of the pyramidal regions created by the planes.

A valid output for this instance would therefore be the 14 points.

0.96716439  0.21970818  0.12775507
-0.84226594 -0.37162037 -0.39049504
0.60145419  0.58557551 -0.54346497
-0.4850864   0.26433373 -0.83355796
-0.69380885 -0.62124759  0.36425364
0.09880807 -0.09275638  0.99077405
0.19131628 -0.97946254 -0.06364919
-0.98100047  0.17593773 -0.08175574
-0.76331009  0.63554838 -0.11591362
0.89104146  0.34728895  0.29229352
-0.61696755 -0.26411912  0.74134482
0.18382632 -0.03398322 -0.98237112
-0.44983933 -0.18244731  0.87427545
0.8391045   0.18244157 -0.51246338

Answer 1

Python, 423 376 353

S=sum
from numpy import*
C=cross
e=1e-15
def P(N,I,M):
 if M:v=M[0];L=[C(n,v)for m in N for n in-m,m if all(dot(N,C(n,v))<e)];[P(N+[u*v],[i for i in I if all(dot(N+[u*v],i)<e)]+L,M[1:])for u in-1,1]if L else P(N,I,M[1:])
 else:
    for u in-1,1:print list(u*S(I))
N=map(array,input())
a,b=N[:2]
c=C(a,b)
for u in-1,1:P([a,u*b],[c,-c,C(a,c)+u*C(b,c)],N[2:])

Reads a list of comma-separated triplets from STDIN and writes a list of triplets to STDOUT. Assumes at least two planes are given.

Example

$ python cone.py <<< [0.44060338,-0.55711491,-0.70391167],[0.15640806,0.89972402,0.40747172],[-0.48916566,0.86280671,0.12759912],[0.25920378,-0.94828262,0.18323068]
[0.50892804716566153, 0.17353382779474077, -1.1968194626314141]
[-0.50892804716566153, -0.17353382779474077, 1.1968194626314141]
[0.1422280331828244, -0.62348456475000202, -0.24583080842150512]
[-0.1422280331828244, 0.62348456475000202, 0.24583080842150512]
[-1.2182327603643608, 0.1230166261435828, -1.5722574943149592]
[1.2182327603643608, -0.1230166261435828, 1.5722574943149592]
[-1.5849327743471977, -0.67400176640115994, -0.62126884010505024]
[1.5849327743471977, 0.67400176640115994, 0.62126884010505024]
[2.0938608215128593, 0.84753559419590085, -0.5755506225263638]
[-2.0938608215128593, -0.84753559419590085, 0.5755506225263638]
[1.7271608075300224, 0.050517201651157972, 0.37543803168354511]
[-1.7271608075300224, -0.050517201651157972, -0.37543803168354511]
[0.36670001398283714, 0.79701839254474283, -0.95098865420990886]
[-0.36670001398283714, -0.79701839254474283, 0.95098865420990886]

How it works

Note that when three planes (whose normals are non-coplanar) partition the space, each of the eight resulting regions is bounded by (subsets of) all three planes. If we added a plane, it would divide some of the regions in two, resulting, again, in a pair of subregions bounded by at least three planes. Any cross-section of such a region is a polygon, whose vertices lie on the lines of intersection of the planes bounding the region. Moreover, since these planes extend to infinity (as planes usually do), the polygon is necessarily convex, the upshot of which is that the region is closed under addition. In particular, if we take one vector along each line of intersection of a pair of bounding planes that lies within the region, the sum of these vectors is also within the region and not on any of the planes (and, hence, not in any other region.)

The program recursively partitions the space by incrementally adding planes and subdividing the currently-explored region, while keeping a track of vectors along the lines of intersection and eventually taking their sum.

This method breaks if only two planes are given, since then there's only one line of intersection. We treat this as a special case and return the sum of a pair of vectors, each of which lies on a different plane, taking the four different sign combinations.