%!PS-Adobe-2.0
%%Creator: Peter Salvi (2011.02.16.)
%%Title: Incremental Delaunay Triangulation
%%BoundingBox: 0 0 500 580
%%EndComments

% Parameters
/large-radius 250 def           % half the bounding box
/vertex-radius 3 def
2 setlinewidth

% User points
/user-points [                  % array of arrays: [ [x1 y1] [x2 y2] ... ]
    [ 0 0 ] [ -0.5 -0.25 ] [ 0.5 -0.25 ]
    [ 0 -0.75 ] [ -0.5 0.75 ] [ 0.5 0.75 ] [ 0 1 ]
] def                           % defined in the [-1 -1] x [1 1] square

% Random points (only applies if there are no user points)
/number-of-points 10 def
/divisions 50 def               % number of grid points to choose from
18 srand                        % random seed

% -----------------------------------------------------------------------------
% ---                        end of user parameters                         ---
% -----------------------------------------------------------------------------

% Math
/v+ {
    /q exch def /p exch def
    [ p 0 get q 0 get add p 1 get q 1 get add ]
} def
/v* {
    /x exch def /p exch def
    [ p 0 get x mul p 1 get x mul ]
} def
/v- { -1 v* v+ } def
/vlength2 { aload pop dup mul exch dup mul add } def
/vlength { vlength2 sqrt } def
/distance { v- vlength } def
/sqr { dup mul } def
/circumcircle {
    aload pop
    dup 0 get /cx exch def 1 get /cy exch def
    dup 0 get /bx exch def 1 get /by exch def
    dup 0 get /ax exch def 1 get /ay exch def
    /d by cy sub ax mul cy ay sub bx mul add ay by sub cx mul add 2 mul def
    [ by cy sub ax sqr ay sqr add mul
      cy ay sub bx sqr by sqr add mul add
      ay by sub cx sqr cy sqr add mul add d div
      cx bx sub ax sqr ay sqr add mul
      ax cx sub bx sqr by sqr add mul add
      bx ax sub cx sqr cy sqr add mul add d div ]
    dup [ ax ay ] distance
} def

% Variables
/random-points user-points length 0 eq def
random-points not { /number-of-points user-points length def } if
/large-triangle [ [ 3 sqrt -2 div -0.5 ] [ 3 sqrt 2 div -0.5 ] [ 0 1 ] ] def
/center [ large-radius large-radius ] def
/points [
    large-triangle { large-radius v* center v+ } forall
    0 1 number-of-points 1 sub {
        random-points
            { pop [
                rand divisions mod divisions 2 div sub divisions 2 div div
                rand divisions mod divisions 2 div sub divisions 2 div div
            ] }
            { user-points exch get }
        ifelse
        large-radius 2 div 2 sqrt div v* center v+
    } for
] def
/triangles [ [ 0 1 2 ] ] def
/all-edges {
    /tri exch def
    [ tri 0 get tri 1 get ] [ tri 1 get tri 2 get ] [ tri 2 get tri 0 get ]
} def
/vertex { points exch get } def
/triangle-points { /tri exch def [ tri { vertex } forall ] } def

% Graphics
/vertex-color { element-type /normal eq { 0 } { 1 } ifelse 0 0 setrgbcolor } def
/triangle-color { element-type /normal eq { 0.9 } { 0.6 } ifelse setgray } def
/normal { /element-type /normal def } def
/current { /element-type /current def } def
/point { vertex-color aload pop vertex-radius 0 360 arc fill } def
/triangle {
    triangle-color
    dup 0 get vertex aload pop moveto
    dup 1 get vertex aload pop lineto
        2 get vertex aload pop lineto
    closepath gsave 0 setgray stroke grestore fill
} def
/circle { vertex-color exch aload pop 3 -1 roll 0 360 arc stroke } def % [x y] r

% Delaunay
/check-triangle {                                    % [i j k] + current-point
    triangle-points circumcircle exch                % r [Ox Oy]
    current-point vertex                             % r [Ox Oy] [Px Py]
    distance le
} def
/check-triangles {
    [ triangles {
        dup check-triangle { pop } if
    } forall ]
} def
/same-triangle {
    /b exch def /a exch def
    a 0 get b 0 get eq
    a 1 get b 1 get eq and
    a 2 get b 2 get eq and
} def
/bad-triangle {                                      % [i j k] + bad-triangles
    /tri exch def
    false bad-triangles {
        tri same-triangle or
    } forall
} def
/same-edge {
    /b exch def /a exch def
    a 0 get b 0 get eq
    a 1 get b 1 get eq and
    a 1 get b 0 get eq
    a 0 get b 1 get eq and or
} def
/count-edge {                                        % edges edge
    /actual exch def
    0 exch {
        actual same-edge { 1 add } if
    } forall
} def
/unique-edges {                                      % edges
    /edges exch def
    [ edges {
        dup edges exch count-edge 1 gt { pop } if
    } forall ]
} def
/walk-polygon {                                      % edges
    /edges exch def
    [ edges 0 get 0 get edges 0 get 1 get
    2 1 edges length { pop
        /next exch def /last exch def last
        edges {
            /actual exch def
            actual 0 get next eq actual 1 get last ne and
            {
                next actual 1 get
            }
            {
                actual 1 get next eq actual 0 get last ne and
                {
                    next actual 0 get
                } if
            } ifelse
        } forall
    } for pop ]
} def
/star-triangles {                                    % polygon point
    /p exch def /poly exch def
    [ 0 1 poly length 1 sub {
        /i exch def
        [ p poly i get poly i 1 add poly length mod get ]
    } for ]
} def
/new-triangles {                                     % triangles + current-point
    /tri exch def
    [ tri { all-edges } forall ]
    unique-edges walk-polygon current-point star-triangles
} def

% -----------------------
% --- Page Generation ---
% -----------------------

% Page 0: only the points
normal
3 1 points length 1 sub {
    vertex point
} for
showpage

% Page 1: points + large triangle
triangles {
    triangle
} forall
points {
    point
} forall
showpage

% Pages 2 - 2n+1: construction phases
3 1 points length 1 sub {
    /current-point exch def
    normal
    triangles {
        triangle
    } forall
    /bad-triangles check-triangles def
    current bad-triangles { triangle } forall
    current bad-triangles { triangle-points circumcircle circle } forall
    0 1 points length 1 sub {
        dup current-point eq { current } { normal } ifelse
        vertex point
    } for
    showpage

    [ triangles {                         % remove bad triangles
        dup bad-triangle { pop } if
    } forall ] /triangles exch def
    
    current
    bad-triangles new-triangles dup {
        triangle
    } forall
    normal
    triangles {
        triangle
    } forall
    0 1 points length 1 sub {
        dup current-point eq { current } { normal } ifelse
        vertex point
    } for
    dup dup length triangles length add array  % new array to hold the triangles
    dup 0 4 -1 roll putinterval                % put the new triangles inside
    dup 3 -1 roll length triangles putinterval % and the old triangles as well
    /triangles exch def
    showpage
} for

% Page 2n+2: the final Delaunay triangulation
normal
triangles {
    /tri exch def    
    tri 0 get 2 gt tri 1 get 2 gt and tri 2 get 2 gt and { tri triangle } if
} forall
3 1 points length 1 sub {
    vertex point
} for
showpage