#! /usr/bin/env python3 # def polygon_test ( ): import matplotlib.pyplot as plt import numpy as np print ( '' ) print ( 'polygon_test():' ) # # Define vertices V of a polygon. # print ( '' ) print ( ' Define polygon vertices V:' ) v = np.array ( [ [0,0], [4,0], [2,2], [1,2], [0,1], [1,1] ] ) v_num = v.shape[0] print ( v ) # # Plot the polygon. # print ( '' ) print ( ' Plot outline of polygon:' ) plt.clf ( ) for i in range ( 0, v_num ): ip1 = ( i + 1 ) % v_num x = [ v[i,0], v[ip1,0] ] y = [ v[i,1], v[ip1,1] ] plt.plot ( x, y, 'k-' ) i = 0 for xy in v: plt.plot ( xy[0], xy[1], 'b.', markersize = 20 ) plt.text ( xy[0], xy[1] + 0.10, str ( i ), fontsize = 20 ) i = i + 1 plt.grid ( True ) plt.axis ( 'equal' ) filename = 'polygon_outline.png' plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + '"' ) plt.show ( ) plt.close ( ) # # Triangulate the polygon. # print ( '' ) print ( ' Triangulate the polygon:' ) tri = np.array ( [ [ 0, 1, 2 ], [2, 3, 5], [3,4,5] ] ) print ( tri ) # # Plot triangulated polygon. # print ( '' ) print ( ' Plot the triangulated polygon:' ) plt.clf ( ) for ti in tri: x = [ v[ti[0],0], v[ti[1],0], v[ti[2],0] ] y = [ v[ti[0],1], v[ti[1],1], v[ti[2],1] ] plt.fill ( x, y ) for ti in tri: for i in range ( 0, 3 ): ip1 = ( i + 1 ) % 3 x = [ v[ti[i],0], v[ti[ip1],0] ] y = [ v[ti[i],1], v[ti[ip1],1] ] plt.plot ( x, y, 'k-' ) i = 0 for xy in v: plt.plot ( xy[0], xy[1], 'b.', markersize = 20 ) plt.text ( xy[0], xy[1] + 0.10, str ( i ), fontsize = 20 ) i = i + 1 plt.grid ( True ) plt.axis ( 'equal' ) filename = 'polygon_filled.png' plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + '"' ) plt.show ( ) plt.close ( ) # # Polygon area by shoelace algorithm. # area = polygon_area_shoelace ( v ) print ( '' ) print ( ' Polygon area by shoelace = ', area ) # # Polygon area by triangulation. # area = polygon_area_triangulated ( v, tri ) print ( '' ) print ( ' Polygon area by triangulation = ', area ) # # Polygon centroid. # centroid = polygon_centroid ( v, tri ) print ( '' ) print ( ' Polygon centroid = ', centroid ) # # Polygon contains point? # print ( '' ) print ( ' polygon_contains_point()?' ) for q in [ [0.0,0.5], [0.75,1.25], [2.5,1.0], [1.5,2.5] ]: contains = polygon_contains_point ( v, tri, q ) print ( ' Polygon does', end = '' ) if ( not contains ): print ( ' NOT', end = '' ) print ( ' contain ', q ) # # Polygon Quad1 # print ( '' ) print ( ' Estimate integral of f(x,y) using 1 point triangle rule' ) f = lambda x, y: x**2 + y**2 Q = polygon_quad1 ( v, tri, f ) print ( Q ) # # Polygon Quad6 # print ( '' ) print ( ' Estimate integral of f(x,y) using 6 point triangle rule' ) f = lambda x, y: x**2 + y**2 Q = polygon_quad6 ( v, tri, f ) print ( Q ) return def line_side ( p0, p1, q ): import numpy as np v0 = p1 - p0 v1 = q - p0 v0xv1 = float ( np.cross ( v0, v1 ) ) return v0xv1 def polygon_area_shoelace ( v ): import numpy as np n = v.shape[0] np.append ( v, v[0] ) area = 0.5 * ( np.sum ( v[0:n-1,0] * v[1:n,1] ) \ - np.sum ( v[1:n,0] * v[0:n-1,1] ) ) return area def polygon_area_triangulated ( v, tri ): import numpy as np area = 0 for ti in tri: t = np.array ( v[ti] ) area = area + triangle_area ( t ) return area def polygon_centroid ( v, tri ): import numpy as np centroid = np.zeros ( 2 ) area = 0.0 for ti in tri: t = np.array ( v[ti] ) t_area = triangle_area ( t ) t_centroid = triangle_centroid ( t ) area = area + t_area centroid = centroid + t_area * t_centroid centroid = centroid / area return centroid def polygon_contains_point ( v, tri, q ): import numpy as np contains = False for ti in tri: t = np.array ( v[ti] ) contains = triangle_contains_point ( t, q ) if ( contains ): break return contains def polygon_quad1 ( v, tri, f ): import numpy as np Q = 0.0 for ti in tri: t = np.array ( v[ti] ) Q = Q + triangle_quad1 ( t, f ) return Q def polygon_quad6 ( v, tri, f ): import numpy as np Q = 0.0 for ti in tri: t = np.array ( v[ti] ) Q = Q + triangle_quad6 ( t, f ) return Q def triangle_area ( t ): import numpy as np s = triangle_sides ( t ) abc = 0.5 * np.sum ( s ) area = np.sqrt ( abc * ( abc - s[0] ) * ( abc - s[1] ) * ( abc - s[2] ) ) return area def triangle_centroid ( t ): import numpy as np centroid = np.sum ( t, axis = 0 ) / 3.0 return centroid def triangle_contains_point ( t, q ): contains = True for i in range ( 0, 3 ): ip1 = ( i + 1 ) % 3 if ( line_side( t[i,:], t[ip1,:], q ) < 0.0 ): contains = False break return contains def triangle_quad1 ( t, f ): A = triangle_area ( t ) xc = ( t[0,0] + t[1,0] + t[2,0] ) / 3.0 yc = ( t[0,1] + t[1,1] + t[2,1] ) / 3.0 Q = A * f ( xc, yc ) return Q def triangle_quad6 ( t, f ): A = triangle_area ( t ) a = 0.816847572980459; b = 0.091576213509771; c = 0.108103018168070; d = 0.445948490915965; u = 0.109951743655322; v = 0.223381589678011; xab = a * t[0,0] + b * t[1,0] + b * t[2,0] xba = b * t[0,0] + a * t[1,0] + b * t[2,0] xbb = b * t[0,0] + b * t[1,0] + a * t[2,0] yab = a * t[0,1] + b * t[1,1] + b * t[2,1] yba = b * t[0,1] + a * t[1,1] + b * t[2,1] ybb = b * t[0,1] + b * t[1,1] + a * t[2,1] xcd = c * t[0,0] + d * t[1,0] + d * t[2,0] xdc = d * t[0,0] + c * t[1,0] + d * t[2,0] xdd = d * t[0,0] + d * t[1,0] + c * t[2,0] ycd = c * t[0,1] + d * t[1,1] + d * t[2,1] ydc = d * t[0,1] + c * t[1,1] + d * t[2,1] ydd = d * t[0,1] + d * t[1,1] + c * t[2,1] Q = 0.0 Q = Q + u * ( f ( xab, yba ) + f ( xba, yab ) + f ( xbb, ybb ) ) Q = Q + v * ( f ( xcd, ydc ) + f ( xdc, ycd ) + f ( xdd, ydd ) ) Q = A * Q return Q def triangle_sides ( t ): import numpy as np s = np.zeros ( 3 ) for i in range ( 0, 3 ): ip1 = ( i + 1 ) % 3 ip2 = ( i + 2 ) % 3 s[i] = np.linalg.norm ( t[ip1,:] - t[ip2,:] ) return s if ( __name__ == "__main__" ): polygon_test ( )