Récupérer deux points à partir d'un triangle et d'un plan

Version imprimable

16/02/2009, 23h56
crismans

Récupérer deux points à partir d'un triangle et d'un plan

Bonjour,

Je possède deux matrices, l'une qui représente des coordonnées de points, matrice notée M la seconde me permettant de définir mes triangles avec des indices de mes points, matrice T

Citation:

M=[8,5,10;8,5,0;5,5,0;5,5,0;5,5,10]
T=[1,2,3;3,5,1]

ps:j'ai bien noté que j'avais deux fois le même point, ce n'est pas une erreur.
Je définis les arrête de mes triangles comme ceci :

Citation:

% A premier sommet, B deuxième sommet
AB=[M(T(:,2),1)-M(T(:,1),1) M(T(:,2),2)-M(T(:,1),2) M(T(:,2),3)-M(T(:,1),3)];
% C troisème sommet
AC=[M(T(:,3),1)-M(T(:,1),1) M(T(:,3),2)-M(T(:,1),2) M(T(:,3),3)-M(T(:,1),3)];

Je définis le vecteur normal de chaque plan à l'aide de la fonction cross.
Je penser procéder comme suit:
Trois points A,B,C
-> 1 plan : ABxAC=N
->n1Xa + n2Ya + n3Za = d
Intersection du plan précédent avec plan z=h où h est une valeur que j'incrémente.
-> droite (D)
-> Intersection de la droite avec AB, on obtient le point D, je vérifie si D appartient au segment [AB], si oui je récupère le point.
-> Intersection de la droite avec AC, on obtient le point F, je vérifie si F appartient au segment [AC], si oui je récupère le point.
Ensuite je trace l'intersection en utilisant la fonction line. avec D mon points de départ et F mon point d'arriver, stockés dans une matrice n lignes 2 colonnes.

Donc j'aimerais récupérer les points D et F pour chaque triangles, mathématiquement avec une feuille et un crayon je sais le faire, sur matlab ça devient beaucoup plus difficiles.
Est ce que mon raisonnement est correct ?

Je vous remercie pour votre aide.

Cordialement
17/02/2009, 09h11
Jerome Briot

Citation:

Envoyé par crismans

mathématiquement avec une feuille et un crayon je sais le faire, sur matlab ça devient beaucoup plus difficiles.
Est ce que mon raisonnement est correct ?

Tu dis que tu sais le faire à la main... mais tu nous demande si ton raisonnement est correcte ? :koi:

De plus il y a exactement la même discussion sur le forum Algorithme > Mathématiques : Problème intersection plan triangles :koi:²

Si la fonction cross réalise bien un produit vectoriel (j'ai pas Matlab sous la main donc je conditionne ma réponse...), alors ton raisonnement est effectivement correct à 2-3 exceptions :

- l'intersection d'un plan avec un triangle peut être le triangle lui-même, l'ensemble vide, ou un segment, selon que ton triangle et ton plan sont (resp.)coplanaires, de normales colinéaires sans être coplanaires ou de normales non colinéaires (je te donne même dans cette phrase les éléments pour faire ton test). Il vaudrait mieux tester ces conditions avant d'essayer de trouver le point d'intersection de AB avec ton plan ! (PS : l'intersection peut aussi être un unique point)

- Ensuite, tu cherches le point d'inetrsection de ton plan avec AB et AC, mais tu conçois bien que tu peux avoir une intersection avec BC ! Donc un troisième système à résoudre !

Donc ton raisonnement sur le papier est bon... quant à la mise en oeuvre en Matlab, vas-y doucement, surtout si c'est ton premier programme (c'est le cas ?).

Pour démarrer, construit une fonction qui détermine l'intersection si elle existe d'un segment et d'un plan. Allez va, un petit coup de pouce, je te file l'en tête, tu complètes et tu reviens au cas où :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 function [intersection, point] = plane_segment_intersection(normal,yintercept,segmenta, segmentb) % Teste si un segment et un plan s'intersecte et renvoi l'intersection % Sortie : % intersection : 'none' signifie pas d'intersection % 'full' signifie que le segment est contenu dan le plan % 'point' signifie que le segment a un point % d'intersection avec le plan % % point : contient les coordonnées du point d'intersection si % intersection == 'point' % Entrées : % normal : contient la normale définissant le plan % yintercept : contient la valeur de d de l'équation du plan % ax+by+cz+d=0 % segmenta : contient le premier point du segment % segmentb : contient le second point du segment % Pseudo-code % Tester les cas particuliers % - intersection 'none' : si les deux points sont du meme coté du plan % - intersection 'full' : si les deux points sont dans le plan % - intersection 'point' : si les deux points sont de part et d'autre du % plan % Cas 'point' (les autres cas sont triviaux) % - résoudre le système (ie poser le système "le point appartient au % segment ET au plan") % Fin
Avec ce début, la suite viendra peut être plus facilement.

17/02/2009, 10h14
crismans

Oui effectivement, c'est mon premier vrai programme, je demandais si mon raisonnement était bon car je ne sais pas si avec matlab on raisonne comme on le fait sur une feuille de papier.

Merci Loic, j'avais bien penser aux autres cas, j'aurais du le mentionner. J'essai de retranscrire ce que j'ai sur ma feuille en suivant l'ordre que vous m'avez donné. Car c'était surtout ça mon problème, comment généraliser ma démarche sous matlab.

Merci, bonne journée

Bonjour,

Voilà d'où en est mon travail, tout d'abord voilà une partie de la matrice sur laquelle je travail :

Code:

1
2
 
Matrice_final = [1 80 50 100; 2 80 50 0;3 50 50 0; 4 50 50 100;5 50 70 0;6 50 70 100;7 80 70 0;8 80 70 100;9 100 100 100;10 0 100 100;11 0 0 100;12 100 0 100;13 100 100 0;14 100 0 0;15 0 0 0;16 0 100 0]

J'ai commencé par traiter la partie où mes triangles sont contenus dans le plan

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
Z=hauteur;                                                              %hauteur du plan de coupe 
    for i=1:length(Matrice_final)-1;                                    %debut de boucle sur la Matrice finale
    A=[Matrice_final(i,2) Matrice_final(i,3) Matrice_final(i,4)];       %premier point du segment
    B=[Matrice_final(i+1,2) Matrice_final(i+1,3) Matrice_final(i+1,4)]; %deuxieme point du segment
    k=(B(3)-A(3))/(Z-A(3)); 
    knu=(B(3)-A(3));
    kde=(Z-A(3));                            %AM=kAB (vecteurs), comme on donne Z, on connait k
 
        if knu==kde                          %A et B sont dans le même plan horizontal
              X(i)=A(1);                     %il y a deux points M donc on doit retrouver deux coordonnées dans la matrice M 
              Y(i)=A(2); 
              M(i,1:2)=[X(i) Y(i)];          %le point A n'est plus dans le plan de coupe 
                                             %mais A et B sont toujours dans le même plan horizontal
 
        else                                              %k~=0 donc A et B sont dans deux plan horizontaux distincts
        X(i)=(1/k)*(B(1)-A(1))+A(1);                                        
        Y(i)=(1/k)*(B(2)-A(2))+A(2);
        M(i,1:2)=[X(i) Y(i)];                             %on place les valeurs de X et Y du point M dans une matrice
        end  
    end
    [M]=M
    [MY]=M(:,2)
    [MX]=M(:,1)
    plot(MX,MY)

Je sais j'avance pas vite, mais c'est le début... ;-)

Bonjour,

J'ai testé autre chose car je n'arrive pas à mettre mon cheminement sous matlab.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
 
% Arêtes
% Deuxième - Premier sommet
AB=[M(T(:,2),1)-M(T(:,1),1) M(T(:,2),2)-M(T(:,1),2) M(T(:,2),3)-M(T(:,1),3)];
% Troisème - Premier sommet
AC=[M(T(:,3),1)-M(T(:,1),1) M(T(:,3),2)-M(T(:,1),2) M(T(:,3),3)-M(T(:,1),3)];
% Troisème - Deuxième sommet
BC=[M(T(:,3),1)-M(T(:,2),1) M(T(:,3),2)-M(T(:,2),2) M(T(:,3),3)-M(T(:,2),3)];
 
Np=[0 0 1] %vecteur normal de mon plan 
d=100 %dans l'équation du plan z+d=0
 
%intersection avec le vecteur AB
%full
for i = 1:1:length(AB); 
    a = 1:3:length(M);
    b = 2:3:length(M);
    c = 3:3:length(M);
    j = 1:1:length(M);
    u1 = M(b,1)-M(a,1);
    u2 = M(b,2)-M(b,2);
    u3 = M(b,3)-M(b,3);
    if dot(AB(i,: ),Np)== 0 %produit scalaire, si =0 alors le vecteur est // au plan
        if M(T(a,1),3)== d %Si le z du point A est égal à d alors le 
                           %vecteur est dans le plan en entier
            A=[M(a,1) M(a,2)]; %essai
            Xs(j)= A(1); %essai
            Ys(j)= A(2); %essai
            V(j,1:2)= [X(j) Y(i)] %essai
        else %none
            resultat='none'
        end
    else %point, produit scalaire différent de 0, il faut vérifier que 
         %le point d'intersection de coordonnées (xi, yi, zi) 
         %appartient au plan et au segment
        zi=-d;
        k=(-d-M(a,3))/u3;
        yi=M(a,2)+u2*k;
        xi=M(a,1)+u1*k;
        Xs(j)=xi; %essai
        Ys(j)=yi; %essai
        V(j,1:2) = [X(j) Y(j)] %essai
    end
end
 
%essai : je voudrais récupérer une matrice deux colonnes contenant tout mes
%points dans l'ordre

Pour chaque triangle faut que je test mes vecteurs AB, AC et BC et que je récupère les données dans une matrice. La j'essai de faire tous les vecteurs AB puis ensuite tous les AC puis les BC. Mais il faudrait peut être que je le fasse pour un triangle et ensuite que je fasse une boucle pour tous mes triangles, comme je les possède dans l'ordre dans ma matrice T ?
Je crois que je suis perdu...
Quelqu'un pourrait il m'aider ?
Merci.

22/02/2009, 18h48
crismans

DUT, j'ai répondu au lien que vous m'avez donné, j'avais déjà vu ce message et j'ai repris son exemple car ce que je cherche à faire sous matlab c'est un petit programme qui permet de faire ça sur les trois vecteurs de chaque triangles. J'ai reformulé la question de bende dans le forum Matlab, car ce n'est pas un problème de math que je rencontre mais plus un problème de "traduction" de ma démarche.
Merci et désolé d'avoir repris le même message.

Essaie de faire comme ceci :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
j = 1;
for i = 1:1:length(AB); 
    ...
    if dot(AB(i,: ),Np)== 0 %produit scalaire, si =0 alors le vecteur est // au plan
        if M(T(a,1),3)== d %Si le z du point A est égal à d alors le 
            ...
            V(j,1:2)= [X(j) Y(i)] %essai
            j = j+1;
        else %none
            resultat='none'
        end
    else 
        V(j,1:2) = [X(j) Y(j)] %essai
        j = j+1;
    end
end

Par contre, je ne sais pas d'où sortent tes variables X et Y ?

22/02/2009, 21h43
crismans

D'accord, je vais essayé ça. Oupss Les variables X et Y sont en fait Xs et Ys.
Merci

J'ai enfin réussi à résoudre mon problème, enfin presque. Il me manque encore deux petite chose mais je n'y arrive pas.

Je remet mon programme pour ceux que ça interresse.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
'0<k<1'
for i = 1:3:length(V);
 
 
    AB = [V(i+1,1)-V(i,1) V(i+1,2)-V(i,2) V(i+1,3)-V(i,3)];
    AC = [V(i+2,1)-V(i,1) V(i+2,2)-V(i,2) V(i+2,3)-V(i,3)];
    BC = [V(i+2,1)-V(i+1,1) V(i+2,2)-V(i+1,2) V(i+2,3)-V(i+1,3)];
 
    if dot(N,AB) == 0
        if V(i,3) == d
            M(i,1:2)=[V(i,1) V(i,2)];
            M(i,1:2)=[V(i+1,1) V(i+1,2)];
        else
        end
    else
        k = (-d-V(i,3))/(V(i+1,3)-V(i,3));
        yi = V(i,2) + k*(V(i+1,2)-V(i,2));
        xi = V(i,1) + k*(V(i+1,1)-V(i,1));
        M(i,1:2)=[xi yi];
    end
 
    if dot(N,AC) == 0;
        if V(i,3) == d;
            M(i,1:2)=[V(i,1) V(i,2)];
            M(i,1:2)=[V(i+2,1) V(i+2,2)];
        else
        end
    else
        k = (-d-V(i,3))/(V(i+2,3)-V(i,3));
        yi = V(i,2) + k*(V(i+2,2)-V(i,2));
        xi = V(i,1) + k*(V(i+2,1)-V(i,1));
        M(i,1:2)=[xi yi];
    end
 
    if dot(N,BC) == 0;
        if V(i+1,3) == d;
            M(i,1:2)=[V(i+1,1) V(i+1,2)];
            M(i,1:2)=[V(i+2,1) V(i+2,2)];
        else
        end
    else
        k = (-d-V(i+1,3))/(V(i+2,3)-V(i+1,3));
        yi = V(i+1,2) + k*(V(i+2,2)-V(i+1,2));
        xi = V(i+1,1) + k*(V(i+2,1)-V(i+1,1));
        M(i,1:2)=[xi yi];
    end
 
end
M=M;
 
idx=M~=0;
MM=M(idx)

Je n'arrive à mettre une condition sur k, je voudrais : 0<k<1, et que si k sort de cette plage, il ne tienne pas compte du calcul.

Quelque chose comme ceci ?

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
    if dot(N,AB) == 0
        if V(i,3) == d
            M(i,1:2)=[V(i,1) V(i,2)];
            M(i,1:2)=[V(i+1,1) V(i+1,2)];
        else
        end
    else
        k = (-d-V(i,3))/(V(i+1,3)-V(i,3));
 
        if k>0 & k<1
            yi = V(i,2) + k*(V(i+1,2)-V(i,2));
            xi = V(i,1) + k*(V(i+1,1)-V(i,1));
            M(i,1:2)=[xi yi];
        end
    end

25/02/2009, 19h39
crismans

Et dire que j'ai réussi à me perdre avec deux conditions If alors qu'il me fallait "juste" un petit &. Merci beaucoup pour votre aide.

Mon programme à encore évoluer quelque peu, j'ai un petit souci sur la comprehension de la fonction plot, en fait j'ai la matrice P suivante
Code:

1 2 3 4 5 6 7 8 9 10 90 %départ segment 1 10 10 %arrivé segment 1 90 10 %départ segment 2 10 10 %arrivé segment 2 90 90 %départ segment 3 10 90 %arrivé segment 3 90 90 %départ segment 4 90 10 %arrivé segment 4
et j'aimerais tracer mes 4 segment?

Code:

plot(P(j,:),P(j+1,:))

Merci

Essaie ceci :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
P = [10    90 %départ segment 1
10    10 %arrivé segment 1
90    10 %départ segment 2
10    10 %arrivé segment 2
90    90 %départ segment 3
10    90 %arrivé segment 3
90    90 %départ segment 4
90    10]; %arrivé segment 4
 
x = reshape(P(:,1),2,[]);
y = reshape(P(:,2),2,[]);
 
figure
line(x,y)
axis equal

27/02/2009, 21h23
crismans

Une nouvelle fois MERCI