"""
Library functions for graham algorithm implementation (find the convex hull of a given list of points).
"""

from copy import deepcopy
import numpy as np


def clean_ROI(image):
    """
    Remove instruments borders from an observation.
    """
    H, J = [], []

    shape = np.array(image.shape)
    row, col = np.indices(shape)

    for i in range(0, shape[0]):
        r = row[i, :][image[i, :] > 0.]
        c = col[i, :][image[i, :] > 0.]
        if len(r) > 1 and len(c) > 1:
            H.append((r[0], c[0]))
            H.append((r[-1], c[-1]))
    H = np.array(H)
    for j in range(0, shape[1]):
        r = row[:, j][image[:, j] > 0.]
        c = col[:, j][image[:, j] > 0.]
        if len(r) > 1 and len(c) > 1:
            J.append((r[0], c[0]))
            J.append((r[-1], c[-1]))
    J = np.array(J)
    xmin = np.min([H[:, 1].min(), J[:, 1].min()])
    xmax = np.max([H[:, 1].max(), J[:, 1].max()])+1
    ymin = np.min([H[:, 0].min(), J[:, 0].min()])
    ymax = np.max([H[:, 0].max(), J[:, 0].max()])+1
    return np.array([xmin, xmax, ymin, ymax])


# Define angle and vectors operations
def vector(A, B):
    """
    Define a vector from the inputed points A, B.
    """
    return (B[0] - A[0], B[1] - A[1])


def determinant(u, v):
    """
    Compute the determinant of 2 vectors.
    --------
    Inputs:
    u, v : vectors (list of length 2)
        Vectors for which the determinant must be computed.
    ----------
    Returns:
    determinant : real
        Determinant of the input vectors.
    """
    return u[0] * v[1] - u[1] * v[0]


def pseudo_angle(A, B, C):
    """
    Define a pseudo-angle from the inputed points A, B, C.
    """
    u = vector(A, B)
    v = vector(A, C)
    return determinant(u, v)


def distance(A, B):
    """
    Compute the euclidian distance betwwen 2 points of the plan.
    ----------
    Inputs:
    A, B : points (list of length 2)
        Points between which the distance must be computed.
    ----------
    Returns:
    distance : real
        Euclidian distance between A, B.
    """
    x, y = vector(A, B)
    return np.sqrt(x ** 2 + y ** 2)


# Define lexicographic and composition order
def lexico(A, B):
    """
    Define the lexicographic order between points A, B.
    A <= B if:
        - y_a < y_b
    or  - y_a = y_b and x_a < x_b
    ----------
    Inputs:
    A, B : points (list of length 2)
        Points compared for the lexicographic order
    ----------
    Returns:
    lexico : boolean
        True if A<=B, False otherwise.
    """
    return (A[1] < B[1]) or (A[1] == B[1] and A[0] <= B[0])


def min_lexico(s):
    """
    Return the minimum of a list of points for the lexicographic order.
    ----------
    Inputs:
    s : list of points
        List of points in which we look for the minimum for the lexicographic
        order
    ----------
    Returns:
    m : point (list of length 2)
        Minimum for the lexicographic order in the input list s.
    """
    m = s[0]
    for x in s:
        if lexico(x, m):
            m = x
    return m


def comp(Omega, A, B):
    """
    Test the order relation between points such that A <= B :
    This is True if :
        - (Omega-A,Omega-B) is a right-hand basis of the plan.
    or  - (Omega-A,Omega-B) are linearly dependent.
    ----------
    Inputs:
    Omega, A, B : points (list of length 2)
        Points compared for the composition order
    ----------
    Returns:
    comp : boolean
        True if A <= B, False otherwise.
    """
    t = pseudo_angle(Omega, A, B)
    if t > 0:
        return True
    elif t < 0:
        return False
    else:
        return distance(Omega, A) <= distance(Omega, B)


# Implement quicksort
def partition(s, left, right, order):
    """
    Take a random element of a list 's' between indexes 'left', 'right' and place it
    at its right spot using relation order 'order'. Return the index at which
    it was placed.
    ----------
    Inputs:
    s : list
        List of elements to be ordered.
    left, right : int
        Index of the first and last elements to be considered.
    order : func: A, B -> bool
        Relation order between 2 elements A, B that returns True if A<=B,
        False otherwise.
    ----------
    Returns:
    index : int
        Index at which have been placed the element chosen by the function.
    """
    i = left - 1
    for j in range(left, right):
        if order(s[j], s[right]):
            i = i + 1
            temp = deepcopy(s[i])
            s[i] = deepcopy(s[j])
            s[j] = deepcopy(temp)
    temp = deepcopy(s[i+1])
    s[i+1] = deepcopy(s[right])
    s[right] = deepcopy(temp)
    return i + 1


def sort_aux(s, left, right, order):
    """
    Sort a list 's' between indexes 'left', 'right' using relation order 'order' by
    dividing it in 2 sub-lists and sorting these.
    """
    if left <= right:
        # Call partition function that gives an index on which the list will be divided
        q = partition(s, left, right, order)
        sort_aux(s, left, q - 1, order)
        sort_aux(s, q + 1, right, order)


def quicksort(s, order):
    """
    Implement the quicksort algorithm on the list 's' using the relation order
    'order'.
    """
    # Call auxiliary sort on whole list.
    sort_aux(s, 0, len(s) - 1, order)


def sort_angles_distances(Omega, s):
    """
    Sort the list of points 's' for the composition order given reference point
    Omega.
    """
    def order(A, B): return comp(Omega, A, B)
    quicksort(s, order)


# Define fuction for stacks (use here python lists with stack operations).
def empty_stack(): return []
def stack(S, A): S.append(A)
def unstack(S): S.pop()
def stack_top(S): return S[-1]
def stack_sub_top(S): return S[-2]


# Alignement handling
def del_aux(s, k, Omega):
    """
    Returne the index of the first element in 's' after the index 'k' that is
    not linearly dependent of (Omega-s[k]), first element not aligned to Omega
    and s[k].
    ----------
    Inputs:
    s : list
        List of points
    k : int
        Index from which to look for aligned points in s.
    Omega : point (list of length 2)
        Reference point in 's' to test alignement.
    ----------
    Returns:
    index : int
        Index of the first element not aligned to Omega and s[k].
    """
    p = s[k]
    j = k + 1
    while (j < len(s)) and (pseudo_angle(Omega, p, s[j]) == 0):
        j = j + 1
    return j


def del_aligned(s, Omega):
    """
    Deprive a list of points 's' from all intermediary points that are aligned
    using Omega as a reference point.
    ----------
    Inputs:
    s : list
        List of points
    Omega : point (list of length 2)
        Point of reference to test alignement.
    -----------
    Returns:
    s1 : list
        List of points where no points are aligned using reference point Omega.
    """
    s1 = [s[0]]
    n = len(s)
    k = 1
    while k < n:
        # Get the index of the last point aligned with (Omega-s[k])
        k = del_aux(s, k, Omega)
        s1.append(s[k - 1])
    return s1


# Implement Graham algorithm
def convex_hull(H):
    """
    Implement Graham algorithm to find the convex hull of a given list of
    points, handling aligned points.
    ----------
    Inputs:
    A : list
        List of points for which the the convex hull must be returned.
    ----------
    S : list
        List of points defining the convex hull of the input list A.
    """
    S = empty_stack()
    Omega = min_lexico(H)
    sort_angles_distances(Omega, H)
    A = del_aligned(H, Omega)
    if len(A) < 3:
        return H
    else:
        stack(S, A[0])
        stack(S, A[1])
        stack(S, A[2])
        for i in range(3, len(A)):
            while pseudo_angle(stack_sub_top(S), stack_top(S), A[i]) <= 0:
                unstack(S)
            stack(S, A[i])
        return S


def image_hull(image, step=5, null_val=0., inside=True):
    """
    Compute the convex hull of a 2D image and return the 4 relevant coordinates
    of the maximum included rectangle (ie. crop image to maximum rectangle).
    ----------
    Inputs:
    image : numpy.ndarray
        2D array of floats corresponding to an image in intensity that need to
        be cropped down.
    step : int, optional
        For images with straight edges, not all lines and columns need to be
        browsed in order to have a good convex hull. Step value determine
        how many row/columns can be jumped. With step=2 every other line will
        be browsed.
        Defaults to 5.
    null_val : float, optional
        Pixel value determining the threshold for what is considered 'outside'
        the image. All border pixels below this value will be taken out.
        Defaults to 0.
    inside : boolean, optional
        If True, the cropped image will be the maximum rectangle included
        inside the image. If False, the cropped image will be the minimum
        rectangle in which the whole image is included.
        Defaults to True.
    ----------
    Returns:
    vertex : numpy.ndarray
        Array containing the vertex of the maximum included rectangle.
    """
    H = []
    shape = np.array(image.shape)
    row, col = np.indices(shape)
    for i in range(0, int(min(shape)/2), step):
        r1, r2 = row[i, :][image[i, :] > null_val], row[-i, :][image[-i, :] > null_val]
        c1, c2 = col[i, :][image[i, :] > null_val], col[-i, :][image[-i, :] > null_val]
        if r1.shape[0] > 1:
            H.append((r1[0], c1[0]))
            H.append((r1[-1], c1[-1]))
        if r2.shape[0] > 1:
            H.append((r2[0], c2[0]))
            H.append((r2[-1], c2[-1]))
    S = np.array(convex_hull(H))

    x_min, y_min = S[:, 0] < S[:, 0].mean(), S[:, 1] < S[:, 1].mean()
    x_max, y_max = S[:, 0] > S[:, 0].mean(), S[:, 1] > S[:, 1].mean()
    # Get the 4 extrema
    # S0 = S[x_min*y_min][np.argmin(S[x_min*y_min][:, 0])]
    # S1 = S[x_min*y_max][np.argmax(S[x_min*y_max][:, 1])]
    # S2 = S[x_max*y_min][np.argmin(S[x_max*y_min][:, 1])]
    # S3 = S[x_max*y_max][np.argmax(S[x_max*y_max][:, 0])]
    S0 = S[x_min*y_min][np.abs(0-S[x_min*y_min].sum(axis=1)).min() == np.abs(0-S[x_min*y_min].sum(axis=1))][0]
    S1 = S[x_min*y_max][np.abs(shape[1]-S[x_min*y_max].sum(axis=1)).min() == np.abs(shape[1]-S[x_min*y_max].sum(axis=1))][0]
    S2 = S[x_max*y_min][np.abs(shape[0]-S[x_max*y_min].sum(axis=1)).min() == np.abs(shape[0]-S[x_max*y_min].sum(axis=1))][0]
    S3 = S[x_max*y_max][np.abs(shape.sum()-S[x_max*y_max].sum(axis=1)).min() == np.abs(shape.sum()-S[x_max*y_max].sum(axis=1))][0]
    # Get the vertex of the biggest included rectangle
    if inside:
        f0 = np.max([S0[0], S1[0]])
        f1 = np.min([S2[0], S3[0]])
        f2 = np.max([S0[1], S2[1]])
        f3 = np.min([S1[1], S3[1]])
    else:
        f0 = np.min([S0[0], S1[0]])
        f1 = np.max([S2[0], S3[0]])
        f2 = np.min([S0[1], S2[1]])
        f3 = np.max([S1[1], S3[1]])

    return np.array([f0, f1, f2, f3]).astype(int)