# -*- coding: utf-8 -*-
"""
Created on Mon Jan 25 15:43:19 2021

@author: Rylei
"""

# https://stackoverflow.com/questions/23073170/calculate-bounding-polygon-of-alpha-shape-from-the-delaunay-triangulation
def alpha_shape(points, triangulation, alpha, only_outer=True):
    """
    Compute the alpha shape (concave hull) of a set of points.
    :param points: np.array of shape (n,2) points.
    :param alpha: alpha value.
    :param only_outer: boolean value to specify if we keep only the outer border
    or also inner edges.
    :return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
    the indices in the points array.
    """
    assert points.shape[0] > 3, "Need at least four points"
    
    def add_edge(edges, i, j):
        """
        Add an edge between the i-th and j-th points,
        if not in the list already
        """
        if (i, j) in edges or (j, i) in edges:
            # already added
            # assert (j, i) in edges, "Can't go twice over same directed edge right?"
            if only_outer:
                # if both neighboring triangles are in shape, it's not a boundary edge
                edges.remove((j, i))
            return
        edges.add((i, j))

    tri = triangulation
    edges = set()
    # Loop over triangles:
    # ia, ib, ic = indices of corner points of the triangle
    for ia, ib, ic in tri.vertices:
        pa = points[ia]
        pb = points[ib]
        pc = points[ic]
        # Computing radius of triangle circumcircle
        # www.mathalino.com/reviewer/derivation-of-formulas/derivation-of-formula-for-radius-of-circumcircle
        a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
        b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
        c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
        s = (a + b + c) / 2.0
        val = s * (s - a) * (s - b) * (s - c)
        if val <=0:
            circum_r = float('nan')
        else:
            area = np.sqrt(s * (s - a) * (s - b) * (s - c))
            circum_r = a * b * c / (4.0 * area)
        if circum_r < alpha:
            add_edge(edges, ia, ib)
            add_edge(edges, ib, ic)
            add_edge(edges, ic, ia)
            
    plt.figure()
    plt.plot(points[:, 0], points[:, 1], '.')
    for i, j in edges:
        plt.plot(points[[i, j], 0], points[[i, j], 1], 'r', linewidth=3.5)
        plt.plot(points[:,0], points[:,1], 'ok')
        plt.triplot(points[:,0], points[:,1], tri.simplices,c='k', linestyle='dotted')



    # plt.show()
    return edges

import numpy as np
points = np.array([[0.1, 0.15],[0.3,1.25],[0.9,-0.5],[0.6,1.15],[0.35,-0.5], [-0.2,0.5], [0, 1.1], [1.3, 0], [1, 1], [0.5,0.7], [1,0.5],[0.2,0.6],[0.5,0.1]])
from scipy.spatial import Delaunay
tri = Delaunay(points)

import matplotlib.pyplot as plt
# plt.triplot(points[:,0], points[:,1], tri.simplices)
# plt.plot(points[:,0], points[:,1], 'ok')
# plt.show()

alpha_shape(points, tri, 0.5)


