LZScene

1
//
2
// This unit is part of the GLScene Engine https://github.com/glscene
3
//
4
{
5
  Class and routines to output isolines.
6

7
   History :  
8
  25/04/15 - PW - Fixed TriangleElevationSegments procedure
9
  06/07/02 - Phil Scadden - Added TContour class and Initialise_Isolining procedure
10
  15/08/01 - Alexander Weidauer - Added CONREC Delphi implementation
11
             (based on Nicholas Yue CONREC.C and  Paul D. Bourke CONREC.F)
12
  15/07/01 - PW - Creation of the unit
13
}
14
unit GLIsolines;
15

16
interface
17

18
uses
19
  System.SysUtils, System.Classes, System.Math,
20
   
21
  GLVectorGeometry, GLVectorLists, GLTypes, GLSpline;
22

23
{$I GLScene.inc}  
24
  
25
type
26
  TVectorL4D = array [0 .. 4] of Single;
27
  TVectorL4I = array [0 .. 4] of Integer;
28
  TCastArray = array [0 .. 2, 0 .. 2, 0 .. 2] of Integer;
29

30
  TGLIsoline2D = array [0 .. 32767] of TGLPoint2D;
31
  PGLIsoline2D = ^TGLIsoline2D;
32

33
  TIsoline = class
34
    NP: integer;
35
    Line: PGLIsoline2D;
36
    constructor Create(LineSize: integer); virtual;
37
    destructor Destroy; override;
38
  end;
39
  PIsoline = ^TIsoline;
40

41
  TIsolineState = (ilsEmpty, ilsCalculating, ilsReady);
42

43
  TIsolines = class
44
    CoordRange: Integer;
45
    LineList: TList;
46
    IsolineState: TIsolineState;
47
    procedure MakeIsolines(var Depths: TGLMatrix; bmSize: Integer;
48
      StartDepth, EndDepth: Single; Interval: Integer);
49
    procedure FreeList;
50
    constructor Create; virtual;
51
    destructor Destroy; override;
52
  end;
53

54
procedure Initialize_Isolining(var DataGrid: TGLMatrix;
55
  NXpoints, NYpoints: integer; CurrentIsoline: Single);
56
procedure Release_Memory_Isolining;
57
function GetNextIsoline(var Isoline: TIsoline): Boolean;
58

59
{ : Defines contouring segments inside a triangle using elevations }
60
procedure TriangleElevationSegments(const p1, p2, p3: TAffineVector;
61
  ElevationDelta: Single; Segments: TAffineVectorList);
62

63
{ : CONREC is a contouring routine for rectangular spaced data or regular 2D grids
64
  It takes each rectangle of adjacent data points and splits it
65
  into 4 triangles after choosing the height at the centre of the rectangle.
66
  For each of the triangles the line segment resulting from the intersection
67
  with each contour plane.
68
  A routine is then called with the starting and stopping coordinates
69
  of the line segment that make up a contour curve and then output these
70
  isolines. See details in http://paulbourke.net/papers/conrec/
71

72
  The input parameters are as follows :
73
  Data -  Scalar field in 2D grid
74
  ilb - lower bound in west - east direction
75
  iub - upper bound in west - east direction
76
  jlb - lower bound in north - south direction
77
  jub upper bound in north - south direction
78
  X - coord. vector for west - east
79
  Y - coord. vector for north - south
80
  nc - number of cut levels
81
  Z - values of cut levels
82
}
83
procedure Conrec(Data: TGLMatrix; ilb, iub, jlb, jub: Integer;
84
  x: TGLVector; y: TGLVector; nc: Integer; z: TGLVector; Isoline: TIsoline);
85

86

87
//----------------------------------------------------------------------
88
//----------------------------------------------------------------------
89
//----------------------------------------------------------------------
90
implementation
91
//----------------------------------------------------------------------
92
//----------------------------------------------------------------------
93
//----------------------------------------------------------------------
94

95
var
96
  ii, jj: integer;
97
  Visited: TGLByteMatrix; // 0 = not visited
98
  // 1 = visited
99
  // if it is a saddle points, then bits 1-4 also encode
100
  // which exit and entry points were used.
101
  Grid: TGLMatrix;
102
  NX, NY: integer;
103
  LineX1, LineY1, LineX2, LineY2: TGLVector;
104

105
function EqAdd(a, b: integer): integer;
106
begin
107
  if a = b then
108
    Result := 1
109
  else
110
    Result := 0;
111
end;
112

113
// Initialize_Isolining
114
//
115
procedure Initialize_Isolining;
116

117
var
118
  i, j: Integer;
119
  maxnp: Integer;
120
begin
121
  ii := 1;
122
  jj := 1;
123
  NX := NXpoints;
124
  NY := NYpoints;
125
  maxnp := NX * NY div 256;
126
  SetLength(Visited, NX, NY);
127
  for i := 0 to NX - 1 do
128
    for j := 0 to NY - 1 do
129
      Visited[i, j] := 0;
130
  SetLength(Grid, NX + 1, NY + 1);
131
  SetLength(LineX1, maxnp);
132
  SetLength(LineY1, maxnp);
133
  SetLength(LineX2, maxnp);
134
  SetLength(LineY2, maxnp);
135
  // Generate a grid of data relative to the Isoline level
136
  for i := 1 to NX do
137
  begin
138
    for j := 1 to NY do
139
    begin
140
      Grid[i][j] := DataGrid[i - 1][j - 1] - CurrentIsoline;
141
      (* Don't want any grid points exactly zero *)
142
      if Grid[i][j] = 0 then
143
      begin
144
        Grid[i][j] := 1E-8;
145
      end;
146
    end;
147
  end;
148
end;
149

150
// Release_Memory_Isolining
151
//
152
procedure Release_Memory_Isolining;
153
begin
154
  SetLength(Visited, 0);
155
  SetLength(Grid, 0);
156
  SetLength(LineX1, 0);
157
  SetLength(LineY1, 0);
158
  SetLength(LineX2, 0);
159
  SetLength(LineY2, 0);
160
end;
161

162
// Cuts
163
//
164
procedure Cuts(const g: TGLMatrix; i, j: integer; var s: array of integer);
165
begin
166
  s[0] := 0;
167
  if g[i][j + 1] * g[i + 1][j + 1] < 0 then
168
  begin
169
    Inc(s[0]);
170
    s[s[0]] := 1;
171
  end;
172
  if g[i + 1][j + 1] * g[i + 1][j] < 0 then
173
  begin
174
    Inc(s[0]);
175
    s[s[0]] := 2;
176
  end;
177
  if g[i + 1][j] * g[i][j] < 0 then
178
  begin
179
    Inc(s[0]);
180
    s[s[0]] := 3;
181
  end;
182
  if g[i][j] * g[i][j + 1] < 0 then
183
  begin
184
    Inc(s[0]);
185
    s[s[0]] := 4;
186
  end;
187
end;
188

189
// Intercept
190
//
191
procedure Intercept(const g: TGLMatrix; i, j, s: Integer;
192
  var x, y: single);
193
begin
194
  case s of
195
    1:
196
      begin
197
        x := abs(g[i][j + 1] / (g[i + 1][j + 1] - g[i][j + 1])) + i;
198
        y := 1 + j;
199
      end;
200
    2:
201
      begin
202
        y := abs(g[i + 1][j] / (g[i + 1][j + 1] - g[i + 1][j])) + j;
203
        x := 1 + i;
204
      end;
205
    3:
206
      begin
207
        x := abs(g[i][j] / (g[i + 1][j] - g[i][j])) + i;
208
        y := j;
209
      end;
210
    4:
211
      begin
212
        y := abs(g[i][j] / (g[i][j + 1] - g[i][j])) + j;
213
        x := i;
214
      end;
215
  end;
216
end;
217

218
// Free_Exit
219
//
220
function Free_Exit(const Visited: TGLByteMatrix;
221
  i, j, NX, NY, Lexit: Integer): Boolean;
222
var
223
  ni, nj: Integer;
224
  entry: Integer;
225
begin
226
  nj := j + EqAdd(lexit, 1) - EqAdd(lexit, 3);
227
  ni := i + EqAdd(lexit, 2) - EqAdd(lexit, 4);
228
  if (ni < 1) or (ni >= NX) or (nj < 1) or (nj >= NY) or (Visited[ni][nj] = 0)
229
  then
230
    Result := True // can always exit on an edge
231
  else
232
  begin
233
    entry := ((lexit + 1) mod 4) + 1;
234
    Result := (((Visited[ni][nj] shr entry) and 1) = 0);
235
  end;
236
end;
237

238
// TraceIsoline
239
//
240
procedure TraceIsoline(i, j, Lexit, NX, NY: Integer; const Grid: TGLMatrix;
241
  const Visited: TGLByteMatrix; var LineX, LineY: TGLVector;
242
  var NP: Integer; var OffGrid: Boolean);
243
var
244
  ni, nj, si, sj: Integer;
245
  p, q: Integer;
246
  s: array [0 .. 5] of Integer;
247
  entry: integer;
248

249
begin
250
  ni := i;
251
  nj := j;
252
  si := i;
253
  sj := j;
254
  NP := 0;
255
  offgrid := False;
256
  Visited[i][j] := 1;
257
  Intercept(Grid, i, j, lexit, LineX[NP], LineY[NP]);
258
  NP := 1;
259
  while True do
260
  begin
261
    nj := nj + EqAdd(lexit, 1) - EqAdd(lexit, 3);
262
    ni := ni + EqAdd(lexit, 2) - EqAdd(lexit, 4);
263
    if (ni < 1) or (ni > NX - 1) or (nj < 1) or (nj > NY - 1) then
264
    begin
265
      offgrid := True;
266
      break;
267
    end;
268
    Visited[ni][nj] := 1;
269
    entry := ((lexit + 1) mod 4) + 1;
270
    Cuts(Grid, ni, nj, s);
271
    // Have come to a new point on the Isoline
272
    lexit := 0;
273
    if (s[0] = 2) then
274
    begin
275
      // If there are two cuts then choose the one that is not the entry
276
      if entry = s[1] then
277
        lexit := s[2]
278
      else
279
        lexit := s[1];
280
    end
281
    else
282
    begin
283
      // If there are four cuts (saddle) then work round from the left
284
      p := (entry mod 4) + 1;
285
      while p <> entry do
286
      begin
287
        for q := 1 to s[0] do
288
        begin
289
          if (s[q] = p) and Free_Exit(Visited, NX, NY, ni, nj, p) then
290
          begin
291
            lexit := p;
292
            break;
293
          end;
294
        end;
295
        // Aim is to find first
296
        if lexit <> 0 then
297
          break;
298
        p := (p mod 4) + 1;
299
      end;
300
      (* exit from cell, going *)
301
      // We found a way out, make a note of way in and way out.
302
      // Need to do this as saddle may be visited twice.
303
      Visited[ni][nj] := (Visited[ni][nj]) or (1 shl entry) or (1 shl lexit);
304
    end;
305
    // clockwise from entry point
306
    Assert(lexit > 0, 'Contour routine caught in a loop');
307
    if (lexit = 0) then
308
      break;
309
    Intercept(Grid, ni, nj, lexit, LineX[NP], LineY[NP]);
310
    Inc(NP);
311
    if (ni = si) and (nj = sj) then
312
      break;
313
  end;
314
  // Have finished loop
315
end;
316

317
{ LineX and LineY are (pointers to) zero-offset vectors, to which
318
  sufficient space has been allocated to store the coordinates of
319
  any feasible Isoline }
320
//
321
function GetNextIsoline(var Isoline: TIsoline): Boolean;
322
var
323
  OffGrid: boolean;
324
  lexit: integer;
325
  np1, np2: integer;
326
  i, j, k: integer;
327
  s: array [0 .. 4] of integer;
328
begin
329
  for i := ii to NX - 1 do
330
  begin
331
    for j := 1 + (jj - 1) * EqAdd(i, ii) to NY - 1 do
332
    begin
333
      if (Visited[i][j] = 0) then
334
      begin
335
        Cuts(Grid, i, j, s);
336
        if s[0] = 2 then
337
        begin
338
          Lexit := s[2];
339
          TraceIsoline(i, j, lexit, NX, NY, Grid, Visited, LineX1, LineY1,
340
            np1, offgrid);
341
          // Follow the Isoline along
342
          if offgrid then
343
          begin
344
            // Go back to start of Isoline and trace in opposite direction
345
            Lexit := s[1];
346
            TraceIsoline(i, j, Lexit, NX, NY, Grid, Visited, LineX2, LineY2,
347
              np2, offgrid);
348
            // Copy both bits of line into Isoline
349
            Isoline := TIsoline.Create(np1 + np2);
350
            for k := 0 to np2 - 1 do
351
            begin
352
              Isoline.Line^[k].x := LineX2[np2 - k - 1];
353
              Isoline.Line^[k].y := LineY2[np2 - k - 1];
354
            end;
355
            for k := 0 to np1 - 1 do
356
            begin
357
              Isoline.Line^[k + np2].x := LineX1[k];
358
              Isoline.Line^[k + np2].y := LineY1[k];
359
            end;
360
          end
361
          else
362
          begin
363
            // Just copy the single Isoline loop into LineX & LineY
364
            Isoline := TIsoline.Create(np1);
365
            for k := 0 to np1 - 1 do
366
            begin
367
              Isoline.Line^[k].x := LineX1[k];
368
              Isoline.Line^[k].y := LineY1[k];
369
            end;
370
          end;
371
          // scale Isoline into true units
372
          { for k:=1 to np do
373
            begin
374
            LineX[k-1]:= xlo+(LineX[k]-1)*(xhi-xlo) / (nx-1);
375
            LineY[k-1]:= ylo+(LineY[k]-1)*(yhi-ylo) / (ny-1);
376
            // LineX and LineY are zero offset vectors
377
            end; }
378
          ii := i;
379
          jj := j;
380
          Result := True;
381
          Exit;
382
        end;
383
      end;
384
    end;
385
  end;
386
  Result := False;
387
end;
388

389
// TriangleElevationSegments
390
//
391
procedure TriangleElevationSegments(const p1, p2, p3: TAffineVector;
392
  ElevationDelta: Single; Segments: TAffineVectorList);
393

394
  function SegmentIntersect(const a, b: TAffineVector; e: Single): Integer;
395
  var
396
    f: Single;
397
  begin
398
    if a.Z < b.Z then
399
    begin
400
      if (e >= a.Z) and (e < b.Z) then
401
      begin
402
        f := (e - a.Z) / (b.Z - a.Z);
403
        Segments.Add(VectorLerp(a, b, f));
404
        Result := 1;
405
      end
406
      else
407
        Result := 0;
408
    end
409
    else if a.Z > b.Z then
410
    begin
411
      if (e > b.Z) and (e <= a.Z) then
412
      begin
413
        f := (e - b.Z) / (a.Z - b.Z);
414
        Segments.Add(VectorLerp(b, a, f));
415
        Result := 1;
416
      end
417
      else
418
        Result := 0;
419
    end
420
    else
421
      Result := 0;
422
  end;
423

424
var
425
  i, n, LowElev, HighElev: Integer;
426
  e: Single;
427

428
begin
429
  LowElev := Round(MinFloat(p1.Z, p2.Z, p3.Z) / ElevationDelta - 0.5);
430
  HighElev := Round(MaxFloat(p1.Z, p2.Z, p3.Z) / ElevationDelta + 0.5);
431
  for i := LowElev to HighElev do
432
  begin
433
    e := i * ElevationDelta + 0.1;
434
    // add a real offset - this avoids all the special cases
435
    n := SegmentIntersect(p1, p2, e);
436
    if n < 2 then
437
      n := n + SegmentIntersect(p2, p3, e);
438
    if n < 2 then
439
      n := n + SegmentIntersect(p3, p1, e);
440
    Assert((n = 2) or (n = 0));
441
  end;
442
end;
443

444
// TIsolines class
445
//
446
constructor TIsolines.Create;
447
begin
448
  inherited;
449
  LineList := TList.Create;
450
  IsolineState := ilsEmpty;
451
end;
452

453
destructor TIsolines.Destroy;
454
begin
455
  FreeList;
456
  inherited;
457
end;
458

459
procedure TIsolines.FreeList;
460
var
461
  i: integer;
462
begin
463
  for i := LineList.Count - 1 downto 0 do
464
  begin
465
    with TIsoline(LineList.Items[i]) do
466
      Free;
467
  end;
468
  LineList.Clear;
469
  IsolineState := ilsEmpty;
470
end;
471

472
procedure TIsolines.MakeIsolines(var Depths: TGLMatrix; bmSize: integer;
473
  StartDepth, EndDepth: Single; Interval: Integer);
474
var
475
  Isoline: TIsoline;
476

477
begin
478
  IsolineState := ilsCalculating;
479
  CoordRange := bmSize;
480
  FreeList;
481
  repeat
482
    Initialize_Isolining(Depths, bmSize, bmSize, StartDepth);
483
    while GetNextIsoline(Isoline) do
484
    begin
485
      LineList.Add(Isoline);
486
    end;
487
    Release_Memory_Isolining;
488
    StartDepth := StartDepth + Interval;
489
  until StartDepth > EndDepth;
490
  IsolineState := ilsReady;
491
end;
492

493
constructor TIsoline.Create(LineSize: integer);
494
begin
495
  inherited Create;
496
  NP := LineSize;
497
  Getmem(Line, NP * 2 * Sizeof(Single));
498
end;
499

500
destructor TIsoline.Destroy;
501
begin
502
  inherited;
503
  if Assigned(Line) then
504
    Freemem(Line);
505
  NP := 0;
506
end;
507

508
// Conrec
509
//
510
procedure Conrec(Data: TGLMatrix; ilb, iub, jlb, jub: Integer;
511
  x: TGLVector; y: TGLVector;  nc: Integer; z: TGLVector; Isoline: TIsoline);
512
// ------------------------------------------------------------------------------
513
const
514
  im: array [0 .. 3] of Integer = (0, 1, 1, 0); // coord. cast array west - east
515
  jm: array [0 .. 3] of Integer = (0, 0, 1, 1);
516
  // coord. cast array north - south
517
  // ------------------------------------------------------------------------------
518
var
519
  m1, m2, m3,
520
  Deside: Integer;
521
  dmin, dmax, x1, x2, y1, y2: Single;
522
  lcnt, i, j, k, m: Integer;
523
  CastTab: TCastArray;
524
  h: TVectorL4D;
525
  sh: TVectorL4I;
526
  xh, yh: TVectorL4D;
527
  temp1, temp2: Single;
528
  r: Byte;
529

530
  // ------- service xsec west east lin. interpol -------------------------------
531
  function Xsec(p1, p2: Integer): Single;
532
  begin
533
    Result := (h[p2] * xh[p1] - h[p1] * xh[p2]) / (h[p2] - h[p1]);
534
  end;
535

536
// ------- service ysec north south lin interpol -------------------------------
537
  function Ysec(p1, p2: Integer): Single;
538
  begin
539
    Result := (h[p2] * yh[p1] - h[p1] * yh[p2]) / (h[p2] - h[p1]);
540
  end;
541

542
begin
543
  // set casting array
544
  CastTab[0, 0, 0] := 0;
545
  CastTab[0, 0, 1] := 0;
546
  CastTab[0, 0, 2] := 8;
547
  CastTab[0, 1, 0] := 0;
548
  CastTab[0, 1, 1] := 2;
549
  CastTab[0, 1, 2] := 5;
550
  CastTab[0, 2, 0] := 7;
551
  CastTab[0, 2, 1] := 6;
552
  CastTab[0, 2, 2] := 9;
553

554
  CastTab[1, 0, 0] := 0;
555
  CastTab[1, 0, 1] := 3;
556
  CastTab[1, 0, 2] := 4;
557
  CastTab[1, 1, 0] := 1;
558
  CastTab[1, 1, 1] := 3;
559
  CastTab[1, 1, 2] := 1;
560
  CastTab[1, 2, 0] := 4;
561
  CastTab[1, 2, 1] := 3;
562
  CastTab[1, 2, 2] := 0;
563

564
  CastTab[2, 0, 0] := 9;
565
  CastTab[2, 0, 1] := 6;
566
  CastTab[2, 0, 2] := 7;
567
  CastTab[2, 1, 0] := 5;
568
  CastTab[2, 1, 1] := 2;
569
  CastTab[2, 1, 2] := 0;
570
  CastTab[2, 2, 0] := 8;
571
  CastTab[2, 2, 1] := 0;
572
  CastTab[2, 2, 2] := 0;
573

574
  // set line counter
575
  lcnt := 0;
576
  // -----------------------------------------------------------------------------
577
  for j := jub - 1 downto jlb do
578
  begin // over all north - south and              +for j
579
    for i := ilb to iub - 1 do
580
    begin // east - west coordinates of datafield    +for i
581
      // set casting bounds from array
582
      temp1 := Min(Data[i, j], Data[i, j + 1]);
583
      temp2 := Min(Data[i + 1, j], Data[i + 1, j + 1]);
584
      dmin := Min(temp1, temp2);
585
      temp1 := Max(Data[i, j], Data[i, j + 1]);
586
      temp2 := Max(Data[i + 1, j], Data[i + 1, j + 1]);
587
      dmax := Max(temp1, temp2);
588
      if (dmax >= z[0]) and (dmin <= z[nc - 1]) then
589
      begin // ask horizontal cut avail.    +If dmin && dmax in z[0] .. z[nc-1]
590
        for k := 0 to nc - 1 do
591
        begin // over all possible cuts ---- +for k
592
          if (z[k] > dmin) and (z[k] <= dmax) then
593
          begin // aks for cut intervall ----- +If z[k] in dmin .. dmax
594
            // -----------------------------------------------------------------------
595
            for m := 4 downto 0 do
596
            begin // deteriening the cut casts and set the ---- +for m
597
              if (m > 0) then
598
              begin // height and coordinate vectors
599
                h[m] := Data[i + im[m - 1], j + jm[m - 1]] - z[k];
600
                xh[m] := x[i + im[m - 1]];
601
                yh[m] := y[j + jm[m - 1]];
602
              end
603
              else
604
              begin
605
                h[0] := (h[1] + h[2] + h[3] + h[4]) / 4;
606
                xh[0] := (x[i] + x[i + 1]) / 2;
607
                yh[0] := (y[j] + y[j + 1]) / 2;
608
              end; // if m>0 then else
609
              if h[m] > 0 then
610
                sh[m] := 1
611
              else If h[m] < 0 then
612
                sh[m] := -1
613
              else
614
                sh[m] := 0;
615
            end; // --- -for m
616

617
            // -----------------------------------------------------------
618
            for m := 1 to 4 do
619
            begin // set directional CastTable
620
              //
621
              // Note: at this stage the relative heights of the corners and the
622
              // centre are in the h array, and the corresponding coordinates are
623
              // in the xh and yh arrays. The centre of the box is indexed by 0
624
              // and the 4 corners by 1 to 4 as shown below.
625
              // Each triangle is then indexed by the parameter m, and the 3
626
              // vertices of each triangle are indexed by parameters m1,m2,and
627
              // m3.
628
              // It is assumed that the centre of the box is always vertex 2
629
              // though this isimportant only when all 3 vertices lie exactly on
630
              // the same contour level, in which case only the side of the box
631
              // is drawn.
632
              //
633
              // AS ANY BODY NOWS IST FROM THE ORIGINAL
634
              //
635
              // vertex 4 +-------------------+ vertex 3
636
              // | \               / |
637
              // |   \    m-3    /   |
638
              // |     \       /     |
639
              // |       \   /       |
640
              // |  m=2    X   m=2   |       the centre is vertex 0
641
              // |       /   \       |
642
              // |     /       \     |
643
              // |   /    m=1    \   |
644
              // | /               \ |
645
              // vertex 1 +-------------------+ vertex 2
646
              //
647
              //
648
              // Scan each triangle in the box
649
              //
650
              m1 := m;
651
              m2 := 0;
652
              if not(m = 4) then
653
                m3 := m + 1
654
              else
655
                m3 := 1;
656
              Deside := CastTab[sh[m1] + 1, sh[m2] + 1, sh[m3] + 1];
657
              if not(Deside = 0) then
658
              begin // ask is there a desition available ---+if if not(Deside=0)
659
                case Deside of
660
                // ---- determin the by desided cast cuts ---- +Case deside;
661
                  1:
662
                    begin
663
                      x1 := xh[m1];
664
                      y1 := yh[m1];
665
                      x2 := xh[m2];
666
                      y2 := yh[m2];
667
                    end;
668
                  2:
669
                    begin
670
                      x1 := xh[m2];
671
                      y1 := yh[m2];
672
                      x2 := xh[m3];
673
                      y2 := yh[m3];
674
                    end;
675
                  3:
676
                    begin
677
                      x1 := xh[m3];
678
                      y1 := yh[m3];
679
                      x2 := xh[m1];
680
                      y2 := yh[m1];
681
                    end;
682
                  4:
683
                    begin
684
                      x1 := xh[m1];
685
                      y1 := yh[m1];
686
                      x2 := xsec(m2, m3);
687
                      y2 := ysec(m2, m3);
688
                    end;
689
                  5:
690
                    begin
691
                      x1 := xh[m2];
692
                      y1 := yh[m2];
693
                      x2 := xsec(m3, m1);
694
                      y2 := ysec(m3, m1);
695
                    end;
696
                  6:
697
                    begin
698
                      x1 := xh[m3];
699
                      y1 := yh[m3];
700
                      x2 := Xsec(m1, m2);
701
                      y2 := Ysec(m1, m2);
702
                    end;
703
                  7:
704
                    begin
705
                      x1 := Xsec(m1, m2);
706
                      y1 := Ysec(m1, m2);
707
                      x2 := Xsec(m2, m3);
708
                      y2 := Ysec(m2, m3);
709
                    end;
710
                  8:
711
                    begin
712
                      x1 := Xsec(m2, m3);
713
                      y1 := Ysec(m2, m3);
714
                      x2 := Xsec(m3, m1);
715
                      y2 := Ysec(m3, m1);
716
                    end;
717
                  9:
718
                    begin
719
                      x1 := Xsec(m3, m1);
720
                      y1 := Ysec(m3, m1);
721
                      x2 := Xsec(m1, m2);
722
                      y2 := Ysec(m1, m2);
723
                    end;
724
                end; // ---  -Case deside;
725
                // -------Output of results ---------------------
726
                GetNextIsoline(Isoline);
727
                //Writeln(Format('%2.2f %2.2f %2.2f %2.2f %2.2f', [z[k], x1, y1, x2, y2]));
728
                // ---------------------------------------------------------
729
              end; // ----  -if Not(deside=0)
730
            end; // ----  -for m
731
          end; // ----  -if z[k] in dmin .. dmax
732
        end; // ----  -for k
733
      end; // ----  -if dmin && dmax in z[0] .. z[nc-1]
734
    end; // ----  -for i
735
  end; // ----  -for j
736
end;
737
// ------ End of ----------------------------------------------------------------
738

739
end.
740