# include <math.h>
# include <stdio.h>
# include <time.h>

# include "ppc_array.h"
# include "ppc_haar.h"
# include "ppc_image.h"
# include "ppc_image_analysis.h"
# include "ppc_xmalloc.h"

/******************************************************************************/

void clip_matrix ( double **a, int m, int n, int M )

/******************************************************************************/
/*
  Purpose:

    clip_matrix() forces matrix entries to be integers in a given range.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    16 September 2023

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double **a: the m x n matrix.

    int m, n: the number of rows and columns.

    int M: the maximum value allowed for the values in the file.

  Output:

    double **a: the clipped matrix.
*/
{
  int i;
  int j;

  for ( i = 0; i < m; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      if ( a[i][j] < 0 )
      {
        a[i][j] = 0;
      }
      else if ( M < a[i][j] )
      {
        a[i][j] = M;
      }
      else
      {
        a[i][j] = lrint ( a[i][j] );
      }
    }
  }
  return;
}
/******************************************************************************/

double matrix_norm_frobenius ( double **a, int m, int n )

/******************************************************************************/
/*
  Purpose:

    matrix_norm_frobenius() evaluates the Frobenius norm of a matrix.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 September 2023

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double **a: the m x n matrix.

    int m, n: the number of rows and columns.

  Output:

    double matrix_norm_frobenius: the Frobenius norm of the matrix.
*/
{
  int i;
  int j;
  double norm;
  double s;

  norm = 0.0;

  for ( i = 0; i < m; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      s = a[i][j];
      norm = norm + s * s;
    }
  }
  norm = sqrt ( norm );

  return norm;
}
/******************************************************************************/

int prune_pgm ( double **a, int m, int n, double rel_err )

/******************************************************************************/
/*
  Purpose:

    prune_pgm() drops relatively small entries from the matrix.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 September 2023

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double a[m,n]: the image to be pruned.

    int m, n: the number of rows and columns in the image.

    double relerr: the relative error tolerance, between 0 and 1.

  Output:

    double a[m,n]: the pruned image.

    int deleted: the number of entries deleted from the image.
*/
{
  double abs_err2;
  int deleted;
  double e;
  int i;
  int j;
  double max;
  double min;
  double norm;
  double s;
  double t;
/*
  Find max and min of matrix entries.
*/
  s = fabs ( a[0][0] );
  max = s;
  min = s;

  for ( i = 0; i < m; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      s = fabs ( a[i][j] );
      if ( max < s )
      {
        max = s;
      }
      if ( s < min )
      {
        min = s;
      }
    }
  }
/*
  Find the Frobenius norm of matrix.
*/
  norm = matrix_norm_frobenius ( a, m, n );

  printf ( "\n" );
  printf ( "  min matrix = %g\n", min );
  printf ( "  max matrix = %g\n", max );
/*
  Determine the cutoff threshhold t.
*/
  abs_err2 = pow ( norm, 2 ) * rel_err;
  printf ( "  norm^2 = %g\n", pow ( norm, 2 ) );
  printf ( "  rel_err = %g\n", rel_err );
  printf ( "  abs_err2 = %g\n", abs_err2 );

  while ( 0.001 * norm <= max - min )
  {
    t = ( min + max ) / 2.0;
    e = truncation_error ( t, a, m, n );
    if ( abs_err2 < e )
    {
      max = t;
    }
    else
    {
      min = t;
    }
    printf ( "min = %g, max = %g, e = %g\n", min, max, e );
  }
/*
  Zero out matrix entries below t.
*/
  deleted = 0;
  for ( i = 0; i < m; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      if ( fabs ( a[i][j] ) < t )
      {
        a[i][j] = 0.0;
        deleted = deleted + 1;
      }
    }
  }
  return deleted;
}
/******************************************************************************/

int prune_ppm ( double **r, double **g, double **b, int m, int n, 
  double rel_err )

/******************************************************************************/
/*
  Purpose:

    prune_ppm() drops relatively small entries from a PPM image

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    23 April 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double r[m,n], g[m,n], b[m,n]: the RGB components of the image to be pruned.

    int m, n: the number of rows and columns in the image.

    double relerr: the relative error tolerance, between 0 and 1.

  Output:

    double r[m,n], g[m,n], b[m,n]: the pruned image components.

    int deleted: the number of entries deleted from the image.
*/
{
  double abs_err2;
  int deleted;
  double e;
  int i;
  int j;
  double max;
  double min;
  double norm;
  double **RGB;
  double s;
  double t;
  double v;

  make_matrix ( RGB, m, n );

  norm = 0.0;

  for ( i = 0; i < m; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      v = r[i][j] * r[i][j] 
        + g[i][j] * g[i][j]
        + b[i][j] * b[i][j];

      norm = norm + v;

      RGB[i][j] = v;
    }
  }
/*
  Find max and min of RGB component vectors, and Frobenius norm of RGB matrix.
*/
  s = fabs ( RGB[0][0] );
  max = s;
  min = s;

  for ( i = 0; i < m; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      s = fabs ( RGB[i][j] );

      if ( max < s )
      {
        max = s;
      }
      if ( s < min )
      {
        min = s;
      }
    }
  }

  printf ( "\n" );
  printf ( "  min matrix = %g\n", min );
  printf ( "  max matrix = %g\n", max );
/*
  Determine the cutoff threshhold t.
*/
  abs_err2 = pow ( norm, 2 ) * rel_err;
  printf ( "  norm^2 = %g\n", pow ( norm, 2 ) );
  printf ( "  rel_err = %g\n", rel_err );
  printf ( "  abs_err2 = %g\n", abs_err2 );

  while ( 0.001 * norm <= max - min )
  {
    t = ( min + max ) / 2.0;
    e = truncation_error ( t, RGB, m, n );
    if ( abs_err2 < e )
    {
      max = t;
    }
    else
    {
      min = t;
    }
    printf ( "min = %g, max = %g, e = %g\n", min, max, e );
  }
/*
  Zero out matrix entries below t.
*/
  deleted = 0;
  for ( i = 0; i < m; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      if ( RGB[i][j] < t )
      {
        r[i][j] = 0.0;
        g[i][j] = 0.0;
        b[i][j] = 0.0;
        deleted = deleted + 1;
      }
    }
  }

  free_matrix ( RGB );

  return deleted;
}
/******************************************************************************/

void reduce_pgm_image ( struct image *img, double rel_err )

/******************************************************************************/
/*
  Purpose:

    reduce_pgm_image() reduces PGM image size using Haar transforms.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 April 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct image *img: the PGM image to be modified.

    double rel_err: the relative error tolerance, 0.0 <= relerr <= 1.0.

  Output:

    struct image *img: the modified image.
*/
{
  int m = img->pam.height;
  int n = img->pam.width;
  int M = img->pam.maxval;
  int zero_count;

  haar_transform_matrix ( img->g, m, n, WT_FWD );
  zero_count = prune_pgm ( img->g, m, n, rel_err );
  haar_transform_matrix ( img->g, m, n, WT_REV );
  clip_matrix ( img->g, m, n, M );

  fprintf ( stderr, "\n" );
  fprintf ( stderr, "reduce_pgm_image():\n" );
  fprintf ( stderr, "  zeroed %d of %d wavelet coefficients\n", zero_count, m*n );
  fprintf ( stderr, "  %d remaining.\n", m * n - zero_count );

  return;
}
/******************************************************************************/

void reduce_ppm_image ( struct image *img, double rel_err )

/******************************************************************************/
/*
  Purpose:

    reduce_ppm_image() reduces PPM image size using Haar transforms.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    23 April 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct image *img: the PPM image to be modified.

    double rel_err: the relative error tolerance, 0.0 <= relerr <= 1.0.

  Output:

    struct image *img: the modified image.
*/
{
  int m = img->pam.height;
  int n = img->pam.width;
  int M = img->pam.maxval;
  int zero_count;

  haar_transform_matrix ( img->r, m, n, WT_FWD );
  haar_transform_matrix ( img->g, m, n, WT_FWD );
  haar_transform_matrix ( img->b, m, n, WT_FWD );

  zero_count = prune_ppm ( img->r, img->g, img->b, m, n, rel_err );

  haar_transform_matrix ( img->r, m, n, WT_REV );
  haar_transform_matrix ( img->g, m, n, WT_REV );
  haar_transform_matrix ( img->b, m, n, WT_REV );

  clip_matrix ( img->r, m, n, M );
  clip_matrix ( img->g, m, n, M );
  clip_matrix ( img->b, m, n, M );

  fprintf ( stderr, "\n" );
  fprintf ( stderr, "reduce_ppm_image():\n" );
  fprintf ( stderr, "  zeroed %d of %d wavelet coefficients\n", zero_count, m*n );
  fprintf ( stderr, "  %d remaining.\n", m * n - zero_count );

  return;
}
/******************************************************************************/

void show_usage ( char *progname )

/******************************************************************************/
/*
  Purpose:

    show_usage() prints a reminder of how to use the program.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 April 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    char *progname: the name of the program.
*/
{
  fprintf ( stderr, "Usage: %s relerr infile outfile\n", progname );
  fprintf ( stderr, "0.0 <= relerr <= 1.0\n" );
  fprintf ( stderr, "1. Reads a PGM or PPM image from infile.\n" );
  fprintf ( stderr, "2. Applies the Haar wavelet transform.\n" );
  fprintf ( stderr, "3. Sets small coefficients to zero.\n" );
  fprintf ( stderr, "4. Reconstructs the image and writes to outfile.\n" );

  return;
}

/******************************************************************************/

double truncation_error ( double t, double **a, int m, int n )

/******************************************************************************/
/*
  Purpose:

    truncation_error() reports the truncation error for a given threshold.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    16 September 2023

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

  Output:
*/
{
  double e;
  int i;
  int j;
  double s;

  e = 0.0;
  for ( i = 0; i < m; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      s = fabs ( a[i][j] );
      if ( s < t )
      {
        e = e + s * s;
      }
    }
  }
  return e;
}