Calcul de correspondance entre 2 géométries 3D

[MCphylyss]

Bonjour,
Je me suis inscrit sur ce forum en désespoir de cause.
Mon problème est très simple en théorie, mais je n'arrive pas à trouver le bon algo.

J'ai deux géométries 3D, qui se ressemblent beaucoup. On pourrait même dire que l'une est un clone de l'autre. La seule différence, c'est que ces géométries n'ont pas le même nombre de points.

Je voudrais simplement calculer pour chaque point de la géométrie A à quel triangle il appartient dans la géométrie B. c'est a dire que pour un point de A, je voudrais connaitre les 3 points de B qui entourent ce point.

Ma solution actuelle est de calculer les 3 points de B les plus proches du point de A. ça marche, mais les résultats sont loin d'être parfait alors qu'ils le seraient si je connaissais le bon algo pour faire ça.

Par avance, merci pour votre aide,
En espérant avoir été suffisamment claire.

[fatal_error]

salut,

A supposer que ton point soit bien "contenu" dans ton triangle (et non un peu en hauteur, au quel cas il faudrait le projeter), pour déterminer s'il est dans le triangle ABC,
tu peux prendre les vecteurs AB et AC, et exprimer AP en fonction de AB et AC (où P est ton point)
AP = aAB+bAC
si 0<=a <=1 et 0<=b<=1
alors ton point est dans le triangle

[MCphylyss]

Bonjour, merci pour ta réponse.

Le point n'est pas réellement contenu dans le triangle. Il en est très proche, mais il y a, au moins par endroits des différences de hauteur. Une projection me parait donc nécessaire, mais je ne sais pas comment la mettre en place.

Mettons que p soit le point de la géométrie A, et a, b, c soient trois points de la géométrie B que j'étudie. Pour savoir si a, b, c sont les 3 points les plus susceptibles d'entourer p, j'ai pensé à étudier les vecteurs (pa), (pb) et (pc), mais honnêtement je ne sais pas quoi faire de ces vecteurs.

Pourrais tu m'aider à réaliser une telle projection (si j'ai bien compris, il faudrait projeter le point p dans le triangle abc, et vérifier si p' appartient au triangle), ou même simplement m'orienter sur un cours / tuto / outil mathématique stp?

Merci.

[fatal_error]

si j'ai bien compris, il faudrait projeter le point p dans le triangle abc, et vérifier si p' appartient au triangle

oui

Tu peux par exemple (ya ptet mieux) considérer
- P' comme tu l'as défini.
- (AB') orthogonal à (AB) dans le plan (ABC) (de facto orthogonal à (P'P) aussi)

travaillons avec des vecteurs unitaires:



1) On commence par décomposer en


on déduit alpha


on déduit (proportionnel à v)


et donc beta (ce qui ne sert pas forcément mais permet de vérifier les calculs...)


2) on exprime en fonction de et


De la même manière que pour alpha et beta on déduit a et b
Il n'y a plus qu'à remonter en remplacant par son expression en et




or

D'où on déduit en fonction de et

[MCphylyss]

Merci!

ça fait un moment que je n'ai pas fait de calcul vectoriel et j'ai un peu de mal avec les notations donc ça me prendra un peu de temps pour me remettre dedans, mais grâce à toi, je pense que j'ai tout ce qu'il faut pour me sortir de là!

Mille merci pour cette explication complète et détaillée! ça fait plusieurs jours que je cherche, et jusque là j'avais aucune piste...

[fatal_error]

Code: Tout sélectionner

var scal = function(u,v){
    return v[0]*u[0]+v[1]*u[1]+v[2]*u[2];
}
var vec = function(u,v){
    return [v[0]-u[0],v[1]-u[1],v[2]-u[2]]
}
var norm = function(u){
    return Math.sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
}
var smul = function(a,u){
    return [a*u[0],a*u[1],a*u[2]]
}

/**
 * assumes A,B,C not linear
 * A,B,C,P 3d point (x,y,z)
 * @return {a:number, b:number} such that av_{AB}+bv_{AC} = AP' where P' is P projected orthogonally onto (ABC)
 */
function decomp(A,B,C,P){
    var AB = vec(A,B);
    var AC = vec(A,C);
    var AP = vec(A,P);
    var u = smul(1/norm(AB), AB);
    var alpha = scal(AC,u);
    var ABprime = vec(smul(alpha,u), AC);
    var v = smul(1/norm(ABprime), ABprime)
    //we can check that beta is ok: AC - au - bv = 0
    var beta = scal(AC, v);
    console.log(vec(smul(beta,v),vec(smul(alpha,u), AC)))   //expects [0,0,0]

    var a = scal(AP,u);
    var b = scal(AP,v);
    return {
        a:(a - alpha*b/norm(ABprime)) / norm(AB),
        b:b/norm(ABprime)
    }
}

var A = [1,1,2]
var B = [7,1,3]
var C = [8,4,1]
var P = [21,1.5,5.2]
var res = decomp(A,B,C,P);
//we can check as well that aAB+bAC == aAP' and P'P.AB==0 and P'P.AC==0
var PprimeP = vec(smul(res.b, vec(A,C)), vec(smul(res.a, vec(A,B)), vec(A,P)))
console.log(scal(PprimeP, vec(A,B))) //expects 0
console.log(scal(PprimeP, vec(A,C))) //expects 0


j'ai pas testé visuellement, mais c'est dans l'idée..

[MCphylyss]

Oh alors ça, c'est vraiment sympa! si si, j'insiste! Je pensais juste demander une piste, mais si tu me fais tout le travail, alors là je prends!

Je testerais ça demain, il faut que je convertisse en python et que je l'intègre dans mon convertisseur de morphologies, je te dirais si ça marche. Mais dans tous les cas, tu me sauves la vie! Encore merci!

[fatal_error]

attention si 0<=a <=1 et 0<=b<=1 est insuffisant

Code: Tout sélectionner

     C
     /\ M2
    / /\
   / /  \
  x------>
    M1
  A      B


si on conserve a devant AB et b devant AC alors
si on prend evidemment on est en dehors du triangle.
en notant a = AM1 et b=M1M2 puis thales on a
(1-AM1)/1 = M1M2/1 <=> a+b=1 pour etre sur la "bordure" de BC

idem il faut considérer
0<=a<=1
0<=b<=1
a+b<=1
pour avoir le point à l'intérieur du triangle

[MCphylyss]

Merci,

J'ai testé le code que tu as posté hier, je pense qu'il y a une erreur, mais je ne trouve pas ou (et c'est peut être moi, quand je l'ai convertit en python). Quoique il en soit, tu m'as fait le plus gros du travail. J'ai pu trouver le reste en suivant quelques tutos, et ça semble fonctionner. Je n'en aurais la certitude que quand j'aurais intégré le code dans mon convertisseur, ça me prendra un peu de temps. Mais ça marche avec tous les exemples simples que j'ai essayé.

Si quelqu'un est intéressé, voici le code python :

Code: Tout sélectionner

import math

def vec(a,b) :return [b[0]-a[0], b[1]-a[1], b[2]-a[2]]
def norm(a) :return math.sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2])
def cross(u, v) :return [u[1]*v[2]-u[2]*v[1], u[2]*v[0]-u[0]*v[2], u[0]*v[1]-u[1]*v[0]]
def dot(u,v) :return v[0]*u[0]+v[1]*u[1]+v[2]*u[2]
def triangle_normal_vec(a, b, c) :return cross(vec(a,b), vec(c,b))

def project_point_on_triangle_plane(p, a, b, c) :
    n = triangle_normal_vec(a, b, c)
    t = (n[0]*a[0]-n[0]*p[0] + n[1]*a[1]-n[1]*p[1] + n[2]*a[2]-n[2]*p[2])/float(n[0]*n[0]+n[1]*n[1]+n[2]*n[2])
    return [p[0]+t*n[0], p[1]+t*n[1], p[2]+t*n[2]]

def is_same_side(p1,p2, a,b) :
    cp1 = cross(vec(a,b), vec(a,p1))
    cp2 = cross(vec(a,b), vec(a,p2))
    if dot(cp1, cp2)>=0 :return True
    else :return False

def is_point_in_triangle(p, a,b,c) :
    p = project_point_on_triangle_plane(p, a, b, c)
    if is_same_side(p,a, b,c) and is_same_side(p,b, a,c) and is_same_side(p,c, a,b) :return True
    else :return False

a, b, c = [1,1,1], [7,1,5], [12,1,1]
p = [7,0,5]

raw_input(is_point_in_triangle(p, a,b,c))

Globalement, la démarche est la suivante :
- recupérer le vecteur normal au triangle (simple produit en croix de deux des vecteurs composant le triangle)
- a partir de là, on peut projeter le point dans le plan du triangle
- vérifier si le point appartient au triangle (pas trop comprit cette partie, mais ça à l'air de marcher)

Encore merci pour ton aide, j'y serais pas arrivé sans toi!