function [ r, diff, energy ] = cvt_iterate ( n, r, ratio )

%*****************************************************************************80
%
%% cvt_iterate() takes one step of the CVT iteration.
%
%  Discussion:
%
%    The routine is given a set of points, called "generators", which
%    define a tessellation of the region into Voronoi cells.  Each point
%    defines a cell.  Each cell, in turn, has a centroid, but it is
%    unlikely that the centroid and the generator coincide.
%
%    Each time this CVT iteration is carried out, an attempt is made
%    to modify the generators in such a way that they are closer and
%    closer to being the centroids of the Voronoi cells they generate.
%
%    A large number of sample points are generated, and the nearest generator
%    is determined.  A count is kept of how many points were nearest to each
%    generator.  Once the sampling is completed, the location of all the
%    generators is adjusted.  This step should decrease the discrepancy
%    between the generators and the centroids.
%
%    The centroidal Voronoi tessellation minimizes the "energy",
%    defined to be the integral, over the region, of the square of
%    the distance between each point in the region and its nearest generator.
%    The sampling technique supplies a discrete estimate of this
%    energy.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 December 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Qiang Du, Vance Faber, and Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review, Volume 41, 1999, pages 637-676.
%
%  Input:
%
%    integer N, the number of Voronoi cells.
%
%    real R(2,N), the Voronoi cell generators.
%
%    integer RATIO, the number of sample points per generator.
%
%  Output:
%
%    real R(NDIM,N), the updated Voronoi cell generators.
%
%    real DIFF, the L2 norm of the difference between the iterates.
%
%    real ENERGY, the discrete "energy", divided by the number 
%    of sample points.
%
  ndim = 2;
%
%  Generate the sampling points S.
%
  sample_num = ratio * n;
  s = rand ( 2, sample_num );
%
%  Find the index N of the nearest cell generator to each sample point S.
%
  nearest = find_closest ( ndim, n, sample_num, s, r );
%
%  Add S to the centroid associated with generator N.
%
  r2 = zeros ( 2, n );
  energy = 0.0;
  count = zeros ( n, 1 );
  for j = 1 : sample_num
    r2(1:ndim,nearest(j)) = r2(1:ndim,nearest(j)) + s(1:ndim,j);
    energy = energy + sum ( ( r(1:ndim,nearest(j)) - s(1:ndim,j) ).^2 );
    count(nearest(j)) = count(nearest(j)) + 1;
  end

  energy = energy / sample_num;
%
%  Estimate the centroids.
%
  for j = 1 : n
    if ( 0 < count(j) )
      r2(1:ndim,j) = r2(1:ndim,j) / count(j);
    end
  end
%
%  Determine the sum of the distances between generators and centroids.
%
  diff = 0.0;
  for j = 1 : n
    diff = diff + sqrt ( sum ( ( r2(1:ndim,j) - r(1:ndim,j) ).^2 ) );
  end
%
%  Replace the generators by the centroids.
%
  r(1:ndim,1:n) = r2(1:ndim,1:n);
%
%  Normalize the discrete energy estimate.
%
  energy = energy / sample_num;

  return
end

