function [ p, t ] = cvt_example_five_points ( )

%*****************************************************************************80
%
%% CVT_EXAMPLE_FIVE_POINTS demonstrates a CVT on five points in the unit square.
%
%  Discussion:
%
%    This simple example carries out an iterative CVT calculation in the 
%    unit square, with a uniform density, on a particular set of 5 points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 June 2016
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Qiang Du, Vance Faber, Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review, 
%    Volume 41, 1999, pages 637-676.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_EXAMPLE_FIVE_POINTS:\n' );
  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, and five points.\n' );
  
  n = 5;

  sample_num = 5000;

  delaunay_display = 0;
%
%  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.
% 
  accumarray_available = 1;
%
%  Clear the figure screen, if already open.
%
  clf
%
%  Specify the five points.
%
  p = [ ...
    0.7,  0.2; ...
    0.4,  0.4; ...
    0.3,  0.3; ...
    0.5,  0.7; ...
    0.8,  0.5 ];

  plot ( p(:,1), p(:,2), 'b.', 'MarkerSize', 30 );
  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', 'k', 'LineWidth', 2 );
  title_string = sprintf ( 'Initial Generators' );
  title ( title_string );
  axis equal
  axis off
  drawnow
  filename = 'cvt_five_00.png';
  print ( '-dpng', filename );

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

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

  it = 0;
  
  while ( 1 )
%
%  Compute the Delaunay triangle information T for the current nodes.
%
    t = delaunay ( p(:,1), p(:,2) );
%
%  Display the Delaunay triangulation, if requested.
%
    if ( delaunay_display )
      subplot ( 1, 2, 2 )
      trimesh ( t, p(:,1), p(:,2), zeros(n,1), 'EdgeColor', 'b' )
      hold on
      plot ( p(:,1), p(:,2), 'b.' );
      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', 'k', 'LineWidth', 2 );
      title_string = sprintf ( 'Delaunay, step %d', it );
      title ( title_string );
      axis equal
      axis off
      view ( 2 )
      hold off
    end
%
%  Display the CVT generators, and the associated Voronoi diagram.
%
    if ( delaunay_display )
      subplot ( 1, 2, 1 )
    end
    
    H = voronoi ( p(:,1), p(:,2), t );
    set ( H(1), 'MarkerSize', 30 )
    set ( H(2), 'LineWidth', 3 )
    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', 'k', 'LineWidth', 2 );
    title_string = sprintf ( 'Voronoi, step %d', it );
    title ( title_string );
    axis equal
    axis off
    axis tight
    drawnow
    filename = filename_inc ( filename );
    print ( '-dpng', filename );
%
%  Generate sample points.  
%
%  These sample points implicitly define the geometry of the region.  
%  If the region is not a unit square, then the range of the sample 
%  data must be changed.
%   
    s = rand ( sample_num, 2 );
%
%  For each sample point, find K, the index of the nearest generator.
%  We do this efficiently by using the Delaunay information with
%  Matlab's DSEARCHN command, rather than a brute force nearest neighbor
%  computation.
% 
    k = dsearchn ( p, t, s );
%
%  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,1) = accumarray ( k, ones(sample_num,1) );
      centroid(1:n,1) = accumarray ( k, s(:,1) );
      centroid(1:n,2) = accumarray ( k, s(:,2) );

    else

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

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

    end
%
%  Replace the generators by the centroids.
%
    pnew(1:n,1) = centroid(1:n,1) ./ count(1:n);
    pnew(1:n,2) = centroid(1:n,2) ./ count(1:n);
%
%  Display the NEW generators, and the OLD Voronoi diagram.
%   
    pause

    clf
    hold on
    H = voronoi ( p(:,1), p(:,2), t );
    set ( H(1), 'MarkerSize', 30 )
    set ( H(2), 'LineWidth', 3 )
    plot ( pnew(:,1), pnew(:,2), 'r.', 'MarkerSize', 30 )
    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', 'k', 'LineWidth', 2 );
    title_string = sprintf ( 'Old and New Generators, step %d', it );
    title ( title_string );
    axis equal
    axis tight
    axis off
    drawnow
    hold off
    filename = filename_inc ( filename );
    print ( '-dpng', filename );

    p = pnew;

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

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

    it = it + 1;
    
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_EXAMPLE_FIVE_POINTS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function filename = filename_inc ( filename )

%*****************************************************************************80
%
%% FILENAME_INC increments a partially numeric filename.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP!         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILENAME, the string to be incremented.
%
%    Output, string FILENAME, the incremented string.
%
  lens = length ( filename );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILENAME_INC - Fatal error!\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILENAME_INC - Fatal error!' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = filename(i);

    if ( '0' <= c && c <= '8' )

      c = c + 1;
      
      filename(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';
      
      filename(i) = c;

    end

  end
%
%  No digits were found.  Return blank.
%
  if ( change == 0 )
    filename = ' ';
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

