/* Last edited on 2003-04-07 02:45:22 by stolfi */

    /* wk->mk = alloc_line_mask(wrad); */
  
double *alloc_line_mask(int wrad)
  { int nw = 2*wrad+1;
    double *mk = (double *)malloc(nw*nw*sizeof(double));
    if (mk == NULL) { pm_error("out of memory for line mask"); }
    return mk;
  }


    double *wt,     /* Weight table. */
    double *mk,     /* Pattern. */
    int wrad        /* Half-width of window, times {NSTEPS}. */

    double *wt,
    int wrad,
    double *prof,
    int prad,
    double *mk

    int hfrac = (hp - col;
    int vfrac = (vp - ;
    int col_min = ceil(hp - wrad/(
    double radius2 = (double)(radius*radius);
            /* Displacements from window center {hp,vp} to pixel center: */
            double fdh = ((double)dh)/FN;
            double fdv = ((double)dv)/FN;
            /* Window weight function: */
            double w = window_weight((fdh*fdh + fdv*fdv)/radius2);
      /* Normalize {prof} to have zero mean: */
    { /* Use weighted mean to account for effect of gaussian window: */
      double sp = 0.0, sw = 0.0;
      for (i = 0; i < nprof; i++) 
        { double t = (double)(i - prad)/(double)prad;
          double w = sqrt(fabs(1.0 - t*t));  /* Extent of circular window. */
          sp += w*prof[i];
          sw += w;
        }
      sp /= sw;
      for (i = 0; i < nprof; i++) { prof[i] -= sp; }
    }
    /* Normalize {prof} to have unit variance: */
    { /* Use weighted mean to simulate circular window: */
      double spp = 0.0, sw = 0.0;
      for (i = 0; i < nprof; i++) 
        { double t = (double)(i - prad)/(double)prad;
          double w = sqrt(fabs(1.0 - t*t));  /* Extent of circular window. */
          spp += w*prof[i]*prof[i];
          sw += w;
        }
      spp = sqrt(spp/sw);
      for (i = 0; i < nprof; i++) { prof[i] /= spp; }
    }
              
  
double *make_ideal_profile(double linewd, double eps, int *pradP);
  /* Creates a profile for a line with nominal width {linewd}. The
    profile is approximately a Gaussian with standard deviation
    {linewd/2}, modified so that it is zero when the Gaussian would be
    less than {eps}. The profile is returned as an array
    {prof[0..2*prad-1]}, sampled every {1/NSTEPS} pixels; where {prad =
    *pradP} is a large enough integer, computed by the procedure.
    The line axis lies between the two elements {prof[0..prad-1]}
    and {prof[0..prad]}. */

/* For the following procedures, the line mask {mk} and the weight
  mask {wk} are assumed to be square of side {2*wrad}. The center of
  either mask (which lies between four elements) is assumed to be be
  the center of the current window. In other words, when the
  processing window is centered at point {(hp/NSTEPS,vp/NSTEPS)}, the
  weight of a pixel {img[v,h]} will be {wt[wrad+dv,wrad+dh]}, and its
  ideal value will be {mk[wrad+dv,wrad+dh]}; where {dh = (h+0.5)*N -
  hp}, {dv = (v+0.5)*N - vp}, and {N = NSTEPS}. */
    
double *make_weight_mask(double radius, double eps, int *wradP);
  /* Returns an array {wt} containing a weight function with
    nominal radius {radius} (in pixels). The half-width {wrad =
    *wradP} is computed by the procedure.
    
    The weight function is approximately a Gaussian bump, with
    standard deviation {radius} when projected along each axis,
    modified so that it is zero when the Gaussian would be less than
    {eps}. */

double *alloc_line_mask(int wrad);
  /* Allocates a square array {mk} with {2*wrad} rows and columns,
    suitable for {compute_line_mask} below.  */

void compute_line_mask
  ( double hn, double vn, 
    double *wt,
    int wrad,
    double prof,
    int prad,
    double *mk
  );
  /* Stores in {*mk} the image of an ideal line with profile {prof},
    swept orthogonally to the vector {(hn,vn)} in the {H,V} coordinate
    system. The center of the {mk} array will lie on the line axis.
    
    The mask is normalized to have zero mean and unit variance
    (weighted by the {wt} function) when it is sub-sampled with stride
    {NSTEPS}, for every possible alignment. */

double hfrac = hp - (double)col, vfrac = vp - (double)row;
double radius2 = (double)(radius*radius);
            
          { /* Find best position of ``horizontal'' line around {hpcur,vpcur}: */
          double best_h_score = -MAXDOUBLE;
          double best_vp = vpcur;
          double best_hp = hpcur;
          double best_ht = htcur;
          for (dt = -maxdt; dt <= +maxdt; dt++) 
            { /* Candidate inclination of ``horizontal'' lines: */
              double ht_c = htcur + tstep * (double)dt;
              /* Unit vector perpendicular to ``horizontal'' lines: */
              double nh_c = sin((ht_c)*PI/180.0);
              double nv_c = cos((ht_c)*PI/180.0);
              for (dv = -maxdp; dv <= +maxdp; dv++)
                { /* Candidate line displacement in pixels: */
                  double disp_c = pstep * (double)dv;
                  /* Candidate center of line segment: */
                  double hp_c = hpcur + disp_c * nh_c;
                  double vp_c = vpcur + disp_c * nv_c;
                  
                  int trash = (debugit = ((dt == 0) & (dv == 0)));
                  double score = line_match(img, hp_c, vp_c, nh_c, nv_c, rad, prof);
                  fprintf(stderr, "  p = (%6.1f %6.1f ) ht = %6.1f  score = %8.3f",
                    hp_c, vp_c, ht_c, score
                  );
                  if (score > best_h_score) 
                    { best_h_score = score; 
                      best_vp = vp_c; 
                      best_hp = hp_c; 
                      best_ht = ht_c;
                      fprintf(stderr, " (improved)");
                    }
                  fprintf(stderr, "\n");
                }
            }
          vpcur = best_vp;
          hpcur = best_hp;
          htcur = best_ht;
        }
        
        { /* Find best position of vertical line: */
          double best_v_score = -MAXDOUBLE;
          double best_hp = hpcur;
          double best_vp = vpcur;
          double best_vt = vtcur;
          fprintf(stderr, "  vert line search...\n");
          for (dt = -maxdt; dt <= +maxdt; dt++)
            { /* Candidate inclination of ``vertical'' lines: */
              double vt_c = vtcur + tstep * (double)dt;
              /* Unit vector perpendicular to ``vertical'' lines: */
              double nh_c = sin((vt_c)*PI/180.0);
              double nv_c = cos((vt_c)*PI/180.0);
              for (dh = -maxdp; dh <= +maxdp; dh++)
                { /* Candidate line displacement in pixels: */
                  double disp_c = pstep * (double)dh;
                  /* Candidate center of line segment: */
                  double hp_c = hpcur + disp_c * nh_c;
                  double vp_c = vpcur + disp_c * nv_c;

                  int trash = (debugit = ((dt == 0) & (dh == 0)));
                  double score = line_match(img, hp_c, vp_c, nh_c, nv_c, rad, prof);
                  fprintf(stderr, "  p = (%6.1f %6.1f ) vt = %6.1f  score = %8.3f",
                    hp_c, vp_c, vt_c, score
                  );
                  if (score > best_v_score) 
                    { best_v_score = score; 
                      best_hp = hp_c; 
                      best_vp = vp_c; 
                      best_vt = vt_c;
                      fprintf(stderr, " (improved)");
                    }
                  fprintf(stderr, "\n");
                }
            }
          hpcur = best_hp;
          vpcur = best_vp;
          vtcur = best_vt;
        }


else if (pm_keymatch(arg, "-tilt", 2))
          { char *rest;
            if (argn+2 >= argc) { pm_usage(usage); }
            o->htilt = strtod(argv[argn+1], &rest); ++argn;
            if (*rest != '\0') { pm_usage(usage); }
            o->vtilt = strtod(argv[argn+1], &rest); ++argn;
            if (*rest != '\0') { pm_usage(usage); }
          }
        

int read_image
  ( int **RR, 
    int **GG, 
    int **BB, 
    char *filename, 
    int *width, 
    int *height, 
    int *maxval
  )
  {
    FILE *fp = fopen(filename, "rb");
    int i=0, c=0;
    int t_width, t_height, t_maxval;
    int n, kind;
    char p;
    int *R=*RR, *G=*GG, *B=*BB;

    if (fp == NULL)
      { printf("%s: cannot open.\n", filename);
        return -1;
      }

    fscanf(fp, "%c", &p);
    if (p != 'P')
      { printf("%s: not a PBM/PGM/PPM file.\n", filename);
        return -2;
      }
    fscanf(fp, "%d", &kind);

    switch (kind)
      {
        case 6:
          skip_lines(fp, 1);
          fscanf(fp, "%d %d %d", &t_width, &t_height, &t_maxval);
          if (width != NULL) *width = t_width;
          if (height != NULL) *height = t_height;
          if (maxval != NULL) *maxval = t_maxval;

          n=t_width*t_height;/*numero de pixels na imagem*/

          if (!(R = (int *)calloc(n, sizeof(int)))) printf("no mem for R.\n");
          if (!(G = (int *)calloc(n, sizeof(int)))) printf("no mem for G.\n");
          if (!(B = (int *)calloc(n, sizeof(int)))) printf("no mem for B.\n");

          *RR=R; *GG=G; *BB=B;

          break;

        default:
          printf("Unsupported file type.\n");
          return -3;
          break;
      }

    i=0;  
    if (R && G && B)
      while(c != EOF)
        {
          if ((c=getc(fp)) != EOF) R[i]= c; else break;
          if ((c=getc(fp)) != EOF) G[i]= c; else break;
          if ((c=getc(fp)) != EOF) B[i]= c; else break;
          i++;
        }
    else
      if (R)
        while(c != EOF)
          {
            if ((c=getc(fp)) != EOF) R[i]= c; else break;
            i++;
          }
      else
        {
          printf("O primeiro vetor deve ser não nulo em read_image.\n");
          return -4;
        }

    fclose(fp);

    if (i  !=  n)
      {
        printf("Talvez o arquivo tenha acabado antes do término da leitura.\n");
        printf("Isso significa que os valores de t_height e t_width estão inconsitentes, ou que o arquivo está truncado.\n");
      }

    return i;
  }

/*SALVA OS VETORES RGB NO ARQUIVO FILENAME(CRIADO AQUI)*/
int write_image(int kind, int *R, int *G, int *B, char *filename, int t_width, int t_height, int t_maxval)
{
  int i;
  int tam = t_width*t_height;
  FILE *fp;
  
  switch (kind)
    {
    case 5:/*PGM*/
      if ((!filename) || (!R))
	{printf("Inconsistência de dados em write_image.\n"); return -2;}
      if (!(fp=fopen(filename, "wb")))
	{printf("O arquivo %s não pode ser criado.\n", filename); return -1;}
      fprintf(fp, "P5\n");
      fprintf(fp, "# Created by Amandio Sena\n");
      fprintf(fp, "%d %d\n%d\n", t_width, t_height, t_maxval);
      for (i=0; i<tam; i++)
	{
	  fputc(R[i], fp);
	}
      fclose(fp);
      
      break;      
      
    case 6:/*PPM*/
      if ((!filename) || (!R) || (!G) || (!B))
	{printf("Inconsistência de dados em write_image.\n"); return -2;}
      if (!(fp=fopen(filename, "wb")))
	{printf("O arquivo %s não pode ser criado.\n", filename); return -1;}
      fprintf(fp, "P6\n");
      fprintf(fp, "# Created by Amandio Sena\n");
      fprintf(fp, "%d %d\n%d\n", t_width, t_height, t_maxval);
      for (i=0; i<tam; i++)
	{
	  fputc(R[i], fp);
	  fputc(G[i], fp);
	  fputc(B[i], fp);
	}
      fclose(fp);
      
      break;
      
    }

  return i;
}


double *detectalinhas2(int nlinhas, int *R, int width, int height)
{
  int i=0, j=0, k=0;/*contadores*/
  double ang=0;
  double *prof, *m, *mm;
  int n=640*180;/*distancia L máxima:640; ang de -90 a 89. */

  prof = curvaideal(0, fator, R, width);
  
  if ((m=(double *)calloc(n,sizeof(double))) == NULL)
    printf("Defeito ao alocar matriz em detectalinhas2.\n");
  for (i=0; i<n; i++)
    { m[i]=0; }

  /*Para cada pixel da imagem grande*/
  for (i=radius; i<(height-radius); i++)
    { fprintf(stderr, ".");
      for (j=radius; j<(width-radius); j++)
        { 
	  double cmax = -1e30, Gv;
	  double mang, mcang, msang;
	  
          for (ang=-90; ang<90; ang+=5)
            { double cang = cos(ang/180.0*pi);
              double sang = sin(ang/180.0*pi);
              double c;
              c=compute_window_profile(R,width,i,j,cang,sang,prof);
              if (c > cmax)
		{
		  cmax = c;
		  mang=ang;
		  mcang=cang;
		  msang=sang;
		}
            }
          /* cmax é a melhor nota para todos os angulos ang neste pixel */
          /* cmax = 0 se o casamento foi perfeito. */
          /* cmax << 0 se o casamento foi ruim. */
          /* Gv = 1.0 se casamento perfeito, 0.0 se péssimo. */
          Gv = (15 - fabs(cmax))/15.0;
          if (Gv < 0.0) { Gv = 0.0; }
          if (Gv > 1.0) { Gv = 1.0; }
	  
	  /*deve pegar uma matriz (o bucket) e somar as notas de cada pixel*/
	  /*esta matriz deve ter o mesmo tamanho da figura*/
	  /*a figura tem 512x384 pixels */
	  /* ang varia de -90 a +89 */
	  /* a distancia maxima e de 640 pixels */
	  m[indice(mang,distancia(i,j,mcang,msang),640)] += Gv;
        }
    }

  ang=0;/*soh aproveitando a variavel, Nao tem nada a ver com ang. Guarda o pior valor encontrado */
  j=k=0;/*contadores auxiliares. */
  for (i=0; i<n; i++)
    {
      if (m[i]>ang)
	{
	  ang=m[i];
	  j++;
	  k=i;/*k guarda o indice do ultimo elemento aproveitado*/
	}
      else
	{
	  if (j<nlinhas)
	    {
	      ang=m[i];
	      j++;
	      k=i;/*k guarda o indice do ultimo elemento aproveitado*/
	    }
	}
    }
  if ((mm=(double *)calloc(2*nlinhas +1, sizeof(double))) == NULL)
    printf("Falha ao alocar mm em detecta2.\n");
  mm[0]=2*nlinhas;
  j=1;
  for (i=0; i<=k; i++)
    {
      if (i==k)
	{
	  mm[j]=linha(i,640);/*ang*/
	  mm[10+j]=coluna(i,640);/*L*/
	  j++;
	}
      else
	if (m[i]>ang)
	  {
	    mm[j]=linha(i,640);
	    mm[10+j]=coluna(i,640);
	    j++;
	  }
    }
  return mm;
}

void detectalinhas(int *R, int width, int height, int *G)
{
  int i=0, j=0;/*contadores*/
  double ang=0;
  double *prof;

  /* Na figura grande(f0a.ppm) o ponto (306,173) é uma linha*/
  /* (129,144) é reta na vertical(angulo 87 graus) na figura f0.ppm */
  /* chutando fator o */
  prof = curvaideal(0, fator, R, width);

  /*Para cada pixel da imagem grande*/
  for (i=radius; i<(height-radius); i++)
    { fprintf(stderr, ".");
      for (j=radius; j<(width-radius); j++)
        { double cmax = -1e30, Gv;
          for (ang=-90; ang<90; ang+=5)
            { double cang = cos(ang/180.0*pi);
              double sang = sin(ang/180.0*pi);
              double c;
              c=compute_window_profile(R,width,i,j,cang,sang,prof);
              if (c > cmax) { cmax = c; }
            }
          if ((i > 100) && (i < 150) && (j > 100) && (j < 150))
            { fprintf(stderr, "(%f)", cmax); }
          /* cmax é a melhor nota para todos os angulos ang neste pixel */
          /* cmax = 0 se o casamento foi perfeito. */
          /* cmax << 0 se o casamento foi ruim. */
          /* Gv = 1.0 se casamento perfeito, 0.0 se péssimo. */
          Gv = (15 - fabs(cmax))/15.0;
          if (Gv < 0.0) { Gv = 0.0; }
          if (Gv > 1.0) { Gv = 1.0; }
          G[indice(i,j,width)] = (int)(255 * Gv + 0.5);
        }
    }
}


int *extrailinhas(int n, int *R, int width, int height)
{
  int i=0, j=0;/*contadores*/
  double ang=0;
  no_t *no=NULL;
  lista_t *lista=NULL;
  int *m;
  double *prof;

  /* Na figura grande(f0a.ppm) o ponto (306,173) é uma linha*/
  /* (129,144) é reta na vertical(angulo 87 graus) na figura f0.ppm */
  /* chutando o fator o */
  prof = curvaideal(0, fator, R, width);

  lista=novalista(TAMMAX);

  /*Para cada pixel da imagem grande*/
  for (i=radius; i<(height-radius); i++)
    {
      for (j=radius; j<(width-radius); j++)
        {
          for (ang=-90; ang<90; ang+=5)
            { double cang = cos(ang/180.0*pi);
              double sang = sin(ang/180.0*pi);
              no=novono();
              no->ang=ang;
              no->v=i;
              no->h=j;
              no->c=compute_window_profile(R,width,i,j,cang,sang,prof);
              insereno(no, lista);
            }
        }
    }
  printf("Inicia o Bucketing.\n");
  m=bucketing(lista, width, height);
  printf("Termina o Bucketing.\n");
  liberalista(lista);
  printf("Antes de extrair picos.\n");
  lista = extraipicos(n,m,width,height);
  printf("Termina extração de picos.\n");
  free(m);
  if ((m=(int *)calloc(2*n, sizeof(int))) == NULL)
    printf("Problemas ao alocar m.\n");
  no=lista->primno;
  for(i=0; i<2*n; i++)
    {
      m[i]=no->v;
      m[2*i]=no->ang;
      no=no->prox;
    }
  printf("Terminou(quase).\n");
  liberalista(lista);
  
  return m;
}
#ifndef linhas_H
#define linhas_H

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

#ifndef PI
#define PI 3.1415926
#endif

/* Raio da janela projetada: */
#define CRAIO 6
#define CRAIO 6

#define TAMMAX 500

#define linha(i,width) (i)/width
#define coluna(i,width) (i)%width
#define indice(l,c,width) ((int)((l)*width+(c)))

/*calcula a distância do ponto (x,y) a uma reta que passa por (0,0) de inclinação ang)*/
  /*eu acho que é assim: */
  /*d = x*sin(ang)-y*cos(ang) */
  /* o stolfi disse que era -x.sen + y.cos */
  /* ele tbm disse que tem que dividir d por um fator o ... */
  /*tô fazendo do jeito do stolfi*/
#define distancia(y,x,cang,sang) ((x)*(cang)-(y)*(sang))


/*um nó representa um pixel(i,j) e a profile de intensidade a partir de um ângulo ang*/
typedef struct strucno_t
{
  double c;
  double ang;
  int v;
  int h;
  /* int width; */
  struct strucno_t *prox;
  struct strucno_t *ant;
} no_t;

/*uma lista contém um ponteiro para o primeiro nó, o tamanho atual da lista e seu tamanho máximo possível(definido arbitrariamente).*/
/*cmin é o pior valor de c(semelhança da profile representada com a profile ideal) encontrado, ou seja, o valor que estáh no último nó.*/
typedef struct struclista_t
{
  int tam;
  int tammax;
  double cmin;
  no_t *primno;
} lista_t;



lista_t *novalista(int tammax);

no_t *novono();

int liberano(no_t *no);

int libno(no_t *no);

int liberalista(lista_t *lista);

int arredonda(double d);

double compute_window_profile(int *R, int width, int v, int h, double cang, double sang, double *prof);

int insereno(no_t *no, lista_t *lista);

int *bucketing(lista_t *lista, int H, int V);

lista_t *extraipicos(int n, int *m, int width, int height);

double *curvaideal(double d, double o, int *R, int width);

void find_crossing(int *R, int width, int height, int ld, int cd, double *pmi, double *pmj);

double *detectalinhas2(int nlinhas, int *R, int width, int height);

void detectalinhas(int *R, int width, int height, int *G);

void normaliza(double *profile);

int *extrailinhas(int n, int *R, int width, int height);


#endif


int skip_lines(FILE *fp, int l)
{
  char c;
  int i=0;

  fscanf(fp, "%c", &c);
  for (i=0; i<l; i++)
    {
      while (c != '\n') fscanf(fp, "%c", &c);
    }

  return l;
}