from sympy.geometry import *

P1 = Point(0, 0)
P2 = Point(3, 4)
P3 = Point(2, -1)
P4 = Point(-1, 5)

S1 = Segment(P1, P2)
S2 = Segment(P3, P4)

Point.is_collinear(P1, P2, P3)

S1.length

S2.midpoint

S1.slope

S1.intersection(S2)

Segment.angle_between(S1, S2)

S1.contains(P3)

L1 = Line(P1, P2)

L2 = L1.perpendicular_line(P3)  # perpendicular line to L1

L2.arbitrary_point()            # parametric equation of L2

L2.equation()                   # algebraic equation of L2

L2.contains(P4)                 # is P4 in L2?

L2.distance(P4)                 # distance from P4 to L2

L1.is_parallel(S2)              # is S2 parallel to L1?

C1 = Circle(P1, 3)
C2 = Circle(P2, P3, P4)

C2.area

C2.radius

C2.equation()

C2.center

C2.circumference

C2.intersection(C1)

C2.intersection(S1)

C2.is_tangent(L2)

C1.tangent_lines(P4)

T = Triangle(P1, P2, P3)

T.area                          # Note it gives a signed area

T.angles

T.sides

T.perimeter

T.is_right()                    # Is T a right triangle?

T.is_equilateral()              # Is T equilateral?

T.is_scalene()                  # Is T scalene?

T.is_isosceles()                # Is T isosceles?

T.altitudes

T.orthocenter                  # Intersection of the altitudes

T.bisectors()                  # Angle bisectors

T.incenter                     # Intersection of angle bisectors

T.incircle

T.inradius

T.medians

T.centroid                   # Intersection of the medians

T.circumcenter               # Intersection of perpendicular bisectors

T.circumcircle

T.circumradius

T.medial

T.intersection(C1)

T.distance(T.circumcenter)

T.is_similar(Triangle(P1, P2, P4))

from sympy import var, pi, sin, cos

var('t', real=True)

Arc = Curve((3*cos(t), 4*sin(t)), (t, 0, 3*pi/4))

T.reflect(L1)

T.rotate(pi/2, P2)

T.translate(5,4)

T.scale(9)

Arc.rotate(pi/2, P3).translate(pi,pi).scale(0.5)

from numpy import array

def read_poly(file_name):
    """
    Simple poly-file reader, that creates a python dictionary 
    with information about vertices, edges and holes.
    It assumes that vertices have no attributes or boundary markers.
    It assumes that edges have no boundary markers.
    No regional attributes or area constraints are parsed.
    """

    output = {'vertices': None, 'holes': None, 'segments': None}

    # open file and store lines in a list
    file = open(file_name, 'r')
    lines = file.readlines()
    file.close()
    lines = [x.strip('\n').split() for x in lines]

    # Store vertices
    vertices= []
    N_vertices, dimension, attr, bdry_markers = [int(x) for x in lines[0]]
    # We assume attr = bdrt_markers = 0
    for k in range(N_vertices):
        label, x, y = [items for items in lines[k+1]]
        vertices.append([float(x), float(y)])
    output['vertices']=array(vertices)

    # Store segments
    segments = []
    N_segments, bdry_markers = [int(x) for x in lines[N_vertices+1]]
    for k in range(N_segments):
        label, pointer_1, pointer_2 = [items for items in lines[N_vertices+k+2]]
        segments.append([int(pointer_1)-1, int(pointer_2)-1])
    output['segments'] = array(segments)

    # Store holes
    N_holes = int(lines[N_segments+N_vertices+2][0])
    holes = []
    for k in range(N_holes):
        label, x, y = [items for items in lines[N_segments + N_vertices + 3 + k]]
        holes.append([float(x), float(y)])

    output['holes'] = array(holes)

    return output

import numpy as np

from scipy.spatial import ConvexHull

import matplotlib.pyplot as plt

%matplotlib inline

lake_superior = read_poly("/Users/francisco.blanco.silva/Dropbox/Documents/Books/Mastering/chapter6/superior.poly")

vertices_ls = lake_superior['vertices']

%time hull = ConvexHull(vertices_ls)

plt.figure(figsize=(18, 18))
plt.xlim(vertices_ls[:,0].min()-0.01, vertices_ls[:,0].max()+0.01)
plt.ylim(vertices_ls[:,1].min()-0.01, vertices_ls[:,1].max()+0.01)
plt.axis('off')
plt.axes().set_aspect('equal')
plt.plot(vertices_ls[:,0], vertices_ls[:,1], 'b.')
for simplex in hull.simplices:
    plt.plot(vertices_ls[simplex, 0], vertices_ls[simplex, 1], 'r-')
    
plt.show()

from mayavi import mlab

points = np.random.rand(320, 3)

hull = ConvexHull(points)

X = hull.points[:, 0]
Y = hull.points[:, 1]
Z = hull.points[:, 2]

mlab.triangular_mesh(X, Y, X, hull.simplices, colormap='gray', 
                     opacity=0.5, representation='wireframe')

mlab.show()

from scipy.spatial import Voronoi, voronoi_plot_2d

vor = Voronoi(vertices_ls)

plt.figure(figsize=(10, 10))
ax = plt.subplot(111, aspect='equal')
voronoi_plot_2d(vor, ax=ax)
plt.xlim( 0.45,  0.50)
plt.ylim(-0.40, -0.35)
plt.show()

vor.point_region

vor.vertices

[vor.vertices[x] for x in vor.regions[vor.point_region[4]]]

vor.ridge_points

from scipy.spatial import Delaunay, delaunay_plot_2d

tri = Delaunay(vertices_ls)

plt.figure(figsize=(18, 18))
ax = plt.subplot(111, aspect='equal')
delaunay_plot_2d(tri, ax=ax)
plt.show()

!pip install triangle

from triangle import triangulate, plot as tplot

cndt = triangulate(lake_superior, 'p')

plt.figure(figsize=(18, 18))
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax, **cndt)
plt.show()

cfdt = triangulate(lake_superior, 'pD')

cncfq20dt = triangulate(lake_superior, 'pq20D')

plt.figure(figsize=(18,18))
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax, **cncfq20dt)
plt.show()

cncfq20adt = triangulate(lake_superior, 'pq20a.001D')

plt.figure(figsize=(18, 18))
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax, **cncfq20adt)
plt.show()

from scipy.spatial import minkowski_distance

X = cncfq20adt['triangles'][:,0]
Y = cncfq20adt['triangles'][:,1]
Z = cncfq20adt['triangles'][:,2]

Xvert = [cncfq20adt['vertices'][x] for x in X]
Yvert = [cncfq20adt['vertices'][y] for y in Y]
Zvert = [cncfq20adt['vertices'][z] for z in Z]

lengthsXY = minkowski_distance(Xvert, Yvert)
lengthsXZ = minkowski_distance(Xvert, Zvert)
lengthsYZ = minkowski_distance(Yvert, Zvert)

from scipy.sparse import lil_matrix
from scipy.sparse.csgraph import shortest_path

nvert = len(cncfq20adt['vertices'])
G = lil_matrix((nvert, nvert))

for k in range(len(X)):
    G[X[k], Y[k]] = G[Y[k], X[k]] = lengthsXY[k]
    G[X[k], Z[k]] = G[Z[k], X[k]] = lengthsXZ[k]
    G[Y[k], Z[k]] = G[Z[k], Y[k]] = lengthsYZ[k]

dist_matrix, pred = shortest_path(G, return_predecessors=True, directed=True, unweighted=False)

index = 370
path = [370]

while index != 197:
    index = pred[197, index]
    path.append(index)
    
Xs = [cncfq20adt['vertices'][x][0] for x in path]
Ys = [cncfq20adt['vertices'][x][1] for x in path]

plt.figure(figsize=(18,18))
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax, **cncfq20adt)
ax.plot(Xs, Ys, '-', linewidth=5, color='blue'); \
plt.show()

P1 = Point(0, 0)
P2 = Point(1, 0)
P3 = Point(-1, 0)
P4 = Point(0, 1)

C = Circle(P2, P3, P4)
T = Triangle(P1, P2, P3)

C.encloses_point(P1)

C.encloses(T)

from matplotlib.path import Path

my_polygon = Path([hull.points[x] for x in hull.vertices])

X = .25 * np.random.randn(100) + .5
Y = .25 * np.random.randn(100) - .5

my_polygon.contains_points([[X[k], Y[k]] for k in range(len(X))])

from scipy.spatial import tsearch

tri = Delaunay(vertices_ls)

points = zip(X, Y)

tsearch(tri, points)

tri.find_simplex(points)

from scipy.spatial import cKDTree

vor  = Voronoi(vertices_ls)
tree = cKDTree(vertices_ls)

tree.query(points)

X = np.linspace( 0.45,  0.50, 256)
Y = np.linspace(-0.40, -0.35, 256)
canvas = np.meshgrid(X, Y)
points = np.c_[canvas[0].ravel(), canvas[1].ravel()]
queries = tree.query(points)[1].reshape(256, 256)

plt.figure(figsize=(18, 18))

ax1 = plt.subplot(121, aspect='equal')
voronoi_plot_2d(vor, ax=ax1)
plt.xlim( 0.45,  0.50)
plt.ylim(-0.40, -0.35)

ax2 = plt.subplot(122, aspect='equal')
plt.gray()
plt.pcolor(X, Y, queries)
plt.plot(vor.points[:,0], vor.points[:,1], 'ro')
plt.xlim( 0.45,  0.50)
plt.ylim(-0.40, -0.35)
plt.show()

points = np.random.rand(320, 2)
range_points = np.random.rand(5, 2)
range_radii = 0.1 * np.random.rand(5)

tree = cKDTree(points)
result = set()

for k in range(5):
    point  = range_points[k]
    radius = range_radii[k]
    partial_query = tree.query_ball_point(point, radius)
    result = result.union(set(partial_query))
    
print result
    
fig = plt.figure(figsize=(9, 9))
plt.axes().set_aspect('equal')
for point in points:
    plt.plot(point[0], point[1], 'ko')
    
for k in range(5):
    point = range_points[k]
    radius = range_radii[k]
    circle = plt.Circle(point, radius, fill=False, color="red", lw=2)
    fig.gca().add_artist(circle)
    
plt.show()

small_superior = lake_superior['vertices'][:9]

tri = Delaunay(small_superior, incremental=True)
vor = Voronoi(small_superior, incremental=True)

plt.figure(figsize=(18,18))

for k in range(4):
    tri.add_points(vertices_ls[10*(k+1):10*(k+2)-1])
    vor.add_points(vertices_ls[10*(k+1):10*(k+2)-1])
    ax1 = plt.subplot(4, 2, 2*k+1, aspect='equal')
    ax1.set_xlim( 0.00,  1.00)
    ax1.set_ylim(-0.70, -0.30)
    delaunay_plot_2d(tri, ax1)
    ax2 = plt.subplot(4, 2, 2*k+2, aspect='equal')
    ax2.set_xlim(0.0, 1.0)
    ax2.set_ylim(-0.70, -0.30)
    voronoi_plot_2d(vor, ax2)
    
plt.show()

import matplotlib.patches as patches

def bezier_parabola(P1, P2, P3):
   return Path([P1, P2, P3],
               [Path.MOVETO, Path.CURVE3, Path.CURVE3])

def bezier_cubic(P1, P2, P3, P4):
    return Path([P1, P2, P3, P4],
                [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4])

def plot_path(path, labels=None):
    Xs, Ys = zip(*path.vertices)
    fig = plt.figure()
    ax  = fig.add_subplot(111, aspect='equal')
    ax.set_xlim(min(Xs)-0.2, max(Xs)+0.2)
    ax.set_ylim(min(Ys)-0.2, max(Ys)+0.2)
    patch = patches.PathPatch(path, facecolor='none', linewidth=2)
    ax.add_patch(patch)
    ax.plot(Xs, Ys, 'o--', color='blue', linewidth=1)
    if labels:
        for k in range(len(labels)):
            ax.text(path.vertices[k][0]-0.1,
                path.vertices[k][1]-0.1, 
                labels[k])
    plt.show()

P1 = (0.0, 0.0)
P2 = (1.0, 1.0)
P3 = (2.0, 0.0)
path_1 = bezier_parabola(P1, P2, P3)
plot_path(path_1, labels=['P1', 'P2', 'P3'])

P4 = (2.0, -1.0)
P5 = (3.0, 0.0)
path_2 = bezier_cubic(P1, P2, P4, P5)
plot_path(path_2, labels=['P1', 'P2', 'P4', 'P5'])

Q1 = P5
Q2 = (4.0, 0.0)
Q3 = (5.0, -1.0)
Q4 = (6.0, 0.0)
path_3 = bezier_cubic(P1, P2, P3, P5)
path_4 = bezier_cubic(Q1, Q2, Q3, Q4)
plot_path(Path.make_compound_path(path_3, path_4))

def rotation(point, angle):
    return (np.cos(angle)*point[0] - np.sin(angle)*point[1], 
            np.sin(angle)*point[0] + np.cos(angle)*point[1])

new_Ps = [rotation(P, np.pi/3) for P in path_3.vertices]
new_Qs = [rotation(Q, np.pi/3) for Q in path_4.vertices]
path_5 = bezier_cubic(*new_Ps)
path_6 = bezier_cubic(*new_Qs)
plot_path(Path.make_compound_path(path_5, path_6))