function cvt_4_movie ( n, sample_num, frame_num )

%*****************************************************************************80
%
%% cvt_4_movie() demonstrates how a CVT can be computed and displayed in MATLAB.
%
%  Discussion:
%
%    This simple example carries out an iterative CVT calculation in the 
%    unit square, with a nonuniform density, which is "heaviest" in the
%    corners.  The initial placement of the generators is random.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 December 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of generators.
%
%    integer SAMPLE_NUM, the number of sample points.
%
%    integer FRAME_NUM, the number of frames to generate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'cvt_4_movie():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Create a movie of a CVT computation.\n' );
  fprintf ( 1, '  This example works in the unit square, but includes\n' );
  fprintf ( 1, '  a nonuniform density function which should drive\n' );
  fprintf ( 1, '  the CVT generators towards the corners of the square.\n' );

  if ( nargin < 1 )
    n = 100;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'cvt_4_movie - 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_4_movie - 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
  
  if ( nargin < 3 )
    frame_num = 20;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'cvt_4_movie - Note:\n' );
    fprintf ( 1, '  No value of FRAME_NUM was supplied.\n' );
    fprintf ( 1, '  FRAME_NUM is the number of "frames" of the movie to make\n' );
    fprintf ( 1, '  (and actually the number of iterations, as well).\n' );
    fprintf ( 1, '  A default value FRAME_NUM = %d will be used.\n', frame_num );
  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, or it is not behaving for you, just set this to 0.
%
  accumarray_available = true;

  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 be simulated.\n' );
  end
%
%  This switch controls whether the Delaunay triangulation is displayed.
%
  delaunay_display = false;
  
  if ( delaunay_display )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'The Delaunay triangulation will be displayed.\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'The Delaunay triangulation will NOT be displayed.\n' );
  end
%
%  Randomize the initial locations of the generators.
%
  p = rand ( n, 2 );
%
%  Initialize the movie.
%
  colormap ( 'jet' );

  moviename = 'cvt_4_movie.avi';
  V = VideoWriter ( moviename, 'Motion JPEG AVI' );
  V.FrameRate = 30;
  open ( V );

  figure ( 1 );

  for frame = 1 : frame_num
%
%  Compute the Delaunay triangle information T for the current nodes.
%
    t = delaunay ( p(:,1), p(:,2) );
%
%  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.
%
    clf ( );
    hold ( 'on' );

    if ( delaunay_display )
      trimesh ( t, p(:,1), p(:,2), zeros(n,1) )
    end
%
%  Show the Voronoi cells and generators.
%
    voronoi ( p(:,1), p(:,2), t );
%
%  Finish the plot.
%
    frame_title = sprintf ( 'cvt\\_4\\_movie: Iteration %d\n', frame );
    title ( frame_title );
    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', 2 );
    axis ( 'equal' );
    axis ( 'tight' );
    view ( 2 );
    drawnow ( );
    hold ( 'off' );
%
%  Extract the graphics information as a frame,
%  and add it to the movie.
%
    frame = getframe ( 1 );
    writeVideo ( V, frame );
%
%  Now we update the CVT.
%  Generate sample points in the unit square.
%    
    s = rand ( sample_num, 2 );
%
%  Apply the density function.
%  The first idea here was to use a cube root function, but 
%  MATLAB will take complex roots, so we have to watch out.
%
    s = 0.5 + 0.5 * ( abs ( 2.0 * s - 1.0 ) ).^(0.3) .* ( 2.0 * s - 1.0 ) ...
      ./ abs ( 2.0 * s - 1.0 );
%
%  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) = 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;
      centroid(1:n,2) = 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.
%
    k = find ( count(1:n,1) ~= 0 );
    p(k,1) = centroid(k,1) ./ count(k,1);
    p(k,2) = centroid(k,2) ./ count(k,1);
    
  end

  close ( V );
  fprintf ( 1, '  Movie saved as "%s"\n', moviename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'cvt_4_movie():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end

