// Traitement numérique du signal
// Projet 3
//
// Redimensionnement d'image
//
// Christophe Boyanique
// Emmanuel Pinard
// Avril 1999
//
//


#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>


// Definition de macros:
#define min(a,b)     (a<b)?(a):(b)
#define max(a,b)     (a>b)?(a):(b)


// Definition de codes de retour du programme
#define NO_ERR       0
#define BAD_PARAMS   1
#define NO_FILE      2
#define BAD_FILE     3
#define READ_ERROR   4
#define WRITE_ERROR  5
#define MEM_ERROR    6

// Definition de la taille de la matrice par defaut
// 1: -1 -> 1 (matrice 3x3)
#define DEF_MATRIX   1


// Definition de types d'images
#define P_UNKNOWN    0
#define P1           1
#define P2           2
#define P3           3


// Definition de type de redimensionnement
#define TYPE_RAW       0
#define TYPE_INTER     1
#define TYPE_SAMPLING  2


// Definition d'une structure 'image' contenant les parametres d'une
// image
typedef struct
{
  char  filename[FILENAME_MAX+2];  // Nom du fichier
  int   w;                         // Largeur
  int   h;                         // Hauteur
  int   type;                      // Type PNM
  int   maxcomp;                   // Valeur max de chaque composante
  void  *adr;                      // Adresse du buffer memoire
  long  len;                       // Taille du buffer memoire
} IMGstr;


// Prototypage des fonctions:
void usage(char *name);
int loadImage(IMGstr *img);
int unloadImage(IMGstr *img);
int redim(int type, int matrix, IMGstr *src, IMGstr *dst);
int read_number(FILE *handle);
int write_number(FILE *handle, int n, int ligne);
void clear(void *adr, size_t len);
double sinc(double x);


// Fonction principale du programme
//
int main(int argc, char **argv)
{
  int     w;                    // Nouvelle largeur
  int     h;                    // Nouvelle hauteur
  IMGstr  src;                  // Image source
  IMGstr  dst;                  // Image de destination
  int     err;                  // Code d'erreur
  int     type;                 // Type de redimensionnement
  int     matrix = DEF_MATRIX;  // Taille de la matrice de reechantillonnage

  if (argc!=6)
  {
    usage(argv[0]);
    return BAD_PARAMS;
  }
  // Type de redimensionnement
  if (argv[1][0] == 'r')
    type = TYPE_RAW;
  else if (argv[1][0] == 'b')
    type = TYPE_INTER;
  else if (argv[1][0] == 's')
  {
    type = TYPE_SAMPLING;
    if (atoi(&argv[1][1]) >= 0)
      matrix = atoi(&argv[1][1]);
  }
  else
  {
    usage(argv[0]);
    return BAD_PARAMS;
  }

  // Preparation de la structure et chargement de l'image source
  strncpy(src.filename, argv[4], FILENAME_MAX);
  err = loadImage(&src);
  if (err != NO_ERR)
    return err;

  // Recuperation de la nouvelle taille entree en ligne de commande
  w = atoi(argv[2]);
  h = atoi(argv[3]);
  if ( (w<1) || (h<1) )
  {
    usage(argv[0]);
    return BAD_PARAMS;
  }

  // Preparation de la structure pour l'image destination
  dst.w    = w;
  dst.h    = h;
  dst.type = src.type;
  strncpy(dst.filename, argv[5], FILENAME_MAX);

  // Redimensionnement de l'image
  err = redim(type, matrix, &src, &dst);

  // Liberation memoire de l'image source
  unloadImage(&src);

  return err;
}


// Fonction usage() affichant la syntaxe du programme
//
void usage(char *name)
{
  printf("Usage:\n");
  printf("  %s r width height input.p[bgp]m output.p[bgp]m\n", name);
  printf("  %s b width height input.p[bgp]m output.p[bgp]m\n", name);
  printf("  %s sx width height input.p[bgp]m output.p[bgp]m\n", name);
  printf("\n");
  printf("Raw (r), Bilinear interpolation (b) or Resampling (s) with\n");
  printf("factor x (1+2*x matrix)\n");
}


// Fonction de chargement d'une image
//
int loadImage(IMGstr *img)
{
  FILE  *ha;             // Handle de fichier
  char  buf[82];         // Buffer de lecture
  int   err;             // Retour d'erreur de fonctions
  long  i;               // compteur boucle for
  char  *p;

  // Type d'image: inconnu (pour le moment)
  img->type = P_UNKNOWN;

  // Ouvre le fichier en lecture
  ha = fopen(img->filename, "r");

  // Test d'erreur: fichier inexistant ?
  if (ha == NULL)
    return NO_FILE;

  // Lecture des deux premiers caracteres du fichier
  err = (int)fread((void *)buf, (size_t)1, (size_t)2, ha);
  if (err == 2 && buf[0]=='P')
  {
    switch (buf[1])
    {
      case '1':
        img->type = P1;  // P1: PBM (monochrome)
        break;
      case '2':
        img->type = P2;  // P2: PGM (niveaux de gris)
        break;
      case '3':
        img->type = P3;  // P3: PPM (true color)
        break;
    }
  }
  if (img->type == P_UNKNOWN)
  {
    fclose(ha);
    return BAD_FILE;
  }

  // Lecture de la largeur de l'image
  img->w = read_number(ha);
  if (img->w<1)
  {
    fclose(ha);
    return BAD_FILE;
  }
  // Lecture de la hauteur de l'image
  img->h = read_number(ha);
  if (img->h<1)
  {
    fclose(ha);
    return BAD_FILE;
  }

  img->maxcomp = 1;
  // Si P[23] alors on lit la valeur max d'une composante
  if ( (img->type==P2) || (img->type==P3) )
  {
    img->maxcomp = read_number(ha);
    if (img->maxcomp != 255)
    {
      fclose(ha);
      return BAD_FILE;
    }
  }

  // Determine la taille de la memoire necessaire suivant le
  // type d'image
  switch (img->type)
  {
    case P1:
    case P2:
      img->len = (long)img->w * (long)img->h;
      break;
    case P3:
      img->len = 3L * (long)img->w * (long)img->h;
      break;
  }

  // Alloue la memoire
  img->adr = malloc(img->len);
  if (img->adr == NULL)
  {
    fclose(ha);
    return MEM_ERROR;
  }

  // Charge le fichier
  p = (char *)img->adr;
  for (i=0; i<img->len; i++)
  {
    err = read_number(ha);
    if (err<0)
    {
      free(img->adr);
      fclose(ha);
      return READ_ERROR;
    }
    ((unsigned char *)img->adr)[i] = err;
  }

  fclose(ha);
  return NO_ERR;

}


// Fonction de liberation memoire pour une image
//
int unloadImage(IMGstr *img)
{
  if (img == NULL)
  {
    return -1;
  }
  free(img->adr);
  return 0;
}


// Les paramètres sont:
// int mid:       seuil de binarisation
// char *input:   nom du fichier d'entrée
// char *output:  nom du fichier de sortie
//
int redim(int type, int matrix, IMGstr *src, IMGstr *dst)
{
  FILE    *ha;           // Handle de fichier
  char    buf[80];       // Buffer de lecture
  int     err;           // Retour d'erreur de fonction
  int     value, nvalue; // Valeur d'un pixel
  long    cnt;           // Compteur de pixel
  int     i, j, k;       // Indice de boucle (image source: x, y, composante)
  int     k1, k2;        // Indice de boucle (image source: x, y, composante)
  int     x, y;          // Coordonnees image source
  int     nc;            // Nombre de composantes de l'image
  double  diff[4];       // Distance pour l'interpolation bilineaire
  double  dx[4],dy[4];   // Coordonnees image source (points voisins, reelles)
  int     ix[4],iy[4];   // Coordonnees image source (points voisins, entieres)
  double  xr, yr;        // Coordonnees image source (reelles)
  int     v[4];          // Valeur des 4 points voisins
  double  V;
  double  sin, ssin;
  double  *adr_sinc;     // Adresse de la table de sinc precalculee


  // Ouverture du fichier de sortie
  ha = fopen(dst->filename, "w");
  if (ha == NULL)
  {
    return WRITE_ERROR;
  }

  // Ecriture du header du fichier de sortie
  // Type de fichier + commentaire
  sprintf(buf,"P%1d\n# CREATOR: C.Boyanique/E.Pinard\n", src->type);
  err = (int)fwrite((const void *)buf, (size_t)1, strlen(buf), ha);
  if (err != strlen(buf))
  {
    fclose(ha);
    return WRITE_ERROR;
  }

  if (type == TYPE_RAW)
    sprintf(buf,"# Raw\n");
  else if (type == TYPE_INTER)
    sprintf(buf,"# Bilinear interpolation\n");
  else if (type == TYPE_SAMPLING)
    sprintf(buf,"# Ideal resampling\n");
  err = (int)fwrite((const void *)buf, (size_t)1, strlen(buf), ha);
  if (err != strlen(buf))
  {
    fclose(ha);
    return WRITE_ERROR;
  }

  // Dimension de l'image
  sprintf(buf,"%3d %3d\n", dst->w, dst->h);
  err = (int)fwrite((const void *)buf, (size_t)1, strlen(buf), ha);
  if (err != strlen(buf))
  {
    fclose(ha);
    return WRITE_ERROR;
  }

  if ( (dst->type==P2) || (dst->type==P3) )
  {
    // Profondeur de chaque composante
    sprintf(buf,"255\n");
    err = (int)fwrite((const void *)buf, (size_t)1, strlen(buf), ha);
    if (err != strlen(buf))
    {
      fclose(ha);
      return WRITE_ERROR;
    }
  }

  // Determination du nombre de composantes
  if (dst->type==P3)
    nc = 3;
  else
    nc = 1;

  // Reservation d'un buffer pour la table de precalcul du sinc
  if (type == TYPE_SAMPLING)
  {
    adr_sinc = (double *)malloc( (2*matrix+1) * (2*matrix+1) * sizeof(double) );
    if (adr_sinc == NULL)
    {
      fclose(ha);
      return MEM_ERROR;
    }
    else
    {
      // Calcul de la matrice:
      for (k1=-matrix; k1<=+matrix; k1++)
        for (k2=-matrix; k2<=+matrix; k2++)
          adr_sinc[(matrix+k1) + (matrix+k2)*(2*matrix+1)] = 
            fabs(sinc(M_PI*sqrt( k1*k1 + k2*k2 )));

      // Calcul de la somme de tous les elements de la matrice:
      V=0;
      for (k1=-matrix; k1<=+matrix; k1++)
        for (k2=-matrix; k2<=+matrix; k2++)
          V += adr_sinc[(matrix+k1) + (matrix+k2)*(2*matrix+1)];

      // Normalisation de la matrice:
      for (k1=-matrix; k1<=+matrix; k1++)
        for (k2=-matrix; k2<=+matrix; k2++)
          adr_sinc[(matrix+k1) + (matrix+k2)*(2*matrix+1)] /= V;
    }
  }

  // Mise a zero du compteur de pixels
  cnt = 0;

  // Balayage des lignes
  for (j=0; j<dst->h; j++)
  {

    // Balayage des colonnes
    for (i=0; i<dst->w; i++)
    {

      // Balayage des composantes
      for (k=0; k<nc; k++)
      {

        // Redimensionnement brut
        if (type == TYPE_RAW)
        {
          // Calcul du point le plus proche:
          x = (int)floor( (double)i * (double)src->w / (double)dst->w );
          y = (int)floor( (double)j * (double)src->h / (double)dst->h );

          x=min(x,src->w-1);
          y=min(y,src->h-1);

          nvalue = (int)( ((unsigned char *)src->adr)[nc*(src->w*y+x)+k] );
        }

        // Interpolation bilineaire
        else if (type == TYPE_INTER)
        {
          // Calcul du point le plus proche:
          x = (int)( -0.5 + (double)i * (double)src->w / (double)dst->w );
          y = (int)( -0.5 + (double)j * (double)src->h / (double)dst->h );

          x = min(x, src->w-1);
          y = min(y, src->h-1);
          x = max(x, 0);
          y = max(y, 0);

          // Calcul du point central (coordonnees reelles):
          xr = (double)( -0.5 + (double)(i) * (double)(src->w) / (double)(dst->w) );
          yr = (double)( -0.5 + (double)(j) * (double)(src->h) / (double)(dst->h) );

          xr = min(xr, (double)src->w-1);
          yr = min(yr, (double)src->h-1);
          xr = max(xr, (double)0);
          yr = max(yr, (double)0);

          // On tombe pile sur un point donc pas d'interpolation:
          if (((double)x == xr) && ((double)y == yr))
          {
            nvalue = (int)( ((unsigned char *)src->adr)[nc*(src->w*y+x)+k] );
          }

          // Interpolation:
          //
          // diff[0] -->  X    X <-- diff[1]
          //                °
          // diff[2] -->  X    X <-- diff[3]
          //
          else
          {
            dx[0] = xr - floor(xr);       ix[0] = (int)floor(xr);
            dy[0] = yr - floor(yr);       iy[0] = (int)floor(yr);

            dx[1] = (double)1 - dx[0];    ix[1] = ix[0] + 1;
            dy[1] = dy[0];                iy[1] = iy[0];

            dx[2] = dx[0];                ix[2] = ix[0];
            dy[2] = (double)1 - dy[0];    iy[2] = iy[0] + 1;

            dx[3] = dx[1];                ix[3] = ix[1];
            dy[3] = dy[2];                iy[3] = iy[2];

            if (dx[0] == 0)
            {
              diff[0] = (double)1 / sqrt( dx[0]*dx[0] + dy[0]*dy[0] );
              diff[2] = (double)1 / sqrt( dx[2]*dx[2] + dy[2]*dy[2] );
              V = diff[0]+diff[2];
              v[0] = ((unsigned char *)src->adr)[nc*(src->w*iy[0]+ix[0])+k];
              v[2] = ((unsigned char *)src->adr)[nc*(src->w*iy[2]+ix[2])+k];
              nvalue = (int)(((double)v[0]*diff[0]+(double)v[2]*diff[2])/V); 
            }
            else if (dy[0] == 0)
            {
              diff[0] = (double)1 / sqrt( dx[0]*dx[0] + dy[0]*dy[0] );
              diff[1] = (double)1 / sqrt( dx[1]*dx[1] + dy[1]*dy[1] );
              V = diff[0]+diff[1];
              v[0] = ((unsigned char *)src->adr)[nc*(src->w*iy[0]+ix[0])+k];
              v[1] = ((unsigned char *)src->adr)[nc*(src->w*iy[1]+ix[1])+k];
              nvalue = (int)(((double)v[0]*diff[0]+(double)v[1]*diff[1])/V); 
            }
            else
            {
              diff[0] = (double)1 / sqrt( dx[0]*dx[0] + dy[0]*dy[0] );
              diff[1] = (double)1 / sqrt( dx[1]*dx[1] + dy[1]*dy[1] );
              diff[2] = (double)1 / sqrt( dx[2]*dx[2] + dy[2]*dy[2] );
              diff[3] = (double)1 / sqrt( dx[3]*dx[3] + dy[3]*dy[3] );
              V = diff[0]+diff[1]+diff[2]+diff[3];
              v[0] = ((unsigned char *)src->adr)[nc*(src->w*iy[0]+ix[0])+k];
              v[1] = ((unsigned char *)src->adr)[nc*(src->w*iy[1]+ix[1])+k];
              v[2] = ((unsigned char *)src->adr)[nc*(src->w*iy[2]+ix[2])+k];
              v[3] = ((unsigned char *)src->adr)[nc*(src->w*iy[3]+ix[3])+k];
              nvalue = (int)( ( (double)v[0]*diff[0]
                              + (double)v[1]*diff[1]
                              + (double)v[2]*diff[2]
                              + (double)v[3]*diff[3] ) / V );
            }
          }
        }

        // Reechantillonnage ideal
        else if (type == TYPE_SAMPLING)
        {
          // Calcul du point le plus proche:
          x = (int)floor( (double)i * (double)src->w / (double)dst->w );
          y = (int)floor( (double)j * (double)src->h / (double)dst->h );

          x = min(x, src->w-1);
          y = min(y, src->h-1);
          x = max(x, 0);
          y = max(y, 0);

          xr = (double)( -0.5 + (double)(i) * (double)(src->w) / (double)(dst->w) );
          yr = (double)( -0.5 + (double)(j) * (double)(src->h) / (double)(dst->h) );

          xr = min(xr, (double)src->w-1);
          yr = min(yr, (double)src->h-1);
          xr = max(xr, (double)0);
          yr = max(yr, (double)0);

          V = 0;
          ssin = 0;
          for (k1=x-matrix; k1<=x+matrix; k1++)
            for (k2=y-matrix; k2<=y+matrix; k2++)
            {
              if (k1>=0 && k1<src->w && k2>=0 && k2<src->h)
              {
                value=(int)(((unsigned char *)src->adr)[nc*(src->w*k2+k1)+k] );
                sin = fabs( sinc( M_PI * sqrt(
                         (xr-(double)k1) * (xr-(double)k1)
                       + (yr-(double)k2) * (yr-(double)k2) )));
                ssin += sin;
                if (sin > 0)
                  V += (double)(value) * sin;
              }
            }
          V /= ssin;
          if ( V < 0)
            nvalue = 0;
          else if (V > 255)
            nvalue = 255;
          else
            nvalue = (int)V;

        }

        // Incremente le compteur et passe a la ligne dans le
        // fichier si on arrive en fin de ligne
        cnt++;
        if (cnt%17L == 0)
          err = write_number(ha, nvalue, 1);
        else
          err = write_number(ha, nvalue, 0);
        if (err != nvalue)
        {
          fclose(ha);
          return WRITE_ERROR;
        }
      }
    }
  }

  // Effacement de la matrice des sinus cardinaux:
  if (type == TYPE_SAMPLING)
    free((void *)adr_sinc);

  // Fermeture du fichier
  fclose(ha);

  return NO_ERR;
}



// Fonction de lecture d'un nombre dans un fichier
//
int read_number(FILE *handle)
{
  int n = 0;                      // Nombre à renvoyer
  int ok = 0;                     // Drapeau d'erreur
  unsigned char c;                // Charactere lu dans le fichier

  // Première boucle destinée à purger tous les caractères inutiles
  // (fin de ligne, espace, commentaires...)
  do
  {
    fread((void *)&c, (size_t)1, (size_t)1, handle);
    if ( c == '#' )               // Cas d'un commentaire:
    {                             // il faut de nouveau boucler jusqu'à
      do                          // la fin de la ligne
      {
        fread((void *)&c, (size_t)1, (size_t)1, handle);
      } while ( c!='\n' && !feof(handle) );
    }
  } while ( ( c<'0' || c>'9') && !feof(handle) );

  // Deuxième boucle lisant le nombre une fois que l'on est bien
  // positionné dans le fichier
  while ( ( c>='0' && c<='9') && !feof(handle) )
  {
    ok = 1;
    n = n*10 + (int)c - (int)'0';
    fread((void *)&c, (size_t)1, (size_t)1, handle);
  }

  if (ok == 1)
    return n;

  return -1;
}



// Fonction d'écriture d'un nombre dans un fichier
//
int write_number(FILE *handle, int n, int ligne)
{
  char buf[6];                    // Buffer d'écriture
  int err;                        // Drapeau d'erreur

  if (n>999)                      // Max: nombre à 3 chiffres
    return -1;

  if (ligne != 0)                 // Retour à la ligne
    sprintf(buf, "%3d\n", n);
  else                            // sinon espace pour prochain nombre
    sprintf(buf, "%3d ", n);

  err = fwrite((const void *)buf, (size_t)1, (size_t)strlen(buf), handle);

  if ( err != strlen(buf) );
    return n;

  return -1;
}



// Efface une zone mémoire
//
void clear(void *adr, size_t len)
{
  char *p = (char *)adr;
  size_t i;

  for (i=0; i<len; i++)
    *p++=0;
}



// Fonction sinus cardinal
//
double sinc(double x)
{
  if (x == 0)
    return 1;

  return (sin(x)/x);
}