MODULE Delaunay;

IMPORT QuadEdge;
FROM QuadEdge IMPORT 
  MakeEdge, Splice,
  Sym, Onext, Oprev, Lnext, Rprev,
  Odata, Ddata, SetOdata, SetDdata;

PROCEDURE Org (e: Arc): REF Site =
  BEGIN
    RETURN Odata(e)
  END Org;
  
PROCEDURE Dest (e: Arc): REF Site =
  BEGIN
    RETURN Ddata(e)
  END Dest;

PROCEDURE Build (READONLY site: ARRAY OF Site): Arc =
  BEGIN
    WITH ss = NEW(REF ARRAY OF REF Site, NUMBER(site))^ DO
      FOR i := FIRST(site) TO LAST(site) DO
        ss[i] := NEW(REF Site); ss[i]^ := site[i]
      END;
      SortSites(ss);
      WITH d = RecDelaunay(ss) DO
        RETURN d.left
      END
    END
  END Build;

PROCEDURE SortSites(VAR site: ARRAY OF REF Site) =
  (* Sorts the "site" array by increasing "x", breaking ties by "y". *)
  VAR step, j: INTEGER; 
      tmp: REF Site;
  BEGIN
    step := NUMBER(site) DIV 2;
    WHILE step > 0 DO 
      FOR i := step TO LAST(site) DO 
        j := i - step;
        WHILE j >= 0 AND 
        ( ( site[j].x > site[j+step].x ) OR 
          ( site[j].x = site[j+step].x  AND site[j].y > site[j+step].y)
        ) DO 
          tmp := site[j];
          site[j] := site[j+step]; 
          site[j+step] := tmp;
          j := j - step
        END;
      END;
      step := step DIV 2
    END;
  END SortSites;

TYPE ArcPair = RECORD left, right: Arc END;

PROCEDURE RecDelaunay(READONLY site: ARRAY OF REF Site): ArcPair =
  BEGIN
    IF NUMBER(site) = 2 THEN
      WITH a = MakeEdge() DO
	SetOdata(a, site[0]);
        SetDdata(a, site[1]);
        RETURN ArcPair{left := a, right := Sym(a)}
      END
    ELSIF NUMBER(site) = 3 THEN
      WITH
        a = MakeEdge(),
	b = MakeEdge(),
	ct = TriSign(site[0]^, site[1]^, site[2]^)
      DO
	Splice(Sym(a), b);
	SetOdata(a, site[0]); SetDdata(a, site[1]);
	SetOdata(b, site[1]); SetDdata(b, site[2]);
	IF ct = 0 THEN 
	  RETURN ArcPair{left := a, right := Sym(b)}
	ELSE
          WITH c = Connect(b, a) DO
	    IF ct > 0 THEN
	      RETURN ArcPair{left := a, right := Sym(b)}
	    ELSE 
	      RETURN ArcPair{left := Sym(c), right := c}
	    END
          END
        END
      END
    ELSE
      WITH
        n2 = NUMBER(site) DIV 2,
        ld = RecDelaunay(SUBARRAY(site, 0, n2)),
        rd = RecDelaunay(SUBARRAY(site, n2, NUMBER(site) - n2))
      DO
	VAR ldo: Arc := ld.left;
            ldi: Arc := ld.right;
            rdi: Arc := rd.left;
            rdo: Arc := rd.right;
            basel, lcand, rcand: Arc;
        BEGIN
          LOOP
            IF LeftOf(Org(rdi), ldi) THEN ldi := Lnext(ldi) 
            ELSIF RightOf(Org(ldi), rdi) THEN rdi := Rprev(rdi)
            ELSE EXIT
            END
          END;

          basel := Connect(Sym(rdi), ldi);
          IF Org(ldi) = Org(ldo) THEN ldo := Sym(basel) END;
          IF Org(rdi) = Org(rdo) THEN rdo := basel END;

          LOOP
            lcand := Onext(Sym(basel));
            IF RightOf(Dest(lcand), basel) THEN
              WHILE 
                InCircle(Dest(basel)^, Org(basel)^, Dest(lcand)^, Dest(Onext(lcand))^)
              DO 
                VAR t := Onext(lcand); BEGIN Disconnect(lcand); lcand := t END
              END;
            END;
            
            rcand := Oprev(basel);
            IF RightOf(Dest(rcand), basel) THEN
              WHILE 
                InCircle(Dest(basel)^, Org(basel)^, Dest(rcand)^, Dest(Oprev(rcand))^)
              DO
                VAR t := Oprev(rcand); BEGIN Disconnect(rcand); rcand := t END;
              END
            END;
            
            IF NOT RightOf(Dest(lcand), basel)
            AND NOT RightOf(Dest(rcand), basel) THEN
              EXIT
            END;

            IF NOT RightOf(Dest(lcand), basel) OR 
              ( RightOf(Dest(rcand), basel) AND
                InCircle(Dest(lcand)^, Org(lcand)^, Org(rcand)^, Dest(rcand)^)
              )
            THEN
              basel := Connect(rcand, Sym(basel))
            ELSE
              basel := Connect(Sym(basel), Sym(lcand));
            END;
          END;
          RETURN ArcPair{left := ldo, right := rdo}
        END
      END
    END
  END RecDelaunay;

PROCEDURE Connect (a, b: Arc): Arc =
  (* Connects Dest(a) to Org(b) across the face Left(a) = Left(b). *)
  VAR e: Arc := MakeEdge();
  BEGIN
    SetOdata(e, Dest(a));
    SetDdata(e, Org(b));
    Splice(e, Lnext(a));
    Splice(Sym(e), b);
    RETURN e
  END Connect;
  
PROCEDURE Disconnect(a: Arc) =
  BEGIN
    IF a # Onext(a) THEN Splice(a, Oprev(a)) END;
    a := Sym(a);
    IF a # Onext(a) THEN Splice(a, Oprev(a)) END;
  END Disconnect;
  
PROCEDURE RightOf(s: REF Site; e: Arc): BOOLEAN =
  (* TRUE iff site "s" lies to the right of the line of edge "e". *)
  BEGIN
    RETURN TriSign(s^, Dest(e)^, Org(e)^) > 0
  END RightOf;

PROCEDURE LeftOf(s: REF Site; e: Arc): BOOLEAN =
  (* TRUE iff site "s" lies to the left of the line of edge "e". *)
  BEGIN
    RETURN TriSign(s^, Dest(e)^, Org(e)^) < 0
  END LeftOf;

PROCEDURE TriSign(READONLY a, b, c: Site): [-1 .. +1] =
  (* Sign of triangle "a b c" *)
  BEGIN
    WITH
      x1 = FLOAT(a.x, LONGREAL), y1 = FLOAT(a.y, LONGREAL),
      x2 = FLOAT(b.x, LONGREAL), y2 = FLOAT(b.y, LONGREAL),
      x3 = FLOAT(c.x, LONGREAL), y3 = FLOAT(c.y, LONGREAL),
      det = (x2*y3-y2*x3) - (x1*y3-y1*x3) + (x1*y2-y1*x2)
    DO
      IF det > 0.0D0 THEN
        RETURN +1
      ELSIF det < 0.0D0 THEN
        RETURN -1
      ELSE
        RETURN 0
      END
    END
  END TriSign;

PROCEDURE InCircle(READONLY a, b, c, d: Site): BOOLEAN =
  (* TRUE iff "d" is inside the circumcircle of "a b c" *)
  BEGIN
    WITH
      x1 = FLOAT(a.x, LONGREAL), y1 = FLOAT(a.y, LONGREAL),
      x2 = FLOAT(b.x, LONGREAL), y2 = FLOAT(b.y, LONGREAL),
      x3 = FLOAT(c.x, LONGREAL), y3 = FLOAT(c.y, LONGREAL),
      x4 = FLOAT(d.x, LONGREAL), y4 = FLOAT(d.y, LONGREAL),
      d1 = ((y4-y1)*(x2-x3)+(x4-x1)*(y2-y3))*((x4-x3)*(x2-x1)-(y4-y3)*(y2-y1)),
      d2 = ((y4-y3)*(x2-x1)+(x4-x3)*(y2-y1))*((x4-x1)*(x2-x3)-(y4-y1)*(y2-y3))
    DO
      RETURN d1 > d2
    END
  END InCircle;

BEGIN
END Delaunay.