SharedResumen Geometría Computacional.sagewsOpen in CoCalc



                                #              FUNCIONES BASICAS

# Calcula la distancia euclidea entre dos puntos
def dist(A,B):
    return sqrt((A[0]-B[0])**2+(A[1]-B[1])**2)

# Evitamos tener que hacer la raiz cuadrada cuando tenemos que comparar distancias
def dist2(A,B):
     return (A[0]-B[0])**2+(A[1]-B[1])**2
# Nos calcula el area signada de los punto A, B y C con esa orientación. Siendo positiva si la orientación sigue el sentido antihorario , negativa en caso contrario y cero si los puntos están alineados.
def sarea(A,B,C):
    return (1/2)*((B[0]-A[0])*(C[1]-A[1])-(C[0]-A[0])*(B[1]-A[1]))
# Nos da el signo del area signada. Siendo -1 si va en sentido horario, 1 en sentido antihorario y 0 si están alineados.
def orientation(A,B,C):
    return sign((1/2)*((B[0]-A[0])*(C[1]-A[1])-(C[0]-A[0])*(B[1]-A[1])))
# Calcula el punto medio de A y B.
def midPoint(A,B):
    return [((A[0]+B[0])/2),((A[1]+B[1])/2)]

#Primero comprobamos que la coordenada x de P esté entre las coordenadas x de los puntos que definen el segmento. 
#Después repetimos el mismo proceso para las coordenadas y.
#Por último comprobamos si están alineados calculando la orientación.
#Devolviendo True si el punto está en el segmento y False si no lo está.

def inSegment(P,s):
        if ((s[0][0]<=P[0]<=s[1][0] or s[1][0]<=P[0]<=s[0][0]) & (s[1][1]<=P[1]<=s[0][1] or s[0][1]<=P[1]<=s[1][1])):
            if orientation(P,s[0],s[1])==0:
                return True
            else:
                return False
        else:
            return False


#Calculamos las orientaciones del conjutno formado por  P y cada par distinto de vértices y si dan a la vez 1 y -1 quiere decir
#que el punto no pertenece al triángulo. Devuelve True si pertenece al triángulo y False en caso contrario.
def inTriangle(P,t):
    L=(orientation(t[0],t[1],P),orientation(t[0],P,t[2]),orientation(P,t[1],t[2]))
    if (1 in L and -1 in L): return False
    return True
# Si se intersecan hay dos casos:
# -Se cortan en un punto que no es extremo de ningún segmento. Para ello hacemos la orientación de los extremos de un            segmento con cada uno de los extremos del otro.Si son distintas, hacemos el mismo paso pero invirtiendo los papeles. Si      vuelven a ser disntintas, quiere decir que se cortan en un punto.
# -Uno de los extremos de un segmento está contenido en el otro. Para ello usamos la función inSegment y vemos si se cortan o    no.

def segmentIntersectionTest(a,b):
    if (orientation(a[0],a[1],b[0])!=orientation(a[0],a[1],b[1])) & (orientation(b[0],b[1],a[0])!=orientation(b[0],b[1],a[1])):
        return True
    else:
        return inSegment(a[0],b) or inSegment(a[1],b) or inSegment(b[0],a) or inSegment(b[1],a)

# Sacamos los vectores de cada recta. Comprobamos si el determinante de esos dos vectores es cero, sino lo es se cortan en un punto. Hallamos el punto resolviendo el sistema Ax=b.
def lineIntersection(r,s):
    A=vector(r[0])
    B=vector(r[1])
    C=vector(s[0])
    D=vector(s[1])
    AB=B-A
    CD=D-C
    if(det(matrix([AB,CD]))==0):
        print("no se cortan en un punto")
        return
    else:
        AA=matrix([[A[1]-B[1],B[0]-A[0]],[C[1]-D[1],D[0]-C[0]]])
        bb=vector([A[1]*B[0]-A[0]*B[1],C[1]*D[0]-C[0]*D[1]])
        return list(AA\bb)

# Calcula el volumen signado de 4 punto en tres dimensiones. Siendo positivo si la orientación es en sentido antihorario, negativa en sentido horario y cero si están alienados.
def svolume(a,b,c,d):
    Ap = [a[0],a[1],a[0]**2+a[1]**2,1]
    Bp = [b[0],b[1],b[0]**2+b[1]**2,1]
    Cp = [c[0],c[1],c[0]**2+c[1]**2,1]
    Dp = [d[0],d[1],d[0]**2+d[1]**2,1]
    M =matrix([Ap,Bp,Cp,Dp])
    return M.det()
#incircle test. Entrada cuatro puntos a,b,c,d y salida 1, -1 o 0 dependiendo de que el punto d sea interior, exterior o este en la frontera del circulo que determinan a,b y c.


def inCircle(a,b,c,d):
    if (sarea(a,b,c) == 0):
        print ("error")
        return
    elif ((sarea(a,b,c)*svolume(a,b,c,d))>0):

        return int(1)
    elif ((sarea(a,b,c)*svolume(a,b,c,d))==0):

        return int(0)
    else:
        return int(-1)

#Devulve el circuncentro de tres puntos.Primero comprobamos que no estén alineados y luego utilizamos la fórmula del circuncentro.
def circumcenter(a,b,c):
    sa=sarea(a,b,c)
    if (sa==0):
        print ("alineados")
        return
    cx=matrix([[1,1,1],[a[1],b[1],c[1]],[a[0]**2+a[1]**2,b[0]**2+b[1]**2,c[0]**2+c[1]**2]]).det()/(-4*sa)
    cy=matrix([[1,1,1],[a[0],b[0],c[0]],[a[0]**2+a[1]**2,b[0]**2+b[1]**2,c[0]**2+c[1]**2]]).det()/(4*sa)
    return [cx,cy]
#                                                ORDENACION GEOMETRICA
#1.- Ordenar una lista de números de menor a mayor.

# Ordena una lista de número de menor a mayor. Empleando la estrategia de divide y vencerás.
def mergeSort(alist):
    
    if len(alist)>1:
        mid = len(alist)//2
        lefthalf = alist[:mid]
        righthalf = alist[mid:]

        mergeSort(lefthalf)
        mergeSort(righthalf)

        i=0
        j=0
        k=0
        while i < len(lefthalf) and j < len(righthalf):
            if lefthalf[i] < righthalf[j]:
                alist[k]=lefthalf[i]
                i=i+1
            else:
                alist[k]=righthalf[j]
                j=j+1
            k=k+1

        while i < len(lefthalf):
            alist[k]=lefthalf[i]
            i=i+1
            k=k+1

        while j < len(righthalf):
            alist[k]=righthalf[j]
            j=j+1
            k=k+1
    return alist
#2.- Ordenar una lista de números de menor a mayor, sabiendo que se encuentra ordenada de mayor a
#menor.
# Como sabemos que está ordenada de mayor a menor, solo hay que darle la vuelta con reverse.
def ejercicio2(L):
    L.reverse()
    return L
#Concatenamos las dos listas y le aplicamos el metodo mergeSort
def ejercicio3(L1,L2):
    L1.extend(L2)
    mergeSort(L1)
    return L1

# Tan solo hay que hacer un sorted de la lista
def ejercicio4(p):
    return sorted(p)
def sortv(p,v):
    
    return

def angularSort(p,c):
    def compare(p,q):
        if p==q:
            return int(0)
        if p[1]>c[1]>=q[1]:
            return int(1)
        elif q[1]>c[1]>=p[1]:
            return int(-1)
        elif orientation(c,p,q)==1:
            return int(1)
        elif orientation(c,p,q)==-1:
            return int(-1)
        elif dist2(c,p)<dist2(c,q):
            return int(1)
        else:
            return int(-1)
    
    q=deepcopy(p)
    q=sorted(q,cmp=compare)  
    q.reverse()
    return q

 
# para visualizar
#   una lista de puntos p:   point(p)
#   una poligonal definida por una lista de puntos p:  line(p)
# se pueden añadir parametros (color, tamaño del punto, grosor de la linea, aspecto,...
point([[random(),random()] for i in range(100)],aspect_ratio=1,color="red",size=100,faceted=100)
 
line([[random(),random()] for i in range(100)],aspect_ratio=1,color="red",thickness=2)
 



# poligonizacion X-monotona
#cogemos los puntos máximo y mínimo de la lista de puntos y a partir de ellos 
#dividimos los puntos en dos listas según el área signada respecto de estos dos puntos.
#Ordenamos ambas listas y unimos.
def polygonization(p):
    U=[]
    L=[]
    #Buscamos el valor max y min de p
    maxim=max(p)
    minim=min(p)
    for i in range(len(p)):
        if sarea(maxim,p[i],minim)>=0: U.append(p[i])
        else: L.append(p[i])
    U= sorted(U)
    L= sorted(L)
    U.reverse()
    return L+U       
    
# poligonizacion estrellada
# cogemos el punto medio de la lista de puntos y aplicamos angularsort

def starPolygonization(p):
    C=midPoint(p[0],p[1])
    return angularSort(p,C)


#Le pasamos un poligono y una recta sobre la que se hará el corte. El resultado es el polígono formado por el corte de la recta sobre el anterior polígono.
# Primero hallamos el area signada para todos los puntos.
# -Si el area signada con la recta es mayor o igual que cero, añadimos el punto. Ahora nos fijamos en el siguiente punto, si    el area signada con la recta es negativo, quiere decir que la recta que une el punto anterior con el actual, corta a la      recta. Calculamos el punto de intersección y le añadimos a la lista de puntos.
# -Sino era mayor, quiere decir que está por debajo. Hallamos el sarea del siguiente y si es positivio, añadimos el punto de interseccion de la recta definida por los puntos y la recta de corte inicial.
def clipping(P,r):
    n=len(P)
    L=[]
    for i in range(n):
        if sarea(r[0],r[1],P[i])>=0:
            L.append(P[i])
            if sarea(r[0],r[1],P[(i+1)%n])<0:
                M=lineIntersection(r,[P[(i+1)%n],P[i]])
                L.append(M)
        elif sarea(r[0],r[1],P[(i+1)%n])>=0:
            N=lineIntersection(r,[P[(i+1)%n],P[i]])
            L.append(N)
    return L

### funcion auxiliar que detecta si un vértice es oreja en tiempo O(n)
# Primero escogemos un punto.
# Después calculamos el area signada del anterior a ese punto, el punto y el siguiente. Si es negativa o iguala cero, no es     oreja.
# Si lo anterior no era cierto, aún queda por comprobar que ningún vértice del polígono esté dentro del triángulo que forman.   Si ninguno de los puntos estaba dentro, es un vértice oreja.
def earTest(P,i):
    n=len(P)
    if sarea(P[i-1],P[i],P[(i+1)%n])<=0:return False
    j=0
    while j<n:
            if j in [(i-1)%n,i,(i+1)%n]:j+=1
            else:
                if inTriangle(P[j],[P[i-1],P[i],P[(i+1)%n]]):return False
                else:
                    j+=1
    return True
# devolver el indice del vertice oreja
# Aplicamos el test anterior a todos los puntos del polígono y si alguno es oreja, imprimimos el indice de la primera oreja     encontrada.
def ear(P):
    n=len(P)
    a=0
    m=0
    while(a<n)and(m==0):
        if (earTest(P,a))==True: m=1
        else:a=a+1
        
    return a
            
 
# Hacemos una copia del poligono
# Buscamos una oreja y metemos en la lista un vector con los vertices de la oreja. A continuación eliminamos el vértice oreja   y repetimos este paso hasta que solo queden 3 vertices(vertices necesarios para formar un triángulo) y los añadimos           también.
def otectomyTriangulation(P):
    p=deepcopy(P)
    T=[]
    n=len(p)
    while n>3:
        e=ear(p)
        T.append([p[e-1],p[e],p[(e+1)%n]])
        del p[e]
        n -= 1 
    T.append(p)
    return T

# solucion 1: con la funcion ear(P)

def diagonal_1(P):
    n=len(P)
    a=ear(P)
    
    return [P[a-1],P[(a+1)%n]]
#Fuerza bruta búsqueda de vértices
# Tenemos que hallar los puntos interiores. Para ello hacemos todas las combinaciones posibles de 3 puntos y vamos probando con los restantes si son interiores al triangulo que forman. Así conseguimos una lista con puntos interiores y los restantes seran los puntos de la envolvente. Elegimos la ordenacion angular para hacer al envolvente.
def convexHullBF1(p):
    n=len(p)
    interiores=[]
    for i in Combinations(n,3):
        for j in range(n):
            if j <> i[0] and j <> i[1] and j <> i[2] and inTriangle(p[j],[p[i[0]],p[i[1]],p[i[2]]]):
                if j not in interiores:
                    interiores.append(j)
    C=[]               
    for i in range(n):
        if i not in interiores:
            C.append(i)
    CH=[p[i] for i in C]
    M=midPoint(CH[0],CH[1])
    CH=angularSort(CH,M)
    return CH
#Fuerza bruta búsqueda de aristas
# Hacemos un bucle con todas las combinaciones de aristas del poligono posibles. Luego cogemos otro punto y hallamos la orientación. Si es distinta de cero, comprobamos si ya teniamos esa orientación, 
def convexHullBF2(P):
    n=len(P)
    E=[]
    for i in Combinations(n,2):
        L=[]
        j=0
        while j < n:
            A = orientation(P[i[0]], P[i[1]], P[j])
            if A == 0:
                j += 1
                continue
            else:
                if A not in L:
                    L.append(A)
                if len(L) > 1:
                    L=[]
                    j += 1
                    break
            j += 1
        if len(L) == 1:
            E.append(i)
    C=[]
    for i in E:
        C.append(i[0])
        C.append(i[1])
    C = sorted(C)
    CC = [C[0]]
    for i in range(1, len(C)):
        if C[i] <> C[i-1]:
            CC.append(C[i])
    CH = starPolygonization([P[i] for i in CC])
    return CH
 
#Algoritmo de Graham
#Este algoritmo utiliza una lista en la que va almacenando y ordenando correctamente los puntos que constituyen los extremos del cierre convexo.
#Primero elegimos un primer punto para almacenar en la lista. A continuación vamos almacenando el resto en la lista, ordenados según el valor del ángulo que forman con el primer punto que escogimos.
#Después recorremos la lista que hemos formado, tomando en cada momento tres puntos: Inicial, Medio y Final.
#Comprobamos la orientación que forman Inicial, Medio y Final (IMF).
#Si es positiva, Medio pasará a ser Inicial, Final pasará a ser Medio y el siguiente elemento de la lista pasará a ser Final.
#Si es negativa, Medio será borrado de la lista, Inicial pasará a ser Medio, el anterior elemento de la lista pasará a ser Inicial, y Final queda inalterado.
def convexHullGraham(p):
    P=deepcopy(p)
    from operator import itemgetter
    M = min(P, key=itemgetter(1))
    P = angularSort(P, M)
    i = 1 
    while P[i] <> M:
        if orientation(P[i-1], P[i], P[i+1])<0:
            del P[i]
            i -= 1
        else:
            i += 1
    return P
#Quickhull 
#Tomamos los puntos más a la izquierda y el más a la derecha (ambos podemos
#confirmar que forman parte de la envolvente). Con recursividad y a partir de estos dos, cogiendo
#puntos de P por encima y por debajo encontramos el cierre convexo.

def dome(p,A,B):
    def comp(a):
        return sarea(B,A,a)
    if len(p)==0:
        return [B]
    M=max(p,key=comp)
    
    L1,L2=[],[]
    for i in p:
        if sarea(A,i,M)>0:
            L1.append(i)
        if sarea(M,i,B)>0:
            L2.append(i)
            
    return dome(L1, A, M) + dome(L2,M,B)+[B]
    
def quickhull(p):
    A,B=max(p), min(p)
    return dome(p,A,B)+dome(p,B,A)

#segunda parte
# funcion que crea un DCEL, para un poligono
def dcel(P):
    n=len(P)
    V=[[P[i],i] for i in range(len(P))]
    e=[[i,n+i,(i-1)%n,(i+1)%n,1]for i in range(n)]+[[(i+1)%n,i,n+(i+1)%n,n+(i-1)%n,0]for i in range(n)]
    f=[n,0]
    return [V,e,f]

# funciones para referirse a los elementos asociados a un elemento del DCEL

# indice del origen de una arista e
def origin(e,D):
    return D[1][e][0]

# coordenadas del origen de la arista e
def originCoords(e,D):
    return D[0][origin(e,D)][0]

# arista gemela de la arista e    
def twin(e,D):
    return D[1][e][1]

# arista previa de la arista e    
def prev(e,D):
    return D[1][e][2]

# arista siguiente de la arista e    
def next(e,D):
    return D[1][e][3] 

# indice de la cara de cuyo borde forma parte la arista e    
def face(e,D):
    return D[1][e][4]               

# indice de una de las aristas del borde de la cara c
def edge(c,D):
    return D[2][c]
        
# funcion para dibujar las aristas de un DCEL

def plotDCEL(D):
    return sum(line([originCoords(i,D),originCoords(twin(i,D),D)],aspect_ratio=1) for i in range(len(D[1])))  
    
# funcion para colorear una cara de un DCEL

def plotFace(c,D,col):

    f=D[2][c]
    C=[f]
    f=next(f,D)
    while f <> C[0]:
        C.append(f)
        f=next(f,D)
    
    P=[originCoords(j,D) for j in C]
    return polygon(P,color=col, alpha=.5)     

# funcion para colorear las caras de un DCEL
    
def colorDCEL(D):
    return sum(plotFace(i,D,(random(),random(),random())) for i in range(1,len(D[2])))
# funcion para dividir una cara del DCEL D por una diagonal
# e1 y e2 son las aristas cuyos orígenes son los extremos de la diagonal que divide la cara

def splitFace(e1,e2,D):
    
    # si no son aristas de la misma cara o si son adyacentes sus origenes no definen una diagonal
    if face(e1,D) <> face(e2,D) or origin(e2,D) == origin(twin(e1,D),D) or origin(e1,D) == origin(twin(e2,D),D):
        print "no diagonal"
        return
    
    nv, ne, nf = len(D[0]), len(D[1]), len(D[2])
    preve1 = prev(e1,D)
    preve2 = prev(e2,D)
    k=face(e1,D)
    
    # añadimos las aristas nuevas
    D[1].append([origin(e1,D),ne+1,preve1,e2,k])
    D[1].append([origin(e2,D),ne,preve2,e1,nf])
    
    # modificamos aristas afectadas
    D[1][preve1][3]=ne
    D[1][e1][2]=ne+1
    D[1][preve2][3]=ne+1
    D[1][e2][2]=ne
    i=e1
    while i<>ne+1:
        D[1][i][4]=nf
        i=next(i,D)
    
    #modificamos la cara afectada
    D[2][k]=ne
    
    # añadimos la nueva cara
    D[2].append(ne+1)
def faceEdges(c,D):
    e=edge(c,D)
    C=[e]
    e=next(e,D)
    while e <> C[0]:
        C.append(e)
        e=next(e,D)
    return C

def faceVertices(c,D):
    A = faceEdges(c,D)
    V = [D[1][A[i]][0] for i in range(len(A))]
    return V

def faceVerticesCoords(c,D):
    A = faceEdges(c,D)
    coords = [originCoords(A[i],D) for i in range(len(A))]
    return A

def faceNeighbours(c,D):
    A = faceEdges(c,D)
    coords = [twin(A[i],D) for i in range(len(A))]
    return A

def vertexEdges(v,D):
    e0=D[0][v][1]
    E=[e0]
    e=twin(prev(e0,D),D)
    while e<>e0:
        E.append(e)
        e=twin(prev(e,D),D)
    return E

def vertexFan(v,D):
    A=vertexEdges(v,D)
    C=[face(A[i],D) for i in range(len(A))]
    return C

# CONVEX HULL
# esta funcion no sirve en general
# pero si para casos en que la cara externa es un poligono estrellado o monotono
def DCELconvexHull(D): #scan de Graham
    E=faceEdges(0,D)
    V=[[originCoords(i,D),i] for i in E]
    m=min(V)[1]
    e=next(m,D)
    while origin(e,D)<>origin(m,D):
        if sarea(originCoords(prev(e,D),D),originCoords(e,D),originCoords(next(e,D),D))>0:
            e1=prev(e,D)
            e2=next(e,D)
            splitFace(e1,e2,D)
            if origin(prev(e2,D),D)<>origin(m,D):
                e=prev(e2,D)
            else:
                e=e2
        else:
            e=next(e,D)
            

            
def triangulation(p):
    n=len(p)
    P=[[i[1],i]for i in p]
    M=min(P)[0]  #punto de menor ordenada
    P=angularSort(p,M)
    D=dcel(P)
    for i in range(3,n):  #triangulacion de la cara acotada
        splitFace(n-i,n-1,D)
    DCELconvexHull(D)  #scan de Graham
    return D         

   
def triangulation(p):
    n=len(p)
    P=[[i[1],i]for i in p]
    M=min(P)[1]  #punto de menor ordenada
    P=angularSort(p,M)
    D=dcel(P)
    DCELconvexHull(D)  #scan de Graham
    for i in range(3,n):  #triangulacion de la cara acotada
        splitFace(n-i,n-1,D)
    
    return D

 


 
def flip(a,D):
    oa=D[1][a][0]
    ga=D[1][a][1]
    ca=D[1][a][4]
    aa=D[1][a][2]
    pa=D[1][a][3]
    cb=D[1][ga][4]
    ab=D[1][ga][2]
    pb=D[1][ga][3]
    oga=D[1][ga][0]
    D[1][a]=[D[1][aa][0],ga,pa,ab,ca]
    D[1][ga]=[D[1][ab][0],a,pb,aa,cb]
    D[1][pa][2]=ab
    D[1][pa][3]=a
    D[1][aa][2]=ga
    D[1][aa][3]=pb
    D[1][aa][4]=cb
    D[1][pb][2]=aa
    D[1][pb][3]=ga
    D[1][ab][2]=a
    D[1][ab][3]=pa
    D[1][ab][4]=ca
    D[2][ca]=a
    D[2][cb]=ga
    D[0][oa][1]=pb
    D[0][oga][1]=pa

    return
#flipable(a,D)
def flipable(a,D):
    if face(a,D)==0 or face(twin(a,D),D)==0:
        return false
    A,B,C,D = originCoords(a,D),originCoords(prev(twin(a,D),D),D),originCoords(twin(a,D),D),originCoords(prev(a,D),D)
    if sarea(A,B,C)>0 and sarea(B,C,D)>0 and sarea(C,D,A)>0 and sarea(D,A,B)>0:
        return True
    return false
#legal(a,D)
def legal(a,D):
    if flipable(a,D)==-1:
        return true
    A,B,C,D = originCoords(a,D),originCoords(prev(twin(a,D),D),D),originCoords(twin(a,D),D),originCoords(prev(a,D),D)
    if inCircle(A,C,D,B)==false:
        return True
    else:
        return False
#legalize(T)

def legalize(t):
    T=deepcopy(t)
    c=1
    vueltas=0
    flips=0
    while c>0:
        vueltas=vueltas+1
        c=0
        for i in range(len(T[1])):
            if legal(i,T)==-1:
                flip(i,T)
                flips=flips+1
                c=c+1
    print "num vueltas: ",vueltas,"num flips: ",flips
    return T

#delone(p)
    k=triangulation(p)
    k=legalize(k)
    return k

def infiniteVertex(e,D):
    A = originCoords(e,D)
    B = originCoords(next(e,D),D)
    v = vector(B)-vector(A)
    v0 = vector([-v[1],v[0]])
    v0 *= 100/v0.norm()
    M = midPoint(A,B)
    C = vector(M)+v0
    return list(C)

def voronoiRegion(i, D):
    F = vertexEdges(i,D)
    R = []
    for i in F:
        if face(i,D) != 0:
            A,B,C = originCoords(prev(i,D),D),originCoords(i,D),originCoords(next(i,D),D)
            R.append(circumcenter(A,B,C))
        else:
            R.append(infiniteVertex(i,D))
            R.append(infiniteVertex(prev(i,D),D))
    return R

def voronoi(D): #2n
    V=[]
    for i in range(len(D[0])):
        r = voronoiRegion(i,D)
        V.append(r)
    return V

# Calcula la distancia euclidea entre dos puntos
def dist(A,B):
    return sqrt((A[0]-B[0])**2+(A[1]-B[1])**2)
# prueba
dist([0,0],[1,1])
sqrt(2)

2. cuadrado de la distancia entre los puntos


# Si una distancia A es mayor que otra B haciendo la raiz cuadrada, sin la raiz cuadrada A seguirá siendo mayor que B sin la raiz, así evitamos tener que hacer la raiz cuadrada cuando tenemos que comparar distancias.
def dist2(A,B):
     return (A[0]-B[0])**2+(A[1]-B[1])**2
# prueba
dist2([0,0],[1,1])
2

3. area signada

# Nos calcula el area signada de los punto A, B y C con esa orientación. Siendo positiva si la orientación sigue el sentido antihorario , negativa en caso contrario y cero si los puntos están alineados.
def sarea(A,B,C):
    return (1/2)*((B[0]-A[0])*(C[1]-A[1])-(C[0]-A[0])*(B[1]-A[1]))
# prueba
sarea([1,1],[0,0],[1,0])
1/2
# prueba
sarea([1,0],[0,0],[1,1])
-1/2
# Nos da el signo del area signada. Siendo -1 si va en sentido horario, 1 en sentido antihorario y 0 si están alineados.
def orientation(A,B,C):
    return sign((1/2)*((B[0]-A[0])*(C[1]-A[1])-(C[0]-A[0])*(B[1]-A[1])))
#pruebas
orientation([0,0],[1,0],[0,-1]), orientation([1,0],[0,0],[-1,0]), orientation([1,1],[0,0],[1,0])
(-1, 0, 1)

5. punto medio de dos puntos

# Calcula el punto medio de A y B.
def midPoint(A,B):
    return [((A[0]+B[0])/2),((A[1]+B[1])/2)]
#prueba
midPoint([2,4],[4,8])
[3, 6]

#Primero comprobamos que la coordenada x de P esté entre las coordenadas x de los puntos que definen el segmento. 
#Después repetimos el mismo proceso para las coordenadas y.
#Por último comprobamos si están alineados calculando la orientación.
#Devolviendo True si el punto está en el segmento y False si no lo está.

def inSegment(P,s):
        if ((s[0][0]<=P[0]<=s[1][0] or s[1][0]<=P[0]<=s[0][0]) & (s[1][1]<=P[1]<=s[0][1] or s[0][1]<=P[1]<=s[1][1])):
            if orientation(P,s[0],s[1])==0:
                return True
            else:
                return False
        else:
            return False


#pruebas
inSegment([1,1],[[0,0],[2,2]]), inSegment([2,2],[[0,0],[2,2]]), inSegment([1.000000000000000000001,1],[[0,0],[2,2]])
(True, True, False)

7. test de punto en triángulo

#Calculamos las orientaciones del conjutno formado por  P y cada par distinto de vértices y si dan a la vez 1 y -1 quiere decir
#que el punto no pertenece al triángulo. Devuelve True si pertenece al triángulo y False en caso contrario.
def inTriangle(P,t):
    L=(orientation(t[0],t[1],P),orientation(t[0],P,t[2]),orientation(P,t[1],t[2]))
    if (1 in L and -1 in L): return False
    return True
#pruebas
inTriangle([0,1],[[-2,0],[2,0],[0,2]]),inTriangle([1,1],[[-2,0],[2,0],[0,2]]),inTriangle([1,1.0000000000000000000000000000000000000000000001],[[-2,0],[2,0],[0,2]])

(True, True, False)
p= [1,1]
t=[[-2,0],[2,0],[0,2]]

g1=line([t[0]]+[t[1]]+[t[2]]+[t[0]],thickness=2,color="red")+point(p,size=50,color="blue")+text(inTriangle([0,1],[[-2,0],[2,0],[0,2]]), (1,2.5))

p= [0,1]
t=[[-2,0],[2,0],[0,2]]

g2=line([t[0]]+[t[1]]+[t[2]]+[t[0]],thickness=2,color="red")+point(p,size=50,color="blue")+text(inTriangle([0,1],[[-2,0],[2,0],[0,2]]), (1,2.5))

P=(random(),random())
t=[(random(),random()),(random(),random()),(random(),random())]

g3=point(P,color="red") + polygon (t)+text(inTriangle(P,t), (1,2.5))

graphics_array([[g1,g2,g3]]).show(figsize=[10,5])


P=(random(),random())
t=[(random(),random()),(random(),random()),(random(),random())]
print inTriangle(P,t)
point(P) + polygon (t)
False

8. test de intersección de dos segmentos

# Si se intersecan hay dos casos:
# -Se cortan en un punto que no es extremo de ningún segmento. Para ello hacemos la orientación de los extremos de un            segmento con cada uno de los extremos del otro.Si son distintas, hacemos el mismo paso pero invirtiendo los papeles. Si      vuelven a ser disntintas, quiere decir que se cortan en un punto.
# -Uno de los extremos de un segmento está contenido en el otro. Para ello usamos la función inSegment y vemos si se cortan o    no.

def segmentIntersectionTest(a,b):
    if (orientation(a[0],a[1],b[0])!=orientation(a[0],a[1],b[1])) & (orientation(b[0],b[1],a[0])!=orientation(b[0],b[1],a[1])):
        return True
    else:
        return inSegment(a[0],b) or inSegment(a[1],b) or inSegment(b[0],a) or inSegment(b[1],a)
t1=segmentIntersectionTest([[0,0],[2,2]],[[2,0],[0,2]])
p11=line([[0,0]]+[[2,2]],color="cyan")+line([[2,0]]+[[0,2]],color="red")
t11= text(t1, (1,2))
p1=p11+t11

t22=segmentIntersectionTest([[0,0],[2,2]],[[1,1],[1.5,1.5]])
p22=line([[0,0]]+[[2,2]],color="cyan")+line([[1,1]]+[[1.5,1.5]],color="red")
t2= text(t22, (1,2))
p2=p22+t2

t33=segmentIntersectionTest([[0,0],[0,2]],[[0,1],[2,1]])
p33=line([[0,0]]+[[0,2]],thickness=3,color="cyan")+line([[0,1]]+[[2,1]],color="red")
t3=text(t33, (1,2))
p3=p33+t3

t44=segmentIntersectionTest([[0,0],[2,2]],[[2,2],[1,0]])
p44=line([[0,0]]+[[2,2]],color="cyan")+line([[2,2]]+[[1,0]],color="red")
t4= text(t44, (1,2))
p4= p44+t4

t55=segmentIntersectionTest([[0,0],[2,2]],[[0,-1],[2,1]])
p55=line([[0,0]]+[[2,2]],color="cyan")+line([[0,-1]]+[[2,1]],color="red")
t5= text(t55, (1,2))
p5= p55+t5

graphics_array([[p1, p2,p3, p4,p5]]).show(figsize=[10,4])

9. punto de intersección de dos rectas

# Sacamos los vectores de cada recta. Comprobamos si el determinante de esos dos vectores es cero, sino lo es se cortan en un punto. Hallamos el punto resolviendo el sistema Ax=b.
def lineIntersection(r,s):
    A=vector(r[0])
    B=vector(r[1])
    C=vector(s[0])
    D=vector(s[1])
    AB=B-A
    CD=D-C
    if(det(matrix([AB,CD]))==0):
     
        return
    else:
        AA=matrix([[A[1]-B[1],B[0]-A[0]],[C[1]-D[1],D[0]-C[0]]])
        bb=vector([A[1]*B[0]-A[0]*B[1],C[1]*D[0]-C[0]*D[1]])
        return list(AA\bb)



# pruebas
punto=lineIntersection([[0,0],[2,2]],[[2,0],[0,2]])
s1=line([[0,0]]+[[2,2]],color="cyan")+line([[2,0]]+[[0,2]],color="gold")+point(punto,color="red")+text(lineIntersection([[0,0],[2,2]],[[2,0],[0,2]]), (1,2))

punto=lineIntersection([[0,0],[40,40]],[[0,1],[40,44.99]])
s2=line([[0,0]]+[[40,40]],color="cyan")+line([[0,1]]+[[40,44.99]],color="gold")+point(punto,color="red")+text(lineIntersection([[0,0],[40,40]],[[0,1],[40,44.99]]), (30,40))

s3=line([[0,0]]+[[2,2]],color="cyan")+line([[0,1]]+[[2,3]],color="gold")+text(lineIntersection([[0,0],[2,2]],[[0,1],[2,3]]), (1,3))
graphics_array([[s1, s2,s3]]).show(figsize=[10,6])
lineIntersection([[0,0],[2,2]],[[1,1],[3,3]])
no se cortan en un punto
r=[[random(),random()],[random(),random()]]
s=[[random(),random()],[random(),random()]]
line(r,aspect_ratio=1)+line(s)+point(lineIntersection(r,s),color="red",size=20)

10. test de inclusión del punto en círculo

# Calcula el volumen signado de 4 punto en tres dimensiones. Siendo positivo si la orientación es en sentido antihorario, negativa en sentido horario y cero si están alienados.
def svolume(a,b,c,d):
    Ap = [a[0],a[1],a[0]**2+a[1]**2,1]
    Bp = [b[0],b[1],b[0]**2+b[1]**2,1]
    Cp = [c[0],c[1],c[0]**2+c[1]**2,1]
    Dp = [d[0],d[1],d[0]**2+d[1]**2,1]
    M =matrix([Ap,Bp,Cp,Dp])
    return M.det()
svolume([1,2],[1,0],[1,4],[5,4])


64
svolume([0,2],[1,2],[2,2],[3,2])
0
#incircle test. Entrada cuatro puntos a,b,c,d y salida 1, -1 o 0 dependiendo de que el punto d sea interior, exterior o este en la frontera del circulo que determinan a,b y c.


def inCircle(a,b,c,d):
    if (sarea(a,b,c) == 0):
        print ("error")
        return
    elif ((sarea(a,b,c)*svolume(a,b,c,d))>0):

        return int(1)
    elif ((sarea(a,b,c)*svolume(a,b,c,d))==0):

        return int(0)
    else:
        return int(-1)


def circulo3puntos(a,b,c):
    centro=circumcenter(a,b,c)
    radio2=dist(centro,a)
    return circle(centro,radio2,color="cyan")+point(a,color="red",size=30)+point(b,color="red",size=30)+point(c,color="red",size=30)
# pruebas

c1=circulo3puntos([0,0],[0,2],[2,0])+point([1,1],color="black",size=30)+text(inCircle([0,0],[0,2],[2,0],[1,1]), (1,4))

c2=circulo3puntos([-1,0],[1,1],[0,1])+point([-1,-1],color="black",size=30)+text(inCircle([-1,0],[1,1],[0,1],[-1,-1]), (1,2))



c3=circulo3puntos([0,0],[0,2],[2,0])+point([5,5],color="black",size=30)+text(inCircle([0,0],[0,2],[2,0],[5,5]), (3,5))
graphics_array([[c1, c2,c3]]).show(figsize=[10,6])



# A continuación comprobamos el caso en el que los tres puntos estan alienados
inCircle([0,0],[0,2],[0,4],[3,3])

error

11. circuncentro del círculo determinado por tres puntos

#Devulve el circuncentro de tres puntos.Primero comprobamos que no estén alineados y luego utilizamos la fórmula del circuncentro.
def circumcenter(a,b,c):
    sa=sarea(a,b,c)
    if (sa==0):
        print ("alineados")
        return
    cx=matrix([[1,1,1],[a[1],b[1],c[1]],[a[0]**2+a[1]**2,b[0]**2+b[1]**2,c[0]**2+c[1]**2]]).det()/(-4*sa)
    cy=matrix([[1,1,1],[a[0],b[0],c[0]],[a[0]**2+a[1]**2,b[0]**2+b[1]**2,c[0]**2+c[1]**2]]).det()/(4*sa)
    return [cx,cy]
# pruebas

c1=circulo3puntos([0,0],[1,0],[1,1])+point(circumcenter([0,0],[1,0],[1,1]),color="green",size=30)
c1

circumcenter([-1,0],[0,1],[1,0])
[0, 0]
circumcenter([0,0],[1,0],[2,0])
alineados

# 






#                                                ORDENACION GEOMETRICA







#1.- Ordenar una lista de números de menor a mayor.
# Ordena una lista de número de menor a mayor. Empleando la estrategia de divide y vencerás.
def mergeSort(alist):
    
    if len(alist)>1:
        mid = len(alist)//2
        lefthalf = alist[:mid]
        righthalf = alist[mid:]

        mergeSort(lefthalf)
        mergeSort(righthalf)

        i=0
        j=0
        k=0
        while i < len(lefthalf) and j < len(righthalf):
            if lefthalf[i] < righthalf[j]:
                alist[k]=lefthalf[i]
                i=i+1
            else:
                alist[k]=righthalf[j]
                j=j+1
            k=k+1

        while i < len(lefthalf):
            alist[k]=lefthalf[i]
            i=i+1
            k=k+1

        while j < len(righthalf):
            alist[k]=righthalf[j]
            j=j+1
            k=k+1
    return alist
L=[0,1,6,8,4,22,15,39,47,2,2,5,9,47]
L1=[0,1,6,8,4,22,15,39,47,2,2,5,9,47]
print(sorted(L1))
print(mergeSort(L))

[0, 1, 2, 2, 4, 5, 6, 8, 9, 15, 22, 39, 47, 47] [0, 1, 2, 2, 4, 5, 6, 8, 9, 15, 22, 39, 47, 47]

#2.- Ordenar una lista de números de menor a mayor, sabiendo que se encuentra ordenada de mayor a
#menor.
# Como sabemos que está ordenada de mayor a menor, solo hay que darle la vuelta con reverse.
def ejercicio2(L):
    L.reverse()
    return L
L=[22,21,20,19,18,17,15,14,12,13,10,9,5,4,1]
print(ejercicio2(L))

[1, 4, 5, 9, 10, 13, 12, 14, 15, 17, 18, 19, 20, 21, 22]


#Concatenamos las dos listas y le aplicamos el metodo mergeSort
def ejercicio3(L1,L2):
    L1.extend(L2)
    mergeSort(L1)
    return L1
L1=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]
L2=[7,8,9,10,11,12,13,14,15,16]
print(ejercicio3(L1,L2))
[1, 2, 3, 4, 5, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 18]

Dimensión 2:
Dada una lista de n puntos del plano por sus coordenadas cartesianas,

4.- Ordenar lexicográficamente según sus coordenadas cartesianas los n puntos.


# Tan solo hay que hacer un sorted de la lista
def ejercicio4(p):
    return sorted(p)
L=[[1,1],[2,2],[81,0],[2,3],[35,0],[5,0],[7,0],[1,0],[9,0],[2,0],[4,0],[3,0]]
print(ejercicio4(L))
[[1, 0], [1, 1], [2, 0], [2, 2], [2, 3], [3, 0], [4, 0], [5, 0], [7, 0], [9, 0], [35, 0], [81, 0]]

5.- Ordenar los puntos según la dirección de un vector dado.

def sortv(p,v):
    
    return

6.- Ordenar los puntos angularmente respecto de un punto dado


def angularSort(p,c):
    def compare(p,q):
        if p==q:
            return int(0)
        if p[1]>c[1]>=q[1]:
            return int(1)
        elif q[1]>c[1]>=p[1]:
            return int(-1)
        elif orientation(c,p,q)==1:
            return int(1)
        elif orientation(c,p,q)==-1:
            return int(-1)
        elif dist2(c,p)<dist2(c,q):
            return int(1)
        else:
            return int(-1)
    
    q=deepcopy(p)
    q=sorted(q,cmp=compare)  
    q.reverse()
    return q



p=[[int(10*random()),int(10*random())]for i in range(50)]
line(p+[p[0]])+point(p[0],color="red")
a=angularSort(p,p[0])

line(a+[a[0]])+point(p[0],color="red")




                                                    POLIGONIZACION
# poligonizacion X-monotona
#cogemos los puntos máximo y mínimo de la lista de puntos y a partir de ellos 
#dividimos los puntos en dos listas según el área signada respecto de estos dos puntos.
#Ordenamos ambas listas y unimos.
def polygonization(p):
    U=[]
    L=[]
    #Buscamos el valor max y min de p
    maxim=max(p)
    minim=min(p)
    for i in range(len(p)):
        if sarea(maxim,p[i],minim)>=0: U.append(p[i])
        else: L.append(p[i])
    U= sorted(U)
    L= sorted(L)
    U.reverse()
    return L+U       
    
# poligonizacion estrellada
# cogemos el punto medio de la lista de puntos y aplicamos angularsort

def starPolygonization(p):
    C=midPoint(p[0],p[1])
    return angularSort(p,C)

p = [[random(),random()] for i in range(50)]
Pol = polygonization(p)
point(Pol,color="red")+line(Pol+[Pol[0]])+point([p[0],p[1]], color ="green", size=20)
p = [[random(),random()] for i in range(50)]
Pol = starPolygonization(p)
point(Pol,color="red")+line(Pol+[Pol[0]])+point([p[0],p[1]], color ="green", size=20)


#Le pasamos un poligono y una recta sobre la que se hará el corte. El resultado es el polígono formado por el corte de la recta sobre el anterior polígono.
# Primero hallamos el area signada para todos los puntos.
# -Si el area signada con la recta es mayor o igual que cero, añadimos el punto. Ahora nos fijamos en el siguiente punto, si    el area signada con la recta es negativo, quiere decir que la recta que une el punto anterior con el actual, corta a la      recta. Calculamos el punto de intersección y le añadimos a la lista de puntos.
# -Sino era mayor, quiere decir que está por debajo. Hallamos el sarea del siguiente y si es positivio, añadimos el punto de interseccion de la recta definida por los puntos y la recta de corte inicial.
def clipping(P,r):
    n=len(P)
    L=[]
    for i in range(n):
        if sarea(r[0],r[1],P[i])>=0:
            L.append(P[i])
            if sarea(r[0],r[1],P[(i+1)%n])<0:
                M=lineIntersection(r,[P[(i+1)%n],P[i]])
                L.append(M)
        elif sarea(r[0],r[1],P[(i+1)%n])>=0:
            N=lineIntersection(r,[P[(i+1)%n],P[i]])
            L.append(N)
    return L

p = [[random(),random()] for i in range(50)]
Pol = starPolygonization(p)
r=[[random(),random()],[random(),random()]]
CP=clipping(Pol,r)
polygon(CP,alpha=.5)+polygon(Pol,alpha=.5,color="cyan")+line(r,thickness=1,color="red")

3. Dado un polígono, obtener una oreja.

(Un vértice de un polígono se llama oreja si con sus dos vértices adyacentes forman un triángulo contenido en el polígono que no contiene a ningún otro vértice del mismo).


Cálculo de una oreja en tiempo lineal

### funcion auxiliar que detecta si un vértice es oreja en tiempo O(n)
# Primero escogemos un punto.
# Después calculamos el area signada del anterior a ese punto, el punto y el siguiente. Si es negativa o iguala cero, no es     oreja.
# Si lo anterior no era cierto, aún queda por comprobar que ningún vértice del polígono esté dentro del triángulo que forman.   Si ninguno de los puntos estaba dentro, es un vértice oreja.
def earTest(P,i):
    n=len(P)
    if sarea(P[i-1],P[i],P[(i+1)%n])<=0:return False
    j=0
    while j<n:
            if j in [(i-1)%n,i,(i+1)%n]:j+=1
            else:
                if inTriangle(P[j],[P[i-1],P[i],P[(i+1)%n]]):return False
                else:
                    j+=1
    return True
# devolver el indice del vertice oreja
# Aplicamos el test anterior a todos los puntos del polígono y si alguno es oreja, imprimimos el indice de la primera oreja     encontrada.
def ear(P):
    n=len(P)
    a=0
    m=0
    while(a<n)and(m==0):
        if (earTest(P,a))==True: m=1
        else:a=a+1
        
    return a
            
 
p = [[random(),random()] for i in range(50)]
P = starPolygonization(p)
i=0
earTest(P,ear(P))
polygon(P,alpha=.5,color="red")+point(P[ear(P)],size=30,color="green")+line([P[ear(P)-1]]+[P[ear(P)]]+[P[(ear(P)+1)%len(P)]]+[P[ear(P)-1]],color="blue")
True


p = [[random(),random()] for i in range(20)]
P = starPolygonization(p)
i=0
earTest(P,i)
polygon(P,alpha=.5,color="red")+point(P[i],size=30,color="green")+line([P[i-1]]+[P[i]]+[P[(i+1)%len(P)]]+[P[i-1]],color="blue")
False

4. Triangulación por otectomías sucesivas:

Dado un polígono, descomponerlo en triángulos (triangularlo) iterando el cálculo de orejas en el polígono que resulta de eliminar una de ellas.

# Hacemos una copia del poligono
# Buscamos una oreja y metemos en la lista un vector con los vertices de la oreja. A continuación eliminamos el vértice oreja   y repetimos este paso hasta que solo queden 3 vertices(vertices necesarios para formar un triángulo) y los añadimos           también.
def otectomyTriangulation(P):
    p=deepcopy(P)
    T=[]
    n=len(p)
    while n>3:
        e=ear(p)
        T.append([p[e-1],p[e],p[(e+1)%n]])
        del p[e]
        n -= 1 
    T.append(p)
    return T

p = [[random(),random()] for i in range(20)]
P = polygonization(p)
T=otectomyTriangulation(P)
sum(polygon(i,color=(random(),random(),random())) for i in T)+line(P+[P[0]])

Pruebas

5. Dado un polígono, hallar una diagonal interna en tiempo lineal.

# solucion 1: con la funcion ear(P)

def diagonal_1(P):
    n=len(P)
    a=ear(P)
    
    return [P[a-1],P[(a+1)%n]]
p = [[random(),random()] for i in range(20)]
Pol = polygonization(p)
b=diagonal_1(Pol)
point(Pol,color="blue")+line([b[0]]+[b[1]],color="red")+line(Pol+[Pol[0]],color="purple")



        


#Fuerza bruta búsqueda de vértices
# Tenemos que hallar los puntos interiores. Para ello hacemos todas las combinaciones posibles de 3 puntos y vamos probando con los restantes si son interiores al triangulo que forman. Así conseguimos una lista con puntos interiores y los restantes seran los puntos de la envolvente. Elegimos la ordenacion angular para hacer al envolvente.
def convexHullBF1(p):
    n=len(p)
    interiores=[]
    for i in Combinations(n,3):
        for j in range(n):
            if j <> i[0] and j <> i[1] and j <> i[2] and inTriangle(p[j],[p[i[0]],p[i[1]],p[i[2]]]):
                if j not in interiores:
                    interiores.append(j)
    C=[]               
    for i in range(n):
        if i not in interiores:
            C.append(i)
    CH=[p[i] for i in C]
    M=midPoint(CH[0],CH[1])
    CH=angularSort(CH,M)
    return CH
p=[[random(),random()] for i in range(100)]
C=convexHullBF1(p)
line(C+[C[0]])+point(p)







#Fuerza bruta búsqueda de aristas
# Hacemos un bucle con todas las combinaciones de aristas del poligono posibles. Luego cogemos otro punto y hallamos la orientación. Si es distinta de cero, comprobamos si ya teniamos esa orientación, 
def convexHullBF2(P):
    n=len(P)
    E=[]
    for i in Combinations(n,2):
        L=[]
        j=0
        while j < n:
            A = orientation(P[i[0]], P[i[1]], P[j])
            if A == 0:
                j += 1
                continue
            else:
                if A not in L:
                    L.append(A)
                if len(L) > 1:
                    L=[]
                    j += 1
                    break
            j += 1
        if len(L) == 1:
            E.append(i)
    C=[]
    for i in E:
        C.append(i[0])
        C.append(i[1])
    C = sorted(C)
    CC = [C[0]]
    for i in range(1, len(C)):
        if C[i] <> C[i-1]:
            CC.append(C[i])
    CH = starPolygonization([P[i] for i in CC])
    return CH
p=[[random(),random()] for i in range(10)]
C=convexHullBF2(p)
line(C+[C[0]]) + point(p)

def jarvis(p):
    longitud=len(p)
    from operator import itemgetter
    referencia= min(p, key=itemgetter(1))
   
    resultado =[referencia]
    pactual=resultado[0]
    i=0
    while i<longitud:
        sa=[]
        for j in range(longitud):
            if(pactual!=p[i])and(pactual!=p[j])and(p[i]!=p[j]):
                sa.append(orientation(pactual,p[i],p[j]))
        if (0 not in sa)and(1 not in sa):
            if (p[i]!=referencia)and(p[i] not in resultado):
                resultado.append(p[i])
                pactual=p[i]
                i=0
            elif (p[i]==referencia):
                return resultado
                 
                 
            else:
                i=i+1
        if 1 in sa:
            i=i+1
    return resultado



p=[[random(),random()] for i in range(100)]
C=jarvis(p)
line(C+[C[0]],color="red") + point(p,color="green",size=30)
#Algoritmo de Graham
#Este algoritmo utiliza una lista en la que va almacenando y ordenando correctamente los puntos que constituyen los extremos del cierre convexo.
#Primero elegimos un primer punto para almacenar en la lista. A continuación vamos almacenando el resto en la lista, ordenados según el valor del ángulo que forman con el primer punto que escogimos.
#Después recorremos la lista que hemos formado, tomando en cada momento tres puntos: Inicial, Medio y Final.
#Comprobamos la orientación que forman Inicial, Medio y Final (IMF).
#Si es positiva, Medio pasará a ser Inicial, Final pasará a ser Medio y el siguiente elemento de la lista pasará a ser Final.
#Si es negativa, Medio será borrado de la lista, Inicial pasará a ser Medio, el anterior elemento de la lista pasará a ser Inicial, y Final queda inalterado.
def convexHullGraham(p):
    P=deepcopy(p)
    from operator import itemgetter
    M = min(P, key=itemgetter(1))
    P = angularSort(P, M)
    i = 1 
    while P[i] <> M:
        if orientation(P[i-1], P[i], P[i+1])<0:
            del P[i]
            i -= 1
        else:
            i += 1
    return P
p=[[random(),random()] for i in range(1000)]
C=convexHullGraham(p)
line(C+[C[0]]) + point(p)
#Quickhull 
#Tomamos los puntos más a la izquierda y el más a la derecha (ambos podemos
#confirmar que forman parte de la envolvente). Con recursividad y a partir de estos dos, cogiendo
#puntos de P por encima y por debajo encontramos el cierre convexo.

def dome(p,A,B):
    def comp(a):
        return sarea(B,A,a)
    if len(p)==0:
        return [B]
    M=max(p,key=comp)
    
    L1,L2=[],[]
    for i in p:
        if sarea(A,i,M)>0:
            L1.append(i)
        if sarea(M,i,B)>0:
            L2.append(i)
            
    return dome(L1, A, M) + dome(L2,M,B)+[B]
    
def quickhull(p):
    A,B=max(p), min(p)
    return dome(p,A,B)+dome(p,B,A)

p=[[random(),random()] for i in range(5000)]
C=quickhull(p)
line(C+[C[0]]) + point(p)

from time import time
tiempos=[]

for i in range(1,100):
    p=[[random(),random()] for a in range(i*100)]
  
    t0=time()
    C=quickhull(p)
    t1=time()
    tiempos.append([a,t1-t0])

tiempos1=[]
for i in range(1,100):
    p=[[random(),random()] for a in range(i*100)]
  
    t0=time()
    C=convexHullGraham(p)
    t1=time()
    tiempos1.append([a,t1-t0])   


line(tiempos,legend_label='quickull',color="red")+line(tiempos1,legend_label='graham')

 

Práctica 4: ESTRUCTURA DCEL


 

Estructura de datos DCEL

doubly connected edge list (Lista de aristas doblemente enlazadas)

Es una estructura de datos que sirve para representar particiones planas.

En el caso que vamos a usar todas las caras son acotadas y simplemente conexas menos una cuyo borde tiene como borde una curva de Jordan.

Cada cara tiene como borde un polígono determinado por aristas. Todos los bordes estan orientados positivamente excepto el de la cara no acotada.

Cada arista forma parte del borde de una sola cara y posee una arista gemela que forma parte del borde de otra cara.

 

Un DCEL D=[V, E, F] es una lista compuesta por tres listas: la de vértices, la de aristas y la de caras.

Los elementos de la lista de vértices V son pares [[x,y],k] formados por las coordenadas [x,y] de los vértices de la partición y el índice k de una de las aristas de las que el vértice es origen.

Los elementos de la lista de aristas E son quíntuplas de índices enteros [origen, gemela, anterior, posterior, cara]

origen es el índice del vértice origen de la arista

gemela es el índice de la arista gemela

anterior y posterior son los índices de las aristas anterior y posterior a la arista en la cara a la que pertenece

cara es el índice de la cara a la que pertenece la arista

La lista de caras contiene, por cada cara, el índice de una de las aristas de su borde.

Por conveniencia, la cara no acotada será siemmpre la cara 0.

# funcion que crea un DCEL, para un poligono
def dcel(P):
    n=len(P)
    V=[[P[i],i] for i in range(len(P))]
    e=[[i,n+i,(i-1)%n,(i+1)%n,1]for i in range(n)]+[[(i+1)%n,i,n+(i+1)%n,n+(i-1)%n,0]for i in range(n)]
    f=[n,0]
    return [V,e,f]

# funciones para referirse a los elementos asociados a un elemento del DCEL

# indice del origen de una arista e
def origin(e,D):
    return D[1][e][0]

# coordenadas del origen de la arista e
def originCoords(e,D):
    return D[0][origin(e,D)][0]

# arista gemela de la arista e    
def twin(e,D):
    return D[1][e][1]

# arista previa de la arista e    
def prev(e,D):
    return D[1][e][2]

# arista siguiente de la arista e    
def next(e,D):
    return D[1][e][3] 

# indice de la cara de cuyo borde forma parte la arista e    
def face(e,D):
    return D[1][e][4]               

# indice de una de las aristas del borde de la cara c
def edge(c,D):
    return D[2][c]
        
# funcion para dibujar las aristas de un DCEL

def plotDCEL(D):
    return sum(line([originCoords(i,D),originCoords(twin(i,D),D)],aspect_ratio=1) for i in range(len(D[1])))  
    
# funcion para colorear una cara de un DCEL

def plotFace(c,D,col):

    f=D[2][c]
    C=[f]
    f=next(f,D)
    while f <> C[0]:
        C.append(f)
        f=next(f,D)
    
    P=[originCoords(j,D) for j in C]
    return polygon(P,color=col, alpha=.5)     

# funcion para colorear las caras de un DCEL
    
def colorDCEL(D):
    return sum(plotFace(i,D,(random(),random(),random())) for i in range(1,len(D[2])))
# funcion para dividir una cara del DCEL D por una diagonal
# e1 y e2 son las aristas cuyos orígenes son los extremos de la diagonal que divide la cara

def splitFace(e1,e2,D):
    
    # si no son aristas de la misma cara o si son adyacentes sus origenes no definen una diagonal
    if face(e1,D) <> face(e2,D) or origin(e2,D) == origin(twin(e1,D),D) or origin(e1,D) == origin(twin(e2,D),D):
        print "no diagonal"
        return
    
    nv, ne, nf = len(D[0]), len(D[1]), len(D[2])
    preve1 = prev(e1,D)
    preve2 = prev(e2,D)
    k=face(e1,D)
    
    # añadimos las aristas nuevas
    D[1].append([origin(e1,D),ne+1,preve1,e2,k])
    D[1].append([origin(e2,D),ne,preve2,e1,nf])
    
    # modificamos aristas afectadas
    D[1][preve1][3]=ne
    D[1][e1][2]=ne+1
    D[1][preve2][3]=ne+1
    D[1][e2][2]=ne
    i=e1
    while i<>ne+1:
        D[1][i][4]=nf
        i=next(i,D)
    
    #modificamos la cara afectada
    D[2][k]=ne
    
    # añadimos la nueva cara
    D[2].append(ne+1)
# prueba
p=[[random(),random()] for i in range(50)]
P=polygonization(p)
D=dcel(P)
splitFace(0,len(p)-2,D)
plotDCEL(D)+colorDCEL(D)

Ejercicios:

  1. Genera un ejemplo con la prueba anterior y revísalo para comprobar si se trata de una partición correcta delpolígono y si el DCEL es correcto.
  2. Emplea la función splitFace para terminar de triangular el polígono del ejemplo anterior.
  3. Define una función faceEdges(c,D) que obtenga la lista de los indices de las aristas del borde de la cara c del DCEL D (en el orden en que aparecen al recorrer el borde de c)
  4. Define una función faceVertices(c,D) que obtenga la lista de los indices de los vertices del borde de la cara c del DCEL D (en el orden en que aparecen al recorrer el borde de c)
  5. Define una función faceVerticesCoords(c,D) que obtenga el poligono borde de la cara c del DCEL D (con la orientación de c)
  6. Define una función faceNeighbors(c,D) que obtenga la lista de los indices de las caras vecinas de la cara c del DCEL D (en el orden en que aparecen al recorrer el borde de c)
  7. Define una función vertexEdges(v,D) que obtenga la lista de los indices de las aristas del DCEL D cuyo origen es v (en orden angular positivo)
  8. Define una función vertexFan(v,D) que obtenga la lista de los indices de las caras del DCEL D que contienen al vértice v en su borde (en orden angular positivo)
  9. Define una función convexHullDCEL(D) que triangule el exterior del borde de la cara externa iterando la función splitFace
  10. Define una función triangulation(p) que obtenga el DCEL de una triangulación de la lista de puntos p.
# prueba
p=[[random(),random()] for i in range(20)]
P=polygonization(p)
D=dcel(P)
splitFace(0,len(p)-2,D)
V= [D[0][j][0] for j in range (len(D[0]))]
plotDCEL(D)+colorDCEL(D)+point(V,color="red",size=30)

V= [D[0][j][0] for j in range (len(D[0]))]
plotDCEL(D)+colorDCEL(D)+point(V,color="red",size=30)

Define una función faceEdges(c,D) que obtenga la lista de los indices de las aristas del borde de la cara c del DCEL D (en el orden en que aparecen al recorrer el borde de c)

def faceEdges(c,D):
    e=edge(c,D)
    C=[e]
    e=next(e,D)
    while e <> C[0]:
        C.append(e)
        e=next(e,D)
    return C

Define una función faceVertices(c,D) que obtenga la lista de los indices de los vertices del borde de la cara c del DCEL D (en el orden en que aparecen al recorrer el borde de c)

def faceVertices(c,D):
    A = faceEdges(c,D)
    V = [D[1][A[i]][0] for i in range(len(A))]
    return V

Define una función faceVerticesCoords(c,D) que obtenga el poligono borde de la cara c del DCEL D (con la orientación de c)

def faceVerticesCoords(c,D):
    A = faceEdges(c,D)
    coords = [originCoords(A[i],D) for i in range(len(A))]
    return A

Define una función faceNeighbours(c,D) que obtenga la lista de los indices de las caras vecinas de la cara c del DCEL D (en el orden en que aparecen al recorrer el borde de c)

def faceNeighbours(c,D):
    A = faceEdges(c,D)
    coords = [twin(A[i],D) for i in range(len(A))]
    return A

Define una función vertexEdges(v,D) que obtenga la lista de los indices de las aristas del DCEL D cuyo origen es v (en orden angular positivo)

def vertexEdges(v,D):
    e0=D[0][v][1]
    E=[e0]
    e=twin(prev(e0,D),D)
    while e<>e0:
        E.append(e)
        e=twin(prev(e,D),D)
    return E

Define una función vertexFan(v,D) que obtenga la lista de los indices de las caras del DCEL D que contienen al vértice v en su borde (en orden angular positivo)

def vertexFan(v,D):
    A=vertexEdges(v,D)
    C=[face(A[i],D) for i in range(len(A))]
    return C

# CONVEX HULL
# esta funcion no sirve en general
# pero si para casos en que la cara externa es un poligono estrellado o monotono
def DCELconvexHull(D): #scan de Graham
    E=faceEdges(0,D)
    V=[[originCoords(i,D),i] for i in E]
    m=min(V)[1]
    e=next(m,D)
    while origin(e,D)<>origin(m,D):
        if sarea(originCoords(prev(e,D),D),originCoords(e,D),originCoords(next(e,D),D))>0:
            e1=prev(e,D)
            e2=next(e,D)
            splitFace(e1,e2,D)
            if origin(prev(e2,D),D)<>origin(m,D):
                e=prev(e2,D)
            else:
                e=e2
        else:
            e=next(e,D)
def triangulation(p):
    n=len(p)
    P=[[i[1],i]for i in p]
    M=min(P)[0]  #punto de menor ordenada
    P=angularSort(p,M)
    D=dcel(P)
    for i in range(3,n):  #triangulacion de la cara acotada
        splitFace(n-i,n-1,D)
    DCELconvexHull(D)  #scan de Graham
    return D         

   
def triangulation(p):
    n=len(p)
    P=[[i[1],i]for i in p]
    M=min(P)[1]  #punto de menor ordenada
    P=angularSort(p,M)
    D=dcel(P)
    DCELconvexHull(D)  #scan de Graham
    for i in range(3,n):  #triangulacion de la cara acotada
        splitFace(n-i,n-1,D)
    
    return D

 

Triangulación de Delaunay

Ejercicios:

  1. Define una función flipable(a,D) que analice si laarista a del DCEL D es flipable
  2. Define una función legal(a,D) que analice si la arista a del DCEL D es legal
  3. Define una función legalize(T) que legalice las aristas de la triangulación T mediante flips
  4. Define la función delone(p) que calcule la triangulación de Delaunay de una lista de puntos p

Puedes emplear la siguiente función para flipar una arista flipable:


def flip(a,D):
    ga=D[1][a][1]
    ca=D[1][a][4]
    aa=D[1][a][2]
    pa=D[1][a][3]
    cb=D[1][ga][4]
    ab=D[1][ga][2]
    pb=D[1][ga][3]
    D[1][a]=[D[1][aa][0],ga,pa,ab,ca]
    D[1][ga]=[D[1][ab][0],a,pb,aa,cb]
    D[1][pa][2]=ab
    D[1][pa][3]=a
    D[1][aa][2]=ga
    D[1][aa][3]=pb
    D[1][aa][4]=cb
    D[1][pb][2]=aa
    D[1][pb][3]=ga
    D[1][ab][2]=a
    D[1][ab][3]=pa
    D[1][ab][4]=ca
    return

def flip(a,D):
    oa=D[1][a][0]
    ga=D[1][a][1]
    ca=D[1][a][4]
    aa=D[1][a][2]
    pa=D[1][a][3]
    cb=D[1][ga][4]
    ab=D[1][ga][2]
    pb=D[1][ga][3]
    oga=D[1][ga][0]
    D[1][a]=[D[1][aa][0],ga,pa,ab,ca]
    D[1][ga]=[D[1][ab][0],a,pb,aa,cb]
    D[1][pa][2]=ab
    D[1][pa][3]=a
    D[1][aa][2]=ga
    D[1][aa][3]=pb
    D[1][aa][4]=cb
    D[1][pb][2]=aa
    D[1][pb][3]=ga
    D[1][ab][2]=a
    D[1][ab][3]=pa
    D[1][ab][4]=ca
    D[2][ca]=a
    D[2][cb]=ga
    D[0][oa][1]=pb
    D[0][oga][1]=pa

    return 
#flipable(a,D)
def flipable(a,D):
    if face(a,D)==0 or face(twin(a,D),D)==0:
        return false
    A,B,C,D = originCoords(a,D),originCoords(prev(twin(a,D),D),D),originCoords(twin(a,D),D),originCoords(prev(a,D),D)
    if sarea(A,B,C)>0 and sarea(B,C,D)>0 and sarea(C,D,A)>0 and sarea(D,A,B)>0:
        return True
    return false
################################
# contador de aristas illegal
################################
def illegal(D):
    j=0
    for i in range(len(D[1])):
        if legal(i,D)==-1:
            j=j+1
    return j/2
################################
# contador de aristas flipables
################################
def flipables(D):
    j=0
    for i in range(len(D[1])):
        if flipable(i,D)==1:
            j=j+1
    return j/2
#legal(a,D)
def legal(a,D):
    if flipable(a,D)==-1:
        return true
    A,B,C,D = originCoords(a,D),originCoords(prev(twin(a,D),D),D),originCoords(twin(a,D),D),originCoords(prev(a,D),D)
    if inCircle(A,C,D,B)==false:
        return True
    else:
        return False

#legalize(T)
def legalize(t):
    T=deepcopy(t)
    c=1
    vueltas=0
    flips=0
    while c>0:
        vueltas=vueltas+1
        c=0
        for i in range(len(T[1])):
            if legal(i,T)==-1:
                flip(i,T)
                flips=flips+1
                c=c+1
    print "num vueltas: ",vueltas,"num flips: ",flips
    return T

def delone(p):
    k=triangulation(p)
    k=legalize(k)
    return k
# Prueba

p=[[random(),random()]for i in range(20)]
T=triangulation(p)

DT=delone(p)

plotDCEL(DT)+colorDCEL(DT)
num vueltas: 1 num flips: 0





 

Diagrama de Voronoi

El diagrama de Voronoi de n puntos del plano es la descomposición del plano en n regiones cada una de las cuales es el lugar geométrico de los puntos del plano más próximos a uno de los puntos dados que al resto.

Cada una de esas regiones se llama región de Voronoi del punto correspondiente.

Demuestra las siguientes propiedades:

  1. El diagrama de Voronoi de n puntos tiene n regiones no vacías.
  2. Las regiones de Voronoi son regiones poligonales convexas.
  3. La región de Voronoi de un punto que sea vértice de la envolvente convexa es no acotada.
  4. La región de Voronoi de un punto que no sea vértice de la envolvente convexa es acotada.
  5. La intersección de las fronteras de dos regiones de Voronoi es:
    1. Vacía
    2. o un punto
    3. o un segmento
    4. o una semirrecta
    5. o una recta (en el caso de solo dos puntos o de más de dos todos  alineados)
  6. La frontera común de las regiones de Voronoi de dos puntos pip_i y pjp_j está contenida en la mediatriz de pip_i y pjp_j.
  7. Los vértices del diagrama de Voronoi equidistan de tres o más puntos generadores del diagrama.
  8. El grafo dual del diagrama de Voronoi de n puntos es el que tiene como vértices los n puntos y dos de ellos están conectados por una arista rectilínea si, y solo si, las regiones de los dos puntos comparten frontera con más de un punto. Si los puntos están en posición general (no degenerada), se verifica que el grafo dual del diagrama de Voronoi es la triangulación de Delaunay.

Vamos a obtener el diagrama de Voronoi a partir de la triangulación de Delaunay.

Ejercicios:

  1. Define una función voronoiRegion(v,D,R) que calcule la región de Voronoi del vértice v, recortada por el rectángulo R, de la triangulación de Delaunay D (definida por un DCEL)
  2. Define una función Voronoi(p,R) que calcule el diagrama de Voronoi de un conjunto de puntos p recortado por el rectángulo R.
def infiniteVertex(e,D):
    A = originCoords(e,D)
    B = originCoords(next(e,D),D)
    v = vector(B)-vector(A)
    v0 = vector([-v[1],v[0]])
    v0 *= 100/v0.norm()
    M = midPoint(A,B)
    C = vector(M)+v0
    return list(C)

def voronoiRegion(i, D):
    F = vertexFan(i,D)
    R = []
    for i in len(F):
        if face(i,D) != 0:
            A,B,C = originCoords(prev(i,D),D),originCoords(i,D),originCoords(next(i,D),D)
            R.append(circumcenter(A,B,C))
        else:
            R.append(infiniteVertex(i,D))
            R.append(infiniteVertex(prev(i,D),D))
    return R

def voronoi(D): #2n
    V=[]
    for i in range(len(D[0])):
        r = voronoiRegion(i,D)
        V.append(r)
    return V
p=[[random(),random()]for i in range(10)]
DT=delone(p)
V = voronoi(DT)

print(V)


num vueltas: 1 num flips: 0
Error in lines 3-3 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1013, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "", line 4, in voronoi File "", line 4, in voronoiRegion TypeError: 'int' object is not iterable

def dentroCirculoGeneradoPorP(p,r,a):
    if(dist(a,p)<=r): return 1
    else: return 0
def voronoiRegion(i, D):
    punto=D[0][i][0]
    listavertices=[]
    for c in len(D[0]):
        listavertices.append(D[a][0])
    listapuntos=[]
    while(len(listapuntos))<3:
        for a in Combinations(listavertices,3):
            if(circumcenter(a[0],a[1],a[2]))==punto:
                listapuntos.append(a[0])
                listapuntos.append(a[0])
                listapuntos.append(a[0])
                radio=dist(punto,a[0])
                for c in range(len(listavertices)):
                    if(listavertices[c]!=a[0])and(listavertices[c]!=a[0])and(listavertices[c]!=a[0])and(dentroCirculoGeneradoPorP(punto,radio,listavertices[c])==1):listapuntos.append(listavertices[c])
                        
        
    for a in range(len(lista)):
        
    F = vertexEdges(i,D)
    R = []
    for i in F:
        if face(i,D) != 0:
            A,B,C = originCoords(prev(i,D),D),originCoords(i,D),originCoords(next(i,D),D)
            R.append(circumcenter(A,B,C))
        else:
            R.append(infiniteVertex(i,D))
            R.append(infiniteVertex(prev(i,D),D))
    return R