function cvt_metric ( n, sample_num, metric, filename )

%*****************************************************************************80
%
%% cvt_metric() computes a CVT in a square with a variable metric function.
%
%  Discussion:
%
%    This program computes a CVT, but assumes that distances are determined
%    by a metric function:
%
%      d(x,y) = sqrt ( (x-y)' * A * (x-y) )
%
%    where A is a symmetric positive definite matrix.  We will NOT assume
%    the matrix A is constant.  Instead, we will assume that it varies
%    spatially.  For simplicity, when computing d(x,y), we will evaluate
%    A at the midpoint (x+y)/2.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 September 2024
%
%  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 generators.
%
%    integer SAMPLE_NUM, the number of sample points.
%
%    A = metric(p): returns the 2x2 metric matrix at point p.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'cvt_metric():\n' );
  fprintf ( 1, '  A simple demonstration of a CVT computation\n' );
  fprintf ( 1, '  (Centroidal Voronoi Tessellation)\n' );
  fprintf ( 1, '  in a square, using a metric matrix.\n' );
  
  if ( nargin < 1 )
    n = 100;
    fprintf ( 1, '\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 );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User specified number of generators = %d\n',  n );
  end

  if ( nargin < 2 )
    sample_num = 1000 * n;
    fprintf ( 1, '\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 );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User specified number of sample points = %d\n', ...
      sample_num );
  end

  if ( nargin < 3 )
    metric = @ metric_04;
  end

  if ( nargin < 4 )
    filename = 'cvt_metric.png';
  end
%
%  This switch is set to 1 (TRUE) if the accumarray() command is available.
%  That speeds up the calculation a lot.  If you don't have the accumarray()
%  command, just set this to 0.
% 
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The accumarray() command can be used for faster\n' );
  fprintf ( 1, '  computation.  This command is not available in\n' );
  fprintf ( 1, '  some versions of MATLAB.  If accumarray() is available,\n' );
  fprintf ( 1, '  simply make sure that the accumarray_available variable\n' );
  fprintf ( 1, '  is set to 1!\n' );
  
  accumarray_available = 1;

  if ( accumarray_available )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The accumarray() command will be used.\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The accumarray() command will NOT be used.\n' );
  end
   
  for it = 0 : 15

    if ( it == 0 )

      p = rand ( n, 2 );

    else
%
%  Generate sample points.  
%
      ps = rand ( sample_num, 2 );
%
%  For each sample point, find K, the nearest generator.
%
      k(1:sample_num,1) = 0;
      d(1:sample_num,1) = Inf;
        
      for g = 1 : n
%
%  For some asine reason, can't say p(g,1:2) here...
%
        p_matlab_idiots = p(g,1:2);
        A = metric ( p_matlab_idiots );
      
        for s = 1 : sample_num
          d2 = ( ps(s,1:2) - p(g,1:2) ) * A * ( ps(s,1:2) - p(g,1:2) )';
          if ( d2 < d(s,1) )
            k(s,1) = g;
            d(s,1) = d2;
          end
        end

      end
%
%  The centroid of the Voronoi region associated with each generator
%  is approximated by the average of the sample points it was closest to.
%
      if ( accumarray_available )

        count(1:n) = accumarray ( k, ones(sample_num,1) );
        centroid(1,1:n) = accumarray ( k, ps(:,1) );
        centroid(2,1:n) = accumarray ( k, ps(:,2) );

      else

        count(1:n) = 0;
        centroid(1,1:n) = 0;
        centroid(2,1:n) = 0;

        for i = 1 : sample_num
          j = k(i);
          count(j) = count(j) + 1;
          centroid(1,j) = centroid(1,j) + ps(i,1);
          centroid(2,j) = centroid(2,j) + ps(i,2);
        end

      end
%
%  Replace the generators by the centroids.
%
      p(1:n,1) = ( centroid(1,1:n) ./ count(1:n) )';
      p(1:n,2) = ( centroid(2,1:n) ./ count(1:n) )';

    end
%
%  Display the CVT generators, and the associated Voronoi diagram.
%
    clf ( );
    scatter ( p(:,1), p(:,2), [], 'b', 'filled' );
    axis ( [ -0.1, 1.1, -0.1, 1.1 ] );
    line ( [ 0.0, 1.0, 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0, 1.0, 0.0 ], ...
      'Color', 'r', 'linewidth', 3 );
    title_string = sprintf ( 'Points, step %d', it );
    title ( title_string );
    axis ( 'equal' );
    drawnow ( );
    
  end
%
%  Save last configuration.
%
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end
