#include "exttrimesh.h"
#include <string.h>
#include <stdlib.h>

const char *input_filename;

double closestPair(List *bl1, List *bl2, Vertex **closest_on_bl1, Vertex **closest_on_bl2)
{
 Node *n, *m;
 Vertex *v,*w;
 double adist, mindist = DBL_MAX;

 FOREACHVVVERTEX(bl1, v, n)
  FOREACHVVVERTEX(bl2, w, m)
   if ((adist = w->squaredDistance(v))<mindist)
   {
    mindist=adist;
    *closest_on_bl1 = v;
    *closest_on_bl2 = w;
   }

 return mindist;
}

/**
 * Joins the closest components, that have boundaries (holes).
 */
bool joinClosestComponents(ExtTriMesh *tin, bool justconnect = false, bool refine = true, bool fair = true) {
    Vertex *v,*w, *gv, *gw;
    Triangle *t, *s;
    Node *n;
    List triList, boundary_loops, *one_loop;
    List **bloops_array;
    int i, j, numloops;

    i=0;
    // delete info of all triangles
    FOREACHVTTRIANGLE((&(tin->T)), t, n) t->info = NULL;
    // initialize info of all triangles with their component number starting by 1.
    FOREACHVTTRIANGLE((&(tin->T)), t, n) {
        if (t->info == NULL) {
            i++;
            triList.appendHead(t);
            t->info = (void *)i;
            while(triList.numels()) {
                t = (Triangle *)triList.popHead();
                if ((s = t->t1()) != NULL && s->info == NULL) {triList.appendHead(s); s->info = (void *)i;}
                if ((s = t->t2()) != NULL && s->info == NULL) {triList.appendHead(s); s->info = (void *)i;}
                if ((s = t->t3()) != NULL && s->info == NULL) {triList.appendHead(s); s->info = (void *)i;}
            }
        }
    }
    // if less then 2 components
    if (i<2) {
        // unset info again
        FOREACHVTTRIANGLE((&(tin->T)), t, n) t->info = NULL;
        JMesh::info("Mesh is a single component. Nothing done.");
        return false;
    }
    // copy triangle component number to the vertices
    FOREACHVTTRIANGLE((&(tin->T)), t, n) {
        t->v1()->info = t->v2()->info = t->v3()->info = t->info;
    }
    // create list boundary loop lists (= lists of connected vertices on a boundary)
    FOREACHVVVERTEX((&(tin->V)), v, n) {
        // find next vertex of an unmarked boundary
        if (!IS_VISITED2(v) && v->isOnBoundary()) {
            w = v;
            one_loop = new List;
            // mark all vertices at this boundary
            do {
                one_loop->appendHead(w);
                MARK_VISIT2(w);
                w = w->nextOnBoundary();
            } while (w != v);
            boundary_loops.appendHead(one_loop);
        }
    }
    FOREACHVVVERTEX((&(tin->V)), v, n) UNMARK_VISIT2(v);

    bloops_array = (List **)boundary_loops.toArray();
    numloops = boundary_loops.numels();

    int numtris = tin->T.numels();
    double adist, mindist=DBL_MAX;

    gv=NULL;
    for (i=0; i<numloops; i++) {
        for (j=0; j<numloops; j++) {
            // if i,j are indices of vertices of different boundary loops, search for the closes pair of vertices and update mindist
            if (((Vertex *)bloops_array[i]->head()->data)->info != ((Vertex *)bloops_array[j]->head()->data)->info) {
                adist = closestPair(bloops_array[i], bloops_array[j], &v, &w);
                if (adist<mindist) {mindist=adist; gv=v; gw=w;}
            }
        }
    }
    if (gv!=NULL) tin->joinBoundaryLoops(gv, gw, justconnect, refine, fair);

    FOREACHVTTRIANGLE((&(tin->T)), t, n) t->info = NULL;
    FOREACHVVVERTEX((&(tin->V)), v, n) v->info = NULL;

    free(bloops_array);
    while ((one_loop=(List *)boundary_loops.popHead())!=NULL) delete one_loop;

    return (gv!=NULL);
}

//#define DISCLAIMER

void usage()
{
 printf("%s v%s - by %s.\n================================================================================\n", JMesh::app_name, JMesh::app_version, JMesh::app_authors);
 printf("USAGE: meshfix <file1> [<file2>] [OPTIONS]\n");
 printf("  Processes file1 and saves the result to <file1>_fixed.off.\n");
 printf("  An optionally passed file2 is merged with the first one.\n");
 printf("OPTIONS:\n");
 printf(" -a <epsilon_angle>  Allowed range: 0 < epsilon_angle < 2, default: 0 (degrees).\n");
 printf(" -j                  Join 2 biggest components if they overlap, remove the rest.\n");
 printf(" -jc                 Join the closest pair of components.\n");
 printf(" -h, --help          Print this help and exit.\n");
 printf(" --shells <n>        Only the n biggest shells are kept.\n");
 printf(" -o <output>         Set the output filename (without extension).\n");
 printf(" -q                  Quiet mode, don't write much to stdout.\n");
 printf(" --remove-handles    Remove all handles of the mesh.\n");
 printf(" -u <steps>          Uniform remeshing of the whole mesh, steps > 0\n");
 printf("   --vertices <n>    Constrain number of vertices to n (only with -u)\n");
 printf(" --no-clean          Don't clean.\n");
 printf(" --smooth <n>        Apply n laplacian smoothing steps.\n");
 printf(" -s, --stl           Result is saved in STL     format instead of OFF.\n");
 printf(" -w, --wrl           Result is saved in VRML1.0 format instead of OFF.\n");
 printf(" --fsmesh            Result is saved in FreeSurfer format instead of OFF.\n");
 printf(" --xshift <d>        Shift x-coordinates of vertices by d when saving output.\n");
 printf("                     Only works with --fsmesh; used to deal with small FreeSurfer glitch\n");
 printf(" --msh               Result is saved in gmsh format for debugging (including vertex and triangle masks)\n");
 printf(" == Cutting, decoupling, dilation ==\n");
 printf(" --cut-outer <d>     Remove triangles of 1st that are outside of the 2nd shell.\n");
 printf(" --cut-inner <d>     Remove triangles of 1st that are inside  of the 2nd shell.\n");
 printf("                     Dilate 2nd by d; Fill holes and keep only 1st afterwards.\n");
 printf(" --decouple-inin <d> Treat 1st file as inner, 2nd file as outer component.\n");
 printf("                     Resolve overlaps by moving inners triangles inwards.\n");
 printf(" --decouple-outin <d> Treat 1st file as outer, 2nd file as inner component.\n");
 printf("                     Resolve overlaps by moving outers triangles inwards.\n");
 printf(" --decouple-outout <d> Treat 1st file as outer, 2nd file as inner component.\n");
 printf("                     Resolve overlaps by moving outers triangles outwards.\n");
 printf("                     Constrain the min distance between the components > d.\n");
 printf(" --fineTuneIn <d> <n> Used to fine-tune the minimal distance between surfaces \n");
 printf("                     A minimal distance d is ensured, and reached in n substeps \n");
 printf("                     When using the surfaces for subsequent volume meshing by gmsh\n");
 printf("                     this step prevent too flat tetrahedra\n");
 printf(" --fineTuneOut <d> <n> Similar to --fineTuneIn, but ensures minimal distance in the other direction\n");
 printf(" --dilate <d>        Dilate the surface by d. d < 0 means shrinking.\n");
 printf(" --intersect         If the mesh contains intersections, return value = 1.\n");
 printf(" --intersect -o fname.msh  If the mesh contains intersections, return value = 1.\n");
 printf("                     In addtion, save mesh with highlighted intersections in Gmsh format\n");
 printf("Accepted input formats are OFF, PLY and STL.\nOther formats are supported only partially.\n");
 printf("See http://jmeshlib.sourceforge.net for details on supported formats.\n");
 printf("\nIf MeshFix is used for research purposes, please cite the following paper:\n");
 printf("M. Attene - A lightweight approach to repairing digitized polygon meshes.\nThe Visual Computer, 2010. (c) Springer.\n");
 exit(0);
}

char *createFilename(const char *iname, const char *subext, const char *newextension, bool stripExt)
{
 static char tname[2048];
 char *oname = (char *)malloc(strlen(iname)+strlen(subext)+strlen(newextension)+1);
 strcpy(tname, iname);
 if(!stripExt)
 {
  int tl = strlen(tname), nl = strlen(newextension);
  stripExt = true;
  for (int n=0; n < nl; n++) stripExt &= tname[tl-nl+n] == newextension[n]; // strip if identical extension found
 }
 if(stripExt) for (int n=strlen(tname)-1; n>0; n--) if (tname[n]=='.') {tname[n] = '\0'; break;}
 sprintf(oname,"%s%s%s",tname,subext,newextension);
 return oname;
}

int main(int argc, char *argv[])
{
 JMesh::init();
 JMesh::app_name = "MeshFix";
 JMesh::app_version = "1.2-alpha";
 JMesh::app_year = "2011";
 JMesh::app_authors = "Marco Attene, Mirko Windhoff, Axel Thielscher";
 JMesh::app_maillist = "attene@ge.imati.cnr.it, mirko.windhoff@tuebingen.mpg.de, axel.thielscher@tuebingen.mpg.de";
 ExtTriMesh tin;

#ifdef DISCLAIMER
 printf("\n*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n");
 printf("This software can be used ONLY with an explicit authorization of the author.\n");
 printf("If you do not have such an authorization, you must delete this software.\n");
 printf("In no event this version of MeshFix can be redistributed.\n");
 printf("\n*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*\n");
#endif

 if (argc < 2) usage();

 float par = 0;
 unsigned numberComponentsToKeep = 1;
 bool joinOverlappingComponents = false;
// float minAllowedDistance = 0;
 bool haveJoinClosestComponents = false;
 int uniformRemeshSteps = 0, numberOfVertices = 0;
 int smoothingSteps = 0, nsteps = -1;
 double cutOuterMinDist = -1, cutInnerMinDist = -1;
 double decoupleOuterOutMinDist = -1, decoupleOuterInMinDist = -1, decoupleInnerInMinDist = -1;
 double dilateDist = 0;
 double fineTuneIn = -1, fineTuneOut = -1;
 bool clean = true;
 bool removeHandles = false;
 bool save_vrml = false;
 bool save_stl = false;
 bool save_msh = false;
 bool save_fsmesh = false;
 float xshift = 0.0; // used by saveFSMESH
 bool haveOutputFile = false;
 bool haveIntersectText = false;
 const char *outputFile;
 if (!strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) usage();
 for (int i=2; i<argc; i++)
 {
  if (!strcmp(argv[i], "-a"))
  {
   if (i<argc-1) par = (float)atof(argv[i+1]); else par = 0;
   if (par < 0) JMesh::error("Epsilon angle must be > 0.\n");
   if (par > 2) JMesh::error("Epsilon angle must be < 2 degrees.\n");
   tin.epsilon_angle = par;
   if (tin.epsilon_angle)
   {
    JMesh::acos_tolerance = asin((M_PI*tin.epsilon_angle)/180.0);
	printf("Fixing asin tolerance to %e\n",JMesh::acos_tolerance);
    i++;
   }
  }
  else if (!strcmp(argv[i], "-h") || !strcmp(argv[i], "--help")) usage();
  else if (!strcmp(argv[i], "--shells")) {
      if (i<argc-1) {
          numberComponentsToKeep = atoi(argv[i+1]);
          if (numberComponentsToKeep < 1)
              JMesh::error("# components to keep must be >= 1.\n");
          else
              i++;
      }
  }
  else if (!strcmp(argv[i], "-w") || !strcmp(argv[i], "--wrl")) save_vrml = true;
  else if (!strcmp(argv[i], "-s") || !strcmp(argv[i], "--stl")) save_stl = true;
  else if (!strcmp(argv[i], "--msh")) save_msh = true;
  else if (!strcmp(argv[i], "--fsmesh")) save_fsmesh = true;
  else if (!strcmp(argv[i], "-j")) joinOverlappingComponents = true;
  else if (!strcmp(argv[i], "-u")) {
      if (i>=argc-1 || (uniformRemeshSteps = atoi(argv[i+1]))<1)
          JMesh::error("# uniform remesh steps must be >= 1.\n");
      i++;
  }
  else if (!strcmp(argv[i], "--xshift")) {
      if (i>=argc-1) JMesh::error("error reading xshift\n");
      xshift = (float) atof(argv[i+1]);
      JMesh::info("xshift set to %f\n",xshift);
      i++;
  }
  else if (!strcmp(argv[i], "--vertices")) {
      if (i>=argc-1 || (numberOfVertices = atoi(argv[i+1]))<1)
          JMesh::error("# of vertices must be >= 0.\n");
      i++;
  }
  else if (!strcmp(argv[i], "--smooth")) {
      if (i>=argc-1 || (smoothingSteps = atoi(argv[i+1]))<1)
          JMesh::error("# smoothing steps must be >= 1.\n");
      i++;
  }
  else if (!strcmp(argv[i], "--cut-outer")) {
      if (i<argc-1) {
          cutOuterMinDist = atof(argv[i+1]);
          if (cutOuterMinDist < 0)
              JMesh::error("cutOuterMinDist must be >= 0.\n");
          else
              i++;
      }
  }
  else if (!strcmp(argv[i], "--cut-inner")) {
      if (i<argc-1) {
          cutInnerMinDist = atof(argv[i+1]);
          if (cutInnerMinDist < 0)
              JMesh::error("cutInnerMinDist must be >= 0.\n");
          else
              i++;
      }
  }
  else if (!strcmp(argv[i], "--decouple-outout")) {
      if (i<argc-1) {
          decoupleOuterOutMinDist = atof(argv[i+1]);
          if (decoupleOuterOutMinDist < 0)
              JMesh::error("decoupleOuterOutMinDist must be >= 0.\n");
          else
              i++;
      }
  }
  else if (!strcmp(argv[i], "--decouple-outin")) {
      if (i<argc-1) {
          decoupleOuterInMinDist = atof(argv[i+1]);
          if (decoupleOuterInMinDist < 0)
              JMesh::error("decoupleOuterInMinDist must be >= 0.\n");
          else
              i++;
      }
  }
  else if (!strcmp(argv[i], "--decouple-inin")) {
      if (i<argc-1) {
          decoupleInnerInMinDist = atof(argv[i+1]);
          if (decoupleInnerInMinDist < 0)
              JMesh::error("decoupleInnerInMinDist must be >= 0.\n");
          else
              i++;
      }
  }
  else if (!strcmp(argv[i], "--decouple-inin")) {
      if (i<argc-1) {
          decoupleInnerInMinDist = atof(argv[i+1]);
          if (decoupleInnerInMinDist < 0)
              JMesh::error("decoupleInMinDist must be >= 0.\n");
          else
              i++;
      }
  }
  else if (!strcmp(argv[i], "--dilate")) {
      if (i<argc-1) {
          dilateDist = atof(argv[i+1]);
          i++;
      }
  }
  else if (!strcmp(argv[i], "--fineTuneIn")) {
      if (i<argc-1) { fineTuneIn = atof(argv[i+1]);  i++; }
      if (fineTuneIn <= 0) JMesh::error("fineTuneIn MinDist must be > 0.\n");

      if (i<argc-1) { nsteps = atoi(argv[i+1]);  i++; }
      if (nsteps <= 0) JMesh::error("number of substeps must be > 0.\n");
  }
  else if (!strcmp(argv[i], "--fineTuneOut")) {
      if (i<argc-1) { fineTuneOut = atof(argv[i+1]);  i++; }
      if (fineTuneOut <= 0) JMesh::error("fineTuneOut MinDist must be > 0.\n");

      if (i<argc-1) { nsteps = atoi(argv[i+1]);  i++; }
      if (nsteps <= 0) JMesh::error("number of substeps must be > 0.\n");
  }
  else if (!strcmp(argv[i], "--remove-handles")) removeHandles = true;
  else if (!strcmp(argv[i], "--intersect")) haveIntersectText = true;
  else if (!strcmp(argv[i], "--no-clean")) clean = false;
  else if (!strcmp(argv[i], "-jc")) haveJoinClosestComponents = true;
  else if (!strcmp(argv[i], "-o")) {
      if (i<argc-1) {
          haveOutputFile = true;
          outputFile = argv[i+1];
          i++;
      }
  }
  else if (!strcmp(argv[i], "-q")) JMesh::quiet = true;
  else if (argv[i][0] == '-') JMesh::warning("%s - Unknown operation.\n",argv[i]);
 }

 printf("meshfix %s\n", argv[1]);
 // The loader performs the conversion to a set of oriented manifolds
 if (tin.load(argv[1]) != 0) JMesh::error("Can't open file '%s'.\n", argv[1]);

 // Join the second input argument if existing
 if (tin.append(argv[2]) == 0)
     JMesh::info("%s was joined.\n", argv[2]);
 input_filename = argv[1];

 // Keep only the biggest components
 tin.removeSmallestComponents( numberComponentsToKeep );

 // Fill holes by taking into account both sampling density and normal field continuity
 tin.fillSmallBoundaries(tin.E.numels(), true, true);

 if (joinOverlappingComponents) {
     tin.removeSmallestComponents(2);
     if(!tin.joinOverlappingComponentPair2()) {
         if(!joinClosestComponents(&tin, false, true, false))
             JMesh::warning("Joining didn't succeed.\n");
     } else numberComponentsToKeep = 1; // for subsequent cleaning
 }
 if (haveJoinClosestComponents)
 {
  printf("\nJoining input components ...\n");
  JMesh::begin_progress();
  while (joinClosestComponents(&tin, false, true, true)) JMesh::report_progress("Num. components: %d       ",tin.shells());
  JMesh::end_progress();
  tin.deselectTriangles();
 }

 if (removeHandles) {
     printf("Removing all handles ...\n");
     if(tin.shells() > 1)
         JMesh::warning("Remove handles works only on single component meshes. Keeping only biggest shell.\n");
     if(!tin.removeHandles()) JMesh::warning("Remove handles didn't succeed.\n");
 }

 if (dilateDist != 0.0) {
     printf("Dilating by %g.\n", dilateDist);
     tin.dilate(dilateDist);
 }
 if (cutOuterMinDist >= 0) {
     printf("Cutting triangles of the first component away, that are outside of the second one; Fill holes.\n");
     if(tin.shells() != 2) JMesh::warning("Incorrect number of components, won't cut. Having %d and should have 2.\n", tin.shells());
     else tin.cutFirstWithSecondComponent(cutOuterMinDist, true);
 }
 if (cutInnerMinDist >= 0) {
     printf("Cutting triangles of the first component away, that are inside of the second one; Fill holes.\n");
     if(tin.shells() != 2) JMesh::warning("Incorrect number of components, won't cut. Having %d and should have 2.\n", tin.shells());
     else tin.cutFirstWithSecondComponent(cutInnerMinDist, false);
 }
 if (decoupleOuterOutMinDist >= 0) {
     if(numberComponentsToKeep == 1) JMesh::warning("Use --shells 2 for decoupling.\n");
     printf("Decoupling first (outer) component from second one (move outwards). Min. distance: %g.\n", decoupleOuterOutMinDist);
     if(tin.shells() != 2) JMesh::warning("Incorrect number of components, won't decouple. Having %d and should have 2.\n", tin.shells());
     else tin.decoupleFirstFromSecondComponent(decoupleOuterOutMinDist, 15, true, true);
     numberComponentsToKeep = 1; // for subsequent cleaning
 } else if(decoupleOuterInMinDist >= 0) {
     if(numberComponentsToKeep == 1) JMesh::warning("Use --shells 2 for decoupling.\n");
     printf("Decoupling first (outer) component from second one (move inwards). Min. distance: %g.\n", decoupleInnerInMinDist);
     if(tin.shells() != 2) JMesh::warning("Incorrect number of components, won't decouple. Having %d and should have 2.\n", tin.shells());
     else tin.decoupleFirstFromSecondComponent(decoupleOuterInMinDist, 15, true, false);
     numberComponentsToKeep = 1; // for subsequent cleaning
 } else if(decoupleInnerInMinDist >= 0) {
     if(numberComponentsToKeep == 1) JMesh::warning("Use --shells 2 for decoupling.\n");
     printf("Decoupling first (inner) component from second one (move inwards). Min. distance: %g.\n", decoupleInnerInMinDist);
     if(tin.shells() != 2) JMesh::warning("Incorrect number of components, won't decouple. Having %d and should have 2.\n", tin.shells());
     else tin.decoupleFirstFromSecondComponent(decoupleInnerInMinDist, 15, false, false);
     numberComponentsToKeep = 1; // for subsequent cleaning
 }

 if (fineTuneIn >=0) {
     printf("Fine-tuning for later volume meshing ...\n");
     if(tin.shells() != 2) JMesh::warning("Incorrect number of components, won't decouple. Having %d and should have 2.\n", tin.shells());
     else tin.fineTune(fineTuneIn, nsteps, true);
 } else if (fineTuneOut >=0) {
     printf("Fine-tuning for later volume meshing ...\n");
     if(tin.shells() != 2) JMesh::warning("Incorrect number of components, won't decouple. Having %d and should have 2.\n", tin.shells());
     else tin.fineTune(fineTuneOut, nsteps, false);
 }


 if (uniformRemeshSteps) {
     printf("Uniform remeshing ...\n");
     tin.uniformRemesh(uniformRemeshSteps, numberOfVertices, tin.E.numels());
 } else if(numberOfVertices) { JMesh::warning("--vertices works only together with -u."); }

 if(smoothingSteps) {
     printf("Smoothing %d steps.\n", smoothingSteps);
     tin.laplacianSmooth(smoothingSteps, 1);
 }
 // Run geometry correction
 if (clean) {
     printf("Cleaning intersections, degeneracies ...\n");
     if (!tin.clean(20, 3, numberComponentsToKeep)) {
      fprintf(stderr,"MeshFix failed!\n");
      fprintf(stderr,"Please try manually using ReMESH: http://remesh.sourceforge.net\n");
      FILE *fp = fopen("meshfix_log.txt","a");
      fprintf(fp,"MeshFix failed on %s\n", input_filename);
      fclose(fp);
     }
 }

 if (haveIntersectText) {
     printf("Testing for intersections ...\n");
     tin.deselectTriangles();

     int its = tin.selectIntersectingTriangles();

     if(haveOutputFile) tin.saveMSH(outputFile);

     if(its) return 0;
     return 1;
 }

 char *fname = createFilename( haveOutputFile ? outputFile : argv[1], haveOutputFile ? "": "_fixed", (save_vrml? ".wrl" : (save_stl? ".stl": (save_msh? ".msh" : (save_fsmesh? ".fsmesh" : ".off")))), !haveOutputFile);
 printf("Saving output mesh to '%s'\n",fname);
 if (save_vrml)
     tin.saveVRML1(fname);
 else if (save_stl)
     tin.saveSTL(fname);
 else if (save_msh)
     tin.saveMSH(fname);
 else if (save_fsmesh)
     tin.saveFSMESH(fname,xshift);
 else
     tin.saveOFF(fname);
 return 0;
}