Point in Polyhedron

[fatalaccidents]

Hello all,

I was hoping to figure out a code that could find out whether a point was inside a group of points. I already noticed that I could create a convex hull from the points with the scipy.spatial.ConvexHull module. What I'm trying to do is then find a way to figure out if certain points are within that convex hull. For 2D point in poly I was using the ray casting method that I found online:

def point_inside_polygon(x,y,poly):

    n = len(poly)
    inside =False

    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xinters:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside

I was wondering if there was a way to adapt this to 3 dimensions, or if there was something prepackaged that could be used. Thank you all for your time!

[Gribouillis]

It seems to me that this question is solvable by linear programming. The function scipy.optimize.linprog() should work with convenient parameters.

[Gribouillis]

Here is the basic idea

#!/usr/bin/env python
# -*-coding: utf8-*-
'''doc
'''
from __future__ import (absolute_import, division,
                        print_function, unicode_literals)

import numpy as np
from scipy.optimize import linprog

def main():
    point = np.array([0.8, 0.12, 0.4])
    n = len(point)
    N = 30
    a = np.random.rand(n, N)
    a = np.insert(a, n, 1.0, axis=0)
    b = np.append(point, 1.0).reshape(n+1,1)
    bounds = [(0.0, 1.0) for i in range(N)]
    c = np.zeros(N)
    c[0] = 1.0
    print(a)
    print(b)
    print(c)
    result = linprog(c, A_eq=a, b_eq=b, bounds=bounds)
    print(result)
    if result.status == 0:
        print('Point is in convex hull')
    else:
        print('Point is not in convex hull')

if __name__ == '__main__':
    main()

""" my output -->
17:31 2015-03-31] python convex.py
[[ 0.26366509  0.92411732  0.9081528   0.88461942  0.68029591  0.42997518
   0.17189674  0.07542255  0.6627228   0.46370079  0.29521928  0.50230071
   0.73392039  0.70076308  0.39031932  0.73670852  0.64815496  0.52975825
   0.78314098  0.60898234  0.63308929  0.32659967  0.02298467  0.17494737
   0.01397131  0.27248117  0.98193604  0.69879033  0.90454457  0.06127812]
 [ 0.16418386  0.05802131  0.07673203  0.50894301  0.45885838  0.13100915
   0.43334997  0.3613589   0.70104947  0.47547909  0.64204601  0.18567562
   0.07257718  0.49943999  0.96490179  0.16901995  0.14288372  0.55158995
   0.12849092  0.4287121   0.33522966  0.32258957  0.54137264  0.24533261
   0.26382428  0.35957173  0.45413843  0.96331716  0.41831803  0.34442463]
 [ 0.77288019  0.04946268  0.00280707  0.07217767  0.62230434  0.97351116
   0.78616099  0.74724063  0.35259706  0.82299753  0.8465832   0.7326628
   0.65494121  0.88904928  0.67263953  0.23434955  0.74143111  0.4536738
   0.81357322  0.7408426   0.52483849  0.44427614  0.75840825  0.04656877
   0.84187546  0.94461097  0.42198365  0.20315269  0.21713841  0.29179429]
 [ 1.          1.          1.          1.          1.          1.          1.
   1.          1.          1.          1.          1.          1.          1.
   1.          1.          1.          1.          1.          1.          1.
   1.          1.          1.          1.          1.          1.          1.
   1.          1.        ]]
[[ 0.8 ]
 [ 0.12]
 [ 0.4 ]
 [ 1.  ]]
[ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
  status: 0
   slack: array([ 1.        ,  1.        ,  0.64624474,  1.        ,  1.        ,
        1.        ,  1.        ,  1.        ,  1.        ,  1.        ,
        1.        ,  1.        ,  0.41893992,  1.        ,  1.        ,
        1.        ,  1.        ,  1.        ,  1.        ,  1.        ,
        1.        ,  1.        ,  1.        ,  1.        ,  1.        ,
        1.        ,  0.97621747,  0.95859787,  1.        ,  1.        ])
 success: True
     fun: -0.0
       x: array([ 0.        ,  0.        ,  0.35375526,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.58106008,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.02378253,  0.04140213,  0.        ,  0.        ])
 message: 'Optimization terminated successfully.'
     nit: 5
Point is in convex hull
"""

In this code, there is a dummy objective function. The simplex method could be halted after first phase when a feasable point has been reached. On the other hand, optimizing the dummy function may be cheaper in most cases than detecting the end of phase 1 with a callback.

Edit: it seems that the algorithm works with line 20 commented, which means that the dummy objective function is zero. The effect of this is to stop at the end of phase 1: for a constant objective, any feasable point is a solution. Also the test at the end should probably be replaced by if result.success .

[Gribouillis]

linprog() needs scipy version >= 0.15.0. In ubuntu, I had to remove package scipy and reinstall scipy with the pip command (as of march 2015).

[snippsat]

because my anaconda distribution doesn't seem to have linprog (?)

You can use pip install or conda install(Scripts folder) to get new stuff into Anaconda.
For this is best to use conda install scikit-learn,then you get all new packages that scikit-learn need.
Look like this when run conda install scikit-learn.

The following packages will be downloaded:   

package                    |            build
    ---------------------------|-----------------
    conda-3.10.0               |           py27_0         207 KB
    conda-env-2.1.3            |           py27_0          54 KB
    nose-1.3.4                 |           py27_1         233 KB
    numpy-1.9.2                |           py27_0        23.2 MB
    requests-2.6.0             |           py27_0         590 KB
    scikit-learn-0.16.0        |       np19py27_0         3.5 MB
    scipy-0.15.1               |       np19py27_0        71.3 MB
    ------------------------------------------------------------
                                           Total:        99.0 MB

As you see you get scipy version 0.15.1.

You can even create a virtual environment,with all stuff you need.
conda create -n test scikit-learn
Create a folder test(virtual environment) that has Python and scikit-learn installed.