From adee794550a258c9307eb0d9e343840abef1bd2f Mon Sep 17 00:00:00 2001
From: Runar <runar.berge@uib.no>
Date: Thu, 2 Feb 2017 14:30:59 +0100
Subject: [PATCH] New function: half_space_int - find points that lie in the
 intersection of half planes

Also with unit tests
---
 .gitignore                   |  2 ++
 half_space_int.py            | 48 ++++++++++++++++++++++++++++++++++++
 tests/test_half_space_int.py | 26 +++++++++++++++++++
 3 files changed, 76 insertions(+)
 create mode 100644 half_space_int.py
 create mode 100644 tests/test_half_space_int.py

diff --git a/.gitignore b/.gitignore
index 0d20b6487c..840e386596 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,3 @@
 *.pyc
+*.py#
+*.py~
\ No newline at end of file
diff --git a/half_space_int.py b/half_space_int.py
new file mode 100644
index 0000000000..23f3f0d15c
--- /dev/null
+++ b/half_space_int.py
@@ -0,0 +1,48 @@
+import numpy as np
+
+
+def half_space_int(n,x0, pts):
+    """
+    Find the points that lie in the intersection of half spaces (3D)
+
+    Parameters
+    ----------
+    n : ndarray
+        This is the normal vectors of the half planes. The normal 
+        vectors is assumed to point out of the half spaces.
+    x0 : ndarray
+        Point on the boundary of the half-spaces. Half space i is given
+        by all points satisfying (x - x0[:,i])*n[:,i]<=0
+    pts : ndarray
+        The points to be tested if they are in the intersection of all
+        half-spaces or not.
+
+    Returns
+    -------
+    out : ndarray
+        A logical array with length equal number of pts.
+
+        out[i] is True if pts[:,i] is in all half-spaces
+
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> n = np.array([[0,1],[1,0],[0,0]])
+    >>> x0 = np.array([[0,-1],[0,0],[0,0]])
+    >>> pts = np.array([[-1,1,4],[2,1,-2],[0,0,0]])
+    >>> out = half_plane_int(n,x0,pts)
+    array([False,  True, False], dtype=bool)
+    """
+    assert n.shape[0] == 3, ' only 3D supported'
+    assert x0.shape[0] == 3, ' only 3D supported'
+    assert pts.shape[0] == 3, ' only 3D supported'
+    assert n.shape[1] == x0.shape[1], 'ther must be the same number of normal vectors as points'
+    n_pts   = pts.shape[1]
+    in_hull = np.zeros(n_pts)
+    x0      = np.repeat(x0[:,:,np.newaxis], n_pts,axis=2)
+    n       = np.repeat( n[:,:,np.newaxis], n_pts,axis=2)
+    for i in range(x0.shape[1]):
+        in_hull += np.sum((pts - x0[:,i,:])*n[:,i,:],axis=0)<=0
+        
+    return in_hull==x0.shape[1]
diff --git a/tests/test_half_space_int.py b/tests/test_half_space_int.py
new file mode 100644
index 0000000000..8887917763
--- /dev/null
+++ b/tests/test_half_space_int.py
@@ -0,0 +1,26 @@
+import unittest
+import numpy as np
+
+from utils import half_space_int
+
+
+def test_one_half_space():
+
+    n  = np.array([[-1],[0],[0]])
+    x0 = np.array([[0],[0],[0]])
+    pts = np.array([[1,-1],[0,0],[0,0]])
+    out = half_space_int.half_space_int(n,x0,pts)
+    assert np.all(out == np.array([True,False]))
+
+def test_two_half_spaces():
+    n  = np.array([[-1,0],[0,-1],[0,0]])
+    x0 = np.array([[0,0],[0,1],[0,0]])
+    pts = np.array([[1,-1,1,0],[2,0,2,0],[0,0,0,0]])
+    out = half_space_int.half_space_int(n,x0,pts)
+    assert np.all(out == np.array([True,False,True,False]))
+
+
+if __name__ == '__main__':
+    test_one_half_space()
+    test_two_half_spaces()
+