Discussion :

[ToTo13]

Bonjour,

tout est dans le titre :
- Je dispose d'un plan (équation, vecteur normal, ...) et d'un cube unitaire (arète de longueur 1) dont je connais la position et qui n'a subit aucune rotation.
- Je souhaiterai savoir comment calculer la surface du polygone formé par l'intersection du plan et du cube.

Trouver les points et les arètes d'intersections entre le plan et les arètes du cube est évident.
Calculer de manière "bourrin" la surface avec ces informations est faisable.


Ce que souhaiterai savoir, c'est s'il existe une manière intélligente d'y parvenir ?
Disons que intuitivement et d'après les caractéristiques mathématiques, il doit y avoir une fonction pouvant apporter la réponse...


Merci par avance...

[PRomu@ld]

En fait, intuitivement je ferais comme ça :

une fois que tu as les points d'intersections du plan et du cube, tu triangularises le polygone obtenu et pour chacun des triangles, tu calcules l'aire (méthode de héron) et tu sommes le tout.

[adaneels]

Comme le polygone obtenu par l'intersection du plan et du cube est quelconque, tu dois calculer les sommets de cette forme.
Donc calculer au maximum 12 points.

Si c'est un calcul d'aire, décomposé en triangles.
Mais si c'est toi qui a laissé ce message, tu le sais déjà :-/

Le mieux que je vois après pour optimiser le code à l'exécution si c'est un calcul d'aire, c'est de différencier les cas : calculer l'aire comme pour un losange pour 4 sommets par exemple.

Je ne vois pas de solution rapide et intelligente de faire cela, ou quoi que ce soit qui réduit le nombre de calculs.

[j.p.mignot]

1- choisir un système d'axe où le plan est la plan xOy et z sa normale
2- calculer dans ce repère les coordonnées du sommet du cube ( rotation + translation )
3- en déduire les points d'intersection des côtés du cube avec le plan xOy
4- utliser Stokes ( somme ( rot(A). ds ) = somme ( A.dl ) )
on choisi le champ A tel que rot(A) = +-k ( par exemple A=-i où A=j où A = 1/2(j-i), ou ... )
d'où S = 0.5 * some ( ydx-xdy) = somme ( ydx) = -somme(xdy) = ...
avec n intercection qui sont des droites on ramène cette intégralle à une somme de différences.
Attention c'est une aire algébrique!

[larnicebafteur]

Le polygone obtenu par l'intersection d'un cube par un plan peut être :

- un triangle ou
- un quadrilatère ou
- un pentagone ou
- un hexagone.

Dans un premier temps, il faut determiner dans quel cas on se situe puis connaitre les coordonnées des sommets de ce polygone.

Ensuite, si ce n'est pas un triangle, décomposer ce polygone en 2, 3 ou 4 triangles selon le cas puis calculer leurs aires.

Malheureusement, je ne pense pas qu'il existe de méthode directe !

[PRomu@ld]

- un triangle ou
- un quadrilatère ou
- un pentagone ou
- un hexagone.

Un point et une droite sont aussi des solutions.

[larnicebafteur]

Un point et une droite sont aussi des solutions.

Oui, mais je ne parlais que des polygones possibles. Pour calculer l'aire d'une droite ou d'un point, c'est immédiat ...

[borisd]

Un point et une droite sont aussi des solutions.

Pas plutôt un point et un segment ?
Ca fait avancer le schmiltruc tout ça

[PRomu@ld]

Pas plutôt un point et un segment ?

Oui, c'est vrai ...

Ca fait avancer le schmiltruc tout ça

Que proposes tu ? En fait pour avoir une solution analytique, il faudrait une équation du cube ( je ne pense pas que ce soit un conique ... )

[borisd]

Je ne propose rien (c'est ma remarque qui ne fait rien avancer, je ne juge pas la tienne).
A part peut-être ça : quand l'intersection est un point ou un segment, l'aire est nulle.

[j.p.mignot]

voici 1 code pascal qui illustre la solution que j'ai proposé ci-dessu.
j'ai vérifié dans quelques cas simples mais je ne garanti pas qu'il soit 'buggless'.
Par ailleur ce n'est pas optimisé.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
 
 
unit cross_area;
 
interface
 
type I2 = record n1,n2 : byte end;
const Paires_liees : array[1..12] of I2 = ( (n1:0;n2:1),(n1:0;n2:2),
                                            (n1:0;n2:4),(n1:7;n2:6),
                                            (n1:7;n2:5),(n1:7;n2:3),
                                            (n1:6;n2:2),(n1:6;n2:4),
                                            (n1:1;n2:3),(n1:1;n2:5),
                                            (n1:4;n2:5),(n1:2;n2:3));
type crd = record x,y,z : double end;    
var Cabc : array[0..7] of crd;
 
 
function Compute_area(Point_plan, Normale  : crd ; var n_cross : byte) : double;
 
implementation
uses math;   // max
procedure Q ( var C : crd; x1,y1,z1 : double);
   begin
   with C do
      begin
      x:=x1;
      y:=y1;
      z:=z1
      end
   end;
procedure Make_Cube;
   begin
   Q(Cabc[0],0,0,0);
   Q(Cabc[1],1,0,0);
   Q(Cabc[2],0,1,0);
   Q(Cabc[3],1,1,0);
   Q(Cabc[4],0,0,1);
   Q(Cabc[5],1,0,1);
   Q(Cabc[6],0,1,1);
   Q(Cabc[7],1,1,1);
   end;
function normalize_normale( var Normale : crd) : boolean;
var r : double;
   begin
   normalize_normale := false;
   r := sqr(Normale.x) + sqr(Normale.y) + sqr(Normale.z);
   if  r < 1e-20 then
      exit;
   r:=1/sqrt(r);
   Normale.x:=Normale.x*r;
   Normale.y:=Normale.y*r;
   Normale.z:=Normale.z*r;
   normalize_normale := true;
   end;
procedure r2 ( var C : crd; Point_Plan : crd );
   begin
   with C do
      begin
      x:=x - Point_Plan.x;
      y:=y - Point_Plan.y;
      z:=z - Point_Plan.z;
      end;
   end;
procedure translat(Point_Plan : crd);
var i : byte;
   begin
   for i:=0 to 7 do r2(Cabc[i], Point_Plan);
   end;
var a11,a12,a13,a21,a22,a23,a31,a32,a33 : double;   
procedure chx( i : byte );
var C1 : crd;
   begin
   C1 :=Cabc[i];
   Cabc[i].x:=C1.x*a11 + C1.y*a12 + C1.z*a13;
   Cabc[i].y:=C1.x*a21 + C1.y*a22 + C1.z*a23;
   Cabc[i].z:=C1.x*a31 + C1.y*a32 + C1.z*a33;
   end;
procedure change_base( Normale : crd );
var T : double;
var xu,xv,xw,yu,yv,yw : double;
    i         : byte;
   begin
   T:=sqr(Normale.x) + sqr(Normale.y);
   if ( T < 1e-35) then exit;   // pas de changement de repere si la normale est deja z
   T:=1/sqrt(T);
   // plan orthogonal => xu + yv + Normale.z = 0
   //  T * (v,-u,0) <-> (1,0,0)
   xu:=Normale.y*T; xv:=-Normale.x*T; xw:=0; // choix arbitraire de l'axe x' dans le plan normal à Normal
   // on en deduit Y' tel que       Y' = z' vect. x'
   yw:=Normale.x*xv - Normale.y*xu;
   yv:=-(Normale.x*xw-Normale.z*xu);
   yu:=Normale.y*xw-Normale.z*xv;
   (*
   transformation  f(xu,yu,Normale.x)= (1,0,0),  f(yu,yv,yw) =0,1,0)  ..
                                      -1
   =>     |  xu   yu   Normale.x  |
      M = |  xv   yv   Normale.y  |              M^-1 = transposee( commat(M))  *1/det(M)
          |  xw   yw   Normale.z  |
   *)
   a11:=yv*Normale.z-yw*Normale.y; a21:=-(xv*Normale.z-Normale.y*xw); a31:=xv*yw-xw*yv;
   a12:=-(yu*Normale.z-Normale.x*yw); a22:=xu*Normale.z-Normale.x*xw; a32:=-(xu*yw-xw*yu);
   a13:=yu*Normale.y-yv*Normale.x; a23:=-(xu*Normale.y-xv*Normale.x); a33:=xu*yv-xv*yu;
   for i:=0 to 7 do chx(i);
   end;
type C24 = array[1..24] of crd;
var Cross : C24;
procedure Selection_des_Paires_croisant_le_paln_Z0( var n_cross : byte );
var u,v,w,p : double; i : byte;
   begin
   n_cross := 0;
   for i:=1 to 12 do with paires_liees[i] do
      begin
      if ( Cabc[n1].z * Cabc[n2].z <=0 ) then // le segment croise le plan
         begin
         u:=Cabc[n2].x-Cabc[n1].x;
         v:=Cabc[n2].y-Cabc[n1].y;
         w:=Cabc[n2].z-Cabc[n1].z;
         if w <> 0 then // normalement toujours vrai
            begin
            p:=-Cabc[n1].z/w;
            inc(n_cross);
            with cross[n_cross] do
               begin
               z:=0;
               x:=Cabc[n1].x + p*u;
               y:=Cabc[n1].y + p*v;
               end;
            end
         else  // w=0 => le seg [n1,n2] est horizontal avec z1*z2 <=0 donc il est dans le plan
            begin // on ajoute donc les 2 points
            inc(n_cross);
            with cross[n_cross] do
               begin
               z:=0;
               x:=Cabc[n1].x;
               y:=Cabc[n1].y;
               end;
            inc(n_cross);
            with cross[n_cross] do
               begin
               z:=0;
               x:=Cabc[n2].x;
               y:=Cabc[n2].y;
               end;
            end
         end;
      end;
   end;    
function dst( i,j : byte ) : double;
   begin
   dst:=sqr(cross[i].x-cross[j].x) + sqr(cross[i].y-cross[j].y) + sqr(cross[i].z-cross[j].z);
   end;
procedure remove_doublons( var n_cross : byte);
var ok : boolean; i,j,k : byte;
   begin
   if ( n_cross <=1) then exit;
   i:=0;
      repeat
      inc(i);
      for j:= i+1 to n_cross do
      if dst(i,j) < 1e-10 then // point i = point j
         begin
         for k:=j to n_cross-1 do cross[k] := cross[k+1];
         dec(n_cross);
         i:=0;
         end;
      ok := i= n_cross-1;
      until ok=true;
   end;
function ydx ( a,b : byte ) double;
   begin
   ydx := (cross[a].y+cross[b].y)*(cross[a].x-cross[b].x)/2;
   end;
function A4( i1,i2,i3,i4: byte ) : double;
   begin
          A4 := abs( ydx(i1,i2) +
                         ydx(i2,i3) +
                         ydx(i3,i4) +
                         ydx(i4,i1)
                  );
   end;                     
function A6( i1,i2,i3,i4,i5,i6: byte ) : double;
   begin
          A6 := abs( ydx(i1,i2) +
                         ydx(i2,i3) +
                         ydx(i3,i4) +
                         ydx(i4,i5) +
                         ydx(i5,i6) +
                         ydx(i6,i1)
                  );
   end;
function compute_area_with_stokes(n_cross : byte) : double;
var A : double; i1,i2,i3,i4,i5,i6 : byte;
   begin
   case n_cross of
      3 : begin
          A := abs( 
                        ydx(1,2) +
                        ydx(2,3) +
                        ydx(3,1)
                  );
          end;
      4 : begin
          A:=  A4(1,2,3,4);
          A:= max( A,A4(1,2,4,3));
          A:= max( A,A4(1,3,2,4));
          A:= max( A,A4(1,3,4,2));
          A:= max( A,A4(1,4,2,3));
          A:= max( A,A4(1,4,3,2));
          end;
      6 :
         begin
         A:=0;
         for i2:= 2 to 6 do
         for i3:= 2 to 6 do if (i3 <> i2 ) then
         for i4:= 2 to 6 do if ((i4 <> i2) and ( i4 <> i3)) then
         for i5:= 2 to 6 do if ((i5 <> i2) and ( i5 <> i3) and ( i5 <> i4)) then
         for i6:= 2 to 6 do if ((i6 <> i2) and ( i6 <> i3) and ( i6 <> i4) and(i6 <> i5)) then
            A:=max(A,A6(1,i2,i3,i4,i5,i6));
         end
      else
         A := -2;  // normalement impossible
      end;
   compute_area_with_stokes := A;
   end;
function Compute_area(Point_plan, Normale  : crd; var n_cross : byte ) : double;
var err : double;
   begin
   err:=-1;
   Make_Cube;
   if normalize_normale(Normale) then
      begin
      err:=0;
      Translat(Point_Plan);  // faire passer le plan par (0,0,0)
      change_base(Normale);  // mettre la normale comme z et choisir un repere xy orthonormé
      Selection_des_Paires_croisant_le_paln_Z0(n_cross); // selection des segment croisant le plan z=0
      remove_doublons(n_cross);   // suppression des points comptes n X
      if ( n_cross <= 2) then
         err:=0    // avec 0,1 ou 2 points, l'aire est nulle
      else
         err:=compute_area_with_stokes(n_cross);   // n_cross = 3,4 ou 6 ici
      end;
   Compute_area:=err;  // aire = -1 si probleme dans norme, -2 si probleme dans aire
   end;
 
end.

[ToTo13]

Bonjour,

merci pour ce morceau de code...

Je vais voir cela...

[j.p.mignot]

Je reste à votre service pour tout commentaire ou discussion.
J'ai écri

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
...
 for i2:= 2 to 6 do
         for i3:= 2 to 6 do if (i3 <> i2 ) then
         for i4:= 2 to 6 do if ((i4 <> i2) and ( i4 <> i3)) then
         for i5:= 2 to 6 do if ((i5 <> i2) and ( i5 <> i3) and ( i5 <> i4)) then
         for i6:= 2 to 6 do if ((i6 <> i2) and ( i6 <> i3) and ( i6 <> i4) and(i6 <> i5)) then
            A:=max(A,A6(1,i2,i3,i4,i5,i6));
...

Ceci me semble "sale" au posible... Tout particulièrement cela n'est pas 'autoadaptatif' avec le nembre d'éléments à permuter! Toutes suggestions 'propre' récurcives pour scanner les différentes permutations est la bien venue!!!

merci!

[ToTo13]

Merci pour cette correction...

[ToTo13]

Bonjour,

le temps et mon directeur ne m'ont pas permis de tester plus avant cette méthode, mais je la garde sous le coude. Qui sait, peut être en parlerai je dans ma thèse pendant ma rédaction, mais là, il faut que l'on avance...


Merci à tous...