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

#define val(i,j) ( img->data[((j)*img->sx+(i))] )

unsigned char inxt[5] = {1,2,0,1,2};


/* definition of data structures */
typedef struct {
  int            sx,sy;
  unsigned char *data;
  char          *name;
} PPMimage;
typedef PPMimage * pPPMimage;

typedef struct {
  int            x,y,val;
} Point;
typedef Point * pPoint;
  
typedef struct {
  int    v[3];
} Tria;
typedef Tria * pTria;

typedef struct {
  int       np,nt,tol;
  char     *name;
  Point    *point;
  Tria     *tria;
} Mesh;
typedef Mesh * pMesh;


/* load PPM image */
int loadPPM(PPMimage *img) {
  FILE      *fp;
  int        i,k,ret,r,g,b,maxval,bitsize;
  char      *ptr,c,typimg,buff[1024],data[256];

  /* search for image */
  ptr = strstr(img->name,".ppm");
  strcpy(data,img->name);
  if ( !ptr ) {
    ptr = strstr(img->name,".pgm");
    if ( !ptr )  strcat(data,".ppm");
  }
  fp = fopen(data,"rb");
  if ( !fp ) {
    fprintf(stderr,"  ## UNABLE TO OPEN FILE %s.\n",data);
    return(0);
  }
  fprintf(stdout,"     LOADING  %s\n",data);

  if ( !fgets(buff,sizeof(buff),fp) ) {
    fprintf(stderr,"  ## INVALID HEADER.\n");
    return(0);
  }

  /* check header file */
  if ( buff[0] != 'P' ) {
    fprintf(stderr,"  ## INVALID IMAGE FORMAT (MUST BE 'PX').\n");
    return(0);
  }
  typimg = buff[1];

  do {
    ret = fscanf(fp,"%s",buff);
    if ( ret == EOF ) break;
    /* check and strip comments */
    if ( buff[0] == '#' )
      do
        c = getc(fp);
      while ( c != '\n' );
    else break;
  }
  while (1);

  /* read columns + lines */
  ret  = sscanf(buff,"%d",&img->sx);
  ret += fscanf(fp,"%d",&img->sy);
  if ( ret != 2 ) {
    fprintf(stderr,"  ## ERROR LOADING IMAGE.\n");
    return(0);
  }
  if ( fscanf(fp,"%d",&maxval) != 1 ) {
    fprintf(stderr,"  ## INVALID IMAGE SIZE.\n");
    return(0);
  }

  /* strip line */
  while ( fgetc(fp) != '\n' );

  /* size based on type */
  bitsize = img->sx * img->sy;
  if ( typimg == '3' || typimg == '6' )
    bitsize *= 3;
  fprintf(stdout,"     SIZE     %d x %d = %d pixels\n",img->sx,img->sy,bitsize);

  /* memory allocation */
  img->data = (unsigned char*)malloc(bitsize*sizeof(unsigned char));
  assert(img->data);

  /* read data file */
  if ( typimg == '2' || typimg == '3' ) {
    for (i=0; i<bitsize; i++) {
      fscanf(fp,"%d",&r);
      img->data[i] = r;
    }
  }
  else {  
    ret = fread(img->data,sizeof(unsigned char),bitsize,fp);
    if ( ret != bitsize ) {
      fprintf(stderr,"  ## ERROR LOADING IMAGE.\n");
      return(0);
    }
  }
  fclose(fp);

  /* convert to grey levels */
  if ( typimg == '3' || typimg == '6' ) {
    fprintf(stdout,"  Converting to grey levels (%d)\n",bitsize);
    for (i=0,k=0; i<bitsize; i+=3,k++) {
      r = (int)img->data[i];
      g = (int)img->data[i+1];
      b = (int)img->data[i+2];
      img->data[k] = (unsigned char)(0.3*r+0.59*g+0.11*b);
    }
    img->data = (unsigned char*)realloc(img->data,sizeof(unsigned char)*bitsize/3+1);
  }

  return(1);
}


int saveMesh(pMesh mesh) {
  FILE     *out;
  int       k,nt;
  char     *ptr,data[512];

  strcpy(data,mesh->name);
  ptr = strstr(data,".ppm");
  if ( ptr )  *ptr = '\0';
  strcat(data,".mesh");
  out = fopen(data,"w");
  if ( !out ) {
    fprintf(stderr,"  ## UNABLE TO SAVE MESH FILE.\n");
    exit(2);
  }
  fprintf(stdout,"  %%%% %s OPENED\n",data);

  fprintf(out,"MeshVersionFormatted 1\n");
  fprintf(out,"\nDimension\n 2\n");

  /* save vertices */
  fprintf(out,"\nVertices \n%d\n",mesh->np);
  for (k=1; k<=mesh->np; k++)
    fprintf(out,"%d %d 0\n",mesh->point[k].x,mesh->point[k].y);

  /* save triangles */
  nt = 0;
  for (k=1; k<=mesh->nt; k++)
    if ( mesh->tria[k].v[0] > 0 )  nt++;
  fprintf(out,"\nTriangles\n %d\n",nt);
  for (k=1; k<=mesh->nt; k++) {
    if ( mesh->tria[k].v[0] > 0 )
      fprintf(out,"%d %d %d 0\n",mesh->tria[k].v[0],mesh->tria[k].v[1],mesh->tria[k].v[2]);
  }
  fprintf(out,"\nEnd\n");
  fclose(out);

  /* save solution */
  strcpy(data,mesh->name);
  ptr = strstr(data,".mesh");
  if ( ptr )  *ptr = '\0';
  strcat(data,".sol");
  out = fopen(data,"w");
  if ( !out ) {
    fprintf(stderr,"  ## UNABLE TO SAVE MESH FILE.\n");
    exit(2);
  }
  fprintf(stdout,"  %%%% %s  OPENED\n",data);

  fprintf(out,"MeshVersionFormatted 1\n");
  fprintf(out,"\nDimension\n 2\n");

  /* save vertices */
  fprintf(out,"\nSolAtVertices \n%d\n",mesh->np);
  fprintf(out,"1 1\n");
  for (k=1; k<=mesh->np; k++)
    fprintf(out,"%d\n",mesh->point[k].val);
  fprintf(out,"\nEnd\n");
  fclose(out);
  
  fprintf(stdout,"     %d VERTICES  %d TRIANGLES\n",mesh->np,nt);
  return(1);
}




int main(int argc,char *argv[]) {   
  PPMimage  img;
  Mesh      mesh;

  if ( argc < 2 ) {
    img.name = (char *)calloc(128,sizeof(char));
    assert(img.name);
    fprintf(stdout,"  -- Image name ?\n");
    fflush(stdin); 
    fscanf(stdin,"%s",img.name);
    fprintf(stdout,"  -- Tolerance ?\n");
    fflush(stdin);
    fscanf(stdin,"%d",&mesh.tol);
  }
  else {
    img.name = argv[1];
  }
  mesh.name = (char*)calloc(512,sizeof(char));
  strcpy(mesh.name,img.name);
  
  /* read image */
  fprintf(stdout,"\n  -- READING FILE\n");
  clock();
  if ( !loadPPM(&img) ) {
    fprintf(stdout,"  ## Exit.\n");
    exit(1);
  }
  fprintf(stdout,"  -- DATA READING COMPLETED. %6.2f\n",clock()/(double)CLOCKS_PER_SEC);

  /* memory allocation */
	mesh.np = img.sx * img.sy;
  mesh.point = (Point *)calloc(mesh.np+2,sizeof(Point));
  assert(mesh.point);
  mesh.nt = 2 * mesh.np;
  mesh.tria = (Tria *)calloc(mesh.nt+2,sizeof(Tria));
  assert(mesh.tria);

  /* Delaunization */
  fprintf(stdout,"\n  -- STEP 2: TRIANGULATING\n");
  clock();

  fprintf(stdout,"  -- STEP 2 COMPLETED.       %6.2f sec.\n",clock()/(double)CLOCKS_PER_SEC);

  /* save mesh result */
  fprintf(stdout,"\n  -- SAVING MESH\n");
  clock();
  if ( !saveMesh(&mesh) ) {
    fprintf(stdout,"  ## Unable to save mesh.\n");
    exit(2);
  }
  fprintf(stdout,"  -- MESH WRITING COMPLETED. %6.2f sec.\n",clock()/(double)CLOCKS_PER_SEC);
  fprintf(stdout,"\n   Normal completion of program.\n");

  free(img.data);
  return(0);
}