Discussion :

[diam's]

voila le problème :

J'ai un plan (de maison, d'immeuble, de champ, ...) sout forme de fichier JPEG

Sur ce plan, je place des points, lesquels définissent des segments de droites permettant de délimiter une surface CLOSE (pièce, appartement, parcelle, etc...), cad le premier point et le dernier point de la poylignes définie par les segments de droite sont les mêmes.

Je définis aussi les échelles (horizontale et verticale) par rapport des élements du plan dont je connais les dimensions.

Là où je sèche, c'est pour calculer la surface : la forme peut être quelconque (ceux qui ont fait de l'arpentage ou qui voient régulièrement certains plans d'architectes me comprendront, il n'y a qu'a voir l'image ci jointe).


J'ai déjà regardé les forums, notamment celui "Algo". Mais très souvent on parle de définir un tableau qu'on va "colorier", de découper en triangles, ou d'utiliser des outils mathématiques que j'ai un peu de mal à comprendre (Pour le découpage en triangle à partir d'un sommet fixe, je vois franchement pas comment faire pour l'image ci dessus). La solution la plus 'envisageable' que j'ai trouvé se trouve là :Code de j.p.mignot (en C++)

Pour mon grand malheur, en plus d'être largué en math, je ne comprends pas le C++ (je comprends en gros comment fonctionne le programme, mais je n'arrive pas à le traduire en Pascal). En fait que représente le * dans

Code :

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

float *DX

S'agit il d'une forme de pointeur ? et que représente de signe de pourcentage dans la partie suivante :

Code :

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

X[(i+1)%Npts]

je suppose qu'il ne s'agit pas d'une division, mais alors, qu'est-ce ?

Donc, quelqu'un pourrait-il m'aider à traduire cet algo en pascal ?
Ou, quelqu'un connait-il une méthode pour calculer la surface d'un élement (comme sur l'image ci-dessus) sans sortir la grosse artillerie (intégrale et autres joyeusetés )?

[sjrd]

En fait que représente le * dans

Code :

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

float *DX

S'agit il d'une forme de pointeur ?

Ce n'est pas une forme de pointeur, c'est un pointeur. DX est un pointeur sur float. En Pascal ça donne :

Code :

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

var DX : ^Single;

et que représente de signe de pourcentage dans la partie suivante :

Code :

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

X[(i+1)%Npts]

je suppose qu'il ne s'agit pas d'une division, mais alors, qu'est-ce ?

C'est stricutement équivalent au mod du Pascal.

Pour le calcul d'aire, je ne connais aucun de ces algorithmes, mais j'associerait le nom "triangularisation à un sommet fixe" à ce type de résolution :

• Tu choisis un des sommets de ta forme. Notons-le O. Les autres sommets seront, dans l'ordre, A, B, C, etc. jusque N. Tu peux tourner dans un sens ou dans l'autre, peu importe.
• Pour chaque paire d'autres sommets consécutifs : A-B, B-C, ..., M-N. On nomme cette paire X-Y.
• Tu calcules l'aire du triangle OXY par la formule de Héron :
  p = demi-périmètre = ( |OX| + |XY| + |YO| ) / 2
  aire = \/¨( p (p-|OX|) (p-|XY|) (p-|YO|) )
• Tu détermines si l'angle XOY est positif ou négatif (tu peux choisir le sens, mais par logique on prendra le sens trigonométrique).
• Si l'angle est positif, alors tu ajoutes l'aire du triangle OXY à l'aire totale de ta forme ; si l'angle est négatif, tu l'en retranches.
• Si à la fin, tu te retrouves avec une aire négative, tu prends sa valeur absolue - c'est juste que tu auras tourné dans le "mauvais" sens.

Si un connaisseur d'algo peut confirmer... Là je viens d'inventer ça sur le vif

[diam's]

Merci Sjrd pour les infos sur les conventions en C++

Par contre, pour ce qui est de la triangulation à sommet fixe, je ne vois pas comment faire. C'est effectivement simple avec une forme géométrique simple, mais avec la forme ci-dessous, c'est plus compliqué. Je m'explique :

Prenons ce dessin, tiré d'un cas réel (et oui, y'a parfois des architectes 'tordus')


Il est aisé, avec ta méthode, de calculer la surface du triangle BOC. Mais il est plus délicat de calculer celle du triangle IOJ. En effet, IOJ n'est pas 'plein' et en plus, une partie de sa surface se recoupe avec celle des triangles QOP, SOT, PON et TOU. A moins que je n'ai mal compris

Alors, j'ai pensé à déplacer le point O à l'intérieur de la forme. Mais quelque soit l'endroit où je le place, je retombe forcément sur le cas 'IOJ ' ci dessus.

N'y aurait-il pas quelqu'un qui aurait développé une fonction/un compo permettant de connaitre la surface d'une forme géométrique fermée dont un connait TOUS les sommets ?

[Pedro]

(et oui, y'a parfois des architectes 'tordus')

En tant que géomètre, je te réponds +100000

Par contre, si c'est bien la surface que tu veux, il existe une méthode pour la calculer à partir des coordonnées

[diam's]

Par contre, si c'est bien la surface que tu veux, il existe une méthode pour la calculer à partir des coordonnées

ça tombe bien, c'est ce que je cherche

par contre, cette méthode, tu l'aurais pas sous le coude par hasard ?

C'est pas un truc du style remplissage de tableau j'espère ? Par ce que, les tableaux, c'est va pas etre facile car :
- si cotation en cm, une case doit représenté 1cm sur 1 cm
- pour un batiment dont la surface projetée peut être inscrite dans un carré de 100 m x 100 m, cela fait 100 000 000 de cases. En admettant que chaque case contienne un booleen (la case est DANS la forme ou AU DEHORS), cela fait 100 000 000 d'octets (un bool est, je crois, codé sur un octet), soit environ 95 mO en mémoire...

Bon d'accord, c'est vrai qu'un bâtiment de 100/100 m, c'est pas fréquent, jusqu'à présent j'ai jamais dépassé le 120 x 50 m, mais quand même...

[LadyWasky]

Va vite voir ici si ça répond à ta question, parce que personnellement, c'est ce que j'utilise
http://www.efg2.com/Lab/Graphics/PolygonArea.htm

[Pedro]

ça tombe bien, c'est ce que je cherche

par contre, cette méthode, tu l'aurais pas sous le coude par hasard ?

Je pense que l'exemple que donne Waskol est le bon

C'est pas un truc du style remplissage de tableau j'espère ? Par ce que, les tableaux, c'est va pas etre facile car :
- si cotation en cm, une case doit représenté 1cm sur 1 cm
- pour un batiment dont la surface projetée peut être inscrite dans un carré de 100 m x 100 m, cela fait 100 000 000 de cases. En admettant que chaque case contienne un booleen (la case est DANS la forme ou AU DEHORS), cela fait 100 000 000 d'octets (un bool est, je crois, codé sur un octet), soit environ 95 mO en mémoire...

Bon d'accord, c'est vrai qu'un bâtiment de 100/100 m, c'est pas fréquent, jusqu'à présent j'ai jamais dépassé le 120 x 50 m, mais quand même...

Non rien de tout cela, ça utilise les coordonnées des extrémités Et puis c'est quoi cette manie qu'ont les architectes de travailler en cm... L'unité SI est le mètre!

[diam's]

Merci Pedro et Waskol

Je viens de lire le post de Waskol, et je venais de trouver ça.

Le fait que Pedro ai dit qu'il était

En tant que géomètre...

et qu'il ai formulé sa réponse dela façon suivante

une méthode pour la calculer à partir des coordonnées

m'avait mis sur la voie.

Les deux méthodes semblent cohérentes à un détail près :
dans le site que j'ai trouvé, l'aire est définie comme suit :


alors que le site donné par Waskol, elle est

avec


Ce qui me gêne c'est qu'on a d'un coté : Xi(Yi+1 - Yi-1)
et de l'autre : XiYi+1 - Xi+1Yi

J'ai donc le doute sur ma solution. Par contre, je n'ai aucun doute sur celle de Waskol notre maitre
J'aimerais juste comprendre la différence mathématique (je taggue le résolu demain ou lundi).

Par contre, ce serait peut-être bien de conserver cette question dans la FAQ, car elle semble être souvent posée, bien que sous des formulations différentes

[Pedro]

J'ai donc le doute sur ma solution. Par contre, je n'ai aucun doute sur celle de Waskol notre maitre

Effectivement si waskol l'utilise, c'est qu'elle doit être juste

[diam's]

ah tiens, Pedro m'a grillé

Mouarf

[LadyWasky]

ah tiens, Pedro m'a grillé

Mouarf

Non , sur ton lien le calcul est faux

D'ailleurs, revoici le bon calcul :
http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea
Ce qui fait du 2-1

Même que tu as un listing en Delphi ici :
http://local.wasp.uwa.edu.au/~pbourk...ea/polygon.pas

Et au moins, sur efg2 (le gars qui tient le site est une pointure en Maths), il a la décence de dire qu'il s'agit du Théorème de Green : au moins, comme ça, en retrouvant le théorème, le calcul est vérifiable

Sinon, faut pas le dire, mais j'était nul en Géométrie au Lycée, je me serais fait jeter comme Géomètere, Pedro est bien meilleur que moi en la matière

[the_clansman]

Bonjour à tous,

Je viens de tester la formule sur la figure suivante :
(1;-1)
(1;2)
(2;2)
(3;1)
(4;2)
(5;1)
(6;2)
(7;1)
(7;-1)

Je la dessine et je trouve A=15 (un rectangle de 6x2 additionné d'un demi-rectangle 6x1)

Les formules me donnent A=20...

1/2 Somme(Abs(Xi*(Yi+1-Yi-1)))
= (1/2)(2+3+0+0+0+0+0+21+14) = 20

1/2 Somme(Abs((Xi*Yi+1)-(Xi+1Yi)))
=(1/2)(2+0+4+2+6+4+8+14+0) = 20

Y'a comme un blême quelque part mais je sais pas où...

Y'a pas quelqu'un qu'aurait testé par hasard ?

James

PS : Merci à tous de me faire prendre la tête sur des problèmes de maths à 1h du mat Meuh non, je plaisante, j'adore ça en fait : mais je suis sûr que ça va m'empêcher de dormir cette histoire ! Mouarf !

[Pedro]

Je viens de dessiner ta surface avec AutoCAD et ça ne correspond pas du tout à ce que tu décris
Par contre, je trouve S=15.5
Vérifie tes données je pense que ça vient de là

[LadyWasky]

Je vous promet de me pencher sur cette question ;p

Edit : Par rapport à ta figure on doit bien trouver 15,5

D'ailleurs, le résultat de ton 2ème calcule devrait être 15,5

Et par rapport au résultat que tu escomptait, ton 2ème point devrait être (1,1), et non pas (1,2)


En plus, le calcul proposé est ABS(SOMME()) et non pas SOMME(ABS()), ce n'est pas pareil

Pourquoi la valeur absolue, tout simplement que si tu parcours les points de ton polygone dans le sens des aiguille d'une montre, tu obtient une valeur négative, et si tu fais le parcours dans le sens trigonométrique, tu obtiens la même valeur, mais positive, bref, les deux formules sont bonnes.
Dans le cas ou tu rajoute la valeur absolu tu obtient bien une aire, dans l'autre cas, tu obtiens l'aire plus une information supplémentaire sur la disposition de tes points (donnée par le signe du résultat)


Les formules me donnent (avec (1,1) comme 2ème point ) :

1/2 Somme( (Xi*(Yi+1-Yi-1)) )
= (1/2)(1+3+0+0+0+0+0-21+14) = -15,5 ==> Donc Aire=15,5 (FAUX)

1/2 Somme( (Xi*Yi+1)-(Xi+1Yi) ))
=(1/2)(2+0-4+2-6+4+8-14-6) = -15 ==> Donc Aire=15 (JUSTE)



En utilisant ton deuxième point (donc (1,2) , on obtient respectivement -16 et -15,5, c'est donc le deuxième résultat qui est juste



Comme on parcours les polygones dans le sens des aiguilles d'une montre, le résultat est négatif, et c'est normal


Voici le code utilisé :

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
 
procedure TForm1.Button1Click(Sender: TObject);
var P:array[0..8] of TPoint;
    i,s1,s2,a1,a2,ip,im:integer;
    r1,r2:double;
begin
  P[0].X:=1;
  P[0].Y:=-1;
 
  P[1].X:=1;
  P[1].Y:=1;
 
  P[2].X:=2;
  P[2].Y:=2;
 
  P[3].X:=3;
  P[3].Y:=1;
 
  P[4].X:=4;
  P[4].Y:=2;
 
  P[5].X:=5;
  P[5].Y:=1;
 
  P[6].X:=6;
  P[6].Y:=2;
 
  P[7].X:=7;
  P[7].Y:=1;
 
  P[8].X:=7;
  P[8].Y:=-1;
 
  s1:=0;
  s2:=0;
  Memo1.Clear;
  for i:=0 to 8 do
  begin
    ip:=(i+1) mod 9; //i+1
    im:=(i-1) mod 9; //i-1
    a1:=P[i].X*(P[ip].Y-P[im].Y); //première formule
    a2:=P[i].X*P[ip].Y-P[ip].X*P[i].Y; //2ème formule
    memo1.lines.Add(inttostr(i)+':'+#9+inttostr(a1)+'   '+inttostr(a2));
    inc(s1,a1); //Somme
    inc(s2,a2); //Somme
  end;
  memo1.lines.Add('________');
  r1:=s1/2; //Résultat1
  r2:=s2/2; //Resultat2
  memo1.lines.Add('r1 = '+FloatToStr(r1));
  memo1.lines.Add('r2 = '+FloatToStr(r2));
end;

[the_clansman]

J'ai enfin trouvé !

Donc, mes erreurs :

- Comme l'a spécifié l'un d'entre vous, il faut bien faire abs(somme(...)) et non pas somme(abs(...)) : mea culpa...
- La figure était bien avec (1,1) en second point comme vous l'avez si bien remarqué...

Maintenant, Waskol, les deux formules sont justes : le premier élément de la première formule est 2 et non pas 1 : en effet il faut prendre, comme spécifié dans la page décrivant la formule Y0=Yn...

Or, dans ton code tu écris

Code :

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

1
2
im:=(i-1) mod 9; //i-1
a1:=P[i].X*(P[ip].Y-P[im].Y); //première formule

Dans le cas où ton i=0 on ne sait pas ce que renvoie P[-1] mais certainement pas P[n] (P[8]) comme il est pourtant mentionné dans la formule...

Les deux formules rendent donc le même résultat :

A=(1/2)Abs(Somme(Xi(Yi+1-Yi-1)))=abs(2+3+0+0+0+0+0-21-14)/2=30/2=15
A=(1/2)Abs(Somme(XiYi+1-Xi+1Yi))=abs(2+0-4+2-6+4-8-14-6)/2=30/2=15

Et maintenant pourquoi rendent-elles le même résultat ?

Et bien oui, chers amis, il n'est plus 1h du matin et j'ai maintenant suffisamment les idées claires pour ne plus faire d'erreurs de copier-coller et confondre abs(somme avec somme(abs

En fait,ces calculs ne sont ni plus ni moins que des calculs de somme d'aires de triangle (le 1/2 factorisé devant la formule est un indice déterminant quand on sait que l'aire d'un triangle se caractérise par 1/2*base*hauteur). Ainsi, dans le cas d'un triangle rectangle, la hauteur n'étant rien d'autre que le côté adjacent à l'angle droit, la formule est d'autant plus simple qu'elle est 1/2*cote1*cote2 où cote1 et cote2 représentent les longueurs des côtés adjacents à l'angle droit dudit triangle.

Donc, on peut penser que ces deux formules ne sont rien d'autres que l'addition d'Aire de triangles calculés à partir de leurs coordonnées.

Maintenant pourquoi la deuxième formule est égale à la première ?

Considérons le triangle ABC, rectangle en B.

Et bien, la première formule consiste à additionner les aires du triangle ABC quand la seconde additionne les aires du triangle CDA, triangle CDA qui n'est ni plus ni moins que le complémentaire du triangle ABC ! Par complémentaire comprenez que ABCD est un rectangle (muni de ces 4 angles à 90degrés).

En d'autres termes, dans un cas, on calcule l'aire d'un demi rectangle et dans l'autre cas, on calcule l'aire de l'autre demi rectangle !

Si vous déssinez un polygone à 3 faces, (1,0) (1,1) (2,0) vous remarquerez sans peine cet état de fait !

Donc, comme dans Jacques Martin tout le monde a gagné car les deux formules sont exactes !

Merci à tous pour votre aide et désolé de vous avoir fait perdre un peu de votre temps par mon étourderie que je mettrai sous le compte de l'heure matinale (...et un peu sur la 1664 ).

A+

James.

[defluc]

J'ai fait une application de calcul de surfaces quelconques il y a des années par la méthode des trapèzes.

C'est enfantin et sans erreur ni approximation tant qu'il n'y a pas de courbes.

[diam's]

@ the_clansman :
Merci pour l'explication

@defluc :
Je connais pa la méthode des trapèzes, mais elle doit pouvoir s'appliquer à mon cas car je n'ai pas de courbe (enfin, si, mais elles sont 'approchée' sous forme de segments)

Je ne pensais pas avoir autant de succès avec ce topic

merci à tous

damien

[LadyWasky]

- Comme l'a spécifié l'un d'entre vous, il faut bien faire abs(somme(...)) et non pas somme(abs(...)) : mea culpa...

Bin.... c'est moi

Maintenant, Waskol, les deux formules sont justes : le premier élément de la première formule est 2 et non pas 1 : en effet il faut prendre, comme spécifié dans la page décrivant la formule Y0=Yn...

Oui, tu as raison, j'ai fais ça en vitesse sans trop regarder...
En fait c'était surtout pour mettre le doigt sur Abs(Somme())

Or, dans ton code tu écris

Code :

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

1
2
im:=(i-1) mod 9; //i-1
a1:=P[i].X*(P[ip].Y-P[im].Y); //première formule

JE vois où tu veux en venir mais...

Dans le cas où ton i=0 on ne sait pas ce que renvoie P[-1] mais certainement pas P[n] (P[8]) comme il est pourtant mentionné dans la formule...

Et toc, et voilà, tu es tombé dans le panneau
si i=0,
alors im=(i-1) mod 9, donc im=8

d'ailleurs, vu les indices de mon tableau (qui vont de 0 à 8), si je passais -1 comme indice à mon tableau (P[-1]), je me prendrais un "Access Violation" des familles, n'est-ce pas ?

Donc rassures toi, dans mon cas, le "truc", c'est qu'en fait il fallait fermer la figure pour avoir Yn=Y0, ce que je n'ai pas fait et ce qui conduit à un résultat faux (vois tu ce que je veux dire ?)





Quoiqu'il en soit, si mon code ne trouve pas le bon résultat, c'est que j'ai du me planter en beauté

[defluc]

D'un point de vue programmation, il faudra que je cherche dans mes archives pour retrouver la source.

Sur le principe, on lit les coordonnées du premier point, P[0], et on les affecte à x1 et y1.

Code :

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

1
2
  x1 := P[0].X;
  y1 := P[0].Y;

Les points suivants sont affectés dans une boucle de 1 au nombre de points.
Dans cette boucle, on affecte les coordonnées des points à x2 et y2.

Code :

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

1
2
3
4
5
6
7
8
9
  For i := 1 to 8 do
  begin
         x2 := P[i].X;
         y2 := P[i].Y;
        Trap := (y1+y2)/2*(x2-x1);  //  La surface du trapèze
         S := S + Trap;                  //  La surface totale
         x1 := x2;
         y1 := y2;
  end;

Fais un dessin et tu verras que l'algorithme s'avère juste.

La seule chose qui, selon mes souvenirs, pose problème, ce sont les valeurs négatives.

Mais il suffit de trouver les x et y minimum et, s'ils sont négatifs, d'effectuer une translation à chaque point en leur appliquant un Delta x et un Delta Y avant calcul.

[defluc]

La méthode des trapèzes confirme ces résultats.
Je la trouve visuellement plus facile à comprendre et donc, à mémoriser (mon application a été créée il y a 25 ans en Quickbasic).
Voici le code

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
procedure TForm1.FormCreate(Sender: TObject);
var P:array[0..8] of TPoint;
    i, x1, y1, x2, y2, DeltaY:integer;
    Surf, Trap:double;
begin
  P[0].X:=1; P[0].Y:=-1;
  P[1].X:=1; P[1].Y:=2;
  P[2].X:=2; P[2].Y:=2;
  P[3].X:=3; P[3].Y:=1;
  P[4].X:=4; P[4].Y:=2;
  P[5].X:=5; P[5].Y:=1;
  P[6].X:=6; P[6].Y:=2;
  P[7].X:=7; P[7].Y:=1;
  P[8].X:=7; P[8].Y:=-1;
 
 { P[0].X:=-2; P[0].Y:=-2;
  P[1].X:=-2; P[1].Y:=0;
  P[2].X:=-1; P[2].Y:=1;
  P[3].X:=0; P[3].Y:=0;
  P[4].X:=1; P[4].Y:=1;
  P[5].X:=2; P[5].Y:=0;
  P[6].X:=3; P[6].Y:=1;
  P[7].X:=4; P[7].Y:=0;
  P[8].X:=4; P[8].Y:=-2;  }
 
  DeltaY := 2147483647;
 
  For i := 0 to 8 do //  Détermination du Minimum pour la translation des Y
     if P[i].Y < DeltaY then DeltaY :=P[i].Y;
 
  if DeltaY > 0 then DeltaY := 0 else DeltaY := Abs(DeltaY);
  Surf := 0;
  Trap := 0;
  x1 := P[0].X;
  y1 := DeltaY + P[0].Y;
  Memo1.Clear;
  memo1.lines.Add('Point 0 : ' + IntToStr(P[0].X) + ' , ' +
                     IntToStr(P[0].Y) +  ' devient  ' + IntToStr(x1) + ' , ' +
                     IntToStr(y1) + '  Trap ' + Floattostr(Trap)+
                     ' Surf : ' + Floattostr(Surf));
  for i:=1 to 8 do
  begin
         x2 := P[i].X;
         y2 := DeltaY + P[i].Y;
         Trap := (y1+y2)/2*(x2-x1);
         Surf := Surf + Trap;
         x1 := x2;
         y1 := y2;
    memo1.lines.Add('Point '+inttostr(i)+' : ' + IntToStr(P[i].X) + ' , ' +
                     IntToStr(P[i].Y) +  ' devient  ' + IntToStr(x1) + ' , ' +
                     IntToStr(y1) + '  Trap ' + Floattostr(Trap)+
                     ' Surf : ' + Floattostr(Surf));
  end;
  memo1.lines.Add('________');
  memo1.lines.Add('Surf = '+FloatToStr(Surf));
end;