Discussion :

[JP CASSOU]

Bjr à vous,

J'ai trouvé une unité Delphi pour la triangulation de Delaunay (merci à son auteur) et que j'ai adapté à mon projet GHTopo.
Je constate un fonctionnement complètement erratique de cette unité:
- En mode graphique, il ne semble pas y avoir de problème (cependant, quelques crashs surviennent)
- En mode piloté, c'est le crash d'emblée ...

A voir ce qui foire

[B]EDIT: Semble résolu. Je laisse cette unité en CC0.[/B]

Ma variante utilise les structures suivantes:

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
 
type TPoint3Df = record
  X:  double;
  Y:  double;
  Z:  double;
end;
type TPoint2Df = record
  X:  double;
  Y:  double;
end;  
type TVertex = record
  ID         : Int64;
  X          : double;
  Y          : double;
  Z          : double;
  NormX      : double;
  NormY      : double;
  NormZ      : double;
  Norme      : double;
end;
type TMNTBoundingBox = record
  C1: TPoint3Df;
  C2: TPoint3Df;
end;      
type TMNTTriangleABC = record
   Numero  :  integer;
   PointA  :  integer;
   PointB  :  integer;
   PointC  :  integer;
   Marked  :  boolean;
   BoundingBox: TMNTBoundingBox;
end;

Les fonctions AfficherMessage() et AfficherMessageErreur() :

Code :

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

1
2
3
4
5
6
7
8
9
 
procedure AfficherMessage(const Msg: string);
begin
  WriteLn(Msg);
end;
procedure AfficherMessageErreur(const Msg: string);
begin
  WriteLn(Msg);
end;

et s'utilise comme suit:

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
 
// Exemple: Soit un tableau de TVertex MonTableauDeVertex déclaré et initialisé préalablement
// avec TTableauDeVertex = array of TVertex;
procedure TriangulerUnTableauDeVertex(const MonTableauDeVertex: TTableauDeVertex);
begin
 MyTriangulationDelaunay := TTriangulationDelaunay.Create;
  try
    // initialiser avec les coordonnées du coin bas gauche et coin haut droit de la zone à mailler
    if (MyTriangulationDelaunay.Initialiser(FCoordsMini, FCoordsMaxi)) then
    begin
      //.... Ajouter les points du maillage ici. La triangulation se fait à la volée
 
      for i := 0 to High(MonTableauDeVertex) do
      begin 
         MyVertex := MonTableauDeVertex [i];
         MyTriangulationDelaunay.AddAPoint(MakeTPoint3Df(MyVertex.X, MyVertex.Y, MyVertex.Z), false);
      end;
 
 
     // récupérer et utiliser le résultat 
 
      NbTriangles := MyTriangulationDelaunay.GetNbTriangles();
      NbVertex    := MyTriangulationDelaunay.GetNbPoints();
      AfficherMessage(format('Triangulation terminée %d vertex, %d triangles', [NbVertex, NbTriangles]));
      if (NbVertex > 0) then
      begin
        FTableVertex.ClearListe();
 
        for i := 1 to NbVertex-1 do
        begin
          QPoint := MyTriangulationDelaunay.GetPoint(i);
          // travail avec le QPoint
          //..........
        end;
      end;
      AfficherMessage('Les triangles');
      if (NbTriangles > 0) then
      begin
        FtableTriangles.ClearListe();
        for i := 0 to NbTriangles - 1 do
        begin
          QTriangle := MyTriangulationDelaunay.GetTriangle(i);
          // ... travail avec le QTriangle
        end;
      end;
    end;
    MyTriangulationDelaunay.Finaliser();
    result := True;
  finally
    MyTriangulationDelaunay.Free;
  end;                                            
end;

L'unité:

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
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
 
unit UnitTriangulationDelaunay;
// 17/01/2019: Fonctionnement correct. TODO: Orientation des triangles
//
{$mode delphi}
{$DEFINE USES_TLISTESIMPLES_POINTS}
{$UNDEF  USES_TLISTESIMPLES_POINTS}
{$DEFINE USES_TLISTESIMPLES_TRIANGLES}
{$UNDEF  USES_TLISTESIMPLES_TRIANGLES}
 
interface
 
uses
  StructuresDonnees,
  Common,
  {$IFDEF USES_TLISTESIMPLES_POINTS}
  UnitListesSimplesWithGeneriques,
  {$ENDIF}
  {$IFDEF USES_TLISTESIMPLES_TRIANGLES}
  UnitListesSimplesWithGeneriques,
  {$ENDIF}
  math,
  Classes, SysUtils;
 
type PPoint3Df = ^TPoint3Df;
type TEquationDroite = record // équation de droite ax+by+c=0
       a,b,c:double;
end;
 
type
  PTriangle=^TTriangle;
  TTriangle = record
            p1,p2,p3: PPoint3Df;       // les trois sommets
            t1,t2,t3: PTriangle;       // les trois triangles ayant comme arretes communes celles
                                       // oposées respectivement à p1 p2 p3
            d1,d2,d3: TEquationDroite; // equations des arretes oposées respectivement à p1 p2 p3
            c       : TPoint2Df;       // centre du cercle circonscrit
            r       : double;          // rayon au carré du cercle circonscrit
           end;
 
{$IFDEF USES_TLISTESIMPLES_POINTS}
type TListePoints = class(TListeSimple<PPoint3Df>)
  private
  public
end;
{$ENDIF}
{$IFDEF USES_TLISTESIMPLES_TRIANGLES}
type TListeTriangles = class(TListeSimple<PTriangle>)
  private
  public
end;
{$ENDIF}
 
const MAX_VERTEX    = 15000;
const MAX_TRIANGLES = 30000;
type
 
{ TTriangulationDelaunay }
 
 TTriangulationDelaunay = class
  private
    {$IFDEF USES_TLISTESIMPLES_POINTS}
    FListePoints    : TListePoints;
    {$ELSE}
    FListePoints    : array[0 .. MAX_VERTEX] of PPoint3Df;
    {$ENDIF}
 
    {$IFDEF USES_TLISTESIMPLES_TRIANGLES}
    FListeTriangles: TListeTriangles;
    {$ELSE}
    FListeTriangles : array[0 .. MAX_TRIANGLES] of PTriangle;
    {$ENDIF}
    FNbPoints       : integer;
    FNbTriangles    : integer;
    FLimitesDuMNT   : TRect2Df;
    FNbErreurs      : integer;
    function AddTriangle(const t: TTriangle; out TriangleAdded: PTriangle): boolean;
    function AddPoint(const QX, QY, QZ: Double; out PointAdded: PPoint3Df): boolean;
    function Algorythm_Delaunay(const x, y, z: double): boolean;
    function DelTriangle(const p: ptriangle; out T: TTriangle): boolean;
    function FindIdxOfPoint(const P: PPoint3Df): integer;
    function FindTriangle(const p: PPoint3Df; out PT: PTriangle): boolean;
    procedure ResetListes();
    function PointIsNotTooNearExisting(const P: TPoint3Df): boolean;
 
 
  public
    function  Initialiser(const QCoinBasGauche, QCoinHautDroit: TPoint3Df): boolean;
    procedure Finaliser();
    function  GetLimitesDuMNT(): TRect2Df;
 
 
    function  GetNbPoints(): integer;
    function  GetNbTriangles(): integer;
 
    procedure AddAPoint(const P: TPoint3Df; const DoMesh: boolean);
    function  GetPoint(const Idx: integer): TPoint3Df;
    function  GetTriangle(const Idx: integer): TMNTTriangleABC;
    function  GetNbErreursTriangulation(): integer;
 
end;
 
implementation
//calcule la distance au carré entre p1 et p2
function sqdistance(const p1, p2: PPoint3Df): double;
begin
 result := (p1.X-p2.X)*(p1.X-p2.X) + (p1.y-p2.y)*(p1.y-p2.y);
end;
 
// crée un triangle
function CreateTriangle(const p1, p2, p3: PPoint3Df; out MyTriangle: TTriangle): boolean;
var
 diviseur: double;
 d1,d2,d3: TEquationDroite;
 c       : TPoint2Df;
begin
 Result := false;
 // calcul de la droite media de [p1;p2]
 d1.a:=-(p1.X-p2.X);
 d1.b:=-(p1.Y-p2.Y);
 d1.c:=-(d1.a*(p1.x+p2.X) / 2) - (d1.b*(p1.y+p2.Y) / 2);
 // calcul de la droite media de [p2;p3]
 d2.a:=-(p2.X-p3.X);
 d2.b:=-(p2.Y-p3.Y);
 d2.c:=-(d2.a*(p2.x+p3.X) / 2) - (d2.b*(p2.y+p3.Y) / 2);
 // calcul de la droite media de [p1;p3]
 d3.a:=-(p1.X-p3.X);
 d3.b:=-(p1.Y-p3.Y);
 d3.c:=-(d3.a*(p1.x+p3.X) / 2) - (d3.b*(p1.y+p3.Y) / 2);
 
 // calcul du point d'inter de d1 et d2
 diviseur := d1.a*d2.b-d2.a*d1.b;
 
  if (IsZero(diviseur)) then Exit(false);   // = 0 si les trois points sont alignés
 
 
 c.X:=(d1.b*d2.c-d2.b*d1.c) / (diviseur);
 c.y:=(d2.a*d1.c-d1.a*d2.c) / (diviseur);
 
 MyTriangle.r:= sqdistance(@c, p1);
 
 MyTriangle.p1 := p1;
 MyTriangle.p2 := p2;
 MyTriangle.p3 := p3;
 
 // calcul de la droite [p1;p2]
 d3.a:=-(p1.y-p2.y);
 d3.b:= (p1.x-p2.x);
 d3.c:=-(d3.a*p1.x) - (d3.b*p1.y);
 // calcul de la droite [p2;p3]
 d1.a:=-(p2.y-p3.y);
 d1.b:= (p2.x-p3.x);
 d1.c:=-(d1.a*p2.x) - (d1.b*p2.y);
 // calcul de la droite [p1;p3]
 d2.a:=-(p1.y-p3.y);
 d2.b:= (p1.x-p3.x);
 d2.c:=-(d2.a*p1.x) - (d2.b*p1.y);
 
 MyTriangle.d1:=d1;
 MyTriangle.d2:=d2;
 MyTriangle.d3:=d3;
 MyTriangle.t1:=nil;
 MyTriangle.t2:=nil;
 MyTriangle.t3:=nil;
 MyTriangle.c:=c;
 
 Result := True;
end;
 
// est-ce que p1 et p2 sont du même coté de d ?
function ptsMemeCotedeD(const p1,p2: PPoint3Df; const d: TEquationDroite):boolean;
begin
 result:=(sign(d.a*p1.x+d.b*p1.y+d.c)*sign(d.a*p2.x+d.b*p2.y+d.c)>0)
end;
 
// p est-il dans t
function PtInTriangle(const p: PPoint3Df; const t: PTriangle):boolean;
begin
 result := ptsMemeCotedeD(p,t.p1,t.d1) and
           ptsMemeCotedeD(p,t.p2,t.d2) and
           ptsMemeCotedeD(p,t.p3,t.d3);
end;
// retourne le sommet de t qui n'est ni a ni b
function Get3rdPoint(const t: PTriangle; const a,b: PPoint3Df; out ThirdPoint: PPoint3Df): boolean;
begin
 result := false;
 if t=nil then exit(false);
 if (( t.p1 = a) and ( t.p2 = b)) or (( t.p1 = b) and ( t.p2 = a)) then ThirdPoint := t.p3 ;
 if (( t.p2 = a) and ( t.p3 = b)) or (( t.p2 = b) and ( t.p3 = a)) then ThirdPoint := t.p1 ;
 if (( t.p3 = a) and ( t.p1 = b)) or (( t.p3 = b) and ( t.p1 = a)) then ThirdPoint := t.p2 ;
 Result := True;
end;
 
//  retourne le triangle voisin de t par l'arete [a,b]
function GetTriangleVoisin(const t: PTriangle; const a,b: PPoint3Df; out TriangleVoisin: ptriangle): Boolean;
begin
 result:=false;
 if ( t = nil) then exit(false);
 if (( t.p1 = a) and ( t.p2 =b)) or (( t.p1 =b) and ( t.p2 = a)) then TriangleVoisin := t.t3;
 if (( t.p2 = a) and ( t.p3 =b)) or (( t.p2 =b) and ( t.p3 = a)) then TriangleVoisin := t.t1;
 if (( t.p3 = a) and ( t.p1 =b)) or (( t.p3 =b) and ( t.p1 = a)) then TriangleVoisin := t.t2;
 Result := True;
end;
 
procedure actualise_voisin(tv, t: PTriangle;const a,b: PPoint3Df);
begin
 if tv=nil then exit;
    if ((tv.p1=a) and (tv.p2=b)) or ((tv.p1=b) and (tv.p2=a)) then tv.t3:=t
    else
  if ((tv.p1=a) and (tv.p3=b)) or ((tv.p1=b) and (tv.p3=a)) then tv.t2:=t
    else tv.t1:=t;
end;
 
 
// inverse les deux triangles d'arrete commune [a,b]
function swaptriangles(t1, t2: PTriangle; const a, b: PPoint3Df): boolean;
var
 p1,p2,p3,p4:PPoint3Df;
 v1,v2,v3,v4:PTriangle;
 foo: Boolean;
begin
 Result := false;
 p1 := a;
 p3 := b;
 foo := Get3rdPoint(t1,a,b, p2);
 foo := Get3rdPoint(t2,a,b, p4);
 
 foo := GetTriangleVoisin(t1,p1,p2, v1);
 foo := GetTriangleVoisin(t1,p3,p2, v2);
 foo := GetTriangleVoisin(t2,p3,p4, v3);
 foo := GetTriangleVoisin(t2,p1,p4, v4);
 
 if (not CreateTriangle(p1,p2,p4, t1^)) then Exit(false);
 if (not CreateTriangle(p3,p2,p4, t2^)) then Exit(false);
 
 
 t1.t1:=t2;
 t1.t2:=v4;
 t1.t3:=v1;
 
 t2.t1:=t1;
 t2.t2:=v3;
 t2.t3:=v2;
 actualise_voisin(v1,t1,p1,p2);
 actualise_voisin(v2,t2,p2,p3);
 actualise_voisin(v3,t2,p3,p4);
 actualise_voisin(v4,t1,p4,p1);
 Result := True;
end;
 
 
 
// applique la 2ème règle de Delaunay, aucun point dans le cercle circonscrit
procedure realigntriangle(t: PTriangle);
var
 t1,t2,t3:ptriangle;
 p1,p2,p3:PPoint3Df;
begin
 if t=nil then exit;
 AfficherMessageErreur('realign');
 try
   Get3rdPoint(t.t1,t.p2,t.p3, p1);
   Get3rdPoint(t.t2,t.p1,t.p3, p2);
   Get3rdPoint(t.t3,t.p1,t.p2, p3);
   begin
 
     // test avec les trois voisins
     if (p1<>nil) and (sqdistance(@t.c,p1)<t.r) then
      begin
       swaptriangles(t,t.t1,t.p2,t.p3);
       realigntriangle(t.t2);
       realigntriangle(t.t3);
      end
     else
     if (p2<>nil) and (sqdistance(@t.c,p2)<t.r) then
      begin
       swaptriangles(t,t.t2,t.p1,t.p3);
       realigntriangle(t.t1);
       realigntriangle(t.t3);
      end
     else
     if (p3<>nil) and (sqdistance(@t.c,p3)<t.r) then
      begin
       swaptriangles(t,t.t3,t.p1,t.p2);
       realigntriangle(t.t1);
       realigntriangle(t.t2);
      end;
   end;
 
 except
   AfficherMessageErreur('Echec en realign');
 end;
end;
 
 
// atribut les triangles t1,t2,t3 à t
procedure defineTriangle(t: PTriangle; const t1,t2,t3: PTriangle);
begin
 t^.t1:=t1; t^.t2:=t2; t^.t3:=t3;
end;
 
///*****************************************************************************
{ TTriangulationDelaunay }
 
function TTriangulationDelaunay.Initialiser(const QCoinBasGauche, QCoinHautDroit: TPoint3Df): boolean;
var
  p1, p2, p3, p4: PPoint3Df;
  dummy: PTriangle;
  QTriangle1, QTriangle2: TTriangle;
 
  {$IFDEF USES_TLISTESIMPLES_TRIANGLES}QTriangle0, QTriangle1: PTriangle;  {$ENDIF}
begin
 AfficherMessageErreur(Format('%s.Initialiser:', [classname]));
 result := false;
 try
    {$IFDEF USES_TLISTESIMPLES_POINTS}
    FListePoints := TListePoints.Create;
    {$ENDIF}
    {$IFDEF USES_TLISTESIMPLES_TRIANGLES}
    FListeTriangles := TListeTriangles.Create;
    {$ENDIF}
    self.ResetListes();
    FLimitesDuMNT.X1 := QCoinBasGauche.X;
    FLimitesDuMNT.Y1 := QCoinBasGauche.Y;
    FLimitesDuMNT.X2 := QCoinHautDroit.X;
    FLimitesDuMNT.Y2 := QCoinHautDroit.Y;
 
    AfficherMessageErreur(Format('-- (%f, %f) -> (%f, %f)', [FLimitesDuMNT.X1, FLimitesDuMNT.Y1, FLimitesDuMNT.X2, FLimitesDuMNT.Y2]));
 
    // crée les quatres points, coins de l'image
    Addpoint(QCoinBasGauche.X, QCoinBasGauche.Y, QCoinBasGauche.Z, p1);
    Addpoint(QCoinHautDroit.X, QCoinBasGauche.Y, QCoinBasGauche.Z, p2);
    Addpoint(QCoinHautDroit.X, QCoinHautDroit.Y, QCoinBasGauche.Z, p3);
    AddPoint(QCoinBasGauche.X, QCoinHautDroit.Y, QCoinBasGauche.Z, p4);
 
    // cre les deux triangles initiaux
    if (createtriangle(p1,p2,p3, QTriangle1)) then addtriangle(QTriangle1, dummy);
    if (createtriangle(p1,p3,p4, QTriangle2)) then addtriangle(QTriangle2, dummy);
 
    // lie les deux premiers triangles
    {$IFDEF USES_TLISTESIMPLES_TRIANGLES}
    QTriangle0 := FListeTriangles.GetElement(0);
    QTriangle1 := FListeTriangles.GetElement(1);
    QTriangle0.t1 := QTriangle1;
    QTriangle1.t3 := QTriangle0;
    FListeTriangles.PutElement(0, QTriangle0);
    FListeTriangles.PutElement(1, QTriangle1);
 
    {$ELSE}
    FListeTriangles[0].t1:=FListeTriangles[1];
    FListeTriangles[1].t3:=FListeTriangles[0];
 
    {$ENDIF}
    result := true;
  except
    AfficherMessageErreur('Echec de démarrage');
  end;
 
end;
procedure TTriangulationDelaunay.Finaliser();
begin
  AfficherMessageErreur(Format('%s.Finaliser:', [classname]));
  try
    self.ResetListes();
 
  finally
    {$IFDEF USES_TLISTESIMPLES_POINTS}
    FListePoints.Free;
    {$ENDIF}
    {$IFDEF USES_TLISTESIMPLES_TRIANGLES}
    FListeTriangles.Free;
    {$ENDIF}
  end;
 
end;
procedure TTriangulationDelaunay.ResetListes();
begin
  {$IFDEF USES_TLISTESIMPLES_POINTS}
  FListePoints.ClearListe();
  {$ELSE}
  //SetLength(FListePoints, 0);
  {$ENDIF}
  {$IFDEF USES_TLISTESIMPLES_TRIANGLES}
  FListeTriangles.ClearListe();
  {$ELSE}
  //SetLength(FListeTriangles, 0);
  {$ENDIF}
 
  FNbPoints    := 0;
  FNbTriangles := 0;
  FNbErreurs   := 0;
end;
// pour éviter des bugs, vérifier si le point candidat est suffisamment éloigné
function TTriangulationDelaunay.PointIsNotTooNearExisting(const P: TPoint3Df): boolean;
const
  DIST_MINI    = 0.5;
  SQ_DIST_MINI = DIST_MINI * DIST_MINI;
var
  i, Nb: Integer;
  PP: TPoint3Df;
  dMin: Double;
  d: ValReal;
begin
  dMin := INFINI;
  result := false;
  {$IFDEF USES_TLISTESIMPLES_POINTS}
  Nb := FListePoints.GetNbElements();
  {$else}
  Nb := Length(FListePoints);
  {$endif};
  if (Nb = 0) then exit(false);
  for i := 0 to Nb-1 do
  begin
    PP := GetPoint(i);
    d  := sqr(P.X - PP.x) + sqr(P.y - PP.y);
    if (d < dMin) then dMin := d;
  end;
  Result := (dmin > SQ_DIST_MINI);
end;
 
 
function TTriangulationDelaunay.GetNbErreursTriangulation(): integer;
begin
 result := FNbErreurs;
end;
 
 
 
function TTriangulationDelaunay.GetLimitesDuMNT(): TRect2Df;
begin
  result := FLimitesDuMNT;
end;
 
function TTriangulationDelaunay.GetNbPoints(): integer;
begin
  {$IFDEF USES_TLISTESIMPLES_POINTS}
  result := FListePoints.GetNbElements();
  {$ELSE}
  result := FNbPoints; //length(FListePoints);
  {$ENDIF}
end;
 
function TTriangulationDelaunay.GetNbTriangles(): integer;
begin
  {$IFDEF USES_TLISTESIMPLES_TRIANGLES}
  result := FListeTriangles.GetNbElements();
  {$ELSE}
  result := FNbTriangles;// length(FListeTriangles);
  {$ENDIF}
end;
 
procedure TTriangulationDelaunay.AddAPoint(const P: TPoint3Df; const DoMesh: boolean);
var
  S: Boolean;
begin
  //s := PointInRectangle(P, FLimitesDuMNT);
  //if (s) then
  Algorythm_Delaunay(P.X, P.Y, P.Z);  // OK ? On ajoute
end;
 
function TTriangulationDelaunay.GetPoint(const Idx: integer): TPoint3Df;
var
  PP: PPoint3Df;
begin
  {$IFDEF USES_TLISTESIMPLES_POINTS}
  PP := FListePoints.GetElement(Idx);
  Result := PP^;
  {$ELSE}
  Result := FListePoints[Idx]^;
  {$ENDIF}
end;
 
function TTriangulationDelaunay.GetTriangle(const Idx: integer): TMNTTriangleABC;
var
  TT: PTriangle;
begin
  {$IFDEF USES_TLISTESIMPLES_TRIANGLES}
  TT := FListeTriangles.GetElement(Idx);
  {$else}
  TT := FListeTriangles[Idx];
  {$endif}
  Result.Numero := Idx;
  Result.PointA := FindIdxOfPoint(TT.p1);
  Result.PointB := FindIdxOfPoint(TT.p2);
  Result.PointC := FindIdxOfPoint(TT.p3);
end;
 
// ajoute un point à listepoint
function TTriangulationDelaunay.AddPoint(const QX, QY, QZ: Double; out PointAdded: PPoint3Df): boolean;
begin
  result := false;
  {$IFDEF USES_TLISTESIMPLES_POINTS}
    New(PointAdded);
    PointAdded^.X := QX;
    PointAdded^.Y := QY;
    PointAdded^.Z := QZ;
    FListePoints.AddElement(PointAdded);
 
  {$else}
  inc(FNbPoints);
  //SetLength(FListePoints, FNbPoints);
  New(FListePoints[FNbPoints-1]);
  FListePoints[FNbPoints-1]^ := MakeTPoint3Df(QX, QY, QZ);
  PointAdded := FListePoints[FNbPoints - 1];
  {$ENDIF}
  result := true;
end;
// ajoute un point à listepoint
function TTriangulationDelaunay.FindIdxOfPoint(const P: PPoint3Df): integer;
var
  Nb, i: Integer;
begin
  result := -1;
  {$IFDEF USES_TLISTESIMPLES_POINTS}
  Nb := FListePoints.GetNbElements();
  if (Nb = 0) then Result := -1;
  for i := 0 to Nb-1 do
  begin
    if (P = FListePoints.GetElement(i)) then Exit(i);
  end;
  {$ELSE}
  Nb := FNbPoints;//length(FListePoints);
  if (Nb = 0) then Result := -1;
  for i := 0 to Nb-1 do
  begin
    if (P = FListePoints[i]) then Exit(i);
  end;
  {$ENDIF}
end;
 
//ajoute un triangle à listetriangle
function TTriangulationDelaunay.AddTriangle(const t: TTriangle; out TriangleAdded: PTriangle): boolean;
begin
 Result := false;
 try
 {$IFDEF USES_TLISTESIMPLES_TRIANGLES}
   New(TriangleAdded);
   TriangleAdded^ := t;
   FListeTriangles.AddElement(TriangleAdded);
 
 {$ELSE}
 
   inc(FNbTriangles);
   //setlength(FListeTriangles, FMaxTriangle);
   new(FListeTriangles[FNbTriangles-1]);
   FListeTriangles[FNbTriangles-1]^:= t;
   TriangleAdded := FListeTriangles[FNbTriangles-1];
 
 
 {$ENDIF}
   result := True;
   AfficherMessageErreur('AddTriangle()' + booltostr(result, 'OK', 'KO'));
 except
   Result := false;
 end;
end;
 
// extract et supprime le triangle n et ses références dans les autres triangles
function TTriangulationDelaunay.DelTriangle(const p: ptriangle; out T: TTriangle): boolean;
var
 i:integer;
 n, Idxfound: integer;
 EWE: PTriangle;
begin
  Result := false;
  Idxfound := -1;
  {$IFDEF USES_TLISTESIMPLES_TRIANGLES}
     // on recherche le triangle dans la liste
     n := FListeTriangles.GetNbElements();
     if (n = 0) then Exit(false);
     for i := 0 to n-1 do
     begin
     if (p = FListeTriangles.GetElement(i)) then
     begin
       Idxfound := i;
       break;
     end;
    end;
    if (Idxfound = -1) then Exit(false);
    T := P^;     // retourne la structure pointé par p
 
    // efface le triangle dans la liste
    FListeTriangles.RemoveElement(Idxfound);
    // remplace les références à p par nil dans le reste de la liste
    n := FListeTriangles.GetNbElements();
    if (n = 0) then Exit(false);
    for i:= 0 to FMaxTriangle-1 do
    begin
      EWE := FListeTriangles.GetElement(i);
      if (EWE.t1 = p) then EWE.t1 := nil;
      if (EWE.t2 = p) then EWE.t2 := nil;
      if (EWE.t3 = p) then EWE.t3 := nil;
      FListeTriangles.PutElement(i, EWE);
    end;
 
  {$ELSE}
 
 
  // on recherche le triangle dans la liste
  n := Length(FListeTriangles);
  for i := 0 to n-1 do
  begin
   if (p = FListeTriangles[i]) then
   begin
     Idxfound := i;
     break;
   end;
  end;
  if (Idxfound = -1) then Exit(false);
  // retourne la structure pointé par p
  T := p^;
  // libere la mémoire
  dispose(FListeTriangles[Idxfound]);
  // efface le triangle dans la liste
  dec(FNbTriangles);
  for i:= Idxfound to FNbTriangles - 1 do FListeTriangles[i]:=FListeTriangles[i+1];
  //SetLength(FListeTriangles, FMaxTriangle);  // setlength() préserve les données du tableau (cf REDIM PRESERVE du Basic)
  // remplace les références à p par nil dans le reste de la liste
  for i:=0 to FNbTriangles-1 do
  begin
   if (FListeTriangles[i].t1 = p) then FListeTriangles[i].t1 := nil;
   if (FListeTriangles[i].t2 = p) then FListeTriangles[i].t2 := nil;
   if (FListeTriangles[i].t3 = p) then FListeTriangles[i].t3 := nil;
  end;
  Result := True;
  {$ENDIF}
end;
 
// cherche et extract le triangle contenant p
function TTriangulationDelaunay.FindTriangle(const p: PPoint3Df; out PT: PTriangle): boolean;
var
 i, Nb :integer;
 t :ptriangle;
begin
 result := false;
 {$IFDEF USES_TLISTESIMPLES_POINTS}
 Nb := FListeTriangles.GetNbElements();
 for i :=0 to Nb - 1 do
 begin
   PT :=  FListeTriangles.GetElement(i);
   if (PtInTriangle(p, PT)) then exit(True);
 end;
 {$ELSE}
 for i :=0 to FNbTriangles-1 do
 begin
   if (PtInTriangle(p, FListeTriangles[i])) then
   begin
     PT:= FListeTriangles[i];
     exit(true);
   end;
 end;
 {$endif}
end;
 
// applique la 1er règle de Delaunay, detruire et construire
function TTriangulationDelaunay.Algorythm_Delaunay(const x, y, z: double): boolean;
var
 ptri:PTriangle;
 oldtri, QTriangle1, QTriangle2, QTriangle3:TTriangle;   // ancien triangle
 t1,t2,t3:PTriangle; // les trois nouveaux triangles
 p1,p2,p3,p4: PPoint3Df; // les trois sommets + le nouveau au centre
 v1,v2,v3:PTriangle; // les trois voisins
 IsOK: boolean;
 
 
 procedure Cuicui(const b: boolean; const Etape: string; out Succeeded: boolean);
 var
   EWE: String;
 begin
   Succeeded := b;
   EWE := Format('%s %s', [Etape, IIF(b, 'OK', 'KO')]);
   AfficherMessageErreur(EWE);
 end;
 
begin
 AfficherMessageErreur(Format('%s.Algorythm_Delaunay: %f, %f, %f', [ClassName, x, y, z]));
 result := false;
 try
   Cuicui(AddPoint(x,y,z, p4), 'Ajout du point', IsOK);
   Cuicui(FindTriangle(p4, ptri), 'Recherche triangle', IsOK);
 
   //supprime le triangle précédent
   Cuicui(DelTriangle(ptri, oldtri), 'Delete this triangle', IsOK);
   //if (not IsOK) then Exit;
 
 
   AfficherMessageErreur('002');
   // récupère les trois sommets et les trois voisins
   p1 := oldtri.p1;
   p2 := oldtri.p2;
   p3 := oldtri.p3;
 
   v1 := oldtri.t1;
   v2 := oldtri.t2;
   v3 := oldtri.t3;
 
   // en construit 3 à la place
   if (CreateTriangle(p4,p2,p3, QTriangle1)) then AddTriangle(QTriangle1, t1);
   if (createtriangle(p4,p3,p1, QTriangle2)) then AddTriangle(QTriangle2, t2);
   if (createtriangle(p4,p1,p2, QTriangle3)) then AddTriangle(QTriangle3, t3);
 
     // attribut les triangles voisins à t1 t2 et t3
     defineTriangle(t1,v1,t2,t3);
     defineTriangle(t2,v2,t3,t1);
     defineTriangle(t3,v3,t1,t2);
 
     // fait pointer les voisins vers les trois triangles
     if v1<>nil then actualise_voisin(v1,t1,t1.p2,t1.p3);
     if v2<>nil then actualise_voisin(v2,t2,t2.p2,t2.p3);
     if v3<>nil then actualise_voisin(v3,t3,t3.p2,t3.p3);
     // réaligne les 3 triangles en appliquant la 2ème règle
 
       realigntriangle(t1);
       realigntriangle(t2);
       realigntriangle(t3);
 
   Result := True;
 except
   FNbErreurs += 1;
   AfficherMessageErreur(Format('Failed %s.Algorythm_Delaunay: %f, %f, %f: Triangle %d', [ClassName, x, y, z, GetNbTriangles()]));
   AfficherMessageErreur(Format('* %.2f < %.2f < %.2f', [FLimitesDuMNT.X1, X, FLimitesDuMNT.X2]));
   AfficherMessageErreur(Format('* %.2f < %.2f < %.2f', [FLimitesDuMNT.Y1, Y, FLimitesDuMNT.Y2]));
 end;
end;
end.

Jeu de données d'exemple:

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
 
44.691947685276    -0.37296687850858    9.000
44.686211870225    -0.37189399490262    5.000
44.684777916462    -0.36631500015164    7.000
44.696106151188    -0.37790214309598    11.000
44.69524577893    -0.372323148345    43.000
44.688792986998    -0.37704383621121    2.000
44.688362800869    -0.37125026473905    5.000
44.694672197425    -0.3650275398245    88.000
44.687502428611    -0.36803161392118    9.000
44.689796754632    -0.37296687850858    5.000
44.687789219364    -0.36760246047879    12.000
44.68563828872    -0.37554179916287    2.000
44.695675965059    -0.35944854507352    53.000
44.68563828872    -0.37146484146024    5.000
44.688076010116    -0.37833129653836    5.000
44.687215637859    -0.37511264572049    4.000
44.69022694076    -0.37811671981717    4.000
44.6918042899    -0.36631500015164    71.000

[BeanzMaster]

Salut pour faire la triangulation 2D tu peux tenter d'adapter cette unité provenant de GLScene

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
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
//
// This unit is part of the GLScene Project, http://glscene.org
//
{
   Classes and methods for 2D triangulation of scatter points. 
  
   History :  
      17/05/15 - PW - Created, based on code from Paul Bourke (http://paulbourke.net/)
         and conversion from VB to Delphi by Dr Steve Evans (steve@lociuk.com)
 
}
 
unit GLTriangulation;
 
interface
 
uses
  Classes,
  Types,
  Dialogs,
  Graphics,
 
  GLVectorGeometry;
 
// Set these as applicable
const
  MaxVertices = 500000;
 
const
  MaxTriangles = 1000000;
 
const
  ExPtTolerance = 0.000001;
 
type
  TDProgressEvent = procedure(State: string; Pos, Max: Integer;
    AlwaysUpdate: Boolean = False) of object;
 
  // Points (Vertices)
type
  DVertex = record
    X: Single;
    Y: Single;
    Z: Single; // Added to save height of terrain.
    U: Single;
    V: Single;
    MatIndex: Integer;
  end;
 
  // Created Triangles, vv# are the vertex pointers
type
  DTriangle = record
    vv0: LongInt;
    vv1: LongInt;
    vv2: LongInt;
    PreCalc: Integer;
    Xc, Yc, R: Single;
  end;
 
type
  TDVertex = array of DVertex;
  TDTriangle = array of DTriangle;
  TDComplete = array of Boolean;
  TDEdges = array of array of LongInt;
 
{  TGLDelaunay2D is a class for Delaunay triangulation of arbitrary points
  Credit to Paul Bourke (http://paulbourke.net/) for the original Fortran 77 Program :))
  Conversion to Visual Basic by EluZioN (EluZioN@casesladder.com)
  Conversion from VB to Delphi6 by Dr Steve Evans (steve@lociuk.com)
  June 2002 Update by Dr Steve Evans (steve@lociuk.com): Heap memory allocation
  added to prevent stack overflow when MaxVertices and MaxTriangles are very large.
  Additional Updates in June 2002: Bug in InCircle function fixed. Radius r := Sqrt(rsqr).
  Check for duplicate points added when inserting new point.
  For speed, all points pre-sorted in x direction using quicksort algorithm and
  triangles flagged when no longer needed. The circumcircle centre and radius of
  the triangles are now stored to improve calculation time.
  You can use this code however you like providing the above credits remain in tact
}
 
type
  TGLDelaunay2D = class
  private
 
    function InCircle(Xp, Yp, X1, Y1, X2, Y2, X3, Y3: Single; var Xc: Single;
      var Yc: Single; var R: Single; j: Integer): Boolean;
    function Triangulate(nvert: Integer): Integer;
  public
 
    Vertex: TDVertex;
    Triangle: TDTriangle;
    HowMany: Integer;
    tPoints: Integer; // Variable for total number of points (vertices)
    OnProgress: TDProgressEvent;
    constructor Create;
    destructor Destroy;override;
    procedure Mesh(sort: Boolean);
    procedure AddPoint(X, Y, Z, U, V: Single; MatIndex: Integer);
    procedure AddPointNoCheck(X, Y, Z, U, V: Single; MatIndex: Integer);
    procedure RemoveLastPoint;
    procedure QuickSort(var A: TDVertex; Low, High: Integer);
  end;
 
//------------------------------------------------------------------------
//------------------------------------------------------------------------
//------------------------------------------------------------------------
implementation
//------------------------------------------------------------------------
//------------------------------------------------------------------------
//------------------------------------------------------------------------
 
constructor TGLDelaunay2D.Create;
begin
  // Initiate total points to 1, using base 0 causes problems in the functions
  inherited;
  tPoints := 1;
  HowMany := 0;
  SetLength(Vertex, MaxVertices);
  SetLength(Triangle, MaxTriangles);
  OnProgress := nil;
end;
 
destructor TGLDelaunay2D.Destroy;
begin
  SetLength(Vertex, 0);
  SetLength(Triangle, 0);
  inherited;
end;
 
function TGLDelaunay2D.InCircle(Xp, Yp, X1, Y1, X2, Y2, X3, Y3: Single;
  var Xc: Single; var Yc: Single; var R: Single; j: Integer): Boolean;
// Return TRUE if the point (xp,yp) lies inside the circumcircle
// made up by points (x1,y1) (x2,y2) (x3,y3)
// The circumcircle centre is returned in (xc,yc) and the radius r
// NOTE: A point on the edge is inside the circumcircle
var
  eps: Single;
  m1: Single;
  m2: Single;
  mx1: Single;
  mx2: Single;
  my1: Single;
  my2: Single;
  dx: Single;
  dy: Single;
  rsqr: Single;
  drsqr: Single;
begin
  eps := 0.000001;
  InCircle := False;
  // Check if xc,yc and r have already been calculated
  if Triangle[j].PreCalc = 1 then
  begin
    Xc := Triangle[j].Xc;
    Yc := Triangle[j].Yc;
    R := Triangle[j].R;
    rsqr := R * R;
    dx := Xp - Xc;
    dy := Yp - Yc;
    drsqr := dx * dx + dy * dy;
  end
  else
  begin
    if (Abs(Y1 - Y2) < eps) and (Abs(Y2 - Y3) < eps) then
    begin
      ShowMessage('INCIRCUM - F - Points are coincident !!');
      Exit;
    end;
    if Abs(Y2 - Y1) < eps then
    begin
      m2 := -(X3 - X2) / (Y3 - Y2);
      mx2 := (X2 + X3) / 2;
      my2 := (Y2 + Y3) / 2;
      Xc := (X2 + X1) / 2;
      Yc := m2 * (Xc - mx2) + my2;
    end
    else if Abs(Y3 - Y2) < eps then
    begin
      m1 := -(X2 - X1) / (Y2 - Y1);
      mx1 := (X1 + X2) / 2;
      my1 := (Y1 + Y2) / 2;
      Xc := (X3 + X2) / 2;
      Yc := m1 * (Xc - mx1) + my1;
    end
    else
    begin
      m1 := -(X2 - X1) / (Y2 - Y1);
      m2 := -(X3 - X2) / (Y3 - Y2);
      mx1 := (X1 + X2) / 2;
      mx2 := (X2 + X3) / 2;
      my1 := (Y1 + Y2) / 2;
      my2 := (Y2 + Y3) / 2;
      if (m1 - m2) <> 0 then // se
      begin
        Xc := (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
        Yc := m1 * (Xc - mx1) + my1;
      end
      else
      begin
        Xc := (X1 + X2 + X3) / 3;
        Yc := (Y1 + Y2 + Y3) / 3;
      end;
    end;
    dx := X2 - Xc;
    dy := Y2 - Yc;
    rsqr := dx * dx + dy * dy;
    R := Sqrt(rsqr);
    dx := Xp - Xc;
    dy := Yp - Yc;
    drsqr := dx * dx + dy * dy;
 
    // store the xc,yc and r for later use
    Triangle[j].PreCalc := 1;
    Triangle[j].Xc := Xc;
    Triangle[j].Yc := Yc;
    Triangle[j].R := R;
  end;
 
  if drsqr <= rsqr then
    InCircle := True;
end;
 
function TGLDelaunay2D.Triangulate(nvert: Integer): Integer;
// Takes as input NVERT vertices in arrays Vertex()
// Returned is a list of NTRI triangular faces in the array
// Triangle(). These triangles are arranged in clockwise order.
var
  Complete: TDComplete;
  Edges: TDEdges;
  Nedge: LongInt;
  // For Super Triangle
  xmin: Single;
  xmax: Single;
  ymin: Single;
  ymax: Single;
  xmid: Single;
  ymid: Single;
  dx: Single;
  dy: Single;
  dmax: Single;
 
  // General Variables
  i: Integer;
  j: Integer;
  k: Integer;
  ntri: Integer;
  Xc: Single;
  Yc: Single;
  R: Single;
  inc: Boolean;
begin
  // Allocate memory
  SetLength(Complete, MaxTriangles);
  SetLength(Edges, 2, MaxTriangles * 3);
 
  // Find the maximum and minimum vertex bounds.
  // This is to allow calculation of the bounding triangle
  xmin := Vertex[1].X;
  ymin := Vertex[1].Y;
  xmax := xmin;
  ymax := ymin;
  for i := 2 to nvert do
  begin
    if Vertex[i].X < xmin then
      xmin := Vertex[i].X;
    if Vertex[i].X > xmax then
      xmax := Vertex[i].X;
    if Vertex[i].Y < ymin then
      ymin := Vertex[i].Y;
    if Vertex[i].Y > ymax then
      ymax := Vertex[i].Y;
  end;
 
  dx := xmax - xmin;
  dy := ymax - ymin;
  if dx > dy then
    dmax := dx
  else
    dmax := dy;
 
  xmid := Trunc((xmax + xmin) / 2);
  ymid := Trunc((ymax + ymin) / 2);
 
  // Set up the supertriangle
  // This is a triangle which encompasses all the sample points.
  // The supertriangle coordinates are added to the end of the
  // vertex list. The supertriangle is the first triangle in
  // the triangle list.
 
  Vertex[nvert + 1].X := (xmid - 2 * dmax);
  Vertex[nvert + 1].Y := (ymid - dmax);
  Vertex[nvert + 2].X := xmid;
  Vertex[nvert + 2].Y := (ymid + 2 * dmax);
  Vertex[nvert + 3].X := (xmid + 2 * dmax);
  Vertex[nvert + 3].Y := (ymid - dmax);
  Triangle[1].vv0 := nvert + 1;
  Triangle[1].vv1 := nvert + 2;
  Triangle[1].vv2 := nvert + 3;
  Triangle[1].PreCalc := 0;
 
  Complete[1] := False;
  ntri := 1;
 
  // Include each point one at a time into the existing mesh
  for i := 1 to nvert do
  begin
    if Assigned(OnProgress) then
      OnProgress('Delaunay triangulation', i - 1, nvert);
    Nedge := 0;
    // Set up the edge buffer.
    // If the point (Vertex(i).x,Vertex(i).y) lies inside the circumcircle then the
    // three edges of that triangle are added to the edge buffer.
    j := 0;
    repeat
      j := j + 1;
      if Complete[j] <> True then
      begin
        inc := InCircle(Vertex[i].X, Vertex[i].Y, Vertex[Triangle[j].vv0].X,
          Vertex[Triangle[j].vv0].Y, Vertex[Triangle[j].vv1].X,
          Vertex[Triangle[j].vv1].Y, Vertex[Triangle[j].vv2].X,
          Vertex[Triangle[j].vv2].Y, Xc, Yc, R, j);
        // Include this if points are sorted by X
        if { usingsort and } ((Xc + R) < Vertex[i].X) then //
          Complete[j] := True //
        else //
          if inc then
          begin
            Edges[0, Nedge + 1] := Triangle[j].vv0;
            Edges[1, Nedge + 1] := Triangle[j].vv1;
            Edges[0, Nedge + 2] := Triangle[j].vv1;
            Edges[1, Nedge + 2] := Triangle[j].vv2;
            Edges[0, Nedge + 3] := Triangle[j].vv2;
            Edges[1, Nedge + 3] := Triangle[j].vv0;
            Nedge := Nedge + 3;
            Triangle[j].vv0 := Triangle[ntri].vv0;
            Triangle[j].vv1 := Triangle[ntri].vv1;
            Triangle[j].vv2 := Triangle[ntri].vv2;
            Triangle[j].PreCalc := Triangle[ntri].PreCalc;
            Triangle[j].Xc := Triangle[ntri].Xc;
            Triangle[j].Yc := Triangle[ntri].Yc;
            Triangle[j].R := Triangle[ntri].R;
            Triangle[ntri].PreCalc := 0;
            Complete[j] := Complete[ntri];
            j := j - 1;
            ntri := ntri - 1;
          end;
      end;
    until j >= ntri;
    // Tag multiple edges
    // Note: if all triangles are specified anticlockwise then all
    // interior edges are opposite pointing in direction.
    for j := 1 to Nedge - 1 do
    begin
      if not(Edges[0, j] = 0) and not(Edges[1, j] = 0) then
      begin
        for k := j + 1 to Nedge do
        begin
          if not(Edges[0, k] = 0) and not(Edges[1, k] = 0) then
          begin
            if Edges[0, j] = Edges[1, k] then
            begin
              if Edges[1, j] = Edges[0, k] then
              begin
                Edges[0, j] := 0;
                Edges[1, j] := 0;
                Edges[0, k] := 0;
                Edges[1, k] := 0;
              end;
            end;
          end;
        end;
      end;
    end;
    // Form new triangles for the current point
    // Skipping over any tagged edges.
    // All edges are arranged in clockwise order.
    for j := 1 to Nedge do
    begin
      if not(Edges[0, j] = 0) and not(Edges[1, j] = 0) then
      begin
        ntri := ntri + 1;
        Triangle[ntri].vv0 := Edges[0, j];
        Triangle[ntri].vv1 := Edges[1, j];
        Triangle[ntri].vv2 := i;
        Triangle[ntri].PreCalc := 0;
        Complete[ntri] := False;
      end;
    end;
  end;
  // Remove triangles with supertriangle vertices
  // These are triangles which have a vertex number greater than NVERT
  i := 0;
  repeat
    i := i + 1;
    if (Triangle[i].vv0 > nvert) or (Triangle[i].vv1 > nvert) or
      (Triangle[i].vv2 > nvert) then
    begin
      Triangle[i].vv0 := Triangle[ntri].vv0;
      Triangle[i].vv1 := Triangle[ntri].vv1;
      Triangle[i].vv2 := Triangle[ntri].vv2;
      i := i - 1;
      ntri := ntri - 1;
    end;
  until i >= ntri;
  Triangulate := ntri;
 
  // Free memory
  SetLength(Complete, 0);
  SetLength(Edges, 2, 0);
end;
 
procedure TGLDelaunay2D.Mesh(sort: Boolean);
begin
  if sort then
    QuickSort(Vertex, 1, tPoints - 1);
  { usingsort:=sort; }
  if tPoints > 3 then
    HowMany := Triangulate(tPoints - 1);
  // 'Returns number of triangles created.
end;
 
procedure TGLDelaunay2D.AddPoint(X, Y, Z, U, V: Single; MatIndex: Integer);
var
  i, AE: Integer;
begin
  // Check for duplicate points
  AE := 0;
  i := 1;
  while i < tPoints do
  begin
    if (Abs(X - Vertex[i].X) < ExPtTolerance) and
      (Abs(Y - Vertex[i].Y) < ExPtTolerance) then
      AE := 1;
    inc(i);
  end;
  if AE = 0 then
  begin
    // Set Vertex coordinates where you clicked the pic box
    Vertex[tPoints].X := X;
    Vertex[tPoints].Y := Y;
    Vertex[tPoints].Z := Z;
    Vertex[tPoints].U := U;
    Vertex[tPoints].V := V;
    Vertex[tPoints].MatIndex := MatIndex;
    // Increment the total number of points
    tPoints := tPoints + 1;
  end;
end;
 
procedure TGLDelaunay2D.AddPointNoCheck(X, Y, Z, U, V: Single; MatIndex: Integer);
begin
  Vertex[tPoints].X := X;
  Vertex[tPoints].Y := Y;
  Vertex[tPoints].Z := Z;
  Vertex[tPoints].U := U;
  Vertex[tPoints].V := V;
  Vertex[tPoints].MatIndex := MatIndex;
  tPoints := tPoints + 1;
end;
 
procedure TGLDelaunay2D.RemoveLastPoint;
begin
  tPoints := tPoints - 1;
end;
 
procedure TGLDelaunay2D.QuickSort(var A: TDVertex; Low, High: Integer);
// Sort all points by x
procedure DoQuickSort(var A: TDVertex; iLo, iHi: Integer);
var
  Lo, Hi: Integer;
  Mid: Single;
  T: DVertex;
begin
  Lo := iLo;
  Hi := iHi;
  Mid := A[(Lo + Hi) div 2].X;
  repeat
    while A[Lo].X < Mid do
      inc(Lo);
    while A[Hi].X > Mid do
      Dec(Hi);
    if Lo <= Hi then
    begin
      T := A[Lo];
      A[Lo] := A[Hi];
      A[Hi] := T;
      inc(Lo);
      Dec(Hi);
    end;
  until Lo > Hi;
  if Hi > iLo then
    DoQuickSort(A, iLo, Hi);
  if Lo < iHi then
    DoQuickSort(A, Lo, iHi);
end;
 
begin
  DoQuickSort(A, Low, High);
end;
 
end.