# Find one point per pyramidal region

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))
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

• Makes no sense to me. How can planes form a cone? Plane are square, cone are round ??? – feersum Oct 3 '14 at 22:53
• This could really do with a better explanation, a picture, and a reference to the foumula. I haven't a clue what you mean by a cone here. Also, I considered that the faces of a tetrahedron if extended infinitely will divide space into a finite closed space (the tetrahedron) plus 14 infinite regions, which matches your formula, but there are other cases (such as the cube, 6 faces + 8 corners + 12 edges = 26 infinite regions) which don't match your formula (2-6+36=32) so the only thing I could guess you might mean seems wrong. Please clarify. – Level River St Oct 4 '14 at 11:55
• @steveverrill I added a picture and also an example valid output for the 4 planes described in the question. If anything is still unclear, please let me know. – user9206 Oct 4 '14 at 19:26
• @feersum Cone has a different meaning in linear algebra as a set of points closed under scaling: en.wikipedia.org/wiki/Cone_(linear_algebra) – xnor Oct 4 '14 at 20:46
• @xnor I changed the word cone so I hope it is less controversial now. I am sad about the close vote however :( – user9206 Oct 4 '14 at 20:50

# Python, 423376 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.

• Thank you for this! I haven't managed to decrypt it yet so any explanation gratefully accepted. In reply to your question, ideally it should work up to about 1000 planes. – user9206 Oct 6 '14 at 14:19
• @Lembik I'll add an explanation later, hopefully. However, I was asking about the minimal number of planes (currently, the program assumes at least two planes.) – Ell Oct 6 '14 at 14:34
• Oh right. Two planes is fine as a minimum. – user9206 Oct 6 '14 at 16:34
• Thank you for the explanation. Would you mind expanding "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." please. In particular, do you only keep track of lines of intersection or actually bounded regions as well? If it is just lines of intersection, let us say we add one new plane and work out all the lines of intersection with the current planes. What do you do now exactly? – user9206 Oct 7 '14 at 19:32
• @Lembik It keeps track of both the bounded region and lines of intersection (LoI). P is where most happen; it takes, in order, the bounding planes, the LoI and the remaining planes. It checks if one of the remaining planes bisects the region by finding its LoI with all the bounding planes and testing if any of them is on the boundary of the region, thus killing two birds with one stone. It then calls itself recursively for each of the new subregions (passing only the LoI of that region,) or for the same region, with the remaining planes. When no planes are left it returns the sum of the LoI. – Ell Oct 7 '14 at 20:41