function svd_gray ( gray_filename )

%*****************************************************************************80
%
%% SVD_GRAY uses SVD to make compressed versions of a gray scale image.
%
%  Discussion:
%
%    The user supplies the name of a file containing an image.
%
%    This is presumably a black and white MxN image.
%
%    The image data is converted to numeric form, and the SVD decomposition 
%    is computed.  It is possible to generate approximate versions of the
%    image for values 1 <= R <= min ( M, N ).
%
%    The program displays versions of the image for various values of R.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 December 2018
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string GRAY_FILENAME, the name of the image file.
%
  timestamp ( );
  fprintf ( 1, 'SVD_GRAY\n' );
  fprintf ( 1, '  Demonstrate SVD image compression on a gray scale image.\n' );
  fprintf ( 1, '  Image consists of M rows and N columns.\n' );
  fprintf ( 1, '  Compute some rank-K SVD approximations.\n' );

  a = imread ( gray_filename );
%
%  Get sizes.
%
  [ m, n ] = size ( a );
%
%  Report the image dimensions.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Image rows M =     %d\n', m );
  fprintf ( 1, '  Image columns N =  %d\n', n );
  fprintf ( 1, '  Image pixels M*N = %d\n', m * n );
%
%  Display the original.
%
  fig_num = 0

  fig_num = fig_num + 1;
  figure ( fig_num )
  imshow ( a )
  title ( 'Original image', 'Fontsize', 16 );
  filename = sprintf ( 'svd_gray_original.png' );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Convert the image to double.
%
  a = double ( a );
%
%  Compute the SVD.
%
  [ u, s, v ] = svd ( a );
%
%  Compute the sum of the singular values.
%
  s_vec = diag ( s );
  s_sum = sum ( s_vec );
%
%  Plot the eigenvalues.
%
  fig_num = fig_num + 1
  figure ( fig_num )
  x = 1 : min ( m, n );
  y = diag ( s );
  plot ( x, y, 'Linewidth', 4 );
  grid
  title ( 'Singular values', 'Fontsize', 16 )
  filename = 'svd_gray_singular_values.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Create and display a sequence of low rank approximations.
%
  r_vec = [ 1, 2, 5, 10, 20, 40, 80 ];
  r_num = length ( r_vec );

  for i = 1 : r_num

    r = r_vec(i);

    if ( min ( m, n ) < r ) then
      break
    end

    fig_num = fig_num + 1;
    figure ( fig_num );
    b = svd_gray_approximate ( m, n, r, u, s, v );
    imshow ( b );
    comp = ( m * r + r + r * n ) / ( m * n );
    psum = sum ( s_vec(1:r) ) / s_sum;
    title ( sprintf ( 'Image approximated with %d singular vectors.', r ), 'Fontsize', 16 )
    fprintf ( 1, '  Compression = %f, Partial sum = %f\n', comp, psum );
    filename = sprintf ( 'svd_gray_approximation_%d.png', r );
    print ( '-dpng', filename );
    fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SVD_GRAY\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
