function [vxx,vy] = voronoi(x,y,arg3,arg4)
%VORONOI Voronoi diagram.
%   VORONOI(X,Y) plots the Voronoi diagram for the points X,Y.
%
%   VORONOI(X,Y,TRI) uses the triangulation TRI instead of
%   computing it via DELAUNAY. 
%
%   H = VORONOI(...,'LineSpec') plots the diagram with color and linestyle
%   specified and returns handles to the line objects created in H.
%
%   [VX,VY] = VORONOI(...) returns the vertices of the Voronoi
%   edges in VX and VY so that plot(VX,VY,'-',X,Y,'.') creates the
%   Voronoi diagram.
%
%   See also DELAUNAY, TRIMESH, TRISURF, DSEARCH, CONVHULL.

%   Clay M. Thompson 7-15-95
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.7 $  $Date: 1997/11/21 23:46:58 $

error(nargchk(2,4,nargin));

if nargin==2,
  tri = delaunay(x,y);
  ls = '';
elseif nargin==3,
  if isstr(arg3),
    tri = delaunay(x,y);
    ls = arg3;
  else
    tri = arg3;
    ls = '';
  end
else
  tri = arg3;
  ls = arg4;
end

n = prod(size(x));
ntri = size(tri,1);
t = (1:ntri)';
T = sparse(tri,tri(:,[3 1 2]),t(:,ones(1,3)),n,n); % Triangle edge if T(i,j) 

E = (T & T').*T; % Voronoi edge if E(i,j) 

[i,j,v] = find(triu(E));
[i,j,vv] = find(triu(E'));
c1 = circle(tri(v,:),x,y);
c2 = circle(tri(vv,:),x,y);

vx = [c1(:,1) c2(:,1)].';
vy = [c1(:,2) c2(:,2)].';

if nargout<2
  if isempty(ls),
    co = get(gcf,'defaultaxescolororder');
    h = plot(vx,vy,'-',x,y,'.','color',co(1,:));
  else
    [l,c,m,msg] = colstyle(ls); error(msg)
    if isempty(m), m = '.'; end
    h = plot(vx,vy,ls,x,y,[c m]);
  end
  if ~ishold,
    view(2), axis([min(x(:)) max(x(:)) min(y(:)) max(y(:))])
  end
  if nargout==1, vxx = h; end
else
  vxx = vx;
end

function c = circle(tri,x,y)
%CIRCLE Return center and radius for circumcircles
%   C = CIRCLE(TRI,X,Y) returns a N-by-3 vector containing [xcenter(:)
%   ycenter(:) radius(:)] for each triangle in TRI.

% Reference: Watson, p32.
x = x(:); y = y(:);

x1 = x(tri(:,1)); x2 = x(tri(:,2)); x3 = x(tri(:,3));
y1 = y(tri(:,1)); y2 = y(tri(:,2)); y3 = y(tri(:,3));

% Set equation for center of each circumcircle: 
%    [a11 a12;a21 a22]*[x;y] = [b1;b2] * 0.5;

a11 = x2-x1; a12 = y2-y1;
a21 = x3-x1; a22 = y3-y1;

b1 = a11 .* (x2+x1) + a12 .* (y2+y1);
b2 = a21 .* (x3+x1) + a22 .* (y3+y1);

% Solve the 2-by-2 equation explicitly
idet = a11.*a22 - a21.*a12;

% Add small random displacement to points that are either the same
% or on a line.
d = find(idet == 0);
if ~isempty(d), % Add small random displacement to points
  delta = sqrt(eps);
  x1(d) = x1(d) + delta*(rand(size(d))-0.5);
  x2(d) = x2(d) + delta*(rand(size(d))-0.5);
  x3(d) = x3(d) + delta*(rand(size(d))-0.5);
  y1(d) = y1(d) + delta*(rand(size(d))-0.5);
  y2(d) = y2(d) + delta*(rand(size(d))-0.5);
  y3(d) = y3(d) + delta*(rand(size(d))-0.5);
  a11 = x2-x1; a12 = y2-y1;
  a21 = x3-x1; a22 = y3-y1;
  b1 = a11 .* (x2+x1) + a12 .* (y2+y1);
  b2 = a21 .* (x3+x1) + a22 .* (y3+y1);
  idet = a11.*a22 - a21.*a12;
end

idet = 0.5 ./ idet;

xcenter = ( a22.*b1 - a12.*b2) .* idet;
ycenter = (-a21.*b1 + a11.*b2) .* idet;

radius = (x1-xcenter).^2 + (y1-ycenter).^2;

c = [xcenter ycenter radius];


