MathgeomGLS

1
unit Vor.Shamos;
2

3
interface
4

5
uses
6
  System.Classes,
7
  Vcl.Graphics,
8
  System.SysUtils,
9
  Vcl.Forms,
10
  Vor.GraphObjects;
11

12
type
13
  TVoronoi = class(TObject)
14
  private
15
    Canvas: TCanvas;
16
    function Iterate(L: TList): TList; // The result is the convex hull !
17
  public
18
    ListPoints: TList; // Global List of Points (and later Lines)
19
    constructor Create(C: TCanvas; L: TList);
20
    procedure ClearLines; // Deletes all Lines from LGlobal
21
    procedure CalcVoronoi(ShowCHull, ShowTriangulation: boolean);
22
    // Builds the Voronoi Diagram into LGlobal
23
  end;
24

25
//=========================================================
26
implementation
27
//========================================================
28

29
uses
30
  fVor2dPick;
31

32
//======================  Service functions ================
33
function GetAvgX(L: TList): extended;
34
var
35
  minx, maxx, x: extended;
36
  z: integer;
37
begin
38
  minx := 66666;
39
  maxx := -66666;
40
  for z := 0 to L.Count - 1 do
41
  begin
42
    if assigned(L.Items[z]) then
43
      if TObject(L.Items[z]) is TGPoint then
44
      begin
45
        x := TGPoint(L.Items[z]).getX;
46
        if x > maxx then
47
          maxx := x;
48
        if x < minx then
49
          minx := x;
50
      end;
51
  end;
52
  result := (maxx + minx) / 2;
53
end;
54

55
// ---------------------------------------------------------------------
56
function Next(i: integer; L: TList): integer;
57
begin
58
  inc(i);
59
  if i = L.Count then
60
    i := 0;
61
  result := i;
62
end;
63

64
function prev(i: integer; L: TList): integer;
65
begin
66
  dec(i);
67
  if i < 0 then
68
    i := L.Count - 1;
69
  result := i;
70
end;
71

72
// ---------------------------------------------------------------------
73
function SortY(i1, i2: pointer): integer;
74
begin
75
  if TGPoint(i1).getY > TGPoint(i2).getY then
76
    result := 1
77
  else
78
    result := -1;
79
end;
80

81
function SortX(i1, i2: pointer): integer;
82
begin
83
  if TGPoint(i1).getX > TGPoint(i2).getX then
84
    result := 1
85
  else
86
    result := -1;
87
end;
88

89
procedure sort(L: TList; cf: TListSortCompare);
90
// performs a bubble sort. TList.sort doesn't work for some reason.
91
var
92
  z, i: integer;
93
  p: pointer;
94
begin
95
  for z := L.Count - 2 downto 0 do
96
    for i := 0 to z do
97
    begin
98
      if cf(L.Items[i], L.Items[i + 1]) < 0 then
99
      begin
100
        p := L.Items[i];
101
        L.Items[i] := L.Items[i + 1];
102
        L.Items[i + 1] := p; // exchange items
103
        TGraphObject(L.Items[i]).ReIndex(false, i);
104
        TGraphObject(L.Items[i + 1]).ReIndex(false, i + 1);
105
      end;
106
    end;
107
end;
108

109
function getPointCount(L: TList): integer;
110
var
111
  z: integer;
112
begin
113
  result := 0;
114
  if L.Count = 0 then
115
    exit;
116
  for z := 0 to L.Count - 1 do
117
    if assigned(L.Items[z]) then
118
      if TObject(L.Items[z]) is TGPoint then
119
        inc(result);
120
end;
121

122
//============= TVoronoi ==================================
123
constructor TVoronoi.Create(C: TCanvas; L: TList);
124
begin
125
  ListPoints := L;
126
  Canvas := C;
127
end;
128

129
function TVoronoi.Iterate(L: TList): TList;
130
// integrates the voronoi diagram from a given point list L into LGlobal, and returns the convex hull
131
var
132
  iterateList1, iterateList2, temp: TList;
133
  z, i, i1, i2, points: integer;
134
  avgX, aex, aey, bex, bey: extended;
135
  ConvexHull1, ConvexHull2: TList;
136
  chain, el, er: TList; // el,er represents the voronoi polygon of pl/pr
137
  Edge: TGLine;
138
  v, pl, pr, ip1, ip2: TGPoint;
139
  mina, minb, maxa, a: extended;
140
  i1maxx, i2minx, i1max, i1min, i2max, i2min: integer;
141
  // indexes for the new convex hull
142

143
label done, abort;
144

145
{$I Service.inc}   // modularized the algorithm to keep it clear
146

147

148
// ---------------------------------------------------------------------------------------------
149
begin
150
  points := getPointCount(L);
151
  if points = 0 then
152
  begin
153
    raise ERangeError.Create('Received 0 point list in iteration!');
154
    exit;
155
  end;
156

157
  result := TList.Create;
158
  if points = 1 then
159
  begin
160
    for z := 0 to L.Count - 1 do
161
      if assigned(L.Items[z]) then
162
        if TObject(L.Items[z]) is TGPoint then
163
        begin
164
          TGPoint(L.Items[z]).MoveToList(result);
165
          exit;
166
        end;
167
    raise ERangeError.Create('Found 0 points instead of 1!');
168
    exit;
169
  end;
170
  if points = 3 then
171
  begin // collinear check, VERY incomplete.
172
    if TGPoint(L.Items[0]).areCollinear(L.Items[1], L.Items[2]) then
173
    begin
174
      if TGPoint(L.Items[0]).getx = TGPoint(L.Items[1]).getx then
175
      begin
176
        Sort(L, @SortY);
177
      end
178
      else if TGPoint(L.Items[0]).gety = TGPoint(L.Items[1]).gety then
179
      begin
180
        Sort(L, @SortX);
181
      end
182
      else
183
        Sort(L, @SortY);
184
      // integrate bisectors
185
      TGPoint(L.Items[0]).bisector(L.Items[1]).MoveToList(ListPoints);
186
      TGPoint(L.Items[1]).bisector(L.Items[2]).MoveToList(ListPoints);
187
      // build chull
188
      TGPoint(L.Items[0]).MoveToList(Result);
189
      TGPoint(L.Items[2]).MoveToList(Result);
190
      exit;
191
    end;
192
  end;
193
  iterateList1 := TList.Create;
194
  iterateList2 := TList.Create;
195
  avgX := getavgX(L);
196

197
  // divide all points into 2 seperate lists according to avgX
198
  for z := 0 to L.Count - 1 do
199
    if assigned(L.Items[z]) then
200
      if TObject(L.Items[z]) is TGPoint then
201
        if TGPoint(L.Items[z]).getx < avgX then
202
          TGPoint(L.Items[z]).MoveToList(iterateList1)
203
        else
204
          TGPoint(L.Items[z]).MoveToList(iterateList2);
205

206
  // for the case that all points have the same x value. we can't iterate an empty list
207
  if iterateList1.Count = 0 then
208
    TGPoint(iterateList2.Items[iterateList2.Count - 1])
209
      .MoveToList(iterateList1);
210
  if iterateList2.Count = 0 then
211
    TGPoint(iterateList1.Items[iterateList1.Count - 1])
212
      .MoveToList(iterateList2);
213
  // now the iteration
214
  ConvexHull1 := Iterate(iterateList1);
215
  ConvexHull2 := Iterate(iterateList2);
216

217
  // the hard part...
218
  chain := TList.Create;
219
  //
220
  // constructing new convex hull according to the Preparata-Hong algorithm
221
  //
222
  maxa := -66666;
223
  i1maxx := -1;
224
  for z := 0 to ConvexHull1.Count - 1 do
225
    if TGPoint(ConvexHull1.Items[z]).getx >= maxa then
226
    begin
227
      i1maxx := z;
228
      maxa := TGPoint(ConvexHull1.Items[z]).getx;
229
    end;
230

231
  mina := 66666;
232
  i2minx := -1;
233
  for z := 0 to ConvexHull2.Count - 1 do
234
    if TGPoint(ConvexHull2.Items[z]).getx <= mina then
235
    begin
236
      i2minx := z;
237
      mina := TGPoint(ConvexHull2.Items[z]).getx;
238
    end;
239
  // find upper supporting line
240
  i1 := i1maxx;
241
  i2 := i2minx;
242
  while (TGPoint(ConvexHull1.Items[i1]).Angle(ConvexHull2.Items[i2]) >
243
    TGPoint(ConvexHull1.Items[i1]).Angle(ConvexHull2.Items[next(i2, ConvexHull2)
244
    ])) or (TGPoint(ConvexHull2.Items[i2]).Angle(ConvexHull1.Items[i1]) <
245
    TGPoint(ConvexHull2.Items[i2]).Angle(ConvexHull1.Items[prev(i1,
246
    ConvexHull1)])) do
247
  begin
248
    while TGPoint(ConvexHull1.Items[i1]).Angle(ConvexHull2.Items[i2]) >
249
      TGPoint(ConvexHull1.Items[i1])
250
      .Angle(ConvexHull2.Items[next(i2, ConvexHull2)]) do
251
      i2 := next(i2, ConvexHull2);
252
    while TGPoint(ConvexHull2.Items[i2]).Angle(ConvexHull1.Items[i1]) <
253
      TGPoint(ConvexHull2.Items[i2])
254
      .Angle(ConvexHull1.Items[prev(i1, ConvexHull1)]) do
255
      i1 := prev(i1, ConvexHull1);
256
  end;
257
  i1max := i1;
258
  i2max := i2;
259
  // find lower supporting line
260
  i1 := i1maxx;
261
  i2 := i2minx;
262
  while (TGPoint(ConvexHull1.Items[i1]).Angle(ConvexHull2.Items[i2]) <
263
    TGPoint(ConvexHull1.Items[i1]).Angle(ConvexHull2.Items[prev(i2, ConvexHull2)
264
    ])) or (TGPoint(ConvexHull2.Items[i2]).Angle(ConvexHull1.Items[i1]) >
265
    TGPoint(ConvexHull2.Items[i2]).Angle(ConvexHull1.Items[next(i1,
266
    ConvexHull1)])) do
267
  begin
268
    while TGPoint(ConvexHull1.Items[i1]).Angle(ConvexHull2.Items[i2]) <
269
      TGPoint(ConvexHull1.Items[i1])
270
      .Angle(ConvexHull2.Items[prev(i2, ConvexHull2)]) do
271
      i2 := prev(i2, ConvexHull2);
272
    while TGPoint(ConvexHull2.Items[i2]).Angle(ConvexHull1.Items[i1]) >
273
      TGPoint(ConvexHull2.Items[i2])
274
      .Angle(ConvexHull1.Items[next(i1, ConvexHull1)]) do
275
      i1 := next(i1, ConvexHull1);
276
  end;
277
  i1min := i1;
278
  i2min := i2;
279
  //
280
  // constructing the dividing chain
281
  //
282
  el := nil;
283
  er := nil;
284
  TGPoint(ConvexHull1.Items[i1max]).bisector(ConvexHull2.Items[i2max])
285
    .MoveToList(chain); // initial ray
286
  Edge := chain.Items[0];
287
  if Edge.ey > 0 then // put v into a "conveniant large ordinate"
288
    v := TGPoint.Create(Edge.p2.getx + 6000 * Edge.ex,
289
      Edge.p2.gety + 6000 * Edge.ey, nil, nil)
290
  else
291
    v := TGPoint.Create(Edge.p2.getx - 6000 * Edge.ex, Edge.p2.gety - 6000 * Edge.ey,
292
      nil, nil);
293
  pl := ConvexHull1.Items[i1max];
294
  pr := ConvexHull2.Items[i2max];
295
  repeat
296
    GetVPoly(pl.GetOrigIndex, el);
297
    i1 := 0;
298
    GetVPoly(pr.GetOrigIndex, er);
299
    i2 := 0;
300
    if (el.Count = 0) and (er.Count = 0) then
301
    // only one bisector, can skip procedure
302
      goto done;
303
    ip1 := nil;
304
    ip2 := nil;
305
    mina := FindMinIPDistance(el, i1, ip1);
306
    minb := FindMinIPDistance(er, i2, ip2);
307
    // if mina=minb then raise ERangeError.Create('Wrong Cut');
308
    if mina = minb then
309
    begin // the bad part. some cut went the wrong way, but to keep it running we skip the procedure.
310
      FormVor2dPick.Label1.Visible := true;
311
      goto abort;
312
    end;
313
    if mina < minb then
314
    begin
315
      aex := TGLine(el.Items[i1]).ex;
316
      aey := TGLine(el.Items[i1]).ey;
317
      bex := Edge.ex;
318
      bey := Edge.ey;
319
      if aey < 0 then
320
        aex := -aex;
321
      if bey < 0 then
322
        bex := -bex;
323
      if ((aex < bex) and (aex >= 0)) or ((aex > bex) and (aex < 0)) then
324
        TGLine(el.Items[i1]).CutLeft(Edge)
325
      else
326
        TGLine(el.Items[i1]).CutRight(Edge);
327
      v.Free;
328
      v := ip1.clone as TGPoint;
329
      if TGLine(el.Items[i1]).BisectorOf[1] = pl.GetOrigIndex then
330
        pl := ListPoints.Items[TGLine(el.Items[i1]).BisectorOf[2]]
331
      else
332
        pl := ListPoints.Items[TGLine(el.Items[i1]).BisectorOf[1]];
333
      pl.bisector(pr).MoveToList(chain);
334
    end
335
    else
336
    begin
337
      aex := TGLine(er.Items[i2]).ex;
338
      aey := TGLine(er.Items[i2]).ey;
339
      bex := Edge.ex;
340
      bey := Edge.ey;
341
      if aey < 0 then
342
        aex := -aex;
343
      if bey < 0 then
344
        bex := -bex;
345
      if ((aex < bex) and (aex >= 0)) or ((aex > bex) and (aex < 0)) then
346
        TGLine(er.Items[i2]).CutRight(Edge)
347
      else
348
        TGLine(er.Items[i2]).CutLeft(Edge);
349
      v.Free;
350
      v := ip2.clone as TGPoint;
351
      if TGLine(er.Items[i2]).BisectorOf[1] = pr.GetOrigIndex then
352
        pr := ListPoints.Items[TGLine(er.Items[i2]).BisectorOf[2]]
353
      else
354
        pr := ListPoints.Items[TGLine(er.Items[i2]).BisectorOf[1]];
355
      pr.bisector(pl).MoveToList(chain);
356
    end;
357
    if ip1 <> nil then
358
      ip1.Free;
359
    if ip2 <> nil then
360
      ip2.Free;
361
    Edge.CutBoth(chain.last);
362

363
  done:
364
    Edge := chain.last;
365

366
  until ((Edge.BisectorOf[1] = TGPoint(ConvexHull1.Items[i1min]).GetOrigIndex) and
367
    (Edge.BisectorOf[2] = TGPoint(ConvexHull2.Items[i2min]).GetOrigIndex)) or
368
    ((Edge.BisectorOf[2] = TGPoint(ConvexHull1.Items[i1min]).GetOrigIndex) and
369
    (Edge.BisectorOf[1] = TGPoint(ConvexHull2.Items[i2min]).GetOrigIndex));
370

371
abort:
372
  DiscardLines; // doesn't work perfect, and is very inefficent
373
  for z := 0 to chain.Count - 1 do // move chain to LGlobal
374
    TGLine(chain.Items[z]).MoveToList(ListPoints);
375

376
  // building new chull
377
  i := i1min;
378
  TGraphObject(ConvexHull1.Items[i]).MoveToList(Result);
379
  while i <> i1max do
380
  begin
381
    i := next(i, ConvexHull1);
382
    TGraphObject(ConvexHull1.Items[i]).MoveToList(Result);
383
  end;
384

385
  i := i2max;
386
  TGraphObject(ConvexHull2.Items[i]).MoveToList(result);
387
  while i <> i2min do
388
  begin
389
    i := next(i, ConvexHull2);
390
    TGraphObject(ConvexHull2.Items[i]).MoveToList(result);
391
  end;
392

393
  // finished
394

395
  chain.Free;
396

397
  ConvexHull1.Free;
398
  ConvexHull2.Free;
399

400
  iterateList2.Free;
401
  iterateList1.Free;
402

403
end;
404

405
procedure TVoronoi.ClearLines; // Deletes all Lines from LGlobal
406
var
407
  z: integer;
408
begin
409
  for z := 0 to ListPoints.Count - 1 do
410
    if assigned(ListPoints.Items[z]) then
411
      if TObject(ListPoints.Items[z]) is TGLine then
412
      begin
413
        TGLine(ListPoints.Items[z]).Free;
414
        ListPoints.Items[z] := nil;
415
      end;
416

417
  ListPoints.Pack;
418
  for z := 0 to ListPoints.Count - 1 do
419
    if TObject(ListPoints.Items[z]) is TGPoint then
420
      TGPoint(ListPoints.Items[z]).ReIndex(true);
421

422
end;
423

424
procedure TVoronoi.CalcVoronoi; // Builds the Voronoi Diagram into ListPoints
425
var
426
  iterateList, chull: TList;
427
  L: TGLine;
428
  z: integer;
429
begin
430
  FormVor2dPick.Label1.Visible := false;
431
  Application.ProcessMessages;
432

433
  // must clone all points out of LGlobal because all points will be moved and killed during iteration
434
  if GetPointCount(ListPoints) < 2 then
435
    exit;
436

437
  IterateList := TList.Create;
438
  for z := 0 to ListPoints.Count - 1 do
439
    if assigned(ListPoints.Items[z]) then
440
      if TObject(ListPoints.Items[z]) is TGPoint then
441
        TGPoint(ListPoints.Items[z]).CloneToList(iterateList);
442

443
  // -------------------------------------------------
444
  chull := Iterate(iterateList);
445
  // -------------------------------------------------
446

447
  if ShowTriangulation then
448
  begin
449
    for z := 0 to ListPoints.Count - 1 do
450
      if TObject(ListPoints.Items[z]) is TGLine then
451
      begin
452
        L := ListPoints.Items[z];
453
        if L.BisectorOf[1] > -1 then
454
        begin
455
          TGLine.Create(ListPoints.Items[L.BisectorOf[1]],
456
                        ListPoints.Items[L.BisectorOf[2]],
457
                        TwoPoint, ListPoints, Canvas).Color.Color := clGreen;
458
        end;
459
      end;
460
  end;
461

462
  if ShowCHull then
463
  begin
464
    for z := 0 to chull.Count - 2 do
465
    begin
466
      TGLine.Create(chull.Items[z], chull.Items[z + 1],
467
        TwoPoint, ListPoints, Canvas).Color.Color := clYellow;
468
    end;
469
    TGLine.Create(chull.Items[chull.Count - 1], chull.Items[0], TwoPoint,
470
      ListPoints, Canvas).Color.Color := clYellow;
471
  end;
472

473
  IterateList.Free;
474
end;
475

476
end.
477