/*************************************************************************\
|* Povray mesh smoothing utility.                                        *|
|* -----------------------------                                         *|
|*                                                                       *|
|* Should compile with almost any c-compiler.                            *|
|* If you are using a UNIX compiler, you probably have to include the    *|
|* math library with -lm like:                                           *|
|*                                                                       *|
|* gcc smooth.c -o smooth -lm                                            *|
|*                                                                       *|
|* See smooth.txt for more information.                                  *|
|*                                                                       *|
|* Send comments, suggestions and bug reports to warp@iki.fi             *|
\*************************************************************************/

/* Const data */

#define LineMax 1024 /* Max length of an input file line */

/*-----------------------------------------------------------------------*/

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

#if (INT_MAX<65536) /* Use more than just 2 bytes for Word */
typedef long Word;
#define IntArg "%li"
#else
typedef int Word;
#define IntArg "%i"
#endif

/* Point list node */
typedef struct P
{ float x,y,z;        /* Point coords */
  float nx,ny,nz;     /* Normal components of this point */
  Word Triangles;     /* Number of triangles using this point */
  struct T* FirstTri; /* First triangle in triangle list using this point */
  struct P* next;
} Point;

/* Triangle list node */
typedef struct T
{ Point *p1,*p2,*p3;  /* Points of the triangle */
  float nx,ny,nz;     /* Normal components of triangle */
  char Reversed;      /* 1 if reversed by CorrectTriangles() else 0 */
  struct T* next;
} Triangle;

Point* PointList=NULL;
Triangle* TriList=NULL;
Triangle* LastTri=NULL;
Word TriCnt,PntCnt;

/* Calculate normal vector for each point */
void CalculatePointNormals(void)
{ Point* Pnt=PointList;
  Triangle* Tri;
  double x,y,z;
  Word Cnt;

  while(Pnt!=NULL)
  { Tri=Pnt->FirstTri;
    x=y=z=0; /* calculate the average triangle normal components to <x,y,z> */
    for(Cnt=0;Cnt<Pnt->Triangles;Cnt++) /* find triangles using this point */
    { while(Tri->p1 != Pnt && Tri->p2 != Pnt && Tri->p3 != Pnt)
        Tri=Tri->next;
      x+=Tri->nx; y+=Tri->ny; z+=Tri->nz; /* sum normal to the average */
      Tri=Tri->next;
    }
    Pnt->nx=x/Pnt->Triangles; /* divide by # of triangles to get */
    Pnt->ny=y/Pnt->Triangles; /* the average */
    Pnt->nz=z/Pnt->Triangles;
    Pnt=Pnt->next;
  }
}

/* If we can't assume that all the triangles are (c)cw, we have to correct
   them to get the right result */
void CorrectAdjacentTriangles(Triangle* Tri)
{ Triangle* List=TriList;
  Point* tmp;
  int i=0;

  Tri->nx=1; /* if nx>0 then we have corrected the triangle, */
             /* if nx<0 we don't */
  while(i<3 && List!=NULL) /* find at most 3 triangles */
  { if(List->nx<0)
    { /* see if the triangles have 2 points in common */
      if((List->p1==Tri->p1 && List->p2==Tri->p2) || /* check if it's */
         (List->p2==Tri->p1 && List->p3==Tri->p2) || /* an adjacent */
         (List->p3==Tri->p1 && List->p1==Tri->p2) || /* triangle and */
         (List->p1==Tri->p2 && List->p2==Tri->p3) || /* wrong sided */
         (List->p2==Tri->p2 && List->p3==Tri->p3) ||
         (List->p3==Tri->p2 && List->p1==Tri->p3) ||
         (List->p1==Tri->p3 && List->p2==Tri->p1) ||
         (List->p2==Tri->p3 && List->p3==Tri->p1) ||
         (List->p3==Tri->p3 && List->p1==Tri->p1))
      { tmp=List->p1; List->p1=List->p2; List->p2=tmp; /* reverse it */
        List->Reversed=1;
        CorrectAdjacentTriangles(List); /* continue recursively */
        i++;
      }
      else
      { if((List->p1==Tri->p2 && List->p2==Tri->p1) || /* check if it's */
           (List->p2==Tri->p2 && List->p3==Tri->p1) || /* an adjacent */
           (List->p3==Tri->p2 && List->p1==Tri->p1) || /* triangle and */
           (List->p1==Tri->p3 && List->p2==Tri->p2) || /* correct */
           (List->p2==Tri->p3 && List->p3==Tri->p2) ||
           (List->p3==Tri->p3 && List->p1==Tri->p2) ||
           (List->p1==Tri->p1 && List->p2==Tri->p3) ||
           (List->p2==Tri->p1 && List->p3==Tri->p3) ||
           (List->p3==Tri->p1 && List->p1==Tri->p3))
        { CorrectAdjacentTriangles(List); /* continue recursively */
          i++;
        }
      }
    }
    List=List->next;
  }
}

void CorrectTriangles(void)
{ Triangle* Tri=TriList;

  while(Tri!=NULL)
  { if(Tri->nx<0) CorrectAdjacentTriangles(Tri);
    Tri=Tri->next;
  }
}

/* Free point list and triangle list */
void FreeLists(void)
{ Point* tmp1;
  Triangle* tmp2;

  while(PointList!=NULL)
  { tmp1=PointList;
    PointList=PointList->next;
    free(tmp1);
  }
  while(TriList!=NULL)
  { tmp2=TriList;
    TriList=TriList->next;
    free(tmp2);
  }
  LastTri=NULL;
}

/* Add a point to the point list */
Point* AddPoint(Point* PList,float p[3],Triangle* Tri)
{ Point* Pnt=(Point*)malloc(sizeof(Point));
  if(Pnt==NULL) { fprintf(stderr,"Not enough memory.\n"); exit(1); }
  Pnt->next=PList;
  Pnt->x=p[0]; Pnt->y=p[1]; Pnt->z=p[2];
  Pnt->Triangles=1;
  Pnt->FirstTri=Tri;
  PntCnt++;
  return Pnt;
}

/* Calculates the normal vector of a triangle */
void TriangleNormal(Triangle* Tri)
{ double x1,y1,z1,x2,y2,z2,nx,ny,nz,len;
  x1= Tri->p3->x - Tri->p1->x; /* cross-product (p3-p1)x(p2-p1) */
  y1= Tri->p3->y - Tri->p1->y;
  z1= Tri->p3->z - Tri->p1->z;
  x2= Tri->p2->x - Tri->p1->x;
  y2= Tri->p2->y - Tri->p1->y;
  z2= Tri->p2->z - Tri->p1->z;
  nx= y1*z2 - y2*z1;
  ny= x2*z1 - x1*z2;
  nz= x1*y2 - x2*y1;
  len=sqrt(nx*nx+ny*ny+nz*nz); /* normalize (and multiply by 10) */
  if(len>0) len=10/len;
  Tri->nx= nx*len;
  Tri->ny= ny*len;
  Tri->nz= nz*len;
}

/* Calculate triangle normals (if we have used CorrectTriangles()) */
void CalculateTriangleNormals(void)
{ Triangle* Tri=TriList;

  while(Tri!=NULL)
  { TriangleNormal(Tri);
    Tri=Tri->next;
  }
}

/* Add a triangle to the end of triangle list */
void AddTriangle(float Points[3][3],char NoAssume)
{ Point* tmp;
  Point* p[3];
  int i;
  Triangle* Tri=(Triangle*)malloc(sizeof(Triangle));
  if(Tri==NULL) { fprintf(stderr,"Not enough memory.\n"); exit(1); }

  for(i=0;i<3;i++)
  { tmp=PointList;
    while(tmp!=NULL) /* search for points that this triangle is using */
    { if(tmp->x==Points[i][0] &&
         tmp->y==Points[i][1] &&
         tmp->z==Points[i][2])
        break;
      tmp=tmp->next;
    }
    if(tmp!=NULL)  /* if point already exists, use it */
    { p[i]=tmp;
      tmp->Triangles++;
    }
    else           /* else create a new point */
    { PointList=AddPoint(PointList,Points[i],Tri);
      p[i]=PointList;
    }
  }

  Tri->p1=p[0]; Tri->p2=p[1]; Tri->p3=p[2];

  if(NoAssume) /* if we can't assume that triangles are defined cw or ccw */
    Tri->nx=-1; /* then set the value for CorrectTriangles() */
  else
    TriangleNormal(Tri); /* else just calculate the normal */

  if(TriList==NULL)
    TriList=LastTri=Tri;
  else
  { LastTri->next=Tri;
    LastTri=Tri;
  }
  Tri->Reversed=0;
  Tri->next=NULL;
  TriCnt++;
}

/* Compare n characters of two strings */
int StrNCmp(char* s1,char* s2,int n)
{ int i;
  for(i=0;i<n && s1[i] && s2[i];i++)
    if(s1[i]!=s2[i]) return 0;
  if(i<n) return 0;
  return 1;
}

/* Scan 's1' for the first occurrence of the substring 's2' */
int StrStr(char* s1,char* s2)
{ int Ind=0;
  while(s1[Ind])
  { if(StrNCmp(&s1[Ind],s2,strlen(s2))) return Ind;
    Ind++;
  }
  return -1;
}

/* Scan 's' for the first occurence of the given character */
int StrChr(char* s,char c)
{ int Ind=0;
  while(s[Ind])
  { if(s[Ind]==c) return Ind;
    Ind++;
  }
  return -1;
}

char* InputLine;
int LinePtr=-1;
long FilePtr;

/* Read one line from file */
char* GetLine(FILE* InFile)
{ FilePtr=ftell(InFile);
  return fgets(InputLine,LineMax,InFile);
}

/* Read data from infile until 'Str' is found and output everything else */
int FindStr(char* Str,FILE* InFile,FILE* OutFile,char Output)
{ int Ptr;
  char c;

  do
  { if(LinePtr==-1)
    { if(GetLine(InFile)==NULL) return 0;
      LinePtr=0;
    }

    Ptr=StrStr(&InputLine[LinePtr],Str);

    if(Ptr==-1)
    { if(Output) fputs(&InputLine[LinePtr],OutFile);
      LinePtr=-1;
    }
    else if(Ptr>0)
    { if(Output)
      { c=InputLine[LinePtr+Ptr]; InputLine[LinePtr+Ptr]=0;
        fputs(&InputLine[LinePtr],OutFile);
        InputLine[LinePtr+Ptr]=c;
      }
      LinePtr+=Ptr;
    }
  } while(LinePtr==-1);

  return 1;
}

/* Read data from infile until 'Chr' is found and output everything else */
int FindChr(char Chr,FILE* InFile,FILE* OutFile,char Output)
{ int Ptr;
  char c;

  do
  { if(LinePtr==-1)
    { if(GetLine(InFile)==NULL)
        return 0;
      LinePtr=0;
    }

    Ptr=StrChr(&InputLine[LinePtr],Chr);

    if(Ptr==-1)
    { if(Output) fputs(&InputLine[LinePtr],OutFile);
      LinePtr=-1;
    }
    else if(Ptr>0)
    { if(Output)
      { c=InputLine[LinePtr+Ptr]; InputLine[LinePtr+Ptr]=0;
        fputs(&InputLine[LinePtr],OutFile);
        InputLine[LinePtr+Ptr]=c;
      }
      LinePtr+=Ptr;
    }
  } while(LinePtr==-1);

  return 1;
}

/* Read data from infile until a digit is found */
int FindDigit(FILE* InFile)
{ while((InputLine[LinePtr]<'0' || InputLine[LinePtr]>'9') &&
        InputLine[LinePtr]!='-' && InputLine[LinePtr]!='.')
  { if(InputLine[LinePtr]==0)
    { if(GetLine(InFile)==NULL)
        return 0;
      LinePtr=-1;
    }
    LinePtr++;
  }
  return 1;
}

/* Strip useless trailing zeros and decimal point from ascii float */
void StripFloat(char* Float)
{ int Ind=0;
  while(Float[Ind]) Ind++;
  Ind--;
  while(Float[Ind]=='0') { Float[Ind]=0; Ind--; }
  if(Float[Ind]=='.') Float[Ind]=0;
}

/* Print float value */
void PrintFloat(float f,FILE* OutFile)
{ char Float[100];
  sprintf(Float,"%f",f);
  StripFloat(Float);
  fputs(Float,OutFile);
}

/* Read triangle list from file */
int ReadTriangles(FILE* InFile,char NoAssume)
{ int Cnt=0,FoundTriangles=0,ReadingTri=0;
  float p[3][3]={{0}};

  TriCnt=PntCnt=0;
  if(FindChr('{',InFile,NULL,0)==0) return 0;

  while(1)
  { switch(InputLine[LinePtr])
    { case 0  : if(GetLine(InFile)==NULL) return 0;
                LinePtr=-1;
                break;
      case '{': Cnt++; break;
      case '}': Cnt--; if(Cnt==0) return FoundTriangles;
                break;
      case '<': if(ReadingTri)
                { if(FindDigit(InFile)==0) return 0;
                  p[ReadingTri-1][0]=atof(&InputLine[LinePtr]);
                  if(FindChr(',',InFile,NULL,0)==0) return 0;
                  if(FindDigit(InFile)==0) return 0;
                  p[ReadingTri-1][1]=atof(&InputLine[LinePtr]);
                  if(FindChr(',',InFile,NULL,0)==0) return 0;
                  if(FindDigit(InFile)==0) return 0;
                  p[ReadingTri-1][2]=atof(&InputLine[LinePtr]);
                }
                break;
      case '>': if(ReadingTri)
                { ReadingTri++;
                  if(ReadingTri==4)
                  { AddTriangle(p,NoAssume);
                    ReadingTri=0;
                  }
                }
                break;
      default : /* See if it's a triangle but not a smooth_triangle */
                if(StrNCmp(&InputLine[LinePtr],"triangle",8) &&
                   (LinePtr==0 || InputLine[LinePtr-1]!='_'))
                { ReadingTri=1; FoundTriangles=1;
                }
    }
    LinePtr++;
  }
}

/* Output smooth triangle mesh */
void PrintMesh(FILE* InFile,FILE* OutFile)
{ int Cnt=0,i;
  Triangle* Tri=TriList;
  Point* pnt[3];

  FindChr('{',InFile,OutFile,1);

  while(1)
  {
    switch(InputLine[LinePtr])
    { case 0  : if(GetLine(InFile)==NULL) return;
                LinePtr=-1;
                break;
      case '{': Cnt++; break;
      case '}': Cnt--; if(Cnt==0) return;
                break;
      default : /* See if it's a triangle but not a smooth_triangle */
                if(StrNCmp(&InputLine[LinePtr],"triangle",8) &&
                   (LinePtr==0 || InputLine[LinePtr-1]!='_'))
                { fputs("smooth_triangle",OutFile);
                  LinePtr+=8; Cnt++;
                  if(Tri->Reversed) { pnt[0]=Tri->p2; pnt[1]=Tri->p1; }
                  else { pnt[0]=Tri->p1; pnt[1]=Tri->p2; }
                  pnt[2]=Tri->p3;
                  for(i=0;i<3;i++)
                  { FindChr('>',InFile,OutFile,1); LinePtr++;
                    fputs(">,<",OutFile); PrintFloat(pnt[i]->nx,OutFile);
                    fputs(",",OutFile); PrintFloat(pnt[i]->ny,OutFile);
                    fputs(",",OutFile); PrintFloat(pnt[i]->nz,OutFile);
                    fputs(">",OutFile);
                  }
                  Tri=Tri->next;
                }
    }
    if(LinePtr!=-1 && InputLine[LinePtr]) fputc(InputLine[LinePtr],OutFile);
    LinePtr++;
  }
}

/* The main smoothing routine */
void Smooth(FILE* InFile,FILE* OutFile,char Verbose,char NoAssume)
{ long Ptr;
  Word Res;
  Point* Pnt;
  Triangle* Tri;

  if((InputLine=(char*)malloc(LineMax))==NULL)
  { fprintf(stderr,"Not enough memory.\n"); exit(1);
  }

  while(FindStr("mesh",InFile,OutFile,1)) /* find a mesh */
  { Ptr=FilePtr+LinePtr;

    if(Verbose) fprintf(stderr,"\nReading mesh...\n");
    Res=ReadTriangles(InFile,NoAssume);

    fseek(InFile,Ptr,SEEK_SET);
    GetLine(InFile); LinePtr=0;

    if(Res)
    { if(Verbose)
        fprintf(stderr,"Found mesh with "IntArg" triangles and "IntArg
                " (different) points.\n",TriCnt,PntCnt);

      if(NoAssume)
      { if(Verbose) fprintf(stderr,"Correcting triangles...\n");
        CorrectTriangles();
        CalculateTriangleNormals();
      }

      if(Verbose) fprintf(stderr,"Calculating normals...\n");
      CalculatePointNormals();

      if(Verbose==2)
      { Pnt=PointList;
        fprintf(stderr,"Points:\n");
        while(Pnt!=NULL)
        { fprintf(stderr,"<%g,%g,%g> used by "IntArg" triangles, "
                  "normal <%g,%g,%g>\n",Pnt->x,Pnt->y,Pnt->z,Pnt->Triangles,
                  Pnt->nx,Pnt->ny,Pnt->nz);
          Pnt=Pnt->next;
        }
        Tri=TriList; Res=1;
        while(Tri!=NULL)
        { fprintf(stderr,"Triangle %i normal <%g,%g,%g>\n",Res++,
                  Tri->nx,Tri->ny,Tri->nz);
          Tri=Tri->next;
        }
      }

      PrintMesh(InFile,OutFile);
      FreeLists();
    }
    else
    { fputs("m",OutFile); LinePtr++;
    }
  }
}

int main(int argc,char* argv[])
{ FILE *InFile,*OutFile=stdout;
  int i=1;
  char c,v=0,n=0;

  while(1)
  { if(argc>i && argv[i][0]=='-')
    { switch(argv[i][1])
      { case 'v': v=1; i++; break;
        case 'e': if(argv[i][2]=='v') { v=2; i++; } break;
        case 'n': n=1; i++; break;
      }
    }
    else break;
  }

  if(argc<i+1)
  { printf(
      "\nUsage: smooth [<options>] <infile> [<outfile>]\n\n"
      "Options:\n"
      "  -v  : Verbose (output information to stderr)\n"
      "  -ev : Extra verbose\n"
      "  -n  : Don't assume that triangles are defined (counter-)clockwise\n"
      "If <outfile> is not specified, the file is outputted to stdout.\n\n"
      "Send comments, suggestions and bug reports to warp@iki.fi\n");
    return 0;
  }

  if((InFile=fopen(argv[i],"r"))==NULL)
  { fprintf(stderr,"Can't open "); perror(argv[i]);
    return 1;
  }
  if(argc>i+1)
  { if(strcmp(argv[i],argv[i+1])==0)
    { fprintf(stderr,"Source and destination files must be different.\n");
      return 1;
    }
    if((OutFile=fopen(argv[i+1],"r"))!=NULL)
    { fclose(OutFile);
      fprintf(stderr,"File '%s' exists. Overwrite? (y/n) ",argv[i+1]);
      scanf("%c",&c);
      if(c!='y') return 0;
    }
    if((OutFile=fopen(argv[i+1],"w"))==NULL)
    { fprintf(stderr,"Can't create "); perror(argv[i+1]);
      return 1;
    }
  }

  Smooth(InFile,OutFile,v,n);

  fclose(InFile); fclose(OutFile);

  return 0;
}