Delaunay-Triangulation

[Bjoerk]

Jens, kannst du mir sagen was der Autor hier macht?

Delphi-Quellcode:

Function TDelaunay.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: PComplete;
  Edges: PEdges;
  Nedge: integer;

  //For Super Triangle
  xmin: double;
  xmax: double;
  ymin: double;
  ymax: double;
  xmid: double;
  ymid: double;
  dx: double;
  dy: double;
  dmax: double;

  //General Variables
  i: integer;
  j: integer;
  k: integer;
  ntri: integer;
  xc: double;
  yc: double;
  r: double;
  inc: Boolean;
begin
  //Allocate memory
  GetMem(Complete, sizeof(Complete^));
  GetMem(Edges, sizeof(Edges^));

  //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
    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 (xc + r) < Vertex[i].X then
          complete[j] := True
        Else
          if inc then
          begin
            Edges^[1, Nedge + 1] := Triangle^[j].vv0;
            Edges^[2, Nedge + 1] := Triangle^[j].vv1;
            Edges^[1, Nedge + 2] := Triangle^[j].vv1;
            Edges^[2, Nedge + 2] := Triangle^[j].vv2;
            Edges^[1, Nedge + 3] := Triangle^[j].vv2;
            Edges^[2, 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^[1, j] = 0) And Not (Edges^[2, j] = 0) then
      begin
        For k := j + 1 To Nedge do
        begin
          if Not (Edges^[1, k] = 0) And Not (Edges^[2, k] = 0) then
          begin
            if Edges^[1, j] = Edges^[2, k] then
            begin
              if Edges^[2, j] = Edges^[1, k] then
              begin
                Edges^[1, j] := 0;
                Edges^[2, j] := 0;
                Edges^[1, k] := 0;
                Edges^[2, 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^[1, j] = 0) And Not (Edges^[2, j] = 0) then
      begin
        ntri := ntri + 1;
        Triangle^[ntri].vv0 := Edges^[1, j];
        Triangle^[ntri].vv1 := Edges^[2, 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
  FreeMem(Complete, sizeof(Complete^));
  FreeMem(Edges, sizeof(Edges^));
end;