In [1]:
import vtk 
import pandas as pd
import numpy as np
import pygslib
from matplotlib import pyplot as plt
import vtk.util.numpy_support as vtknumpy
%matplotlib inline



In [47]:
def inside_polygon(points,polygon_pts = None, polygon = None, normal = None, bounds = None):
    
    result = np.zeros(len(points), dtype = int)
        
    if polygon is None:
        assert polygon_pts is not None
        
        polygon = vtk.vtkPolygon()
        
        for i in polygon_pts:
            polygon.GetPoints().InsertNextPoint(i)
        
    # compute normal 
    if normal is None:
        normal=[0.0,0.0,0.0]
        polygon.ComputeNormal(polygon.GetPoints(), normal)
    # compute bounds
    if bounds is None:
        bounds = [0.,0.,0.,0.,0.,0.];
        polygon.GetPoints().GetBounds(bounds)
    
    for i in range(len(points)):
        result[i] = polygon.PointInPolygon(points[i], 
                       polygon.GetPoints().GetNumberOfPoints(), 
                       vtknumpy.vtk_to_numpy(polygon.GetPoints().GetData()).ravel(), 
                       bounds, 
                       normal) 
        
    return result, polygon, normal, bounds


def polygon_area(polygon_pts = None, polygon = None):
    
    
    if polygon is None:
        assert polygon_pts is not None
        
        polygon = vtk.vtkPolygon()
        
        for i in polygon_pts:
            polygon.GetPoints().InsertNextPoint(i)
    
    # compute normal
    normal=[0.0,0.0,0.0]
    polygon.ComputeNormal(polygon.GetPoints(), normal)    
    
    # compute area
    return polygon.ComputeArea(polygon.GetPoints(),
                        polygon.GetPoints().GetNumberOfPoints(),
                        list(range(polygon.GetPoints().GetNumberOfPoints())),
                        normal)

In [51]:
polygon_pt = [ [0., 0., 0.],
               [1., 0., 0.],
               [1., .9, 0.],
               [0., 1., 0.]]

testIn = [[0.5, 0.5, 0.0]]
testOut = [[2.0, 0.5, 0.0]]

In [52]:
result, polygon, normal, bounds = inside_polygon(points = testIn, 
                                                 polygon_pts = polygon_pt, 
                                                 polygon = None,
                                                 normal = None, 
                                                 bounds = None)

In [53]:
print(result)
print(normal)
print(bounds)


[1]
[0.0, 0.0, 1.0]
[0.0, 1.0, 0.0, 1.0, 0.0, 0.0]

In [54]:
#call with precalculted polygons
result, polygon, normal, bounds = inside_polygon(points = testIn, 
                                                 polygon_pts = None, 
                                                 polygon = polygon,
                                                 normal = None, 
                                                 bounds = None)

print(result)
print(normal)
print(bounds)

result, polygon, normal, bounds = inside_polygon(points = testIn, 
                                                 polygon_pts = None, 
                                                 polygon = polygon,
                                                 normal = normal, 
                                                 bounds = None)

print(result)
print(normal)
print(bounds)

result, polygon, normal, bounds = inside_polygon(points = testIn, 
                                                 polygon_pts = None, 
                                                 polygon = polygon,
                                                 normal = normal, 
                                                 bounds = bounds)

print(result)
print(normal)
print(bounds)


[1]
[0.0, 0.0, 1.0]
[0.0, 1.0, 0.0, 1.0, 0.0, 0.0]
[1]
[0.0, 0.0, 1.0]
[0.0, 1.0, 0.0, 1.0, 0.0, 0.0]
[1]
[0.0, 0.0, 1.0]
[0.0, 1.0, 0.0, 1.0, 0.0, 0.0]

In [55]:
# area of a polygon defined by array of floats
polygon_area(polygon_pts = polygon_pt, polygon = None)


Out[55]:
0.95

In [56]:
# area of a polygon defined by vtkPolygon
polygon_area(polygon_pts = None, polygon = polygon)


Out[56]:
0.95

In [ ]: