function [ p, t ] = cvt_p08_uniform ( n, sample_num )

%*****************************************************************************80
%
%% cvt_p08_uniform() demonstrates how a CVT can be computed and displayed in MATLAB.
%
%  Discussion:
%
%    This simple example carries out an iterative CVT calculation in the 
%    Persson and Strang's pie shaped region, with a uniform density.  
%
%    The initial placement of the generators is random.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of generators.
%
%    integer SAMPLE_num, the number of sample points.
%
%  Output:
%
%    real P(N,2), the location of the generators.
%
%    integer T(NT,3), information defining the Delaunay
%    triangulation of the generators.  NT is the number of triangles,
%    which varies depending on the arrangement of the generators.
%
  if ( nargin < 1 )
    n = 100;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'cvt_p08_uniform() - Note:\n' );
    fprintf ( 1, '  No value of N was supplied.\n' );
    fprintf ( 1, '  N is the number of generators.\n' );
    fprintf ( 1, '  A default value N = %d will be used.\n', n );
  end

  if ( nargin < 2 )
    sample_num = 1000 * n;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CVT_P08_UNIFORM - Note:\n' );
    fprintf ( 1, '  No value of SAMPLE_num was supplied.\n' );
    fprintf ( 1, '  SAMPLE_num is the number of sample points.\n' );
    fprintf ( 1, '  A default value SAMPLE_num = %d will be used.\n', sample_num );
  end
%
%  Randomize the initial locations of the generators.
%
  p = p08_sample ( 2, n );
  p = p';
 
  while ( true )
%
%  Compute the Delaunay triangle information T for the current nodes.
%
    t = delaunayn ( p );
%
%  Display a graph of the generators.
%
%  This was a little tricky.  The TRIMESH command must come before
%  the VORONOI command in order for both items to appear.
%
%  There are obvious discrepancies between the Delaunay and Voronoi
%  graphs which suggest that there are mistakes in the Voronoi diagram
%  as computed by MATLAB's VORONOI command!
%
    pmid = ( p(t(:,1),:) + p(t(:,2),:) + p(t(:,3),:) ) / 3; 
    nt = size ( t, 1 );
    t = t( p08_inside ( 2, nt, pmid' ) == 1, : );
    
    trimesh ( t, p(:,1), p(:,2), zeros(n,1) );
    view(2);
    axis ( 'equal' );
    axis ( 'off' );
    drawnow
%
%  Generate sample points.
%    
    s = p08_sample ( 2, sample_num );
    xs = s(1,1:sample_num);
    ys = s(2,1:sample_num);
%
%  For each sample point, find K, the index of the nearest generator.
%  We do this efficiently by using the Delaunay information.
%  
    k(1:sample_num,1) = dsearch ( p(:,1), p(:,2), t, xs, ys );
%
%  The centroid of the Voronoi region associated with each generator
%  is approximated by the average of the sample points it was closest to.
%
%  The ACCUMARRAY command accumulates (adds) all the appropriate entries.
%    
    count(1:n) = accumarray ( k, ones(sample_num,1) );
    centroid(1,1:n) = accumarray ( k, xs );
    centroid(2,1:n) = accumarray ( k, ys );
%
%  Replace the generators by the centroids.
%
    count(1:n) = max ( count(1:n), 1 );
    p(1:n,1) = centroid(1,1:n) ./ count(1:n);
    p(1:n,2) = centroid(2,1:n) ./ count(1:n);

    string = input ( 'RETURN, or Q to quit: ', 's' );

    if ( string == 'q' | string == 'Q' )
      break
    end
    
  end

  return
end
