Point le plus proche d'un ensemble de droite

[Benjamin]

Salut,

J'aimerai avoir vos lumières sur un problème dont je ne sais pas si il existe une solution matricielle/analytique autre que de faire une optimisation numérique.

J'ai un ensemble de droites dans l'espace (3D), connu par un point et un vecteur directeur. Je cherche le point qui minimise la distance globale à l'ensemble de ces droites (RMS).

Je sais que la distance d'un point à une droite est ||BA ^ u|| (j'ai u unitaire) mais j'ai pas réussi à créer une système matriciel qu'on pourrait résoudre à partir de là.
J'ai tenté de gruger en disant qu'une droite est également l'intersection de deux plans ce qui me donne un système de deux équations pour chaque droite et là j'ai plus ou moins trouvé un système à résoudre mais ça ne me donne pas toujours le bon point (si tant est que ça me donne le bon point quand celui que je trouve est cohérent et pas très loin de la vérité...).

Bref, merci pour vos éclairages J'ai un peu cherché sur le net et pas trouver. Pourtant, ça ne me semblait pas si novateur comme demande ^^

[GaBuZoMeu]

Bonjour,

Qu'appelles-tu "distance globale" ??? La somme des distances à chacune des droites ? Ou alors quoi ?

[Benjamin]

Minimisation du RMS de l'ensemble des distances de ce point à trouver et de chacune des droites.

Si vous voulez plus de contexte, c'est pour trouver le point de convergence d'un système optique qui n'a pas un stigmatisme rigoureux. L'ensemble des rayons converge vers une zone et non un point.

Une autre façon qui me conviendrait (même si c'est pas exactement la même chose), c'est de trouver le centre de la sphère de rayon minimal qui est traversée par toutes les droites (c'est-à-dire trouver le point qui minimise la distance à la droite la plus éloignée).

[GaBuZoMeu]

Peux tu décoder RMS ? Revue Militaire Suisse ? Reims Management School ? Richard M. Stallman ? République des Moluques du Sud ? Radio Morbihan Sud ? Root Mean Square ?

Tu veux minimiser la somme des carrés des distances aux droites ? Minimiser la plus grande distance aux droites ?

[Sake]

Peux tu décoder RMS ? Revue Militaire Suisse ? Reims Management School ? Richard M. Stallman ? République des Moluques du Sud ? Radio Morbihan Sud ? Root Mean Square ?

Tu veux minimiser la somme des carrés des distances aux droites ? Minimiser la plus grande distance aux droites ?

Je pense qu'il veut minimiser la République des Moluques du Sud

[fatal_error]

hi,

si ||AM^u|| donne la distance entre le point M et la droite dont un point est A et et u son vecteur directeur,

si RMS correspond à https://en.wikipedia.org/wiki/Root_mean_square/rmse/ (root mean square)

alors je présume que tu cherches M qui minimise


si je dis pas de bêtises, si tu developpes le produit vectoriel
t'obtiens un vecteur de la forme
(
a1_i y_M + b1_i z_M + c1_i,
a2_i x_M + b2_i z_M + c2_i,
a3_i y_M + b3_i x_M + c3_i,
)

du coup la racine carrée de la norme saute avec ()^2
et il reste juste à développer les carrés respectifs
a1_i y_M + b1_i z_M + c1_i,
a2_i x_M + b2_i z_M + c2_i,
a3_i y_M + b3_i x_M + c3_i,

comme sqrt est strictement croissante tu peux la faire sauter (par rapport à l'optimisation), idem pour 1/n
et il reste juste à trouver M tq la somme des carrés de
a1_i y_M + b1_i z_M + c1_i,
a2_i x_M + b2_i z_M + c2_i,
a3_i y_M + b3_i x_M + c3_i,
est minimale

en développant les carrés, puis en groupant les termes en x_M, ...
en posant x == x_M, y=y_M, z=z_M,
on cherche alors à minimiser
Q(x,y,z) = ax^2 + bxy + cz^2 + dxy + exz + fyz

(cf wiki, forme quadratique)

apres je sais pas, point critiques ou réduction sous forme de carrés..

[Benjamin]

Peux tu décoder RMS ? Revue Militaire Suisse ? Reims Management School ? Richard M. Stallman ? République des Moluques du Sud ? Radio Morbihan Sud ? Root Mean Square ?

Tu veux minimiser la somme des carrés des distances aux droites ? Minimiser la plus grande distance aux droites ?

Root Mean Square dans le cas présent ^^
Au départ, je cherche bien à minimiser la somme des carrés des distances aux droites. Mais je sais que dans mon cas pratique, c'est quasiment la même solution que la minimisation de la plus grande distance aux droites. Donc si mathématiquement, il y en a un qui se résout plus facilement que l'autre, c'est au choix

[GaBuZoMeu]

Minimiser la somme des carrés des distances, pas de problème. Ça revient à résoudre un système linéaire à trois inconnues. Je fais le code en SageMath.

Code: Tout sélectionner

x,y,z=var('x,y,z')
Vars=[x,y,z]

# Calcul du carré de la distance du point [x,y,z]
# à une droite donnée comme le couple d'un point
# et d'un vecteur directeur
def cddroite(droite) :
    point=vector(droite[0])
    dir=vector(droite[1])
    M=vector(Vars)
    crossprod= (M-point).cross_product(dir)
    return (crossprod.dot_product(crossprod))/(dir.dot_product(dir))

# Calcul de la somme des carrés des distances
# du point [x,y,z] aux droites d'une liste de droites
def scddroites(Droites) :
    return add(cddroite(droite) for droite in Droites)

# Minimisation de la somme des carrés des distances
# Retourne le point ou le minimum est atteint et la
# racine carrée de la moyenne des carrés des distances
# depuis ce point
def Minscddroites(Droites) :
    Somme=scddroites(Droites)
    Jac=[diff(Somme,var) for var in Vars]
    Crit=solve(Jac,Vars,solution_dict=True)[0]
    MCD=(Somme.subs([var == Crit[var] for var in Vars]))/len(Droites)
    return Crit, sqrt(MCD)

Un exemple d'exécution :

Code: Tout sélectionner

Minscddroites([[[2,3,4],[-1,2,1]],
               [[0,0,0],[1,1,1]],
               [[5,4,3],[2,1,0]],
               [[1,3,1],[2,0,-1]]])

réponse :

Code: Tout sélectionner

({z: 8879/3635, y: 9026/3635, x: 7187/3635}, 3*sqrt(1073/7270))

Après, vu le contexte, tu souhaiterais peut-être avoir toute un famille de droites dépendant de paramètres et pas seulement un ensemble fini ?

[Benjamin]

Merci à vous deux, je regarde ça.
J'avoue que je ne vois pas trop à quoi correspond ton système d'équation linéaire à 3 inconnus (D'ailleurs, fatal_error obtient une forme quadratique). Les équations linéaires que j'avais obtenues, c'est quand j'avais dit qu'une droite était égale à l'intersection de deux plans. Et je m'étais rendu compte que j'avais pas la bonne solution en remplaçant la distance à la droite par deux distances aux deux plans. Bref, j'ai pas mes notes sur moi là, je regarderai plus tard ^^

Je cherche bien un unique point sinon. J'ai déjà mon ensemble de droites qui correspondent en fait à mes rayons optiques Je veux juste trouver où elles convergent le mieux pour bien placer mon capteur.

[GaBuZoMeu]

La somme des carrés des distances de [x,y,z] à un nombre fini de droites est un polynôme du second degré en x,y,z. Le minimum de ce polynôme sera atteint à son point critique, qui est la solution du système donné par l'annulation des trois dérivées partielles par rapport à x,y,z. C'est un système linéaire de trois équations en x,y,z. Je commente un peu le code de la dernière procédure :

Code: Tout sélectionner

def Minscddroites(Droites) :
#calcul de la somme des carrés des distances
    Somme=scddroites(Droites)
#liste des dérivées partielles / x,y,z
    Jac=[diff(Somme,var) for var in Vars]
#résolution du système donné
#par les dérivées partielles
    Crit=solve(Jac,Vars,solution_dict=True)[0]
    MCD=(Somme.subs([var == Crit[var] for var in Vars]))/len(Droites)
    return Crit, sqrt(MCD)

Ensuite, ma question ne concerne pas l'unicité du point à trouver, mais le fait que tu n'as pas en fait un nombre fini de rayons, mais une famille infinie de rayons. Non, me trompé-je ?

PS. Ce n'est pas une forme quadratique qu'on obtient comme somme des carrés des distances, mais bien un polynôme du second degré (avec un terme constant et des termes du premier degré). Minimiser une forme quadratique positive, ça serait immédiat : le minimum est à l'origine.

[Benjamin]

Ok, merci, c'est plus clair
Non, j'ai un ensemble fini de rayon. Pas besoin d'une infinité. Je fais des tests de convergence ensuite, pour voir quelle densité j'ai besoin.

[GaBuZoMeu]

On pourrait aussi faire des calculs avec une famille complète de rayons (infinie), mais il faudrait pour cela plus de précisions.

[Benjamin]

Merci, ça ira
En fait, j'étais tellement obnubilé par mon envie de trouver un système linéaire Ax=B surdéterminé que j'ai complètement oublié le programme de seconde où on trouve l'extrema d'une fonction avec sa dérivée.........

[GaBuZoMeu]

Mais finalement, c'est bien un système linéaire qu'on résoud !

Et pour une personne imaginaire qui n'aurait pas compris (ça, c'est une habitude) pourquoi c'est linéaire : c'est parce que la dérivée d'un polynôme du second degré est un polynôme du premier degré.

[fatal_error]

hi,

tjs de la forme
Q(x,y,z) = ax^2 + bxy + cz^2 + dxy + exz + fyz + gx+hy+jz+k
(j'ai oublié les monomes et constante en developpant)

pour moi c'est une forme quadratique mais bon... apparemment il faut que le polynome soit homogène de degré 2 (et les monomes gx,..., k sont pas de degré 2). Donc Q est juste un polynome de degré 2 j'imagine...

le système se réécrit
Q(x,y,z) = tXAX - 2tBX + c

et la dérivée par rapport à X donne
Q'(X) = 2AX - 2B

d'où en cherchant le min, AX = B

pour un u donné [a,b,c] et son point A (A_i)
on a [a,b,c]^A_i = [i,j,k] (une valeur qu'on connait)

||u^(M-A)||^2 = ([a,b,c] ^ [x,y,z] - [a,b,c] ^ A_i )^2 =

(bz - cy - i)^2 + (cx - az - j)^2 + (ay - bx - k)^2

en dev:
(c^2+b^2) x^2 + (a^2+c^2)y^2 + (b^2+a^2)z^2
+ (-2ab) xy
+ (-2bc) yz
+ (-2ac) xz
+ 2(kb-jc) x
+ 2(ic-ak) y
+ 2(aj-ib) z
+ i^2+j^2+k^2

on peut donc écrire

et

il reste plus qu'à sommer pour tous les points

Avec l'exemple fourni par gbzm

Code: Tout sélectionner

Minscddroites([[[2,3,4],[-1,2,1]],
               [[0,0,0],[1,1,1]],
               [[5,4,3],[2,1,0]],
               [[1,3,1],[2,0,-1]]])


on obtient

Code: Tout sélectionner

A = [
   1.90000  -0.40000   0.23333
  -0.40000   2.80000  -0.66667
   0.23333  -0.66667   3.30000
]

et

B = [   
3.3333
4.5333
6.8667
]



ci-dessous implem octave, la version "AX=B" et l'optim avec sqp

Code: Tout sélectionner

1;
M = [[[2,3,4],[-1,2,1]],
               [[0,0,0],[1,1,1]],
               [[5,4,3],[2,1,0]],
               [[1,3,1],[2,0,-1]]];
for i = 1:4
    M(i,4:end) = M(i,4:end)/norm(M(i,4:end));
end

function [A,B] = add(v)
    A_i = v(1:3);
    u = v(4:end);
    a = u(1);
    b = u(2);
    c = u(3);
    tmp = cross(u, A_i);
    i = tmp(1);
    j = tmp(2);
    k = tmp(3);
    A = [
        [c^2+b^2, -a*b   , -a*c],
        [-a*b   , a^2+c^2, -b*c],
        [-a*c   , -b*c   , b^2+a^2]
    ];
    B = [
        j*c - k*b; a*k-i*c; i*b-a*j
    ];
endfunction
A = zeros(3,3);
B = zeros(3,1);
for i=1: 4
    [A_i, B_i] = add(M(i,:));
    A += A_i;
    B += B_i;
end
X = inv(A)*B


function obj = phi (M,X)
    obj = 0;
    for i = 1:4
        A = M(i,1:3)';
        u = M(i,4:end)';
        obj += sum(cross(u, X-A).^2);
    end
endfunction

x0 = [-1.8; 1.7; 1.9];

[x, obj, info, iter, nf, lambda] = sqp (x0, @(X) phi(M,X));
x_sqp = x


trouvent bien le même objectif:
X=[
1.9772
2.4831
2.4426
]

edit: lien
https://rextester.com/GXW80205

[Benjamin]

Salut,

Je trouve la même chose dans Matlab, c'est cool
Merci pour votre aide.

Code: Tout sélectionner

% Coordonnées des points
x=[2 0;5 1];
y=[3 0;4 3];
z=[4 0;3 1];
% Vecteur directeur (x,y,z) sur la 3ème dimension
u=zeros(2,2,3);
u(:,:,1)=[-1 1;2 2];
u(:,:,2)=[2 1;1 0];
u(:,:,3)=[1 1;0 -1];
u=u./repmat(sqrt(sum(u.^2,3)),[1 1 3]); % Normalisation

%
Ayz=sum(sum(u(:,:,2).^2+u(:,:,3).^2));
Axy=sum(sum(u(:,:,1).^2+u(:,:,2).^2));
Axz=sum(sum(u(:,:,1).^2+u(:,:,3).^2));
Uxy=sum(sum(u(:,:,1).*u(:,:,2)));
Uxz=sum(sum(u(:,:,1).*u(:,:,3)));
Uyz=sum(sum(u(:,:,2).*u(:,:,3)));
Cx=sum(sum(x.*(u(:,:,2).^2+u(:,:,3).^2)-y.*u(:,:,1).*u(:,:,2)-z.*u(:,:,1).*u(:,:,3)));
Cy=sum(sum(y.*(u(:,:,1).^2+u(:,:,3).^2)-x.*u(:,:,2).*u(:,:,1)-z.*u(:,:,2).*u(:,:,3)));
Cz=sum(sum(z.*(u(:,:,1).^2+u(:,:,2).^2)-x.*u(:,:,3).*u(:,:,1)-y.*u(:,:,3).*u(:,:,2)));
A=[Ayz -Uxy -Uxz ; -Uxy Axz -Uyz ; -Uxz -Uyz Axy]

A =

    1.9000   -0.4000    0.2333
   -0.4000    2.8000   -0.6667
    0.2333   -0.6667    3.3000

B=[Cx;Cy;Cz]

B =

    3.3333
    4.5333
    6.8667

A\B

ans =

    1.9772
    2.4831
    2.4426