Transformée de distances - Python

Bonjour.

La distance dans le plan est classique (distance euclidienne) et si tu as les coordonnées des points de S1 et S2, un algorithme élémentaire répond à la question (test point par point). Comme cet algorithme est polynomial en les données, il peut convenir si tu n'as pas trop de points.

Cordialement.

Je travaille ici sur les indices des tableaux/matrices des points des différentes surfaces; en Français cela se traduit par :

• point $P_1$ de $S_1$ pointe le point $P'_1$ de $S_2$
• point $P_1$ de $S_1$ pointe le point $P'_2$ de $S_2$
• ...
• puis
• point $P_2$ de $S_1$ pointe le point $P'_1$ de $S_2$
• et ainsi de suite

Pour créer les blocs d'indices, j'utilise le produit de Kronecker (np.kron); le code ci-dessous montre qu'avec 2 points sur $S_1$ et 5 points sur $S_2$, on obtient bien 2 vecteurs indice 1 et indice 2 de même taille (ce sont les différentes combinaisons possibles).

import numpy as np

n1 = 2  # nombre de points de la surface 1
n2 = 5  # nombre de points de la surface 2

i1 = np.arange(0,n1, dtype = np.int) # [0 1 2 ... n1]
i2 = np.arange(0,n2, dtype = np.int) # [0 1 2 ... n2]
 
j1 = np.ones(n2, dtype = np.int) # vecteurs unitaires qui servent au produit de Kronecker
j2 = np.ones(n1, dtype = np.int)

indice1 = np.kron(i1,j1) # vecteur d'indices pour les points de S1
indice2 = np.kron(j2,i2) # idem pour S2
print("indice 1 = {}".format(indice1))
print("indice 2 = {}".format(indice2))

Comme montré dans ce fil, il ne reste plus qu'à utiliser ces indices pour calculer un vecteur "distance" et d'en rechercher ensuite les $\min$ et $\max$.

Aucune boucle n'a été utilisée B-)
Paul

nice_boycott a écrit:

Je vais de ce pas essayer de mettre en place cela et voir si c'est bien plus efficace que des > boucles imbriquées

Assurément (par expérience) et donc si tu trouves mieux sans utiliser des artifices comme Numba par exemple, ça m’intéresse

il ne s'applique pas qu'à un problème 1D (voir mon exemple avec 2 surfaces en 3D); ne pas oublier qu'on travaille sur des tableaux/matrices de type $A_{i,j}$ où $i$,$j$ sont tes indices: par exemple $A[1,1]$ en Python pointe sur $A_{2,2}$.

En vectorisation, on ne travaille pas sur 1 indice, mais un "bloc" d'indices qu'on appelle "vecteur" !


Je ne saurais trop te conseiller d'aller voir le fil précédemment cité où tu trouveras des exemples

Comme il parle de la distance de Hausdorff…

nice_boycott si, mais c’est vrai que son message n’est pas très clair.

Alors ce n’est pas la distance de Hausdorff mais le max des max.

Déjà, dans un but d’économie de moyens, la racine carrée est inutile puisque tu veux juste comparer les distances, donc leurs carrés.
Pour le reste, je ne connais pas assez numpy.

En supposant que Points_S1 et Points_S2 sont les tableaux avec les coordonnées des points,alors ce bout de code devrait fonctionner; attention, dans le cas présent, les $X$ sont dans la $2^{ième}$ colonne (voir [:,1]) et les $Y$ dans la $3^{ième}$ (voir [:,2]).

Les termes "indice" regroupent en fait les différentes combinaisons (on fonctionne par "ligne" comme on le voit sur l'image précédente) et $k$ balaie les vecteurs indice1 et indice2 de la première à la dernière ligne.

J'espère que ça te laisse entrevoir les possibilités de la vectorisation (j'utilise sensiblement toujours le même raisonnement, et je suis loin d'en maîtriser toutes les subtilités)

Paul

n1 = np.shape(Points_S2)[0]
n2 = np.shape(Points_S1)[0]


i1 = np.arange(0,n1, dtype = np.int)
i2 = np.arange(0,n2, dtype = np.int)
 
j1 = np.ones(n2, dtype = np.int)
j2 = np.ones(n1, dtype = np.int)

indice1 = np.kron(i1,j1)
indice2 = np.kron(j2,i2)

NbreCombinaisons = np.shape(indice1)[0]
k = np.arange(0,NbreCombinaisons, dtype = np.int)
# 
Normes = np.sqrt( (Points_S1[indice2[k], 1] - Points_S2[indice1[k], 1])**2 + \
                  (Points_S1[indice2[k], 2] - Points_S2[indice1[k], 2])**2  )

MinDistance = np.min(Normes)
MaxDistance = np.max(Normes)

comment peux-tu avoir des Nan avec des termes au carré et aucune division? tu es sûr de tes données? tu pointes bien vers les coordonnées?

Le code ci-dessous donne plus de 42 million de combinaisons sans NaN

import numpy as np
import time

n1 = 20_099
n2 = 2_103

Points_S1 = np.random.random( (n1,3))
Points_S2 = np.random.random( (n2,3))


t0 = time.time()
n1 = np.shape(Points_S1)[0]
n2 = np.shape(Points_S2)[0]


i1 = np.arange(0,n1, dtype = np.int)
i2 = np.arange(0,n2, dtype = np.int)
 
j1 = np.ones(n2, dtype = np.int)
j2 = np.ones(n1, dtype = np.int)

indice1 = np.kron(i1,j1)
indice2 = np.kron(j2,i2)

NbreCombinaisons = np.shape(indice1)[0]
k = np.arange(0,NbreCombinaisons, dtype = np.int)
 
Normes = np.sqrt( (Points_S1[indice1[k], 1] - Points_S2[indice2[k], 1])**2 + \
                  (Points_S1[indice1[k], 2] - Points_S2[indice2[k], 2])**2  )
    
MinDistance = np.min(Normes)
MaxDistance = np.max(Normes)

print("Min = {}".format(MinDistance))
print("Max = {}".format(MaxDistance))
print("Nombre de combinaisons = {}".format(np.shape(k)[0]))

t1 = time.time()
T1 = t1 - t0
print("Durée : {}".format(T1))
del j1; del j2; del k; del Normes

Min = 8.506638984615577e-05
Max = 1.3916897858871962
Nombre de combinaisons = 42268197
Durée : 8.254205465316772

voir dernier apport; quant au code avec A, il fonctionne chez moi (Anaconda et Spyder)

je dois trouver le point de $S_1$ le plus éloigné de tous les points de $S_2$ et vice versa

je cherche un moyen de trouver la distance la plus grande qui soit entre deux points issus respectivement de deux nuages de points dans le plan

Ce n'est pas pareil. Première chose : spécifier clairement le problème.
Pour ne pas être uniquement négatif, une piste qu'il vaudrait peut-être le coup d'explorer pour éviter le calcul de $|S_1|\times |S_2|$ distances.
Ça ne coûte pas très cher de subdiviser la région du plan occupée par les points en petits carrés indexés par des couples d'entiers $(i,j)$ et de ranger les points de $S_1$ et ceux de $S_2$ dans ces carrés. C'est en gros en $n\log(n)$ si les deux paquets de points sont de taille $n$.
Supposons que c'est bien $\max_{x\in S_1,y\in S_2} d(x,y)$ que l'on cherche, Si on s'intéresse alors aux distances d'un $x_0\in S_1$ aux $y\in S_2$ il est inutile de s'embêter à calculer les $d(x_0,y)$ pour $y$ dans un carré proche de celui de $x_0$. La proximité entre le carré $(i_0,j_0)$ et le carré $(i,j)$ peut s'évaluer par $|i-i_0|+|j-j_0|$. Si on a déjà parcouru un bout de la liste des points de $S_1$ et qu'on a un record partiel de distance quand on arrive à $x_0$, je me dis que ça évite un gros paquet de calculs (surtout si on a commencé d'éplucher la liste des points de $S_1$ par ceux qui sont le plus en bas à gauche et la liste des points de $S_2$ par ceux qui sont le plus en haut à droite).
Ça demande à être beaucoup retravaillé, mais ça me semble une piste valable.

C'est simple : on trie la liste des abscisses de l'ensemble de points, ce qui permet de donner le n° de colonne dans la grille de carrés, et on trie la liste des ordonnées, ce qui permet de donner le n° de ligne dans la grille de carrés. Tu connais la complexité des bons algorithmes de tri.

ConvexHull (enveloppe convexe) et pas ConvexHill (colline convexe). ;-)

À part ça, oui : soit $C$ est l'enveloppe convexe du nuage de points $S$ et $x$ un point quelconque. Si $y_0\in S$ réalise le maximum des $d(x,y)$ pour $y\in S$, alors $y_0$ est un sommet de $C$. En effet, $S$ est entièrement contenu dans le disque fermé de centre $x$ et de rayon $d(x,y_0)$, donc $C$ aussi puisque le disque est convexe ; et comme le disque est strictement convexe, $y_0$ est forcément un sommet de $C$.

Le calcul de l'enveloppe convexe de $n$ points dans le plan se fait en $O(n\log(n))$. Est-ce qu'on économise beaucoup en faisant ça ? Ça dépend de la forme du nuage de points. Si par exemple les points sont répartis sur un cercle, pas tellement.