function cvt_2_movie ( n, sample_num, frame_num )

%*****************************************************************************80
%
%% cvt_2_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 uniform density.  The initial placement of the
%    generators is random.
%
%    Using the REBIRTHS option, the animation can be extended for great
%    periods of time, as it will continue to evolve.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    10 December 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of generators.
%    A typical value might be 100.
%
%    integer SAMPLE_NUM, the number of sample points.
%    A typical value might be 100 * N.
%
%    integer FRAME_NUM, the number of frames to generate.
%    Defalt is FRAME_NUM = 100.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'cvt_2_movie():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  cvt_2_movie() creates a movie of a CVT generation\n' );
  fprintf ( 1, '  in the unit square.  For excitement, we randomly\n' );
  fprintf ( 1, '  shift one generator on each iteration.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of generators N = %d.\n', n );
  fprintf ( 1, '  The number of sample points SAMPLE_NUM = %d.\n', sample_num );
  fprintf ( 1, '  The number of "frames" FRAME_NUM = %d.\n', frame_num );
%
%  The switch ACCUMARRAY_AVAILABLE is set to 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 switch to false.
%
  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
%
%  The DELAUNAY_DISPLAY 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(1:n,1:2) = 0.1 * rand ( n, 2 );
  p(1:n,1) = p(1:n,1) + 0.25;
  p(1:n,2) = p(1:n,2) + 0.5;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Generators initialized near (0.25,0.5) in unit square.\n' );
%
%  Initialize the movie.
%
  colormap ( 'jet' );
  moviename = 'cvt_2_movie.avi';
  V = VideoWriter ( moviename, 'Motion JPEG AVI' );
  V.FrameRate = 30;
  open ( V );
  figure ( 1 );
%
%  Generate the frames.
%
  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\\_2\\_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' );
    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.
%
%  Pick one generator at random, and relocate it.
%
    rebirth = randi ( [ 1, n ], 1 );
    p(rebirth,1:2) = rand ( 1, 2 );
%
%  Generate sample points.
%    
    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.
%  
    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;
      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_2_movie():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end

