;;; We're going to try to use enough data abstraction that we can change the
;;; format of how we store data without changing the code that uses it

;;; We will define a triangulation, which will be a decomposition of a polygon
;;; into triangles such that no corner of one triangle is in the interior of
;;; an edge of another triangle.
;;; The triangulation will consist of vertices. edges, and triangles; these
;;; object will reference each other.
;;; The data structures were designed (very) roughly as follows:
;;; 1.  If the number of objects is known a priori, then they are itemized
;;;     (vertex-1, vertex-2 in an edge, etc.)
;;; 2.  If the number of objects is not known a priori, then they are
;;;     contained in some type of container object, usually a list or a vector.
;;;     Programs *in this file* can use the knowledge of what type the
;;;     container object is; programs in other files may not, and must use
;;;     the generic conversion routines *->list and *<-list before using map,
;;;     for-each, etc.
;;; Other design choices could have been made.  For example, we could define
;;; Map and For-each that work on either lists or vectors (with Map returning
;;; a structure of the same type), and Map-with-index and For-each-with-index
;;; that takes a function (f i x y z ...) of n+1 arguments, the first of which
;;; is the index, etc.  (Sometimes we need this, and it obviates the need to
;;; use vectors or the (iota n) function, which returns the values from 0
;;; to n-1.)  If you want to see which of these methods is clearer or faster,
;;; go right ahead.

;;; Note that all vertices, edges, and triangles are unique, so we can use
;;; eq? (and memq, etc.) to compare them.

(define-class Vertex Object
  ((= point)                        ; the location of the vertex
   (= index                         ; each vertex has a unique index
      initializer: (lambda () #f))
   (= edges 
      initializer: (lambda () '())) ; the edges coming into this vertex
   (= fine                          ; data at finer level (not well-defined)
      initializer: (lambda () '()))
   ))

;; Vertex-edges is a vector.

(define (Vertex-for-each-edge proc v)
  (vector-for-each-1 proc (Vertex-edges v)))

(define (Vertex-edges<-list l)
  (list->vector l))

(define (Vertex-edges->list e)
  (vector->list e))

(define (Vertex-edges-length v)
  (vector-length (Vertex-edges v)))

(define (Vertex-contains-edge? v e)
  (declare (not inline) (fixnum))
  (let ((edges (Vertex-edges v)))
    (let loop ((i (- (vector-length edges) 1)))
      (and (>= i 0)
	   (or (eq? e (vector-ref edges i))
	       (loop (- i 1)))))))

(define-class Boundary-vertex Vertex
  ((= boundary-index
      initializer: (lambda () #f))
   (= boundary-edges
      initializer: (lambda () '()))))

(define (Boundary-vertex-for-each-boundary-edge proc v)
  (let ((boundary-edges (Boundary-vertex-boundary-edges v)))
    (declare (fixnum))
    (if (= (vector-length boundary-edges) 2)
	(begin
	  (proc (vector-ref boundary-edges 0))
	  (proc (vector-ref boundary-edges 1))))))

(define (Boundary-vertex-boundary-edges<-list l)
  (list->vector l))

(define (Boundary-vertex-boundary-edges->list e)
  (vector->list e))

(define (Boundary-vertex-boundary-edges-length v)
  (vector-length (Boundary-vertex-boundary-edges v)))

(define (Boundary-vertex-contains-boundary-edge? v e)
  (declare (not inline) (fixnum))
  (let ((edges (Boundary-vertex-boundary-edges v)))
    (let loop ((i (- (vector-length edges) 1)))
      (and (>= i 0)
	   (or (eq? e (vector-ref edges i))
	       (loop (- i 1)))))))

(define-class Edge Object
  ((= vertex-1)                      ; the two vertices
   (= vertex-2)
   (= index                          ; each edge has a unique index
      initializer: (lambda () #f))
   (= fine                           ; data at finer level (not well defined)
      initializer: (lambda () '()))
   ))

(define-class Boundary-edge Edge
  ((= boundary-index maybe-uninitialized:)
   (= triangle maybe-uninitialized:)))

(define (Edge-for-each-vertex proc e)
  (with-access 
   e (Edge vertex-1 vertex-2)
   (proc vertex-1)
   (proc vertex-2)))

(define (Edge-contains-vertex? e v)
  (with-access 
   e (Edge vertex-1 vertex-2)
   (or (eq? vertex-1 v)
       (eq? vertex-2 v))))

(define (Edge-length e)
  (with-access 
   e (Edge vertex-1 vertex-2)
   (Point-L2-norm 
    (Point-subtract (Vertex-point vertex-1)
		    (Vertex-point vertex-2)))))

;;; in triangles, edge-k goes from vertex-k to vertex-k+1 
;;; (point-4 being point-1)

(define-class Triangle Object
  ((= vertex-1)                    ; counter-clockwise
   (= vertex-2)
   (= vertex-3)
   (= edge-1)                      ; counter-clockwise
   (= edge-2)
   (= edge-3)
   (= index                        ; each triangle has a unique index
      initializer: (lambda () #f))
   (= fine                         ; data at finer level (not well defined)
      initializer: (lambda () '()))
   ))

(define (Triangle-for-each-vertex proc t)
  (with-access 
   t (Triangle vertex-1 vertex-2 vertex-3)
   (proc vertex-1)
   (proc vertex-2)
   (proc vertex-3)))

(define (Triangle-for-each-edge proc t)
  (with-access 
   t (Triangle edge-1 edge-2 edge-3)
   (proc edge-1)
   (proc edge-2)
   (proc edge-3)))

(define (Triangle-contains-vertex? t v)
  (with-access 
   t (Triangle vertex-1 vertex-2 vertex-3)
   (or (eq? vertex-1 v)
       (eq? vertex-2 v)
       (eq? vertex-3 v))))

(define (Triangle-contains-edge? t e)
  (with-access 
   t (Triangle edge-1 edge-2 edge-3)
   (or (eq? edge-1 e)
       (eq? edge-2 e)
       (eq? edge-3 e))))

(define (Triangle-area t)
  (with-access
;;; half the determinant of the associated affine transform from
;;; (0,0)->p1, (1,0)->p2, (0,1)->p3.
   t (Triangle vertex-1 vertex-2 vertex-3)
   (let ((point-1 (Vertex-point vertex-1))
	 (point-2 (Vertex-point vertex-2))
	 (point-3 (Vertex-point vertex-3)))
     (let ((x1 (Point-x point-1)) (y1 (Point-y point-1))
	   (x2 (Point-x point-2)) (y2 (Point-y point-2))
	   (x3 (Point-x point-3)) (y3 (Point-y point-3)))
       (declare (flonum))
       (abs (* 0.5 (- (* (- x2 x1) 
			 (- y3 y1)) 
		      (* (- x3 x1) 
			 (- y2 y1)))))))))

;;; For now, vertices, edges, and triangles in a triangulation are vectors.
;;; I rely on being able to randomly access in constant time the i'th
;;; vertex, edge, triangle, etc.

(define-class Triangulation Object
  ((= vertices)
   (= edges)
   (= triangles)
   (= has-polygonal-boundary?)
   (= boundary-vertices maybe-uninitialized:)
   (= boundary-edges maybe-uninitialized:)
   ;; an f64vector that contains min-x, max-x, min-y, and max-y (in that order)
   ;; of the coordinates of the points in the vertices of the triangulation
   (= limits initializer: (lambda () #f))
   ;; data at coarser level (not well defined)
   (= coarse initializer: (lambda () '()))
   ;; data at finer level (not well defined)
   (= fine initializer: (lambda () '()))
   ))

;;; define accessors and setters for limits vector in Triangulation

(define (make-limits #!optional (x-max -inf.0) (x-min +inf.0) (y-max -inf.0) (y-min +inf.0))
  (f64vector x-min x-max y-min y-max))

(define (x-min limits)
  (f64vector-ref limits 0))

(define (x-min-set! limits val)
  (f64vector-set! limits 0 val))

(define (x-max limits)
  (f64vector-ref limits 1))

(define (x-max-set! limits val)
  (f64vector-set! limits 1 val))

(define (y-min limits)
  (f64vector-ref limits 2))

(define (y-min-set! limits val)
  (f64vector-set! limits 2 val))

(define (y-max limits)
  (f64vector-ref limits 3))

(define (y-max-set! limits val)
  (f64vector-set! limits 3 val))

(define (Triangulation-Hilbert-curve-sort! t)
  (with-access
   t (Triangulation limits vertices edges triangles)
   (let ((x-scale (FLOAT (/ (- (x-max limits) (x-min limits)))))
	 (y-scale (FLOAT (/ (- (y-max limits) (y-min limits)))))
	 (x-min (x-min limits))
	 (y-min (y-min limits)))
     (let ((scale-x (lambda (x) (FLOAT (* x-scale (- x x-min)))))
	   (scale-y (lambda (y) (FLOAT (* y-scale (- y y-min)))))
	   (digits (+ 1 (inexact->exact ; in generic arithmetic
			 (ceiling (/ (log (exact->inexact 
					   (Triangulation-vertices-length t)))
				     (log 4.0)))))))
       
       ;; Find the index of a point on the Hilbert curve of order kmax.
       
       (let ((Hilbert-curve-index 
	      (lambda (x-coor y-coor kmax)
		(define (loop x y i k)
		  (if (FIX (< k kmax))
		      (if (FLOAT (< x 0.5))
			  (if (FLOAT (< y 0.5)) 
			      (loop (FLOAT (* 2.0 y))
				    (FLOAT (* 2.0 x))
				    (FIX (* 4 i))
				    (FIX (+ k 1)))
			      (loop (FLOAT (* 2.0 x))
				    (FLOAT (- (* 2.0 y) 1.0))
				    (FIX (+ 1 (* 4 i)))
				    (FIX (+ k 1))))
			  (if (FLOAT (< y 0.5))
			      (loop (FLOAT (- 1.0 (* 2.0 y)))
				    (FLOAT (* 2.0 (- 1.0 x)))
				    (FIX (+ 3 (* 4 i)))
				    (FIX (+ k 1)))
			      (loop (FLOAT (- (* 2.0 x) 1.0))
				    (FLOAT (- (* 2.0 y) 1.0))
				    (FIX (+ 2 (* 4 i)))
				    (FIX (+ k 1)))))
		      i))
		(loop x-coor y-coor 0 0))))
	 
	 ;; sort vertices according to the Hilbert curve
	 
	 (set! vertices (list->vector (map1 (lambda (x) (cdr x))
					    (sort! (map1 (lambda (v) 
							     (with-access
							      v (Vertex point index)
							      (cons (let ((sort-index (Hilbert-curve-index (Point-x point)
													   (Point-y point)
													   digits)))
								      (set! index sort-index)
								      sort-index)
								    v)))
							   (vector->list vertices))
						      (lambda (x y)
							(FIX (< (car x) (car y))))))))
	 
	 ;; sort edges according to the Hilbert curve of (Vertex-point (Edge-vertex-1 e))
	 
	 (set! edges (list->vector (map1 (lambda (x) (cdr x))
					 (sort! (map1 (lambda (e) 
							  (cons (Vertex-index (Edge-vertex-1 e))
								e))
							(vector->list edges))
						   (lambda (x y)
						     (FIX (< (car x) (car y))))))))
	 
	 ;; sort triangles according to the Hilbert curve of (Vertex-point (Triangle-vertex-1 t))
	 
	 (set! triangles (list->vector (map1 (lambda (x) (cdr x))
					     (sort! (map1 (lambda (t) 
							      (cons (Vertex-index (Triangle-vertex-1 t))
								    t))
							    (vector->list triangles))
						       (lambda (x y)
							 (FIX (< (car x) (car y)))))))))))))

(define (Triangulation-z-curve-sort! t)
  (with-access
   t (Triangulation limits vertices edges triangles)
   (let ((x-scale (FLOAT (/ (- (x-max limits) (x-min limits)))))
	 (y-scale (FLOAT (/ (- (y-max limits) (y-min limits)))))
	 (x-min (x-min limits))
	 (y-min (y-min limits)))
     (let ((scale-x (lambda (x) (FLOAT (* x-scale (- x x-min)))))
	   (scale-y (lambda (y) (FLOAT (* y-scale (- y y-min)))))
	   (digits (+ 1 (inexact->exact ; in generic arithmetic
			 (ceiling (/ (log (exact->inexact 
					   (Triangulation-vertices-length t)))
				     (log 4.0)))))))
       
       ;; Find the index of a point on the z-curve of order kmax.
       
       (let ((z-curve-index 
	      (lambda (x-coor y-coor kmax)
		(define (loop x y i k)
		  (if (FIX (< k kmax))
		      (if (FLOAT (< x 0.5))
			  (if (FLOAT (< y 0.5)) 
			      (loop (FLOAT (* 2.0 x))
				    (FLOAT (* 2.0 y))
				    (FIX (* 4 i))
				    (FIX (+ k 1)))
			      (loop (FLOAT (* 2.0 x))
				    (FLOAT (- (* 2.0 y) 1.0))
				    (FIX (+ 1 (* 4 i)))
				    (FIX (+ k 1))))
			  (if (FLOAT (< y 0.5))
			      (loop (FLOAT (- (* 2.0 x) 1.0))
				    (FLOAT (* 2.0 y))
				    (FIX (+ 2 (* 4 i)))
				    (FIX (+ k 1)))
			      (loop (FLOAT (- (* 2.0 x) 1.0))
				    (FLOAT (- (* 2.0 y) 1.0))
				    (FIX (+ 3 (* 4 i)))
				    (FIX (+ k 1)))))
		      i))
		(loop x-coor y-coor 0 0))))
	 
	 ;; sort vertices according to the z-curve
	 
	 (set! vertices (list->vector (map1 (lambda (x) (cdr x))
					    (sort! (map1 (lambda (v) 
							     (with-access
							      v (Vertex point index)
							      (cons (let ((sort-index (z-curve-index (Point-x point)
												     (Point-y point)
												     digits)))
								      (set! index sort-index)
								      sort-index)
								    v)))
							   (vector->list vertices))
						      (lambda (x y)
							(FIX (< (car x) (car y))))))))
	 
	 ;; sort edges according to the z-curve of (Vertex-point (Edge-vertex-1 e))
	 
	 (set! edges (list->vector (map1 (lambda (x) (cdr x))
					 (sort! (map1 (lambda (e) 
							  (cons (Vertex-index (Edge-vertex-1 e))
								e))
							(vector->list edges))
						   (lambda (x y)
						     (FIX (< (car x) (car y))))))))
	 
	 ;; sort triangles according to the z-curve of (Vertex-point (Triangle-vertex-1 t))
	 
	 (set! triangles (list->vector (map1 (lambda (x) (cdr x))
					     (sort! (map1 (lambda (t) 
							      (cons (Vertex-index (Triangle-vertex-1 t))
								    t))
							    (vector->list triangles))
						       (lambda (x y)
							 (FIX (< (car x) (car y)))))))))))))

(define (Triangulation-Gray-code-curve-sort! t)
  (with-access
   t (Triangulation limits vertices edges triangles)
   (let ((x-scale (FLOAT (/ (- (x-max limits) (x-min limits)))))
	 (y-scale (FLOAT (/ (- (y-max limits) (y-min limits)))))
	 (x-min (x-min limits))
	 (y-min (y-min limits)))
     (let ((scale-x (lambda (x) (FLOAT (* x-scale (- x x-min)))))
	   (scale-y (lambda (y) (FLOAT (* y-scale (- y y-min)))))
	   (digits (+ 1 (inexact->exact ; in generic arithmetic
			 (ceiling (/ (log (exact->inexact 
					   (Triangulation-vertices-length t)))
				     (log 4.0)))))))
       
       ;; Find the index of a point on the Gray code curve of order kmax.
       
       (let ((Gray-code-curve-index 
	      (lambda (x-coor y-coor kmax)
		(define (loop x y i k)
		  (if (FIX (< k kmax))
		      (if (FLOAT (< x 0.5))
			  (if (FLOAT (< y 0.5)) 
			      (loop (FLOAT (* 2.0 x))
				    (FLOAT (* 2.0 y))
				    (FIX (* 4 i))
				    (FIX (+ k 1)))
			      (loop (FLOAT (- 1.0 (* 2.0 x)))
				    (FLOAT (* 2.0 (- 1.0 y)))
				    (FIX (+ 1 (* 4 i)))
				    (FIX (+ k 1))))
			  (if (FLOAT (< y 0.5))
			      (loop (FLOAT (- (* 2.0 x) 1.0))
				    (FLOAT (* 2.0 y))
				    (FIX (+ 3 (* 4 i)))
				    (FIX (+ k 1)))
			      (loop (FLOAT (* 2.0 (- 1.0 x)))
				    (FLOAT (* 2.0 (- 1.0 y)))
				    (FIX (+ 2 (* 4 i)))
				    (FIX (+ k 1)))))
		      i))
		(loop x-coor y-coor 0 0))))
	 
	 ;; sort vertices according to the Gray code curve
	 
	 (set! vertices (list->vector (map1 (lambda (x) (cdr x))
					    (sort! (map1 (lambda (v) 
							     (with-access
							      v (Vertex point index)
							      (cons (let ((sort-index (Gray-code-curve-index (Point-x point)
													     (Point-y point)
													     digits)))
								      (set! index sort-index)
								      sort-index)
								    v)))
							   (vector->list vertices))
						      (lambda (x y)
							(FIX (< (car x) (car y))))))))
	 
	 ;; sort edges according to the Gray code curve of (Vertex-point (Edge-vertex-1 e))
	 
	 (set! edges (list->vector (map1 (lambda (x) (cdr x))
					 (sort! (map1 (lambda (e) 
							  (cons (Vertex-index (Edge-vertex-1 e))
								e))
							(vector->list edges))
						   (lambda (x y)
						     (FIX (< (car x) (car y))))))))
	 
	 ;; sort triangles according to the Gray code curve of (Vertex-point (Triangle-vertex-1 t))
	 
	 (set! triangles (list->vector (map1 (lambda (x) (cdr x))
					     (sort! (map1 (lambda (t) 
							      (cons (Vertex-index (Triangle-vertex-1 t))
								    t))
							    (vector->list triangles))
						       (lambda (x y)
							 (FIX (< (car x) (car y)))))))))))))

(define (Triangulation-Sierpinski-curve-sort! t)
  (with-access
   t (Triangulation limits vertices edges triangles)
   (let ((x-scale (FLOAT (/ (- (x-max limits) (x-min limits)))))
	 (y-scale (FLOAT (/ (- (y-max limits) (y-min limits)))))
	 (x-min (x-min limits))
	 (y-min (y-min limits)))
     (let ((scale-x (lambda (x) (FLOAT (* x-scale (- x x-min)))))
	   (scale-y (lambda (y) (FLOAT (* y-scale (- y y-min)))))
	   (digits (+ 1 (inexact->exact ; in generic arithmetic
			 (ceiling (/ (log (exact->inexact 
					   (Triangulation-edges-length t)))
				     (log 2.0)))))))
       
       ;; Find the index of a point in a space-filling-curve order with k-max binary digits.
       
       (let ((curve-index 
	      (lambda (x-coor y-coor kmax)
		(define (loop x y i k kmax-1)
		  (if (FIX (< k kmax-1))
		      (if (FLOAT (> (+ x y) 1.0))
			  (if (FLOAT (> x 0.5))
			      (loop (FLOAT (- (* 2.0 x) 1.0))
				    (FLOAT (- (* 2.0 y) 1.0))
				    (FIX (+ (* 4 i) 3))
				    (FIX (+ k 2))
				    kmax-1)
			      (loop (FLOAT (* 2.0 (- 1.0 y)))
				    (FLOAT (* 2.0 x))
				    (FIX (+ (* 4 i) 2))
				    (FIX (+ k 2))
				    kmax-1))
			  (if (FLOAT (> y 0.5))
			      (loop (FLOAT (- (* 2.0 y) 1.0))
				    (FLOAT (- 1.0 (* 2.0 x)))
				    (FIX (+ (* 4 i) 1))
				    (FIX (+ k 2))
				    kmax-1)
			      (loop (FLOAT (* 2.0 x))
				    (FLOAT (* 2.0 y))
				    (FIX (* 4 i))
				    (FIX (+ k 2))
				    kmax-1)))
		      (if (FIX (= k kmax-1))
			  (if (FLOAT (> (+ x y) 1.0))
			      (FIX (+ (* 2 i) 1))
			      (FIX (* 2 i)))
			  i)))
		(let ((x-coor (scale-x x-coor))
		      (y-coor (scale-y y-coor)))
		  (if (FLOAT (> x-coor y-coor))
		      (loop (FLOAT (- 1.0 x-coor))
			    (FLOAT (- 1.0 y-coor))
			    1
			    1
			    (FIX (- kmax 1)))
		      (loop x-coor
			    y-coor
			    0
			    1
			    (FIX (- kmax 1))))))))
	 
	 ;; sort vertices according to the space-filling curve
	 
	 (set! vertices (list->vector (map1 (lambda (x) (cdr x))
					    (sort! (map1 (lambda (v) 
							     (with-access
							      v (Vertex point index)
							      (cons (let ((sort-index (curve-index (Point-x point)
												   (Point-y point)
												   digits)))
								      (set! index sort-index)
								      sort-index)
								    v)))
							   (vector->list vertices))
						      (lambda (x y)
							(FIX (< (car x) (car y))))))))
	 
	 ;; sort edges according to the space-filling curve of (Vertex-point (Edge-vertex-1 e))
	 
	 (set! edges (list->vector (map1 (lambda (x) (cdr x))
					 (sort! (map1 (lambda (e) 
							  (cons (Vertex-index (Edge-vertex-1 e))
								e))
							(vector->list edges))
						   (lambda (x y)
						     (FIX (< (car x) (car y))))))))
	 
	 ;; sort triangles according to the space-filling curve of (Vertex-point (Triangle-vertex-1 t))
	 
	 (set! triangles (list->vector (map1 (lambda (x) (cdr x))
					     (sort! (map1 (lambda (t) 
							      (cons (Vertex-index (Triangle-vertex-1 t))
								    t))
							    (vector->list triangles))
						       (lambda (x y)
							 (FIX (< (car x) (car y)))))))))))))

;; The following routine sorts vertices, edges, and triangles in lexicographical order, i.e., if
;; first on the y coordinates, then, if they are equal, on the x coordinates.

(define (Triangulation-lexicographical-sort! t)
  
  (sort! (Triangulation-vertices t)
	 (lambda (v1 v2)
	   (let ((point-1 (Vertex-point v1))
		 (point-2 (Vertex-point v2)))
	     (FLOAT (or (< (Point-y point-1) (Point-y point-2))
			(and (= (Point-y point-1) (Point-y point-2))
			     (< (Point-x point-1) (Point-x point-2))))))))

  (sort! (Triangulation-edges t)
	 (lambda (e1 e2)
	   (let ((point-1 (Vertex-point (Edge-vertex-1 e1)))
		 (point-2 (Vertex-point (Edge-vertex-1 e2))))
	     (FLOAT (or (< (Point-y point-1) (Point-y point-2))
			(and (= (Point-y point-1) (Point-y point-2))
			     (< (Point-x point-1) (Point-x point-2))))))))

   (sort! (Triangulation-triangles t)
	 (lambda (t1 t2)
	   (let ((point-1 (Vertex-point (Triangle-vertex-1 t1)))
		 (point-2 (Vertex-point (Triangle-vertex-1 t2))))
	     (FLOAT (or (< (Point-y point-1) (Point-y point-2))
			(and (= (Point-y point-1) (Point-y point-2))
			     (< (Point-x point-1) (Point-x point-2)))))))))

;; This a little bit of magic for Gambit-C.  When you declare block, you are
;; saying that any object that is not redefined in this file will never be
;; redefined.  Well, here we're making sure that Triangulation-sort! is
;; redefined in this file, so later in other files or in the REPL, we can
;; set it to whatever we want and this code will use the new value.

(define Triangulation-sort! #f)

(set! Triangulation-sort! Triangulation-Sierpinski-curve-sort!)

(define-method (initialize! (t Triangulation))
  
  ;; calculate the minimum and maximum extent of the vertices in t if
  ;; it has not been already initialized.
  
  (if (not (Triangulation-limits t))
      (let ((limits (make-limits)))
	(declare (flonum))
	(Triangulation-for-each-vertex
	 (lambda (v)
	   (with-access
	    v (Vertex point)
	    (let* ((point point)
		   (x (Point-x point))
		   (y (Point-y point)))
	      (x-min-set! limits (min (x-min limits) x))
	      (x-max-set! limits (max (x-max limits) x))
	      (y-min-set! limits (min (y-min limits) y))
	      (y-max-set! limits (max (y-max limits) y)))))
	 t)
	(if (or (= (x-max limits) (x-min limits))
		(= (y-max limits) (y-min limits)))
	    (error "Degenerate triangulation"))
	(Triangulation-limits-set! t limits)))
  
  (Triangulation-sort! t)
  
  ;; assign indices to vertices, edges, and triangles to match their order
  
  (vector-for-each-with-index-1
   (lambda (i v)
     (Vertex-index-set! v i))
   (Triangulation-vertices t))
  
  (vector-for-each-with-index-1
   (lambda (i e)
     (Edge-index-set! e i))
   (Triangulation-edges t))
  
  (vector-for-each-with-index-1
   (lambda (i t)
     (Triangle-index-set! t i))
   (Triangulation-triangles t))
  
  ;; initialize the boundary-vertices fields
  
  (initialize-field-value!
   t
   (list->vector
    (let ((vertices (Triangulation-vertices t)))
      (declare (fixnum))
      (do ((i 
	    (- (vector-length vertices) 1) 
	    (- i 1))
	   (result 
	    '() 
	    (let ((v (vector-ref vertices i)))
	      (if (Boundary-vertex? v)
		  (cons v result)
		  result))))
	  ((< i 0) result))))
   'boundary-vertices)
  
  ;; initialize the Boundary-vertex-boundary-indices 
  
  (vector-for-each-with-index-1
   (lambda (i v)
     (Boundary-vertex-boundary-index-set! v i))
   (Triangulation-boundary-vertices t))
  
  ;; initialize the boundary-edges fields
  
  (initialize-field-value!
   t
   (list->vector
    (let ((edges (Triangulation-edges t)))
      (declare (fixnum))
      (do ((i 
	    (- (vector-length edges) 1) 
	    (- i 1))
	   (result 
	    '() 
	    (let ((e (vector-ref edges i)))
	      (if (Boundary-edge? e)
		  (cons e result)
		  result))))
	  ((< i 0) result))))
   'boundary-edges)
  
  ;; initialize the Boundary-edge-boundary-indices 
  
  (vector-for-each-with-index-1
   (lambda (i e)
     (Boundary-edge-boundary-index-set! e i))
   (Triangulation-boundary-edges t))

  ;; add the triangle field to all boundary edges

  (Triangulation-for-each-triangle
   (lambda (t)
     (Triangle-for-each-edge
      (lambda (e)
	(if (Boundary-edge? e)
	    (Boundary-edge-triangle-set! e t)))
      t))
   t)
  
  ;; call the super-class's initialize method
  
  (call-next-method))

(define (Triangulation-vertices<-list l)
  (list->vector l))

(define (Triangulation-vertices->list v)
  (vector->list v))

(define (Triangulation-vertices-length t)
  (vector-length (Triangulation-vertices t)))

(define (Triangulation-boundary-vertices<-list l)
  (list->vector l))

(define (Triangulation-boundary-vertices->list v)
  (vector->list v))

(define (Triangulation-boundary-vertices-length t)
  (vector-length (Triangulation-boundary-vertices t)))

(define (Triangulation-edges<-list l)
  (list->vector l))

(define (Triangulation-edges->list e)
  (vector->list e))

(define (Triangulation-edges-length t)
  (vector-length (Triangulation-edges t)))

(define (Triangulation-boundary-edges<-list l)
  (list->vector l))

(define (Triangulation-boundary-edges->list e)
  (vector->list e))

(define (Triangulation-boundary-edges-length t)
  (vector-length (Triangulation-edges t)))

(define (Triangulation-triangles<-list l)
  (list->vector l))

(define (Triangulation-triangles->list t)
  (vector->list t))

(define (Triangulation-triangles-length t)
  (vector-length (Triangulation-triangles t)))

;;; See the definitions of vector-for-each-1 and vector-for-each-with-index-1
;;; in utilities.scm

(define (Triangulation-for-each-vertex proc t)
  (vector-for-each-1 proc (Triangulation-vertices t)))

(define (Triangulation-for-each-edge proc t)
  (vector-for-each-1 proc (Triangulation-edges t)))

(define (Triangulation-for-each-triangle proc t)
  (vector-for-each-1 proc (Triangulation-triangles t)))

(define (Triangulation-for-each-boundary-vertex proc t)
  (vector-for-each-1 proc (Triangulation-boundary-vertices t)))

(define (Triangulation-for-each-boundary-edge proc t)
  (vector-for-each-1 proc (Triangulation-boundary-edges t)))

(define-class Curved-edge Boundary-edge
  ((= parametrization)))        ; a mapping from [0,1] to the edge points

(define-class Triangle-with-curved-side Triangle
  ((= transformation)))         ; a mapping from the standard triangle with
                                ; vertices (0,0), (1,0), and (0,1) to the
                                ; triangle with curved edge.