class DELAUNAY
        -- Assumed points are given in left-to-right
        -- order.  Suboptimal, since the adjacency lists
	-- are repeatedly searched completely.

creation
  make

feature
  count : INTEGER is
    do
      Result := points.count
    end

  edge_count : INTEGER is
    local
      i : INTEGER
    do
      from
        i := 1
      until
        i > count
      loop
        Result := Result + neighbours.item(i).count
        i := i+1
      end
      Result := Result // 2
    end

  process ( lo, hi: INTEGER ) is
    require
      1 <= lo and lo <= hi and hi <= count
    local
      mid : INTEGER
    do
      if hi = lo+1 then
	add_edge ( lo, hi )
      elseif hi > lo then
        mid := (hi+lo)//2
        process (lo, mid)
        process (mid+1, hi)
        merge (lo,mid,hi)
      end
      triangulation_calculated := true
      io.put_string ("After process " + lo.out + ", "+ hi.out + "%N")
      print_out
    end

  edges : GBN_DLIST [ GBN_PAIR [ INTEGER, INTEGER ] ] is
    require
      precalculated : triangulation_calculated
    local
      i : INTEGER
      pair : GBN_PAIR [ INTEGER, INTEGER ]
      it : GBN_ITERATOR [ INTEGER ]
    do
      !! Result.make
      from
        i := 1
      until
        i > neighbours.count
      loop
        from
	  it := neighbours.item(i).iterator
	until
	  it.finished
	loop
	  !! pair.make ( i, it.item )
	  Result.add_last ( pair )
	  it.forth
	end
	i:= i+1
      end
    end

  print_out is
    local
      i : INTEGER
      it : GBN_DLIST_ITERATOR [ INTEGER ]
    do
      io.put_integer ( count )
      io.put_new_line
      from
        i := 1
      until
        i > count
      loop
        io.put_integer ( i )
        io.put_character ( ' ' )
        io.put_double ( points.item (i) . x )
        io.put_character ( ' ' )
        io.put_double ( points.item (i) . y )
        io.put_new_line
        i := i+1
      end
      io.put_integer ( 2 * edge_count )
        -- edges given twice, to retain the
        -- plane embedding information, since the
        -- edges are sorted anticlockwise around eac
        -- point.
      io.put_new_line
      from
        i := 1
      until
        i > count
      loop
        from
          it := neighbours.item (i) . iterator
        until
          it.finished
        loop
          io.put_integer ( i )
          io.put_character (' ')
          io.put_integer ( it.item )
          io.put_new_line
          it.forth
        end
        i := i+1
      end
    end

  triangulation_calculated: BOOLEAN

  convex_hull : GBN_DLIST [ INTEGER ] is
    require
      precalculated : triangulation_calculated
      local
        i, j, k: INTEGER
	finished: BOOLEAN
    do
      !! Result . make
      if count = 1 then
        Result.add_last ( 1 )
      else
        from
	  i := 1
	  j := second_on_hull
	until
	  finished
	loop
	  Result.add_last (i)
	  if j = 1 then
	    finished := true
	  else
	    k := next_on_hull (i,j)
	    i := j
	    j := k
	  end
	end
      end
    end

feature {NONE}

  second_on_hull : INTEGER is
    require
      two_points : count > 1
    local
      i,j: INTEGER
      it : GBN_ITERATOR [ INTEGER ]
    do
      from
        it := neighbours.item(1).iterator
	i := it.item
	it.forth
      until
        it.finished
      loop
        j := it.item
	if points.item(j).y <= points.item(1).y and then
	   points.item(j).x < points.item(i).x
	then
	  j := i
	end
	it.forth
      end
      Result := i
    end

  next_on_hull (i,j : INTEGER) : INTEGER is
    require
      triangulation_calculated : triangulation_calculated
    local
      list : GBN_DLIST [ INTEGER ]
      it : GBN_DLIST_ITERATOR [ INTEGER ]
    do
      list := neighbours.item (j)
      if list . count = 1 then
        Result := i
      else
        from
	  it := list.iterator
	until
	  it . finished
	loop
	  if it.item = i then
	    Result :=
	      list.item ( list.cyclic_succ ( it . place ) )
	    it . stop
	  else
	    it . forth
	  end
	end
      end
    end


  an_upper_tangent: INTEGER is 1
  a_lower_tangent: INTEGER is 2
  a_double_tangent: INTEGER is 3
  not_a_tangent: INTEGER is 4
  remove_an_edge : INTEGER is 5
  add_an_edge : INTEGER is 6

  find_pred ( j : INTEGER; i: INTEGER ) : 
    GBN_PAIR [ GBN_DLIST_PLACE [ INTEGER ], BOOLEAN ] is
          -- predecessor of point[j] among nbrs of point[i]
        -- in anticlockwise order.  Second component
        -- says whether the angle into which pt fits
        -- is a reflex angle.
    local
      it : GBN_DLIST_ITERATOR [ INTEGER ]
      list : GBN_DLIST [ INTEGER ]
      pt, pt1, pt2, pt3: POINT
      reflex_angle : BOOLEAN
    do
      pt := points.item (j)
      list := neighbours.item (i)
      pt1 := points.item (i)
      if list.count = 0 then
        !! Result.make ( Void, true )
	io.put_string ("neighbours["+i.out+"] empty%N")
      elseif list.count = 1 then
        !! Result.make ( list.first_place, true )
      else
        from
          it := list.iterator
        until
          it.finished
        loop
          pt2 := points.item (it.item)
          pt3 := points.item (list.item(list.cyclic_succ (it.place)))
          reflex_angle :=
              not pt3.left_of ( pt1, pt2 )
          if reflex_angle then
            if pt.left_of ( pt1, pt2 )
               or else
               pt.left_of ( pt3, pt1 )
            then
              !! Result.make ( it.place, true )
            end
          elseif pt.left_of (pt1, pt2) and then
             pt.left_of (pt3, pt1)
          then
            !! Result.make ( it.place, false )
          end
          if Result /= Void then
            it.stop
          else
            it.forth
          end
        end
      end
      io.put_string("find pred "+j.out+" among nbrs "+i.out+": ")
      if Result.item_1 = Void then
        io.put_string ("empty list, reflex angle")
      else
        io.put_string ( neighbours.item(i).item(Result.item_1).out )
	if Result.item_2 then
	  io.put_string ( " --- reflex angle%N")
	else
	  io.put_string ( " --- non-reflex angle%N")
	end
      end
    end

  find_succ ( j: INTEGER; i: INTEGER ) : 
    GBN_PAIR [ GBN_DLIST_PLACE [ INTEGER ], BOOLEAN ] is
    	-- successor of point[j] in list[i]
    do
      Result := find_pred ( j, i )
      Result.put_1 ( neighbours.item (i).cyclic_succ ( Result.item_1 ) )
      io.put_string("find succ "+j.out+" among nbrs "+i.out+": ")
      if Result.item_1 = Void then
        io.put_string ("empty list, reflex angle")
      else
        io.put_string ( neighbours.item(i).item(Result.item_1).out )
	if Result.item_2 then
	  io.put_string ( " --- reflex angle%N")
	else
	  io.put_string ( " --- non-reflex angle%N")
	end
      end
    end

  find_place ( i: INTEGER; list : GBN_DLIST [ INTEGER ] ):
      GBN_DLIST_PLACE [ INTEGER ] is
    do
      from
        Result := list.first_place
      until
        Result = Void or else list.item( Result ) = i
      loop
         Result  := list.succ (  Result  )
      end
    end

  remove_edge ( i, j : INTEGER ) is
    local
      i_in_j, j_in_i : GBN_DLIST_PLACE [ INTEGER ]
    do
      io.put_string("remove edge "+i.out+" "+j.out+"%N")
      j_in_i := find_place ( j, neighbours.item (i) )
      i_in_j := find_place ( i, neighbours.item (j) )
      if j_in_i /= Void then
        neighbours.item(i).remove ( j_in_i )
      end
      if i_in_j /= Void then
        neighbours.item(j).remove ( i_in_j )
      end
    end

  add_edge ( i: INTEGER; j: INTEGER ) is
    local
      pair : GBN_PAIR [ GBN_DLIST_PLACE [ INTEGER ], BOOLEAN ]
      list : GBN_DLIST [ INTEGER ]
      pred_i, succ_j : GBN_DLIST_PLACE [ INTEGER ]
    do
      io.put_string("add edge "+i.out+" to "+j.out+"%N")
      list := neighbours.item (i)
      if not list.is_empty then
        pair := find_pred ( j, i )
	pred_i := pair.item_1
      end

      list := neighbours.item (j)
      if not list.is_empty then
        pair := find_succ ( i, j )
	succ_j := pair.item_1
      end

      neighbours.item(i).add_after ( j, pred_i )
      neighbours.item(j).add_before ( i, succ_j )
    end

  announce_test_tangency ( j: INTEGER; i: INTEGER ) is
    local
      it : GBN_DLIST_ITERATOR [ INTEGER ]
    do
      io.put_string("test tangency "+j.out+" "+i.out+", nbrs")
      from
        it := neighbours.item (i) . iterator
      until
        it . finished
      loop
        io.put_string(", "+it.item.out)
	it.forth
      end
      io.put_string (": ")
    end
        
  test_tangency ( j: INTEGER; i: INTEGER ) : INTEGER is
  	-- test tangency of point[j] among nbrs[i]
    local
      list: GBN_DLIST [INTEGER]
      pair : GBN_PAIR [ GBN_DLIST_PLACE [INTEGER], BOOLEAN ]
      pt, pt1, pt2, pt3: POINT
      reflex_angle : BOOLEAN
    do
      announce_test_tangency (j, i)
      list := neighbours . item (i)
      pt := points.item (j)
      if list.is_empty then
        Result := a_double_tangent
      else
        pair := find_pred ( j, i )
        reflex_angle := pair.item_2
	
	-- io.put_string("preceding pt in nbrs("+i.out+"): "+
        --   list.item ( pair.item_1 ) . out + "%N" ) 

        if not reflex_angle then
          Result := not_a_tangent
        else
          pt1 := points.item (i)
          pt2 := points.item (list.item ( pair.item_1 ))   
          pt3 := points.item (list.item ( list.cyclic_succ (pair.item_1) ) )
          if pt2.left_of ( pt1, pt ) then
            if pt3.left_of ( pt1, pt ) then
	      if pt1.x < pt.x then
                Result := a_lower_tangent
	      else
                Result := an_upper_tangent
	      end
            else
              Result := not_a_tangent
            end
          else
            if not pt3.left_of ( pt1, pt ) then
	      if pt1.x < pt.x then
                Result := an_upper_tangent
	      else
                Result := a_lower_tangent
	      end
            else
              Result := not_a_tangent
            end
          end
        end
	inspect Result
	  when an_upper_tangent then io.put_string("upper tangent%N")
	  when a_lower_tangent then io.put_string("lower tangent%N")
	  when a_double_tangent then io.put_string("double tangent%N")
	  when not_a_tangent then io.put_string("not a tangent%N")
	end
      end
    ensure
      Result >= an_upper_tangent and Result <= not_a_tangent
    end

  upper_common_tangent ( lo, mid, hi: INTEGER ) is
    require
      1 <= lo and lo < mid and mid <= hi and hi <= count
    local
      ok_left, ok_right : BOOLEAN
      i, j: INTEGER
      tangency : INTEGER
      list: GBN_DLIST [INTEGER]
      pair : GBN_PAIR [ GBN_DLIST_PLACE [INTEGER], BOOLEAN ]
    do
      from
        i := mid
        if mid=lo then
          ok_left := true
        end
        j := mid+1
        if mid+1 = hi then
          ok_right := true
        end
      until
        ok_left and ok_right
      loop
        tangency := test_tangency ( j, i )
        ok_left := tangency = an_upper_tangent
	        or tangency = a_double_tangent
        if ok_left then
          tangency := test_tangency ( i, j )
          ok_right := tangency = an_upper_tangent
	        or tangency = a_double_tangent
          if not ok_right then
            list := neighbours.item (j)
            pair := find_pred ( i, j )
            j := list.item ( pair.item_1 )
          end
        else
          list := neighbours.item (i)
          pair := find_succ ( j, i )
          i := list.item ( pair.item_1 )
        end
      end -- loop
      left_upper_tangent := i; right_upper_tangent := j
      io.put_string("upper common tangent("+lo.out+", "+
        mid.out+", "+hi.out+"): "+i.out+", "+j.out+"%N")
    end

  is_lower_tangent ( j: INTEGER; i: INTEGER ) : BOOLEAN is
  	-- is points[j] a lower tangent among nbrs[i]?
    local
      test : INTEGER
    do
      if neighbours.item(i).is_empty then
        Result := true
      else
        test := test_tangency ( j, i )
        Result := test = a_lower_tangent
		or test = a_double_tangent
      end
      io.put_string("is lower tangent "+j.out+" at "+
        i.out +": "+Result.out+"%N")
    end

  is_upper_tangent ( j: INTEGER; i: INTEGER ) : BOOLEAN is
  	-- is points[j] an upper tangent among nbrs[i]?
    local
      test : INTEGER
    do
      if neighbours.item(i).is_empty then
        Result := true
      else
        test := test_tangency ( j, i )
        Result := test = an_upper_tangent
		or test = a_double_tangent
      end
      io.put_string("is upper tangent "+j.out+" at "+
        i.out +": "+Result.out+"%N")
    end

-- merge procedure. lo to mid are delaunayed, ditto mid+1 to hi.
-- begins with i,j upper common tangent.
-- repeatedly tests if i,j lower common tangent.
-- if not,  prepare to form next i,j
--    if i,j lower tangent on left then next_i := i
--    else next_i is next neighbour of i in clockwise order
--    such that the circle i,j,next_i contains no other
--    neighbours of i.
--
--    if i,j lower tangent on right then next_j := j
--    else next_j is next neighbour of j in anticlockwise order
--    such that circle i,j,next_j contains no other neighbours of
--    j.
--
--    If next_i /= i and next_j /= j then test whether
--    next_j is in circle i,j,next_i.  If so, then set
--    next_i to i.  Else set next_j to j.
--
--    Now replace i by next_i and j by next_j and continue.

  merge ( lo, mid, hi : INTEGER ) is
    local
      i, j,  next_i,  next_j: INTEGER
      pt1, pt2, pt3, pt4: POINT
      i_pred, j_succ : GBN_DLIST_PLACE [ INTEGER ]
      pair : GBN_PAIR [ GBN_DLIST_PLACE [ INTEGER ], BOOLEAN ]
      list: GBN_DLIST [ INTEGER ]
      it : GBN_DLIST_ITERATOR [ INTEGER ]
      action : GBN_TRIPLE [ INTEGER, INTEGER, INTEGER ]
      actions : GBN_DLIST [ GBN_TRIPLE [ INTEGER, INTEGER, INTEGER ] ]
      action_it : GBN_DLIST_ITERATOR [ GBN_TRIPLE [ INTEGER, INTEGER, INTEGER ] ]
      ok_on_left, ok_on_right : BOOLEAN
      i_1, i_2: INTEGER
    do
      io.put_string ("merge "+lo.out+" "+mid.out+" "+hi.out+"%N")
      upper_common_tangent ( lo, mid, hi )
      !! actions.make
      from
        i := left_upper_tangent
        j := right_upper_tangent
      until
        ok_on_left and ok_on_right
      loop
        io.put_string("merge: i="+i.out+", j="+j.out+"%N")
	
	pair := find_pred ( j, i )
	if  not pair.item_2 then
	  i_1 := neighbours.item(i).item(pair.item_1)
	  i_2 := neighbours.item(i).item(
	  	neighbours.item(i).cyclic_succ (pair.item_1))
          !!action.make ( remove_an_edge, i_1, i_2 )
	  actions.add_last ( action )
	end
	
	pair := find_pred ( i, j )
	if  not pair.item_2 then
	  i_1 := neighbours.item(j).item(pair.item_1)
	  i_2 := neighbours.item(j).item(
	  	neighbours.item(j).cyclic_succ (pair.item_1))
          !!action.make ( remove_an_edge, i_1, i_2 )
	  actions.add_last ( action )
	end

	!! action.make ( add_an_edge, i, j )
	actions.add_last ( action )

        ok_on_left := is_lower_tangent ( j, i )
        if ok_on_left then 
	  next_i := i
	end
        ok_on_right := is_lower_tangent ( i, j )
        if ok_on_right then 
	  next_j := j
	end

        if not ok_on_left then
          pt1 := points.item (i)
          pt2 := points.item (j)
          list := neighbours.item (i)
          pair := find_pred ( j, i )
          i_pred := pair.item_1
	  next_i := list.item ( i_pred )
          pt3 := points.item ( next_i )
          from 
            it := list.reverse_cyclic_iterator
              ( pair.item_1, list.cyclic_succ ( i_pred ) )
          until
            it.finished
          loop
            pt4 := points.item ( it.item )
            -- if pt4.y >= pt2.y then
            --   it.stop
            -- elseif 
	    if pt4.side_of_circle (pt1, pt2, pt3) = -1 then
              pt3 := pt4
              next_i := it.item
              it.forth
            else
              it.forth
            end
          end
	end

        if not ok_on_right then
          pt1 := points.item (j)
          pt2 := points.item (i)
          list := neighbours.item (j)
          pair := find_succ ( i, j )
          j_succ := pair.item_1
          next_j := list.item ( j_succ )
          pt3 := points.item ( next_j )
          from 
            it := list.cyclic_iterator
              ( pair.item_1, list.cyclic_pred ( j_succ ) )
          until
            it.finished
          loop
            pt4 := points.item ( it.item )
            -- if pt4.y >= pt2.y then
            --   it.stop
            -- elseif 
	    if pt4.side_of_circle (pt1, pt2, pt3) = -1 then
              pt3 := pt4
              next_j := it.item
              it.forth
            else
              it.forth
            end
          end
        end -- if not ok_on_right

        if next_i /= i and next_j /= j then
          pt3 := points.item ( next_i )
          pt4 := points.item ( next_j )
          if pt3.side_of_circle (pt1, pt2, pt4) = -1 then
            i := next_i
          else
            j := next_j
          end
        elseif next_i /= i then
          i := next_i
        else
          j := next_j
        end

      end

      from
        action_it := actions.iterator
      until
        action_it . finished
      loop
        i := action_it . item . item_2
        j := action_it . item . item_3
	if action_it . item .  item_1 = remove_an_edge then
	  remove_edge ( i, j )
	else
	  add_edge ( i, j )
	end
	action_it . forth
      end
    end


  points : ARRAY [ POINT ]
  neighbours : ARRAY [ GBN_DLIST [ INTEGER ] ]
  left_upper_tangent, right_upper_tangent : INTEGER
  left_lower_tangent, right_lower_tangent : INTEGER

  make ( ary : ARRAY [ POINT ] ) is
    local
      i: INTEGER
      list : GBN_DLIST [ INTEGER ]
    do
      !! points . make ( 1, ary.count )
      !! neighbours . make ( 1, points.count )
      from
        i := 1
      until
        i > points.count
      loop
        !! list.make
        neighbours.put ( list, i )
        points.put ( ary.item (i), i )
        i := i+1
      end
    end

end -- class DELAUNAY