#! /usr/bin/env python
#
def svd_image ( ):

#*****************************************************************************80
#
## svd_image uses SVD to carry out PCA on an image.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    24 June 2019
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np
  import platform
  from PIL import Image

  print ( '' )
  print ( 'svd_image:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Read image data as a matrix A.' )
  print ( '  Compute SVD of A.' )
  print ( '  Plot the singular values, and cumulative singular values.' )
  print ( '  Plot 10 and 20 component approximations to A.' )
#
#  Read the image data as a uint8 array.
#
  I = Image.open ( 'casablanca.png' )
#
#  Convert to grayscale mode.
#
  I = I.convert ( 'LA' )
#
#  Display the image.
#
  plt.imshow ( I )
  plt.title ( 'Original image:' )
  plt.show ( )
#
#  Make a real number copy of A and get the dimensions.
#
  A = np.array ( list ( I.getdata(band=0)), float )
  A.shape = ( I.size[1], I.size[0] )
  m, n = A.shape
#
#  Display the numeric version of A.
#
  plt.imshow ( A, cmap = 'gray' )
  plt.title ( 'Numeric version of image' )
  plt.show ( )
#
#  Subtract off column averages.
#
  u0 = np.ones ( m )
  v0 = np.mean ( A, axis = 0 )
  A0 = np.outer ( u0, v0 )
#
#  Compute the SVD of A - A0.
#
  U, s, V = np.linalg.svd ( A - A0 )
  S = np.zeros ( A.shape )
  for i in range ( 0, min ( m, n ) ):
    S[i,i] = s[i]
#
#  Plot the first 25 squared singular values.
#
  x = np.arange ( 25 )
  plt.plot ( x, s[0:25]**2, 'bo-', linewidth = 2 )
  plt.grid ( True )
  plt.xlabel ( '<-- I: index -->' )
  plt.ylabel ( '<-- Square of I-th Singular Value -->' )
  plt.title ( 'Squares of singular values for image columns' )
  filename = 'svd_image_1.png'
  plt.savefig ( filename )
  print ( '  Graphics saved in "%s"' % ( filename ) )
  plt.show ( )
  plt.clf ( )
#
#  Display the relative sum of the first K squared singular values.
#
  x = np.arange ( 26 )
  y = [ 0.0 ]
  y = np.append ( y, np.cumsum ( s[0:25]**2 ) )
  y = y / np.sum ( s**2 )
  plt.plot ( x, y, 'ro-', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<-- K: Number of components -->' )
  plt.ylabel ( '<-- Amount of variance explained -->' )
  plt.title ( 'Information proportion in first K components' )
  filename = 'svd_image_2.png'
  plt.savefig ( filename )
  print ( '  Graphics saved in "%s"' % ( filename ) )
  plt.show ( )
  plt.clf ( )
#
#  Display image approximation with 5 components
#
  r = 5
  A5 = A0 + np.matmul ( np.matmul ( U[:,0:r], S[0:r,0:r] ), V[0:r,:] )

  plt.imshow ( A5, cmap = 'gray' )
  plt.title ( 'Reconstruction using %d components' % ( r ) )
  filename = 'svd_image_3.png'
  plt.savefig ( filename )
  print ( '  Graphics saved in "%s"' % ( filename ) )
  plt.show ( )
  plt.clf ( )
#
#  Display image approximation with 10 components
#
  r = 10
  A10 = A0 + np.matmul ( np.matmul ( U[:,0:r], S[0:r,0:r] ), V[0:r,:] )

  plt.imshow ( A10, cmap = 'gray' )
  plt.title ( 'Reconstruction using %d components' % ( r ) )
  filename = 'svd_image_4.png'
  plt.savefig ( filename )
  print ( '  Graphics saved in "%s"' % ( filename ) )
  plt.show ( )
  plt.clf ( )
#
#  Display image approximation with 20 components
#
  r = 20
  A20 = A0 + np.matmul ( np.matmul ( U[:,0:r], S[0:r,0:r] ), V[0:r,:] )

  plt.imshow ( A20, cmap = 'gray' )
  plt.title ( 'Reconstruction using %d components' % ( r ) )
  filename = 'svd_image_5.png'
  plt.savefig ( filename )
  print ( '  Graphics saved in "%s"' % ( filename ) )
  plt.show ( )
  plt.clf ( )
#
#  Display image approximation with 40 components
#
  r = 40
  A40 = A0 + np.matmul ( np.matmul ( U[:,0:r], S[0:r,0:r] ), V[0:r,:] )

  plt.imshow ( A40, cmap = 'gray' )
  plt.title ( 'Reconstruction using %d components' % ( r ) )
  filename = 'svd_image_6.png'
  plt.savefig ( filename )
  print ( '  Graphics saved in "%s"' % ( filename ) )
  plt.show ( )
  plt.clf ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'svd_image:' )
  print ( '  Normal end of execution.' )

  return

if ( __name__ == '__main__' ):
  svd_image ( )