function cvt_basic ( )

%*****************************************************************************80
%
%% cvt_basic() computes a CVT in a square, with uniform density.
%
%  Discussion:
%
%    This is a "basic" version of cvt_square_uniform(), in which all the
%    choices have been made for the user in advance.  This code may be
%    used to get familiar with the computation of a Centroidal Voronoi
%    Tesselation (CVT).
%
%    This simple example carries out an iterative procedure for estimating
%    the proper locations of N points in a square that will satisfy the
%    Centroidal Voronoi Tessellation (CVT) condition, namely, that they
%    are the centroids of their Voronoi regions.
%
%    The initial placement of the generators is random.
%
%    The desired density of points is uniform.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 April 2022
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Qiang Du, Vance Faber, Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review, 
%    Volume 41, 1999, pages 637-676.
%
%  Local:
%
%    integer N, the number of generators.
%
%    integer IT_NUM, the number of iterations.
%
%    integer SAMPLE_NUM, the number of sample points.
%
%    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.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'cvt_basic():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  A simple demonstration of a CVT computation\n' );
  fprintf ( 1, '  (Centroidal Voronoi Tessellation)\n' );
  fprintf ( 1, '  in a square, with a uniform density.\n' );
  
  n = 100;
  it_num = 50;
  sample_num = 1000 * n;
%
%  Randomize the initial locations of the generators in the unit square.
%  If another region is used, then this initialization should be changed.
%
  p = square_uniform ( n );
  
  for it = 0 : it_num
%
%  Compute the Delaunay triangle information T for the current generators.
%
    t = delaunay ( p(:,1), p(:,2) );
%
%  Generate sample points.  
%
    s = square_uniform ( sample_num );
%
%  Find K, the index of the nearest generator.
% 
    k = dsearchn ( p, t, s );
%
%  Compute centroid of each generator's region.
%
    count(1:n,1) = accumarray ( k, ones(sample_num,1) );
    centroid(1:n,1) = accumarray ( k, s(:,1) );
    centroid(1:n,2) = accumarray ( k, s(:,2) );
%
%  Replace the generators by the centroids.
%
    p(1:n,1) = centroid(1:n,1) ./ count(1:n,1);
    p(1:n,2) = centroid(1:n,2) ./ count(1:n,1);

  end
%
%  Plot the CVT points.
%
  clf ( );
  plot ( p(:,1), p(:,2), '.', 'markersize', 25 );
  axis ( [ -0.0, 1.0, -0.0, 1.0 ] );
  line ( [ 0.0, 1.0, 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0, 1.0, 0.0 ], ...
    'Color', 'r', 'LineWidth', 2 );
  title_string = sprintf ( 'CVT generators' );
  title ( title_string );
  axis ( 'equal' );
  axis ( 'tight' );
  drawnow ( );
  filename = 'cvt_basic.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  close ( );
%
%  Add the "voronoi()" command to see the boundaries around the points.
%
  clf ( );
  hold ( 'on' );
    voronoi ( p(:,1), p(:,2), t );
    plot ( p(:,1), p(:,2), '.', 'markersize', 25 );
    axis ( [ -0.0, 1.0, -0.0, 1.0 ] );
    line ( [ 0.0, 1.0, 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0, 1.0, 0.0 ], ...
      'Color', 'r', 'LineWidth', 2 );
  hold ( 'off' );
  xlim ( [ 0.0, 1.0 ] );
  ylim ( [ 0.0, 1.0 ] );
  title_string = sprintf ( 'CVT generators + Voronoi boundaries' );
  title ( title_string );
  axis ( 'equal' );
% axis ( 'tight' );
  drawnow ( );
  filename = 'cvt_basic_voronoi.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  close ( );

  return
end
function p = square_uniform ( n )

%*****************************************************************************80
%
%% square_uniform() returns sample points from the unit square.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 April 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N: the number of points desired.
%
%  Output:
%
%    real P(N,2): N random points uniformly distributed in the unit square.
%
  p = rand ( n, 2 );

  return
end
    
