#!/usr/bin/env ruby
#
#  ported from p bourke's triangulate.c
#  http://astronomy.swin.edu.au/~pbourke/terrain/triangulate/triangulate.c
#
#  C to Java:
#  fjenett, 20th february 2005, offenbach-germany.
#  contact: http://www.florianjenett.de/
#
#  Java to Ruby:
#  Gregory Seidman, 2006-11-28, Washington, DC, USA
#  contact: gmail account gseidman
#  
#  to view the output: http://processing.org/
#
#  usage: (random vertices) ruby triangulate.rb [num vertices]
#         (benchmark)       ruby triangulate.rb bm

module Delaunay
  extend self
  ITRIANGLE = Struct.new(:p1, :p2, :p3, :complete)
  class IEDGE < Struct.new(:p1, :p2)
    def ==(o)
      (p1 == o.p1 && p2 == o.p2) || (p1 == o.p2 && p2 == o.p1)
    end

    def valid?
      p1 and p2
    end

    def reset!
      self.p1 = self.p2 = nil
    end
  end
  Coord = Struct.new(:x, :y)
  EPSILON = 0.000001

  # Triangulation subroutine
  # Takes as input vertices in array verts
  # Returned is a list of triangular faces in the array tris
  # These triangles are arranged in a consistent clockwise order.

  def triangulate(verts)
    edges = []
    tris = []
    # sort by X coord
    verts = verts.sort_by {|p|p.x}
    # center/radius cache used by circum_circle
    cc_cache = {}

    # Find the maximum and minimum vertex bounds.
    # This is to allow calculation of the bounding triangle
    xmin = verts[0].x
    ymin = verts[0].y
    xmax = xmin
    ymax = ymin
    verts.each { |p|
      xmin = p.x if (p.x < xmin)
      xmax = p.x if (p.x > xmax)
      ymin = p.y if (p.y < ymin)
      ymax = p.y if (p.y > ymax)
    }
    dx = xmax - xmin
    dy = ymax - ymin
    dmax = (dx > dy) ? dx : dy
    xmid = (xmax + xmin) / 2.0
    ymid = (ymax + ymin) / 2.0

    # 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.
    nv = verts.size
    verts << Coord.new(xmid - 2.0 * dmax, ymid - dmax)
    verts << Coord.new(xmid, ymid + 2.0 * dmax)
    verts << Coord.new(xmid + 2.0 * dmax, ymid - dmax)
    tris << ITRIANGLE.new(nv, nv+1, nv+2)


    # Include each point one at a time into the existing mesh
    (0...verts.size).each { |i|
      p = verts[i]
      edges.clear

      # Set up the edge buffer.
      # If the point (xp,yp) lies inside the circumcircle then the
      # three edges of that triangle are added to the edge buffer
      # and that triangle is removed.
      j = 0
      while j < tris.size
        unless tris[j].complete
          p1 = verts[tris[j].p1]
          p2 = verts[tris[j].p2]
          p3 = verts[tris[j].p3]
          inside,xc,yc,r = circum_circle(p, p1, p2, p3, cc_cache)
          if (xc + r) < p.x
            tris[j].complete = true
          end
          if inside
            edges << IEDGE.new(tris[j].p1, tris[j].p2)
            edges << IEDGE.new(tris[j].p2, tris[j].p3)
            edges << IEDGE.new(tris[j].p3, tris[j].p1)
            tri = tris.pop
            tris[j] = tri unless j >= tris.size
            j -= 1
          end
        end
        j += 1
      end #while j

      # Tag multiple edges
      # Note: if all triangles are specified anticlockwise then all
      # interior edges are opposite pointing in direction.
      j = 0
      while j < edges.size - 1
        k = j+1
        while k < edges.size
          if (edges[j] == edges[k])
            edges[j].reset!
            edges[k].reset!
          end
          k += 1
        end #while k
        j += 1
      end #while j

      # Form new triangles for the current point
      # Skipping over any tagged edges.
      # All edges are arranged in clockwise order.
      j = 0
      while j < edges.size
        tris << ITRIANGLE.new(edges[j].p1, edges[j].p2, i) if edges[j].valid?
        j += 1
      end #while j
    } #each i

    # Remove supertriangle vertices
    verts[-3..-1] = nil
    nv = verts.size

    # Remove triangles with supertriangle vertices
    # These are triangles which have a vertex number greater than nv
    i = 0
    while i < tris.size
      if (tris[i].p1 >= nv ||
          tris[i].p2 >= nv ||
          tris[i].p3 >= nv)
        tri = tris.pop
        tris[i] = tri unless i >= tris.size
        i -= 1
      end
      i += 1
    end #while i

    [ verts, tris ]
  end #triangulate

  private

  # Return TRUE if a point p is inside the circumcircle made up of the
  # points p1, p2, p3
  # The circumcircle center is returned in (xc,yc) and the radius r
  # The return value is an array [ inside, xc, yc, r ]
  # Takes an optional cache hash to use for radius/center caching
  # NOTE: A point on the edge is inside the circumcircle
  def circum_circle(p, p1, p2, p3, cache = nil)
    dx,dy,rsqr,drsqr = []
    cached = cache && cache[[p1, p2, p3]]
    xc, yc, r = []
    rsqr = 0

    if cached
      xc, yc, r = cached
      rsqr = r*r
    else
      # Check for coincident points
      if (p1.y-p2.y).abs < EPSILON && (p2.y-p3.y).abs < EPSILON
#        puts("CircumCircle: Points are coincident.")
        return [ false, 0, 0, 0 ]
      end

      if (p2.y-p1.y).abs < EPSILON
        m2 = - (p3.x-p2.x) / (p3.y-p2.y)
        mx2 = (p2.x + p3.x) * 0.5
        my2 = (p2.y + p3.y) * 0.5
        xc = (p2.x + p1.x) * 0.5
        yc = m2 * (xc - mx2) + my2
      elsif (p3.y-p2.y).abs < EPSILON
        m1 = - (p2.x-p1.x) / (p2.y-p1.y)
        mx1 = (p1.x + p2.x) * 0.5
        my1 = (p1.y + p2.y) * 0.5
        xc = (p3.x + p2.x) * 0.5
        yc = m1 * (xc - mx1) + my1
      else
        m1 = - (p2.x-p1.x) / (p2.y-p1.y)
        m2 = - (p3.x-p2.x) / (p3.y-p2.y)
        mx1 = (p1.x + p2.x) * 0.5
        mx2 = (p2.x + p3.x) * 0.5
        my1 = (p1.y + p2.y) * 0.5
        my2 = (p2.y + p3.y) * 0.5
        xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2)
        yc = m1 * (xc - mx1) + my1
      end

      dx = p2.x - xc
      dy = p2.y - yc
      rsqr = dx*dx + dy*dy
      r = Math.sqrt(rsqr)
      cache[[p1, p2, p3]] = [ xc, yc, r ] if cache
    end

    dx = p.x - xc
    dy = p.y - yc
    drsqr = dx*dx + dy*dy

    [ (drsqr <= rsqr), xc, yc, r ]
  end #circum_circle

end #Delaunay

if __FILE__ == $PROGRAM_NAME
  require 'benchmark'

  def main
    if ARGV[0] == 'bm'
      bm
    else
      nv = ARGV[0].to_i
      nv = 20 if (nv <= 0 || nv > 1000)
      Delaunay.output_random(nv)
    end
  end

  def bm
    Benchmark.bm(5) { |b|
      (1..10).each { |i| i *= 100
        b.report("#{i}\t") { 10.times { Delaunay.run_random(i) } }
      }
    }
  end

  module Delaunay

    def run_random(nv)
      points = (0...nv).map { |i| Coord.new(i*4.0, 400.0 * rand) }
      triangulate(points)
    end

    def output(points, tris)
      puts "void setup() { size(800, 800); }\nvoid draw() {"

      puts "\tscale(2);\n\tstrokeWeight(0.5);\n\tnoFill();\n\tbeginShape(TRIANGLES);"
      puts tris.map{ |t|
        tri_verts = [ points[t.p1],points[t.p2],points[t.p3] ]
        tri_verts.map! { |p| "\tvertex(#{p.x}, #{p.y});" }
      }.flatten!.join("\n")
      puts "\tendShape();\n\trectMode(CENTER);\n\tfill(0);"
      puts points.map{ |p| "\trect(#{p.x}, #{p.y}, 3, 3);" }.join("\n")
      puts "}"
    end

    def output_random(nv)
      output(*run_random(nv))
    end

  end

  main
end