(use chicken.random chicken.time)

(define-record-type point-2d
  (make-point-2d x y)
  point-2d?
  (x point-2d-x point-2d-x-set!)
  (y point-2d-y point-2d-y-set!))

(define-record-type circle
  (make-circle center radius radius-squared diameter)
  circle?
  (center circle-center)
  (radius circle-radius)
  (radius-squared circle-radius-squared)
  (diameter circle-diameter))

(define (create-circle center radius)
  (make-circle center radius (* radius radius) (* radius 2.0)))

(define-record-type triangle
  (make-triangle point-1 point-2 point-3 index-1 index-2 index-3 
                 centroid circumcircle)
  triangle?
  (point-1 triangle-point-1)
  (point-2 triangle-point-2)
  (point-3 triangle-point-3)
  (index-1 triangle-index-1)
  (index-2 triangle-index-2)
  (index-3 triangle-index-3)
  (centroid triangle-centroid) 
  (circumcircle triangle-circumcircle))

;; =========== Utilities ===============================

(define (filter pred lst)
  (foldr (lambda (x r) (if (pred x) (cons x r) r)) '() lst))

(define frand
  (let ((big (expt 2 16)))
    (lambda ()
      (+ 0.0 (/ (random big) big)))))

(define (random-point-2d)
  (make-point-2d (frand) (frand)))

(define (random-point-array-2d cnt)
  (let ((res (make-vector cnt)))
    (let loop ((i 0))
       (if (= i cnt)
           res
           (begin 
             (vector-set! res i (random-point-2d))
             (loop (+ i 1)))))))

(define (sort-by-x point-seq)
  (sort point-seq (lambda (p1 p2) (< (point-2d-x p1) (point-2d-x p2)))))

(define (sort-by-y point-seq)
  (sort point-seq (lambda (p1 p2) (< (point-2d-y p1) (point-2d-y p2)))))

(define (get-min-max point-vec)
  (let ((min-x (point-2d-x (vector-ref point-vec 0)))
        (min-y (point-2d-y (vector-ref point-vec 0)))
        (max-x (point-2d-x (vector-ref point-vec 0)))
        (max-y (point-2d-y (vector-ref point-vec 0)))
        (n (vector-length point-vec)))
     (let loop ((i 1))
       (if (>= i n)
           `((min-x ,min-x) (min-y ,min-y) (max-x ,max-x) (max-y ,max-y)) 
           (let ((x (point-2d-x (vector-ref point-vec i)))
                 (y (point-2d-y (vector-ref point-vec i))))
              (set! min-x (min x min-x))
              (set! max-x (max x max-x))
              (set! min-y (min y min-y))
              (set! max-y (max y max-y))
              (loop (+ i 1)))))))

(define (get-bounding-triangle point-vec idx)
  (let* ((res (get-min-max point-vec))
         (min-x (cadr (assq 'min-x res)))
         (min-y (cadr (assq 'min-y res)))
         (max-x (cadr (assq 'max-x res)))
         (max-y (cadr (assq 'max-y res)))
         (fact 2.0)
         (dx (* fact (- max-x min-x)))
         (dy (* fact (- max-y min-y))) 
         (pt1 (make-point-2d (- min-x dx) (- min-y (* dy 3.0))))
         (pt2 (make-point-2d (- min-x dx) (+ max-y dx)))
         (pt3 (make-point-2d (+ max-x (* 3.0 dx)) (+ max-y dx)))
         (new-points (list->vector (append  (vector->list point-vec) (list pt1 pt2 pt3)))))
   (list  new-points (create-triangle (+ idx 1) (+ idx 2) (+ idx 3) new-points))))

(define (indices-and-coordinates tri)
  (let ((i1 (triangle-index-1 tri))
        (i2 (triangle-index-2 tri))
        (i3 (triangle-index-3 tri))
        (p1 (triangle-point-1 tri))
        (p2 (triangle-point-2 tri))
        (p3 (triangle-point-3 tri)))
     (let ((x1 (point-2d-x p1))
           (y1 (point-2d-y p1))
           (x2 (point-2d-x p2))
           (y2 (point-2d-y p2))
           (x3 (point-2d-x p3))
           (y3 (point-2d-y p3)))
        (list (list i1 x1 y1) (list i2 x2 y2) (list i3 x3 y3)))))

;; ==========================================================================
;; Meat & Potatoes

(define (distance-sqr v0 v1)
  (let ((dx (- (point-2d-x v0) (point-2d-x v1)))
        (dy (- (point-2d-y v0) (point-2d-y v1))))
     (+ (* dx dx) (* dy dy))))

(define (distance v0 v1)
  (sqrt (distance-sqr v0 v1)))

(define (centroid p0 p1 p2)
  (make-point-2d (/ (+ (point-2d-x p0) (point-2d-x p1) (point-2d-x p2)) 3.0)
                 (/ (+ (point-2d-y p0) (point-2d-y p1) (point-2d-y p2)) 3.0)))

(define (point-2d-add v0 v1)
  (make-point-2d (+ (point-2d-x v0) (point-2d-x v1))
                 (+ (point-2d-y v0) (point-2d-y v1))))

(define (point-2d-sub v0 v1)
  (make-point-2d (- (point-2d-x v0) (point-2d-x v1))
                 (- (point-2d-y v0) (point-2d-y v1))))

(define (circumcircle v0 v1 v2)
  (let* ((v1-v0 (point-2d-sub v1 v0))
         (v2-v0 (point-2d-sub v2 v0))
         (v2-v1 (point-2d-sub v2 v1))
         (v1+v0 (point-2d-add v1 v0))
         (v2+v0 (point-2d-add v2 v0))
         (a (point-2d-x v1-v0)) 
         (b (point-2d-y v1-v0)) 
         (c (point-2d-x v2-v0)) 
         (d (point-2d-y v2-v0)) 
         (e (+ (* a (point-2d-x v1+v0))
               (* b (point-2d-y v1+v0))))
         (f (+ (* c (point-2d-x v2+v0))
               (* d (point-2d-y v2+v0))))
         (g (* 2.0 (- (* a (point-2d-y v2-v1))
                      (* b (point-2d-x v2-v1)))))
         (colinear? (< (abs g) 1.0e-15))
         (cx 0.0) 
         (cy 0.0)
         (dx 0.0)
         (dy 0.0))
    (if colinear?
        (let ((minx (min (point-2d-x v0) (point-2d-x v1) (point-2d-x v2)))
              (miny (min (point-2d-y v0) (point-2d-y v1) (point-2d-y v2)))
              (maxx (max (point-2d-x v0) (point-2d-x v1) (point-2d-x v2)))
              (maxy (max (point-2d-y v0) (point-2d-y v1) (point-2d-y v2))))
          (set! cx (* 0.5 (+ minx maxx)))
          (set! cy (* 0.5 (+ miny maxy)))
          (set! dx (- cx minx))
          (set! dy (- cy miny)))
        (begin
          (set! cx (/ (- (* d e) (* b f)) g))
          (set! cy (/ (- (* a f) (* c e)) g))
          (set! dx (- cx (point-2d-x v0)))
          (set! dy (- cy (point-2d-y v0)))))
    (let* ((r2 (+ (* dx dx) (* dy dy)))
           (r (sqrt r2)))
       (make-circle (make-point-2d cx cy) r r2 (* 2.0 r)))))
                   
(define (create-triangle idx0 idx1 idx2 point-vec)
   (let ((p0 (vector-ref point-vec idx0))
         (p1 (vector-ref point-vec idx1))
         (p2 (vector-ref point-vec idx2)))
     (make-triangle p0 p1 p2 idx0 idx1 idx2 
                    (centroid p0 p1 p2) (circumcircle p0 p1 p2))))

(define (in-circumcircle? tri pt)
  (let* ((cc (triangle-circumcircle tri))
         (c (circle-center cc))
         (d2 (distance-sqr c pt)))
    (<= d2 (circle-radius-squared cc))))

(define (edge=? a b)
  (or (and (= (car a) (car b))
           (= (cadr a) (cadr b)))
      (and (= (car a) (cadr b))
           (= (cadr a) (car b)))))

(define (unique-edge? edges a)
  (let ((n (length (filter (lambda (b) (edge=? a b)) edges))))
    (< n 2)))


(define (has-shared-vertices? tri1 tri2)
  (let ((ia1 (triangle-index-1  tri1))
        (ia2 (triangle-index-2  tri1))
        (ia3 (triangle-index-3  tri1))
        (ib1 (triangle-index-1  tri2))
        (ib2 (triangle-index-2  tri2))
        (ib3 (triangle-index-3  tri2)))
    (or (= ib1 ia1) (= ib1 ia2) (= ib1 ia3)
        (= ib2 ia1) (= ib2 ia2) (= ib2 ia3)
        (= ib3 ia1) (= ib3 ia2) (= ib3 ia3))))

(define (add-vertex iv triangles points)
  (let ((unaffected-tris (list))
        (edges (list)))
     
     (for-each (lambda (tri)
                  (if (not (in-circumcircle? tri (vector-ref points iv)))
                      (set! unaffected-tris (cons tri unaffected-tris)) 
                      (let* ((i1 (triangle-index-1 tri))
                             (i2 (triangle-index-2 tri))
                             (i3 (triangle-index-3 tri))
                             (e1 (list i1 i2))
                             (e2 (list i2 i3))
                             (e3 (list i3 i1)))
                        (set! edges (append (list e1 e2 e3) edges)))))
                triangles)

    (set! edges (filter (lambda (x) (unique-edge? edges x)) edges))

    (append unaffected-tris
            (map (lambda (x)
                    (let ((ii0 (car x))
                          (ii1 (cadr x)))  
                       (create-triangle ii0 ii1 iv points)))
                 edges))))

;;; triaungulate: procedure (triangulate points) --> (Listof Triangles)
;;;   - points: a vector of point-2d records
;;;   - return value: a list of triangle records
;;;                   (representing a delaunay triangulation of the points)
;;;
(define (triangulate points)
  (let* ((res (get-bounding-triangle points (- (vector-length points) 1)))
         (n (vector-length points))
         (super-triangle (cadr res))
         (points1 (car res))
         (triangles (list super-triangle)))
    (let pt-loop ((iv 0))       
      (if (>= iv n)
          triangles
          (begin
            (set! triangles  
                  (add-vertex iv triangles points1))
            (pt-loop (+ iv 1)))))

    (filter (lambda (t) (not (has-shared-vertices? super-triangle t)))
            triangles)))

(define (display-edge e points)
  (let ((i0 (car e))
        (i1 (cadr e)))
    (for-each (lambda (x) (display x (current-error-port)))
              (list (point-2d-x (vector-ref points i0)) "  " 
                    (point-2d-y (vector-ref points i0)) #\newline
                    (point-2d-x (vector-ref points i1)) "  " 
                    (point-2d-y (vector-ref points i1)) #\newline
                    #\newline))))

(define (display-tris tris)
(for-each (lambda (t)
       (for-each (lambda (pt)
                   (for-each display (list (point-2d-x pt) "  " 
                                           (point-2d-y pt) #\newline)))
                         (list (triangle-point-1 t) (triangle-point-2 t) 
                               (triangle-point-3 t) (triangle-point-1 t)))
               (newline))
            tris))

(define (go)
  (randomize)
  (let* ((points (random-point-array-2d (string->number (list-ref (argv) 1))))
         (triangles (triangulate points)))
     (display-tris triangles)))


(go)