function vqplot2d(C,R)
% vqplot2d(C,R)  Plot Voronoi regions of VQ codebook in 2-D
%     C defines a codebook.  R is a four element vector of the limits 
%     of the plot.
% 2001-03-19 dpwe@ee.columbia.edu

if nargin < 2
  R = [-10 10 -10 10];
end

% Maximum number of line segments to search
maxnls = 100;

xmin = R(1);
xmax = R(2);
ymin = R(3);
ymax = R(4);

[ncw,ndim] = size(C);

for cw = 1:ncw
  
  x = C(cw,:);
  
  % Plot the centroid
  hold on
  plot(x(1), x(2), '.r');
  hold off
  
  % Sort remaining codewords by distance
  d = sum( (((ones(ncw,1)*x) - C).^2)' );
  [dummy, ix] = sort(d);
  
  % Figure line segments for each one
  
  % How many line segments?
  nls = min(maxnls, ncw-1);

  mm = zeros(nls, ndim);
  vv = zeros(nls, ndim);
  
  for ocw = 1:nls
    
    y = C(ix(ocw+1),:);
    % Dividing line as midpoint + dy/dx
    mp = (x+y)/2;
    vp = x-y;
    vp = [vp(2), -vp(1)];
    
    mm(ocw,:) = mp;
    vv(ocw,:) = vp;

  end
  
  % For each line segment, figure intersects with all others & plot
  
  for ls = 1:nls
  
    mp = mm(ls,:);
    vp = vv(ls,:);
    
    % Boundaries of figure space
    if (vp(1) > 0)
      amin = (xmin - mp(1)) / vp(1);
      amax = (xmax - mp(1)) / vp(1);
    else
      amax = (xmin - mp(1)) / vp(1);
      amin = (xmax - mp(1)) / vp(1);
    end
    if (vp(2) > 0)
      amin = max( amin, (ymin - mp(2)) / vp(2) );
      amax = min( amax, (ymax - mp(2)) / vp(2) );
    else
      amax = min( amax, (ymin - mp(2)) / vp(2) );
      amin = max( amin, (ymax - mp(2)) / vp(2) );
    end
    
    % find intersects with other midpoint lines
    ls2 = 1;
    while (amin < amax) & (ls2 <= nls)
  
      if ls2 ~= ls

	mo = mm(ls2,:);
	vo = vv(ls2,:);

	aa = inv([vp',vo'])*((mp-mo)');
      
	aa = -aa(1);
	
	% Which side is it on?
	if  (mo - mp)*(vp') > 0
	  if aa < amax
	    amax = aa;
	  end
	else
	  if aa > amin
	    amin = aa;
	  end
	end
    
      end
      
      ls2 = ls2+1;
    end
        
    % Add the line segment
    if amin < amax
      p1 = mp+amin*vp;
      p2 = mp+amax*vp;
      hold on
      plot([p1(1) p2(1)], [p1(2) p2(2)], '-g');
      hold off
    end
  
  end
  
end