Delaunay Triangulation

breed [points point]
undirected-link-breed [edges edge]
breed [triangles triangle]
globals [numpoints]
points-own [value super?]
triangles-own [Ux Uy rad2 bad? vertices sides]
edges-own [keep?]

to setup
  ca
  set numpoints 20; adjust if interface expanded to include more points
  let xlist (list 0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20)
  let ylist (list 0 y1 y2 y3 y4 y5 y6 y7 y8 y9 y10 y11 y12 y13 y14 y15 y16 y17 y18 y19 y20)
  let vallist (list 0 val1 val2 val3 val4 val5 val6 val7 val8 val9 val10 val11 val12 val13 val14 val15 val16 val17 val18 val19 val20)
  let onlist (list false point1 point2 point3 point4 point5 point6 point7 point8 point9 point10 point11 point12 point13 point14 point15 point16 point17 point18 point19 point20)
  if max xlist > max-pxcor - 1 or max ylist > max-pycor - 1 or min xlist < min-pxcor + 1 or min ylist < min-pycor + 1
     [user-message "One or more coordinate values are outside the screen range.  Please check your inputs and setup again"
      stop]
  create-points numpoints + 1 [ if who != 0 [set shape "circle" set size .6 set color red]
                      setxy item who xlist item who ylist
                      set value item who vallist
                      set label who set label-color yellow
                      set hidden? not item who onlist
                      set super? false
                      face patch 0 0
  ]
  ask points with [not hidden?] [if distance min-one-of other points with [not hidden?] [distance myself] = 0
    [user-message "Two or more active points have identical coordinates.  Please halt, revise your inputs and setup again" stop]]
  create-points 1 [set hidden? true set super? true setxy (min-pxcor + 1) (min-pycor + 1)] ; who numpoints + 1
  create-points 1 [set hidden? true set super? true setxy (min-pxcor + 1) (max-pycor - 1)] ; who numpoints + 2
  create-points 1 [set hidden? true set super? true setxy (max-pxcor - 1) (max-pycor - 1)] ; who numpoints + 3
  create-points 1 [set hidden? true set super? true setxy (max-pxcor - 1) (min-pycor + 1)] ; who numpoints + 4
end 

to triangulate
  setup
  ifelse count points with [not hidden?] < 3 [user-message "You will need at least three active points to triangulate"]
  [

  make-triangle point (numpoints + 1) point (numpoints + 2) point (numpoints + 3) ; super-triangle 1
  make-triangle point (numpoints + 3) point (numpoints + 4) point (numpoints + 1) ; super-triangle 2
  let triset (sort-on [xcor] points with [not hidden?])
  foreach triset ; for each point in the pointlist
    [p ->
     ask triangles [if inside-circum? p self [set bad? true]] ; for each triangle in the triangulation flag points within circumcircle

     let edgelist []
     ask triangles with [bad?] ; for each new bad triangle
      [
        if Pause-between-increments? [ask link-set sides [set color yellow] ask p [set size 2] user-message "Next step?" ask p [set size 1]]
        ask link-set sides [set edgelist ifelse-value member? self edgelist [remove self edgelist][lput self edgelist]]]
     clean-up
     foreach edgelist [ed -> make-triangle [end1] of ed [end2] of ed p];create new triangles

     ]
    ask triangles with [((member? point (numpoints + 1) vertices) or (member? point (numpoints + 2) vertices)
                      or (member? point (numpoints + 3) vertices) or (member? point (numpoints + 4) vertices))][set bad? true]
    if Pause-between-increments? [ user-message "Complete triangulation?"]
    clean-up
    ask links with [color = red][die]
   ]
end 

to make-triangle [p1 p2 p3]
  let me 0
  create-triangles 1
   [set hidden? true
    set bad? false
    set vertices turtle-set (list p1 p2 p3)
    let Ax [xcor] of p1
    let Ay [ycor] of p1
    let Bx [xcor] of p2
    let By [ycor] of p2
    let Cx [xcor] of p3
    let Cy [ycor] of p3
    let midx (max-pxcor + min-pxcor) / 2
    let midy (max-pycor + min-pycor) / 2
    if [super?] of p1 [set Ax midx - 10000 * (midx - Ax) set Ay midy - 10000 * (midy - Ay)] ; super vertices to "infinity"
    if [super?] of p2 [set Bx midx - 10000 * (midx - Bx) set By midy - 10000 * (midy - By)]
    if [super?] of p3 [set Cx midx - 10000 * (midx - Cx) set Cy midy - 10000 * (midy - Cy)]
    let d1x Bx - Ax
    let d1y Ay - By
    let d2x Cx - Bx
    let d2y By - Cy
    let d3x Ax - Cx
    let d3y Cy - Ay
    let A2 (Ax ^ 2 + Ay ^ 2)
    let B2 (Bx ^ 2 + By ^ 2)
    let C2 (Cx ^ 2 + Cy ^ 2)
    let D 2 * (Ax * d2y + Bx * d3y + Cx * d1y)
    set Ux (A2 * d2y + B2 * d3y + C2 * d1y) / D
    set Uy (A2 * d2x + B2 * d3x + C2 * d1x) / D
    set rad2 ((Ax - Ux) ^ 2 + (Ay - Uy) ^ 2)
    set me who
   ]
   ask [vertices] of triangle me [create-edges-with other [vertices] of triangle me [set hidden? false]]
   let vwho [who] of [vertices] of triangle me
   ask triangle me [set sides (list edge item 0 vwho item 1 vwho edge item 1 vwho item 2 vwho edge item 2 vwho item 0 vwho)]
end 

to-report inside-circum? [pt tri]
 let r1sq ([xcor] of pt - [Ux] of tri ) ^ 2 + ([ycor] of pt - [Uy] of tri ) ^ 2
 report r1sq < [rad2] of tri
end 

to clean-up
  ask edges [set keep? false]
  ask triangles [ifelse bad? [die][ask link-set sides [set keep? true set color green]]];remove bad triangles
  ask edges with [not keep?][set color red] ;remove unused edges
end