voronoi.py

import numpy as np
import math
import random as rng
import matplotlib.pyplot as plt

# np.array sans tous les np.array
def mat(*coord):
    try:
        iter(coord[0])
        return np.array([mat(*c) for c in coord])
    except TypeError:
        return np.array(coord)

def produit_scalaire(u, v):
    return sum(u*v)

# Distance entre A et B
def dist(A, B):
    return np.sqrt(dist_sq(A, B))
    
# Distance entre A et B au carré
def dist_sq(A, B):
    u = B-A
    return produit_scalaire(u, u)

# voisin "classique" appartenant au nuage
class Voisin:
    def __init__(self, point):
        self.point = point
    
    def __eq__(self, other):
        if isinstance(other, Voisin):
            return (self.point == other.point).all()
        return False
    
    def __hash__(self):
        return hash(tuple(self.point))
    
    def __str__(self):
        return "V(" + str(self.point) + ")"
        
    # Revoie toujours le point encapsulé
    def get_point(self, P):
        return self.point
    
    def is_border(self):
        return False
    
# bord voisin dans une fenetre (1,1)
class Bord(Voisin):
    def __init__(self, b):
        self.b = b
    
    def __eq__(self, other):
        if isinstance(other, Bord):
            return self.b == other.b
        return False
    
    def __hash__(self):
        return self.b
    
    def __str__(self):
        return "B" + str(self.b)
        
    # Renvoie le symétrique du point P par rapport à ce bord
    def get_point(self, P):
        x,y = P
        if self.b == 0:
            return mat(-x,y)
        elif self.b == 1:
            return mat(x,-y)
        elif self.b == 2:
            return mat(2-x,y)
        elif self.b == 3:
            return mat(x, 2-y)
        else:
            raise Exception("index de bord non défini")
            
    def is_border(self):
        return True
            

# suppression des doublons successifs
def minify(L):
    i = 0
    while i+1 < len(L):
        if L[i] == L[i+1]:
            del L[i]
        else:
            i += 1
            
# Slice circulaire
def super_slice(L, a, b):
    l = len(L)
    a,b = a%l,b%l
    if a < b:
        return L[a:b]
    else:
        return L[a:] + L[:b]
        
            
# Construction des sommets de la cellule de centre P
def mk_cell(P, vois):
    S = []
    
    AB = vois[-1]-P
    I = (vois[-1]+P)/2
    
    ap,bp = AB
    cp = I[0]*ap + I[1]*bp
    
    for i in range(len(vois)):
        AB = vois[i]-P
        I = (vois[i]+P)/2
        
        a,b = AB
        c = I[0]*a + I[1]*b
        
        S.append(mat( (c*bp-cp*b)/(a*bp-ap*b) , (a*cp-c*ap)/(a*bp-ap*b) ))
        
        ap,bp,cp = a,b,c
        
    return S
     
# Construction de nuage aléatoire de taille n   
def random_N(n):
    return [mat(rng.random(), rng.random()) for i in range(n)]

# Nuage de points
#N = random_N(int(input("Taille du nuage")))
B0 = Bord(0)
B1 = Bord(1)
B2 = Bord(2)
B3 = Bord(3)

# Constructioin du diagramme de Voronoï du nuage N (liste de np.array)
# résultat sous forme de dict() tuple -> liste de Voisins
def voronoi(N):
    # Diagramme initial
    Ntemp = [N[0]]
    vor = dict()
    vor[tuple(Ntemp[0])] = [B0,B1,B2,B3]
    
    #draw(vor)
    
    for j in range(1, len(N)):
        Pj = N[j]
        
        # Détermination de Pi
        Pi = Ntemp[0]
        dPi = dist_sq(Pi,Pj)
        for P in Ntemp:
            dP = dist_sq(P, Pj)
            if dP < dPi:
                Pi,dPi = P,dP
            
        vois_Pj = []
        
        print("--- Ajout du point P", j, "=", Pj, "---")
        
        study(Pj, vois_Pj, Pi, vor, set())
        minify(vois_Pj)
        
        print("Ses voisins:", vois_Pj)
        
        Ntemp.append(Pj)
        vor[tuple(Pj)] = vois_Pj
        
        draw(vor)
        
    return vor

# Etude d'une cellule de centre Pk en vue de l'ajout de Pj
def study(Pj, vois_Pj, Pk, vor, studied):
    if tuple(Pk) in studied:
        return
    
    print("Etude de", Pk)
    
    I = (Pj + Pk)/2
            
    vois_Pk = vor[tuple(Pk)]
    # Sommets de la cellule de centre Pk
    Sk = mk_cell(Pk, [v.get_point(Pk) for v in vois_Pk])
            
    # Détermination des index de voisins entrant/sortant
    v_entrant = None
    v_sortant = None
            
    for k in range(len(vois_Pk)):
        A,B = Sk[k-1],Sk[k]
                
        #print("a,b", A,B)
                
        PjPk_IA = produit_scalaire(Pk-Pj, A-I)
        PjPk_IB = produit_scalaire(Pk-Pj, B-I)
        
        if PjPk_IA < 0 and PjPk_IB > 0:
            v_sortant = k-1
        elif PjPk_IA > 0 and PjPk_IB < 0:
            v_entrant = k-1
    
    if v_entrant == None or v_sortant == None:
        return
        #raise Exception("vois e/s", v_entrant, v_sortant)
        
    
        
    vois_Pj.append(Voisin(Pk))
        
    vor[tuple(Pk)] = super_slice(vois_Pk, v_sortant, v_entrant+1) + [Voisin(Pj)]
    
    V = super_slice(vois_Pk, v_entrant, v_sortant+1)
    

    
    studied.add(tuple(Pk))
    
    for V in super_slice(vois_Pk, v_entrant, v_sortant+1):
        print("maybe", V)
        if V.is_border():
            if not V in studied:
                print("ajout du bord", V)
                vois_Pj.append(V)
                studied.add(V)
        else:
            study(Pj, vois_Pj, V.get_point(Pj), vor, studied)

# Affichage d'un diagramme de Voronoï en 2D
def draw(vor):
    plt.figure()
    
    for P,vois in vor.items():
        plt.scatter(P[0], P[1], c='red')
        
        S = mk_cell(P, [v.get_point(P) for v in vois])
        for i in range(len(S)):
            A,B = S[i-1],S[i]
            
            
            plt.scatter(A[0], A[1], c='green')
            plt.scatter(B[0], B[1], c='green')
            
            plt.plot([A[0], B[0]], [A[1],B[1]], color='blue')
            
    plt.show()