/* 
 * PvDim.c, Copyright (C) 1999 - 2002 Jean-Paul Davalan.
 *
 * 29900 Concarneau (France) le : Sun Mar 21 09:57:19 GMT+0200 1999
 *
 * $Modified: dim nov 10 22:49:23 CET 2002
 *
 * PvDim is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 *
 * PvDim is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with PvDim; see the file COPYING.  If not, write to
 * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 * ######### Jean-Paul Davalan <jeanpaul.davalan@wanadoo.fr>
 * Jean-Paul Davalan <jpdvl@wanadoo.fr>
 *
 * gcc -O3 -o pvdim pvdim.c -lm
 *
 */

#include <ctype.h>
#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <unistd.h>
#include <getopt.h>

#ifndef PI
#define PI M_PI
#endif

#define VERSION "0.1"
#define MAXRECURS 8		// iterations
#define KMAX      6

#define max(U,V) ((U>V)?U:V)
#define min(U,V) ((U<V)?U:V)

#define signe(x) ((x < 0 )? -1.0 : ((x == 0)? 0.0 : 1.0 ))

char    Version[] = "0.0",	//
        Appname[] = "PvDim", Date[] = "Sun Mar 21 09:57:19 GMT+0200 1999";	//


char   *InName = NULL,		//
       *OutName = NULL;

/* Gestion des durees */
struct tm *tm;
time_t  date_begin, date_now, date_end;
char   *my_date();
long    cpu_begin, cpu_now, cpu_end;
long    cputime();
char   *cpu_past(long d);
/* Options */

int     opt_stat = 0, opt_verbose = 0;

static char *helparr[] =
{
   "Options :",
   "  --help           -h        aide",
   "  --version        -v        version",
   "  --verbose=#      -V #      0 (quiet) .. 2",
   "  --quiet          -q        verbose = 0 (défaut)",
   "  --copyright      -c",
   "  --input=file     -i file   (default: stdin)",
   "  --output=file    -o file   (default: stdout)",
   "  --xu=#   --ux=#   1ere coordonnee du vecteur u",
   "  --yu=#   --uy=#   2eme    -       -    -     u",
   "  --xv=#   --vx=#   1ere coordonnee du vecteur v",
   "  --yv=#   --vy=#   2eme    -       -    -     v",
   "  --x=#           -x #       ux ou xu (uy = yu = 0)",
   "  --y=#           -y #       vy -  yv (vx = xv = 0)",
   "  --angle=#       -a #       recherche de de 0 a  angle",
   "  --initial=#                angle initial (défaut 0)",
   "  --final=#                        final   (défaut 180)",
   "  --pas=#                    pas (défaut 0.001)",
   "  --step=#        -s #       idem : pas (angulaire du parcours)",
   "  --precision=#   -p #       precision (défaut 1)",
   "  --expansion=#   -e #       rapport maxi des aires",
   "  --klimit=#      -k #       nb maxi |k1| ou |k2|, k1*u + k2*v",
   "                             (défaut 6)",
   "  --recursion=#   -r #       profondeur de recherche",
   "                             (défaut 8) (fract. cont.)",
   "  --statistics    -d         statistiques durees (defaut non)",
   0
};

static const struct option long_options[] =
{
   {"help", no_argument, 0, 'h'},
   {"version", no_argument, 0, 'v'},
   {"copyright", no_argument, 0, 'c'},
   {"input", required_argument, 0, 'o'},
   {"output", required_argument, 0, 'i'},
   {"xu", required_argument, 0, 0},	/* 5 */
   {"ux", required_argument, 0, 0},
   {"yu", required_argument, 0, 0},
   {"uy", required_argument, 0, 0},
   {"xv", required_argument, 0, 0},
   {"vx", required_argument, 0, 0},	/* 10 */
   {"yv", required_argument, 0, 0},
   {"vy", required_argument, 0, 0},
   {"pas", required_argument, 0, 0},
   {"step", required_argument, 0, 's'},
   {"initial", required_argument, 0, 0},
   {"final", required_argument, 0, 0},
   {"angle", required_argument, 0, 'a'},
   {"precision", required_argument, 0, 'p'},
   {"expansion", required_argument, 0, 'e'},
   {"verbose", optional_argument, 0, 'V'},
   {"quiet", no_argument, 0, 'q'},
   {"recursion", required_argument, 0, 'r'},
   {"statistics", no_argument, 0, 'd'},
   {"kmax", required_argument, 0, 'k'},
   {0, 0, 0}
};

int     recursion = MAXRECURS, klimit = 5;

double  X0 = 0,			//
        Y0 = 0,			//
        X1 = 0,			//
        Y1 = 0,			//
        angstartdgr = 0,	//
        angenddgr = 180,	//
        angstart,		//
        angend,			//
        angledgr = 0,		//
        angle,			//
        precision = 1,		//
        expansion = 20,		//
        step = 0.001;		//

void    options(),		/* */
        help(),			/* */
        copyright(),		/* */
        version(),		/* */
        init_rand(),		/* */
        begin_application(),	/* */
        init(),			/* */
        main_run(),		/* */
        end_application(),	/* */
        affiche();		/* */

double  aire(),			/* */
        my_rand();		/* */

int     fcont();

#define Mark(X) fprintf(stderr,"%s\n",X)
/*---------------------------------------------------------------*/

int    main(int argc, char *argv[])
{
  options(argc, argv);
  begin_application(argc, argv);
  init();
  main_run(argc, argv);
  end_application(argc, argv);
  return 0;
}
/*---------------------------------------------------------------*/

void    options(int argc, char *argv[])
{
  int     c, longoptindex = 0, errflag = 0;
  if (argc == 1) {
    help(argc, argv);
    exit(2);
  }
  while ((c = getopt_long_only(argc, argv,
	   "hvcqdi:o:s:a:p:e:x:y:r:V::k:",
	   long_options, &longoptindex)) != EOF) {
    switch (c) {
      case 0:
	switch (longoptindex) {
	  case 5:
	  case 6:
	    X0 = atof(optarg);
	    break;
	  case 7:
	  case 8:
	    Y0 = atof(optarg);
	    break;
	  case 9:
	  case 10:
	    X1 = atof(optarg);
	    break;
	  case 11:
	  case 12:
	    Y1 = atof(optarg);
	    break;
	  case 13:
	  case 14:
	    step = atof(optarg);
	    break;
	  case 15:
	    angstartdgr = atof(optarg);
	    break;
	  case 16:
	    angenddgr = atof(optarg);
	    break;
	  default:
	    break;
	}
	break;
      case 'h':
	help(argc, argv);
	exit(1);
	break;
      case 'v':
	version(argc, argv);
	break;
      case 'c':
	copyright(argc, argv);
	exit(1);
	break;
      case 'i':
	if (optarg && (InName = (char *)malloc((strlen(optarg) + 1) * sizeof(char)))) {
	  *InName = '\x0';
	  InName = strcpy(InName, optarg);
	}
	break;
      case 'o':
	if (optarg && (OutName = (char *)malloc((strlen(optarg) + 1) * sizeof(char)))) {
	  *OutName = '\x0';
	  OutName = strcpy(OutName, optarg);
	}
	break;
      case 's':
	step = atof(optarg);
	break;
      case 'a':
	angledgr = atof(optarg);
	break;
      case 'p':
	precision = atof(optarg);
	break;
      case 'e':
	expansion = atof(optarg);
	break;
      case 'x':
	X0 = atof(optarg);
	break;
      case 'y':
	Y1 = atof(optarg);
	break;
      case 'q':
	opt_verbose = 0;
	break;
      case 'V':
	opt_verbose = 1;
	if (optarg) {
	  opt_verbose = atoi(optarg);
	  if (opt_verbose <= 0 || opt_verbose >= 3)
	    opt_verbose = 1;
	}
	break;
      case 'd':
	opt_stat = 1;
	break;
      case 'r':
	if (optarg)
	  recursion = atoi(optarg);
	break;
      case 'k':
	if (optarg)
	  klimit = atoi(optarg);
	break;
      default:
	errflag++;
	break;
    }
  }
  if (errflag) {
    help(argc, argv);
    exit(2);
  }
}
/*---------------------------------------------------------------*/
void    init(void)
{
/* decomposition simplifiable */
  if (angledgr != 0) {
    if (angstartdgr == 0 && angenddgr == 180) {
      angenddgr = angstartdgr + angledgr;
    } else if (angstartdgr != 0 && angenddgr == 180) {
      angenddgr = angstartdgr + angledgr;
    } else if (angstartdgr != 0 && angenddgr != 180) {
      angledgr = angenddgr - angstartdgr;
    } else {
      angledgr = angenddgr - angstartdgr;
    }
  } else {
    angledgr = angenddgr - angstartdgr;
  }
  if (angledgr == 0) {
    angenddgr = angstartdgr + 180;
    angledgr = 180;
  }
  if (expansion < 1) {
    fprintf(stderr, "Erreur, l'expansion doit etre >= 1\n");
    exit(1);
  }
  if (X0 == 0 && Y0 == 0) {
    fprintf(stderr, "Erreur, le 1er vecteur est nul\n");
    exit(1);
  }
  if (X1 == 0 && Y1 == 0) {
    fprintf(stderr, "Erreur, le 2eme vecteur est nul\n");
    exit(1);
  }
  if (X0 == X1 && Y0 == Y1) {
    fprintf(stderr, "Erreur, les deux vecteurs sont egaux\n");
    exit(1);
  }
}
/*---------------------------------------------------------------*/
double  aire(double u0, double u1, double v0, double v1, double ang,
   double *lw, double *lh, double *err)
{
  int     n;
  double  a0, b0, a1, b1, eh, ew;
  double  co, si, t, x0, y0, x1, y1, ax0, ay0, ax1, ay1, prec = precision;
  t = ang * PI / 180;
  co = cos(t);
  si = sin(t);
  x0 = u0 * co - u1 * si;	/* abscisse 1er vecteur après rotation de +angle ° */
  y0 = u0 * si + u1 * co;
  x1 = v0 * co - v1 * si;
  y1 = v0 * si + v1 * co;
  ax0 = fabs(x0);
  ay0 = fabs(y0);
  ax1 = fabs(x1);
  ay1 = fabs(y1);
  n = fcont(&a0, &b0, -y1, y0, prec);	/* a0 fois le 1er + b0 fois le 2ème vecteur horiz. */
  n = fcont(&a1, &b1, -x1, x0, prec);

  *lw = fabs(a0 * x0 + b0 * x1);
  *lh = fabs(a1 * y0 + b1 * y1);

  eh = fabs(a0 * y0 + b0 * y1);
  ew = fabs(a1 * x0 + b1 * x1);
  *err = max(eh, ew);
  if (*lw < 1 || *lh < 1)
    return -1;
  return (*lw) * (*lh);
}
/*---------------------------------------------------------------*/

int     fcont(double *ap, double *bp, double u, double v, double w)
{
  double  a, b, s, r, a0 = 0, b0 = 1, a1 = 1, b1 = 0, asign, bsign;
  int     k = 0;
  float   x;
/* on doit obtenir ap/bp = u/v */
  asign = signe(u);
  bsign = signe(v);
  u = fabs(u);
  v = fabs(v);
  w = fabs(w);
  w = (w == 0) ? 1 : w;
  if (u == 0 || v == 0) {
    *ap = signe(u);
    *bp = signe(v);
    return 0;
  } else {
    x = u / v;
    r = x;
    while (1) {			/* for(;k<recursion;k++){ */
      s = floor(r);
      r -= s;
      r = (r == 0) ? 0 : 1 / r;
      a = s * a1 + a0;
      b = s * b1 + b0;
      a0 = a1, a1 = a, b0 = b1, b1 = b;
      if (fabs(a * v - b * u) < w || r == 0) {
	*ap = asign * a;
	*bp = bsign * b;
	return k;
      }
  /* Inutile de continuer trop longtemps a, b seraient trop grands.  */
  /* Si on arrive à ce test, le résultat sera  probablement            */
  /* inutilisable. (Ce test n'est pas effectué si u et v sont entiers). */
      else if (a > klimit * u || b > klimit * v) {
	*ap = asign * a0;	/* on récupère la valeur précédente */
	*bp = bsign * b0;
	return (k - 1);
      }
      if (k >= recursion)	/* eliminer si for... */
	return k;
      k++;
    }
    return k;
  }
}
/*---------------------------------------------------------------*/

void    affiche(double Xa, double Ya, double Xb, double Yb, double alp,
   double ai, double wm, double hm, double erreur)
{
  double  a, r, coef = expansion;
  a = fabs(X0 * Y1 - X1 * Y0);
  r = wm * hm / a;
  if (r >= 1 && r <= coef)
    printf(
       "a: %.5f, "		// angle
        "p: %.0f, "		// aire (pixels ?)
        "r: %.1f, "		// quotient des aires
        "l: %.2f, "		// largeur finale
        "h: %.2f, "		// hauteur
        "e: %.4f\n"		// erreur (pixels h ou w)
       ,alp, wm * hm, r, wm, hm, erreur);
}
/*---------------------------------------------------------------*/

long    cputime(void)
{
  return (((long)clock() * 1000) / CLOCKS_PER_SEC);
}
/*---------------------------------------------------------------*/

void    help(int argc, char *argv[])
{
  char  **pt;
  fprintf(stderr, "usage: %s options\n", argv[0]);
  for (pt = helparr; *pt; pt++)
    fprintf(stderr, "%s\n", *pt);
  fprintf(stderr,
//     "exemple:     %s -x 115  -y 65 -e 30\n"
     "exemple :    %s -x 132  -y 94 -e 30\n"
     "--------\n"
     "  Le programme sélectionne les rotations les plus favorables, compte tenu\n"
     "  de la précision et surtout de l'expansion indiquées (1 et 20 par défaut).\n"
     " \n"
//     "  La sortie sur stdout est une suite de lignes comme celle-ci:\n"
     "  a: 35.10200, p: 37224, r: 3.0, l: 162.05, h: 229.71, e: 0.9997\n"
//     "             ang 30.126, aire 29120, largeur 129.5, hauteur 224.9\n"
//     "  la largeur et la hauteur utilisables sont environ 130 et 225 pour un angle\n"
//     "  de rotation de l'image de 30.126 degrés.\n"
     "  la largeur et la hauteur utilisables sont environ 162 et 230 pour un angle\n"
     "  de rotation de l'image de 35.102 degrés.\n"
     "  La rotation est à introduire dans le fichier postscript en écrivant:\n"
//     "             30.126 rotate\n"
     "             35.102 rotate\n"
     "  On peut ensuite obtenir une image aux bonnes dimensions en faisant:\n"
     "  gs -q -dSAFER -dNOPAUSE -sDEVICE=ppm -g162x230 -sOutputFile=-\n"
//     "  gs -q -dSAFER -dNOPAUSE -sDEVICE=ppm -g130x225 -sOutputFile=-\n"
     "                                              - -quit < fic.ps > fic.ppm\n\n"
     "  On change de type de fichier en faisant :\n"
     "             ppmtogif fic.ppm > fic.gif ; rm fic.ppm\n\v"
     "  Il est prudent d'effectuer des essais avec des angles \"alpha\" différents\n"
     "  pour choisirl'image qui convient le mieux.\n"
     "  Des valeurs proches de \"alpha\" (à bien moins de 1° près) peuvent conduire à\n"
     "  des tailles très différentes.\n"
     "  On peut aussi obtenir une image au format gif:\n"
     "            pstoimg -gif -out fic.gif fic.ps\n"
     "  (par exemple le pstoimg fourni avec htmltolatex)\n"
     "  Ou en prenant l'image complète dans xv on peut sélectionner une partie \n"
     "  rectangulaire aux dimensions souhaitées (fonction crop), vérifier le bon\n"
     "  raccordement (Root: center tiled) effectuer des transformations\n"
     "  (Emboss ...) (saturation ...)\n"
     "  La qualité du résultat se jugera à la non visibilité des raccordements\n"
     "   d'images.\n"
     " \n"
     "Bugs: \n"
     "-----\n"
     "        Très peu testé, mal corrigé, des modifications sont\n"
     "          commencées mais ne sont pas terminées.\n"
     "        Ce programme est une reprise rapide d'un programme plus\n"
     "          simple. Une refonte totale est envisagée.\n"
     "        klimit entre autres ne semble pas avoir d'effet.\n"
     "        Utilisation des fractions continues, à revoir.\n"
     "        Les valeurs optimales sont elles obtenues ?\n"
     " \n"
     "Ne pas utiliser les résultats de ce programme,\n"
     "---------------------------------------------\n"
     "  (sauf peut-être si cette utilisation ne peut avoir aucune conséquence\n"
     "  fâcheuse, aussi minime soit-elle et uniquement si les résultats sont\n"
     "  vérifiés.)\n"
     ,argv[0]);
  exit(0);
}
/*---------------------------------------------------------------*/

void    copyright(int argc, char *argv[])
{
  fprintf(stdout,
     " Application %s (PvDim), Copyright (C) 1998 Jean-Paul Davalan.\n\n"
     " 29900 Concarneau (France) le : Sun Mar 21 09:57:19 GMT+0200 1999\n"
     " dim nov 10 22:49:23 CET 2002\n\n"
     " %s (PvDim) is free software; you can redistribute it and/or modify\n"
     " it under the terms of the GNU General Public License as published by\n"
     " the Free Software Foundation; either version 2, or (at your option)\n"
     " any later version.\n\n"
     " %s (PvDim) is distributed in the hope that it will be useful,\n"
     " but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
     " MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n"
     " GNU General Public License for more details.\n"
     " You should have received a copy of the GNU General Public License\n"
     " along with %s (PvDim); see the file COPYING.  If not, write to\n"
     " the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.\n\n"
//     " Jean-Paul Davalan <jeanpaul.davalan@wanadoo.fr>\n"
     " Jean-Paul Davalan <jpdvl@wanadoo.fr>\n"
     ,argv[0], argv[0], argv[0], argv[0]);
}
/*---------------------------------------------------------------*/

void    main_run(int argc, char *argv[])
{
  double  _aire, aire0, aire1, aire2, airemin, airemax, anglemin, alpha;
  double  w0, w1, w2, h0, h1, h2, e0, e1, e2, wmin, hmin, err;
  _aire = fabs(X0 * Y1 - X1 * Y0);
  airemax = airemin = _aire * expansion * 10;

  aire0 = aire(X0, Y0, X1, Y1, angstartdgr - 2 * step, &w0, &h0, &e0);
  if (aire0 > 0) {
    airemin = aire0;
    airemax = aire0;
  }
  err = 1000 * precision;
  aire1 = aire(X0, Y0, X1, Y1, angstartdgr - step, &w1, &h1, &e1);
  for (alpha = angstartdgr; alpha <= angenddgr; alpha += step) {
    aire2 = aire(X0, Y0, X1, Y1, alpha, &w2, &h2, &e2);
    if (aire2 > 0 && e2 < precision && aire2 < airemin) {
      anglemin = alpha;
      airemin = aire2;
      wmin = w2;
      hmin = h2;
      err = e2;
    }
    if (e1 < precision
       &&
       aire2 > 0
       && (
	  (
	     (aire1 <= 0 || aire1 < aire0)
	     &&
	     (aire2 <= 0 || aire1 <= aire2)
	  )
	  ||
	  (
	     (aire0 <= 0 || aire1 <= aire0)
	     &&
	     (aire2 <= 0 || aire1 < aire2)
	  )
       )
       ) {
      affiche(X0, Y0, X1, Y1, alpha - step, aire1, w1, h1, e1);
    }
    aire0 = aire1;
    aire1 = aire2;
    w0 = w1;
    h0 = h1;
    w1 = w2;
    h1 = h2;
    e0 = e1;
    e1 = e2;
  }
  fprintf(stdout, "Optimum : \n");
  affiche(X0, Y0, X1, Y1, anglemin, airemin, wmin, hmin, err);
}
/*---------------------------------------------------------------*/

char   *strtolcase(char *s)
{
  char   *p;
  for (p = s; *p != '\x0'; p++)
    if (isalpha(*p))
      *p = tolower(*p);
  return s;
}
/*---------------------------------------------------------------*/

char   *strtoucase(char *s)
{
  char   *p;
  for (p = s; *p != '\x0'; p++)
    if (isalpha(*p))
      *p = toupper(*p);
  return s;
}
/*---------------------------------------------------------------*/

void    version(int argc, char *argv[])
{
  if (strcasecmp(argv[0], Appname) == 0)
    fprintf(stdout,
       "%s version %s of %s\n",
       strtoucase(Appname), Version, Date);
  else
    fprintf(stdout,
       "%s(%s) version %s of %s\n",
       strtoucase(Appname), argv[0], Version, Date);
  exit(0);
}
/*---------------------------------------------------------------*/

void    begin_application(int argc, char *argv[])
{
  char   *s;
  if (opt_stat) {
    date_begin = time(NULL);
    cpu_begin = cputime();
    fprintf(stdout, "%s running\n", strtoucase(Appname));
    s = my_date(date_begin);
    fprintf(stdout, "starting : %s\n", s);
    free(s);
  }
/* initialisation dans init */
}
/*---------------------------------------------------------------*/

void    end_application(int argc, char *argv[])
{
  char   *s;
  long    cpudiff;
  if (opt_stat) {
    cpudiff = (cpu_end = cputime()) - cpu_begin;
    s = cpu_past(cpudiff);
    fprintf(stdout, "cpu time : %s\n", s);
    fprintf(stdout, "%s done\n", strtoucase(Appname));
  }
}
/*---------------------------------------------------------------*/

#define TIMETOSTRING(X,Y) if(d > X) { \
       if(X > 0) { \
         tmp = d / (X);\
         d = d % (X);\
       } else \
         tmp = d;\
       chn[0] = '\x0';\
       sprintf(chn,"%ld %s ",tmp, Y);\
       strcat(str,chn);\
     }

char   *cpu_past(long d)
{
  long    tmp;
  char    chn[100];
  char   *str = (char *)malloc(100 * sizeof(char));
  *str = '\x0';
  if (d < 10) {
    strcat(str, "0 ms");
  } else {
    TIMETOSTRING(86400000, "d");
    TIMETOSTRING(3600000, "h");
    TIMETOSTRING(60000, "min");
    TIMETOSTRING(1000, "s");
    TIMETOSTRING(0, "ms");
  }
  str = (char *)realloc(str, (strlen(str) + 1));
  return (str);
}
/*---------------------------------------------------------------*/

char   *my_date(time_t t)
{
  struct tm *tm;
  char   *str, chn[100];
  chn[0] = '\x0';
  tm = gmtime(&t);
  strftime(chn, 60, "%A %B %d - %H h %M min %S s GMT - %Y", tm);
  str = (char *)malloc((strlen(chn) + 1) * sizeof(char));
  str = strcpy(str, chn);
  return str;
}
/*---------------------------------------------------------------*/

void    init_rand(void)
{
  srand((unsigned int)(time(NULL) * getpid()));
}
/*---------------------------------------------------------------
  retourne une valeur entre 0 et 1
*/

double  my_rand(void)
{
  return ((double)rand()) / (RAND_MAX + 1.0);
}
/*---------------------------------------------------------------*/