#! /usr/bin/env python3
#
def hull_sample_test ( ):

  from scipy.spatial import ConvexHull
  from scipy.spatial import Delaunay
  import matplotlib.pyplot as plt
  import numpy as np

  points = np.loadtxt ( 'i18.txt' )
#
#  Compute the hull.
#
  hull = ConvexHull ( points )
#
#  List the boundary points.
#
  print ( hull.vertices )

  boundary = points[hull.vertices]
  print ( boundary )

  boundary_wrap = np.vstack ( ( boundary, boundary[0,:] ) )
  print ( boundary_wrap )
#
#  Display the convex hull.
#
  plt.clf ( )
  plt.plot ( boundary_wrap[:,0], boundary_wrap[:,1], 'r-' )
  plt.show ( )
#
#  Triangulate the convex hull.
#
  tri = Delaunay ( boundary )
  print ( tri )

  plt.clf ( )
  plt.triplot ( boundary[:,0], boundary[:,1], tri.simplices )
  plt.plot ( boundary[:,0], boundary[:,1], 'o' )
  filename = 'i18_convex_hull_triangulated.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
#
#  Compute the area of each triangle.
#
  tri_num = len ( tri.simplices )
  area = np.zeros ( tri_num )

  i = 0
  for simplex in tri.simplices:
    print ( boundary[simplex] )
    area[i] = triangle_signed_area ( boundary[simplex] )
    i = i + 1

  print ( area )
#
#  Compute the intervals in [0,1] that correspond to each triangle's
#  relative area.
#
  hull_cdf = np.cumsum ( area )
  area_sum = np.sum ( area )
  hull_cdf = hull_cdf / area_sum
  print ( hull_cdf )
#
#  Pick a triangle i uniformly at random based on area.
#
  r = np.random.random ( )
  for i in range ( 0, tri_num ):
    if ( r <= hull_cdf[i] ):
      break
  print ( '  Random triangle is ', i )
#
#  Pick a point in triangle i at random.
#
  simplex = tri.simplices[i]
  p = triangle_sample ( boundary[simplex] )
  print ( '  Random point is ', p )
#
#  Plot 100 points at random.
#
  plt.clf ( )
  plt.triplot ( boundary[:,0], boundary[:,1], tri.simplices )
  plt.plot ( boundary[:,0], boundary[:,1], 'o' )

  for k in range ( 0, 100 ):
    r = np.random.random ( )
    for i in range ( 0, tri_num ):
      if ( r <= hull_cdf[i] ):
        break
    simplex = tri.simplices[i]
    p = triangle_sample ( boundary[simplex] )
    plt.plot ( p[0], p[1], 'bo' )

  filename = 'i18_convex_hull_sampled.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )

  plt.show ( )

  return

def triangle_sample ( t ):
  import numpy as np
  r1 = np.random.rand ( )
  r2 = np.random.rand ( )
  alpha = 1 - np.sqrt ( r1 )
  beta = np.sqrt( r1 ) * r2
  gamma = np.sqrt ( r1 ) * ( 1.0 - r2 )
  p = alpha * t[0,:] + beta * t[1,:] + gamma * t[2,:]
  return p

def triangle_signed_area ( t ):
  import numpy as np
  value = np.cross ( t[1,:] - t[0,:], t[2,:] - t[0,:] )
  value = 0.5 * float ( value )
  return value

if ( __name__ == "__main__" ):
  hull_sample_test ( )