#include <stdlib.h>
#include <string>
#include <iostream>
#include <fstream>
#include <dirent.h>

#include <float.h>
#include <math.h>
#include <vector>

#include "graphics/internal.h"

#include "Modeler.hxx"

using namespace std;

/*
  This file holds all of the supporting structs and functions 
  for the "loadModel" function available in Modeler.hxx.
  contact: danield@andrew.cmu.edu
*/


const float MINWIDTH = 2*CatomSim::catom_radius;

//for use by STL set
struct comparePoint3D {
  bool operator()(const Point3D&p1, const Point3D&p2) const {
    
    if(!floatEq(p1.getX(), p2.getX())){
      if(floatStrictLess(p2.getX(),p1.getX())){
	return true;
      }
      else if(floatStrictLess(p1.getX(),p2.getX())){
	return false;
      }
    }
    
    if(!floatEq(p1.getY(), p2.getY())){
      if(floatStrictLess(p2.getY(),p1.getY())){
	return true;
      }
      else if(floatStrictLess(p1.getY(),p2.getY())){
	return false;
      }
    }

    if(!floatEq(p1.getZ(), p2.getZ())){
      if(floatStrictLess(p2.getZ(),p1.getZ())){
	return true;
      }
      else if(floatStrictLess(p1.getZ(),p2.getZ())){
	return false;
      }
    }

    return false;
  }
};


//used for STL set
struct comparePoint3Dptr {
  bool operator()(const Point3D* p1, const Point3D* p2) const {
    return comparePoint3D()(*p1, *p2);
  }
};


//a struct which holds a pair of Point3Ds
//used here to represent edges in 3d space
//also contains a hash function for hash_set<Pt3DPair>
struct Pt3DPair{
  Point3D start;
  Point3D end;
};


//for use by STL set
struct comparePt3DPair{
  bool operator()(const Pt3DPair&p1, const Pt3DPair&p2) const {
    if(comparePoint3D()(p1.start,p2.start)) return true;
    if(comparePoint3D()(p2.start,p1.start)) return false;
    if(comparePoint3D()(p1.end, p2.end)) return true;
    if(comparePoint3D()(p2.end, p1.end)) return false;
    return false;
  }
};


//These hold a model's data
struct modelData{
  list<vector<Point3D*> > mesh;
  list<Point3D*> points;
  Point3D largest, smallest;
  bool STLMesh;
};

//loading objs and stls into mesh data objects 
static modelData loadOBJ(string filename);
static modelData loadSTL(string filename);

//function to rescale and tanslate mesh data
//to points between (1,1,1) and (size+1,size+1,size+1)
static modelData resizeMesh(modelData model, int size);

//functions used to create a shell or a solid catom model
//from the loaded mesh data
static set<Point3D, comparePoint3D> catomizeShell(modelData model, string lattice, float tolerance, set<Point3D, comparePoint3D>, string message);
static set<Point3D, comparePoint3D> catomizeSolid(modelData model, string lattice, float tolerance, int size);

//checks if the mesh is "leaky" before creating a solid model
static bool checkMesh(modelData model);

//various 3D math helper functions
static double distancePointToLine(Point3D point, double line[6]);
static bool pointOnPolygon(Point3D point, double eqn[6], vector<Point3D> vertices, double tolerance);
static bool pointInEdgeset(Point3D point, set<Pt3DPair, comparePt3DPair> edgeset, double tolerance);


list<Point3D> loadModel(string filename, int size, string lattice, float tolerance, bool solid){
  
  struct modelData model;
  set<Point3D, comparePoint3D> catomSet;
  list<Point3D> result;
  
  //read data from file, filling points and faces
  if(filename.find(".stl")!=string::npos){
    model = loadSTL(filename);
    model.STLMesh = true;
  }
  else if(filename.find(".obj")!=string::npos){
    model = loadOBJ(filename); 
    model.STLMesh = false;
  } 

  bool watertight;
  if(solid) watertight = checkMesh(model); //no math has been done to these floats!!

  //resize the mesh to fit in a cube from (1, 1, 1) to (size + 1, size + 1, size + 1)
  model = resizeMesh(model, size);

  if(solid){
    //make sure the mesh has no "holes"
    if (watertight){  
      catomSet = catomizeSolid(model, lattice, tolerance, size);
    }
    else{ 
      display_gtk_message("Mesh check found leaks in the model. Solid fill could not be completed.");
      set<Point3D, comparePoint3D> empty;
      catomSet = empty;
    }
  }
  else{
    catomSet = catomizeShell(model, lattice, MINWIDTH, catomSet, "Creating shell (first pass)");
    if(model.STLMesh){catomSet = catomizeShell(model, lattice, tolerance, catomSet, "Creating shell (second pass)");}
  }

  //free model's memory
  list<Point3D*>::iterator ptiter;
  for(ptiter = model.points.begin(); ptiter!=model.points.end(); ptiter++){
    delete (*ptiter);
  }

  set<Point3D, comparePoint3D>::iterator citer;
  for(citer = catomSet.begin(); citer != catomSet.end(); citer++) {
    Point3D temp(citer->getX(), citer->getY(), citer->getZ());
    result.push_back(temp);	
  }

  return result;
}


static modelData loadOBJ(string filename){  
   
  struct modelData model;

  //open input obj file
  ifstream obj(filename.c_str(), ios::in);
  long filesize;
  
  {
    long begin,end;
    begin = obj.tellg();
    obj.seekg (0, ios::end);
    end = obj.tellg();
    filesize = end-begin;
    obj.seekg (0, ios::beg);
  }

  //prepare object data structures
  model.smallest.setX(DBL_MAX);
  model.smallest.setY(DBL_MAX);
  model.smallest.setZ(DBL_MAX);
  model.largest.setX(-DBL_MAX);
  model.largest.setY(-DBL_MAX);
  model.largest.setZ(-DBL_MAX);
  
  //a temporary vector to hold point data
  vector<Point3D*> points;

  //get data from the file one chunk at a time
  string chunk;
  
  //progress bar
  int bar = 0;
  update_widget("catomizeBar", "GtkProgressBar", "Loading .obj file");

  while (obj>>chunk) {		
    
    if(chunk=="v") {
      //if we find a vertex, create a point3D to hold the point
      double x, y, z;
      obj >> x;
      obj >> y;
      obj >> z;
      
      //update largests and smallests
      if (x < model.smallest.getX()) model.smallest.setX(x);				
      if (y < model.smallest.getY()) model.smallest.setY(y);
      if (z < model.smallest.getZ()) model.smallest.setZ(z);
      
      if (x > model.largest.getX()) model.largest.setX(x);				
      if (y > model.largest.getY()) model.largest.setY(y);
      if (z > model.largest.getZ()) model.largest.setZ(z);
      
      //add point to vector of points
      Point3D *point = new Point3D(x, y, z);			
      points.push_back(point);
      model.points.push_back(point);

    }
    else if (chunk=="f") {
      //if we find a face, create a vector to hold the face
      vector<Point3D*> face;
      
      //loop through the pieces of the chunk until we encounter a non-number (not a digit or '-')
      while(obj >> chunk) {
	if (*chunk.c_str() != '-' && (*chunk.c_str() < '0' || *chunk.c_str() > '9')) break;
	//we have another vertex to add; its index is the first part of the '/'-separated piece
	
	//run through the string to get this index
	const char* digit = chunk.c_str();
        int index = 0;
	int sgn = 1;
	if(*digit == '-') {
	  sgn = -1;
	  digit++;
	}
	while (*digit >= '0' && *digit <='9') {
	  index = (index * 10) + (*digit - '0');
	  digit++;
	}
	index*=sgn;
	
	//sometimes, negative indices are used
	//these are relative to the current position
	if(index<0) {
	  index += model.mesh.size();
	}
	
	face.push_back(points[index-1]);
      }//done with this chunk
     
      model.mesh.push_back(face);
      
      //since we end with a non-vertex chunk, it must be returned
      for(const char* c = chunk.c_str(); *c!='\0'; c++) obj.unget(); 		
      
      
    }//any other kind of chunk we will ignore
    
    
    //progress bar
    while(int(obj.tellg()/(filesize/100+1))>bar+1) {
      bar++;
      update_widget("catomizeBar", "GtkProgressBar", ((float)obj.tellg())/((float)filesize));
    }
    
    
  }//done with all chunks
  
  update_widget("catomizeBar", "GtkProgressBar", 0);
  //we are now done getting mesh data
  obj.close();

  return model;
}


static modelData loadSTL(string filename){  
  
  modelData model;
  
  //temporary data structure to hold point data
  set<Point3D*, comparePoint3Dptr> points;

  //open input stl file
  ifstream stl(filename.c_str(), ios::in);
  long filesize;
  
  {
    long begin,end;
    begin = stl.tellg();
    stl.seekg (0, ios::end);
    end = stl.tellg();
    filesize = end-begin;
    stl.seekg (0, ios::beg);
  }

  //prepare object data structures
  model.smallest.setX(DBL_MAX);
  model.smallest.setY(DBL_MAX);
  model.smallest.setZ(DBL_MAX);
  model.largest.setX(-DBL_MAX);
  model.largest.setY(-DBL_MAX);
  model.largest.setZ(-DBL_MAX);
  
  string test;
  stl>>test;
  if(test == "solid"){ //the text file format

    //get data from the file one chunk at a time
    string chunk;
    
    //stl files are a list of triangular faces
    vector<Point3D*> face;
    
    //progress bar
    int bar = 0;
    update_widget("catomizeBar", "GtkProgressBar", "Loading .stl file");
    
    while (stl>>chunk) {		
      if(chunk=="vertex") {
	//if we find a vertex, create a point3D to hold the point
	double x, y, z;
	stl >> x;
	stl >> y;
	stl >> z;
	
	//update largests and smallests
	if (x < model.smallest.getX()) model.smallest.setX(x);				
	if (y < model.smallest.getY()) model.smallest.setY(y);
	if (z < model.smallest.getZ()) model.smallest.setZ(z);
	
	if (x > model.largest.getX()) model.largest.setX(x);				
	if (y > model.largest.getY()) model.largest.setY(y);
	if (z > model.largest.getZ()) model.largest.setZ(z);
	
	//add point to set of points
	Point3D* point = new Point3D(x, y, z);
	
	if(points.count(point)==0){
	  points.insert(point);
	  model.points.push_back(point);
	  face.push_back(point);
	}
	else{
	  point = *(points.find(point));
	  face.push_back(point);
	}

	if(face.size()==3){
	  model.mesh.push_back(face);
	  face.clear();
	}
      }      
      
      //progress bar
      while(int(stl.tellg()/(filesize/100+1))>bar+1) {
	bar++;
	update_widget("catomizeBar", "GtkProgressBar", ((float)stl.tellg())/((float)filesize));
      }
      
    }//done with all chunks
    
    update_widget("catomizeBar", "GtkProgressBar", 0);
    //we are now done getting mesh data

    stl.close();
  }
  else{//binary file format 

    stl.close();
    //open input stl file
    ifstream stlbin(filename.c_str(), ios::in | ios::binary);
    
    //skip the 80-byte header
    char header[80];
    stlbin.read(header, 80);

    //read the 4-byte facet count
    unsigned long int faceCount;
    stlbin.read((char*)(&faceCount), 4);

    //stl files are a list of triangular faces
    vector<Point3D*> face;
    
    //progress bar
    int bar = 0;
    update_widget("catomizeBar", "GtkProgressBar", "Loading .stl file");
    
    for(unsigned long int i=0; i<faceCount; i++) {
      
      //read the normal data for the face
      char normal[12];
      stlbin.read(normal, 12);
      
      //then read 3 points per face
      for(int v=0;v<3;v++){
	
	float tx=0;
	stlbin.read((char*)(&tx), 4);
	float ty=0;
	stlbin.read((char*)(&ty), 4);
	float tz=0;
	stlbin.read((char*)(&tz), 4);

	double x = (double)tx;
	double y = (double)ty;
	double z = (double)tz;

	//update largests and smallests
	if (x < model.smallest.getX()) model.smallest.setX(x);				
	if (y < model.smallest.getY()) model.smallest.setY(y);
	if (z < model.smallest.getZ()) model.smallest.setZ(z);
	
	if (x > model.largest.getX()) model.largest.setX(x);				
	if (y > model.largest.getY()) model.largest.setY(y);
	if (z > model.largest.getZ()) model.largest.setZ(z);

	//add point to set of points
	Point3D* point = new Point3D(x, y, z);
	
	if(points.count(point)==0){
	  points.insert(point);
	  model.points.push_back(point);
	  face.push_back(point);
	}
	else{
	  point = *(points.find(point));
	  face.push_back(point);
	}
	
      }
      //push completed face into faces
      model.mesh.push_back(face);
      face.clear();
      
      //get rid of 2 padding bytes
      char padding[2];
      stlbin.read(padding, 2);
      
      //progress bar
      while(int(stlbin.tellg()/(filesize/100+1))>bar+1) {
	bar++;
	update_widget("catomizeBar", "GtkProgressBar", ((float)stlbin.tellg())/((float)filesize));
      }
    }
    //done with all faces
  
    update_widget("catomizeBar", "GtkProgressBar", 0);
    //we are now done getting mesh data
  }

  return model;
}


static modelData resizeMesh(modelData model, int size){

  //find lengths of axes
  Point3D axes(model.largest.getX() - model.smallest.getX(), model.largest.getY() - model.smallest.getY(), model.largest.getZ() - model.smallest.getZ());
  
  //find proper scaling factor
  double factor = 1;
  if(axes.getX() >= axes.getY() && axes.getX() >= axes.getZ() && axes.getX() != 0) factor = size / axes.getX();
  if(axes.getY() >= axes.getX() && axes.getY() >= axes.getZ() && axes.getY() != 0) factor = size / axes.getY();
  if(axes.getZ() >= axes.getX() && axes.getZ() >= axes.getY() && axes.getZ() != 0) factor = size / axes.getZ();

  //scale and move all points so that:
  // the least coordinates in the model become 1
  // greatest become (size + 1)
  // aspect ratio is preserved
  
  {
    Point3D minposs(1.0f,1.0f,1.0f);
    int i=0;
    
    {
      list<Point3D*>::iterator ptiter;
      for(ptiter = model.points.begin(); ptiter != model.points.end(); ptiter++) {
	*(*ptiter) -= model.smallest;
	*(*ptiter) *= factor;
	*(*ptiter) += minposs;
	double temp = (*ptiter)->getZ();
	(*ptiter)->setZ((*ptiter)->getY());
	(*ptiter)->setY(temp);
	i++;
      }	
    }
    
    model.largest -= model.smallest;
    model.largest *= factor;
    model.largest += minposs;
    double temp = model.largest.getZ();
    model.largest.setZ(model.largest.getY());
    model.largest.setY(temp);
    model.smallest.setData(minposs);
  }

  return model;
}


static set<Point3D, comparePoint3D> catomizeShell(modelData model, string lattice, float tolerance, 
						  set<Point3D, comparePoint3D> presetCatoms, string message)
{
  set<Point3D, comparePoint3D> catoms = presetCatoms;;
  set<Point3D, comparePoint3D> failed;
  set<Point3D, comparePoint3D> painted;
  
  //progress bar
  int pos=0;
  int bar=0;
  int modelmeshsize = model.mesh.size();
  update_widget("catomizeBar", "GtkProgressBar", message);  

  /*
    what we know:
    vertices holds an indexed list of vertices in x,y,z coordinates
    faces holds n-tuples of indices of vertices which make up the face
    what we want to do:
    "fill" each face with catoms
    more precisely:
    given a face, whether it has one, two, or more points,
    we wish to paint in every catom which
    1. lines up with the lattice
    2. lies within some small distance of the plane of the face
    consider that the face may be a point, plane, or line!
    1 point: just approximate point in lattice
    2 points: find line equation
    3 points: find plane equation
    3. lies "within the area created by the points and lines of the face"
    consider that the face may be a point, plane, or line!
    1 point: just approximate point in lattice
    2 points: parameterize line, project point onto it, check that parameter is in range
    3 points: Project point onto plane, drop a coordinate, use point-in polygon test
    OR: take advantage of convex-polygon restriction of OBJ 
    4. has "support" by touching at least 1 catom (this is just guaranteed by the way the algorithm works)
    technique:
    store the count and locations of each vertex in list
    go into cases:
    case 1 point: set equation to point
    case 2 points: set equation to line
    case 3 points: set equation to plane
    choose a point which satisfies 1. 2. and 3. by averaging the locations of all n points and adjusting to lattice
    this seed is added to the back of the queue
    go into cases:
    case 1 point: pull point off, fill point, done
    case 2 or 3 points: pull point off, fill point
    go through all of point's lattice neighbors
    make sure each satisfies 2 and 3 (1 and 4 by construction)
    add to queue
    loop till empty

    this section of the program makes use of the funtions approx and neighbors. approx takes a point's coordinates and
    returns a point on the cubic lattice, making it a legal point for catom placement. neighbors takes a legal lattice
    point and returns a list of coordinates of lattice points "in contact with" this lattice point, or the list of
    coordinates of all possible catoms which share kissing points with the input catom.
  */
  
  //first we loop through to each face
  list<vector<Point3D*> >::iterator fiter;
  for(fiter = model.mesh.begin(); fiter != model.mesh.end(); fiter++) {
    painted.clear();
    failed.clear();
    
    vector<Point3D*> face = *fiter;
    vector<Point3D> vertices;		
    
    Point3D average(0,0,0);
    
    //copy face points, calc average (the average should be in the plane)
    vector<Point3D*>::iterator iiter;
    for(iiter = face.begin(); iiter != face.end(); iiter++) {
      Point3D point;
      point.setData(**iiter);
      
      vertices.push_back(point);
      
      average += point;
    }
    
    int faceSize = vertices.size();
    average /= faceSize;
    //determine type and equation for face
    double eqn[6];
    
  CLASSIFY:
    switch(faceSize) {
    case 1:		//a single-point face
      //this equation is pretty simple
      eqn[0] = vertices[0].getX();
      eqn[1] = vertices[0].getY();				
      eqn[2] = vertices[0].getZ();								
      break;
      
    case 2:		//a linear face
      //this equation will be parametric, where 0<=t<=1:
      //the line is given by all (t*eqn[0]+eqn[1], t*eqn[2]+eqn[3], t*eqn[4]+eqn[5])
      eqn[0] = vertices[1].getX() - vertices[0].getX();
      eqn[1] = vertices[0].getX();
      eqn[2] = vertices[1].getY() - vertices[0].getY();
      eqn[3] = vertices[0].getY();
      eqn[4] = vertices[1].getZ() - vertices[0].getZ();
      eqn[5] = vertices[0].getZ();
      break;
      
    default:	//a polygonal face
      //to be sure that our 3 points make up a face, we must show they are not colinear
      //this normal is consistent with the right-hand rule (in STLs, this normal should face out)
      Point3D normal = (vertices[1] - vertices[0]).cross(vertices[2] - vertices[1]);
      
      //if necessary, we will eliminate the middle point from the colinear points and go back to classify again
      if(floatEq(0, normal.getR())) {
	vector<Point3D>::iterator temp = vertices.begin();
	temp++;
	vertices.erase(temp);
	goto CLASSIFY;
      }
      
      //otherwise we can construct the equation ax + by + cz +d = 0
      eqn[0] = normal.getX();
      eqn[1] = normal.getY();
      eqn[2] = normal.getZ();
      eqn[3] = -normal.dot(vertices[0]);		
      break;
    }
    
    Point3D point;
    list<Point3D> points2paint;
    
    //add the average, moved onto the lattice, as a starting point 
    points2paint.push_back(snapToLattice(average, lattice));
    
    //loop as long as there are points left to be painted
    while(!points2paint.empty()) {
      
      point = points2paint.back();
      points2paint.pop_back();
      
      if(point.getZ()>0) {
	catoms.insert(point);
	painted.insert(point);
      }
      else {
	failed.insert(point);
      }
      
      list<Point3D> neighborlist = latticeNeighbors(point, lattice);
      list<Point3D>::iterator niter;
      for(niter = neighborlist.begin(); niter != neighborlist.end(); niter++) {
	point = *niter;
	if((painted.find(point) != painted.end()) || (failed.find(point) != failed.end())) {
	  continue;
	}		
	bool inModel = false;
	
	switch(faceSize) {
	case 1:		//a single-point face
	  {
	    cout<<"point"<<endl;
	    Point3D goal(eqn[0], eqn[1], eqn[2]);
	    if(point.distanceFrom(goal) <= MINWIDTH) {
	      inModel = true;
	    }
	  }
	  break;
	  
	case 2://a linear face
	  {
	    cout<<"line"<<endl;
	    if(distancePointToLine(point, eqn) <= MINWIDTH) {
	      inModel = true;
	    }
	    break;
	  }
	default: //a polygonal face 
	  {
	    //to find the projection, we find the multiple of the normal vector which would place the point on the plane
	    double t = (-eqn[0]*point.getX() - eqn[1]*point.getY() - eqn[2]*point.getZ() - eqn[3])/(pow(eqn[0], 2) + pow(eqn[1], 2) + pow(eqn[2], 2));
	    
	    Point3D projection(point.getX() + eqn[0]*t, point.getY() + eqn[1]*t, point.getZ() + eqn[2]*t);
	    
	    double distance = fabs(point.distanceFrom(projection));
	    
	    bool alreadyIn = (catoms.find(point) != catoms.end());      
	    
	    if (distance < MINWIDTH || (distance <= tolerance && t<0 && !alreadyIn)) {
	      if (pointOnPolygon(projection, eqn, vertices, tolerance)) {
		inModel = true;
	      }
	    }
	  }					
	  break;
	}
	if(inModel) {
	  points2paint.push_back(point);
	}
	else {
	  failed.insert(point);
	}
      }//end for
      
    }//end while
    
    
    //progress bar		
    pos++;
    while(int(pos/(modelmeshsize/100+1))>bar+1) {
      bar++;
      update_widget("catomizeBar", "GtkProgressBar", ((float)pos)/((float)modelmeshsize));
    }	
  }
  update_widget("catomizeBar", "GtkProgressBar", 0);
  update_widget("catomizeBar", "GtkProgressBar", "");
  
  return catoms;
}


//NOTE: this function assumes we are following the STL format conventions:
//Specifically, the mesh is made of triangles and the mesh is closed.
static set<Point3D, comparePoint3D> catomizeSolid(modelData model, string lattice, float tolerance, int size)
{
  set<Point3D, comparePoint3D> catoms;

  //progress bar
  int bar=0;
  update_widget("catomizeBar", "GtkProgressBar", "Creating solid model");  

  float drawHeight = firstSlice(lattice);

  //start slicing loop
  do {
    //cout<<"height "<<drawHeight<<endl;
        
    //this will hold pairs of points
    //which make up edges in the slice
    list<Pt3DPair> edges;

    //loop through all faces
    //cout<<"face loop"<<endl;
    list<vector<Point3D*> >::iterator fiter;
    for(fiter = model.mesh.begin(); fiter != model.mesh.end(); fiter++) {    
      vector<Point3D*> face = *fiter;
      if(face.size() != 3){
	display_gtk_message("Non-triangular polygon found. Solid fill could not be completed.");
	return catoms;
      }
      else{
	//determine if the face is cut, and if so determine where
	//a face is cut if there are two transitions from above to below (or vice versa)
	//since these are triangles, cuts are as follows:
	//format: p1-p2 is a cut if followed by (p3)
	//over-under(*)
	//under-over(*)
	//over-on(under)
	//under-on(over)
	//over-on(on)
	//under-on(on)
	//on-over(on)
	//on-under(on)
	//an edge is cut if:
	//1: it crosses from under to over or vice versa
	//2: it goes from over to on and the next is not over (the same with unders)
	//3: it goes from on to over (or under) and is followed by on

	//Two cutting points (making one boundary) could be found on
	//two separate triangles iff they shared two points directly on the plane.
	//This problem is taken care of later by the edgeset, which guarantees uniqueness

	//go through the points and record which side each is on
	int side[3];
	{
	  vector<Point3D*>::iterator iiter;
	  int i=0;
	  for(iiter = face.begin(); iiter != face.end(); iiter++) {
	    Point3D point;
	    point.setData(**iiter);	
	    side[i] = (point.getZ()>drawHeight)?1:-1;
	    if (floatEq(drawHeight, point.getZ())) {side[i] = 0;}
	    i++;
	  }
	}

	//go through the lines and check if each one crosses
	for(int i=0; i<3; i++){
	  //cout<<side[i]<<" "<<side[(i+1)%3]<<" "<<side[(i+2)%3]<<" ";
	  if(
	     (side[i]!=0 && side[(i+1)%3]*side[i]==-1) //case 1 
	     ||
	     (side[i]!=0 && side[(i+1)%3]==0 && side[(i+2)%3]*side[i]!=1) //case 2
	     ||
	     (side[i]==0 && side[(i+1)%3]!=0 && side[(i+2)%3]==0) //case 3
	     ){
	    //if it does, add an edge to edges
	    
	    struct Pt3DPair newEdge;
	    newEdge.start.setData(*face[i]);
	    newEdge.end.setData(*face[(i+1)%3]);
	    edges.push_back(newEdge);
	    
	    //cout<<"cut";
	  }
	  //cout<<endl;
	}
      }
      //cout<<endl;
    }

    //cout<<"cut loop"<<endl;
    //loop through all of the edges we've listed
    //find the actual cutting point for each edge
    //each pair of cutting points goes into a Pt3DPair
    set<Pt3DPair, comparePt3DPair> edgeset;
    list<Pt3DPair>::iterator ei;
    {
      int i=0;
      Pt3DPair newEdge;
      for(ei = edges.begin(); ei != edges.end(); ei++){
	//retrieve coords of the 2 points
	Point3D p1(ei->start);
	Point3D p2(ei->end);
	
	//these two lines form the parametric equation
	//p1 + (p2-p1)*t, 0<=t<=1
	//we must find t where z = height in this equation
	//height = z1 + (z2-z1)*t
	//t = (height-z1)/(z2-z1)
	float t = (drawHeight-p1.getZ())/(p2.getZ()-p1.getZ());
	
	//now putting t back in gets us x and y
	float x = p1.getX()+(p2.getX()-p1.getX())*t;
	float y = p1.getY()+(p2.getY()-p1.getY())*t;
	Point3D intersection(x,y,drawHeight);
	
	if(i==0){
	  newEdge.start = intersection;
	}
	else{
	  newEdge.end = intersection;
	  if(comparePoint3D()(newEdge.start, newEdge.end)){
	    Point3D temp = newEdge.start;
	    newEdge.start = newEdge.end;
	    newEdge.end = temp;
	  }

	  edgeset.insert(newEdge);
	}
	i=-i+1;
      }
    }


    set<Point3D, comparePoint3D> failed;
    set<Point3D, comparePoint3D> painted;

    //loop through all of the points in the list
    //cout<<"fill loop"<<endl;
    set<Pt3DPair, comparePt3DPair>::iterator epi;
    for(epi = edgeset.begin(); epi != edgeset.end(); epi++){
      //use each as the start for a flood fill
      //flood fill uses point in edgeset (similar to point in polygon, but different) to make sure points are in
      //if they work and are not already in catoms, add them to painted and catoms
      //otherwise, add them to failed

      list<Point3D> points2paint;
      list<Point3D> frameworkPoints;
      
      frameworkPoints.push_back(snapToLattice(epi->start, lattice));
      frameworkPoints.push_back(snapToLattice(epi->end, lattice));

      //loop as long as there are points left to be painted
      while(!(points2paint.empty() && frameworkPoints.empty())) {
	
	Point3D point;
	if(frameworkPoints.empty()){
	  point = points2paint.back();
	  points2paint.pop_back();
	  catoms.insert(point);
	  painted.insert(point);
	}
	else {
	  point = frameworkPoints.back();
	  frameworkPoints.pop_back();
	}
      
	list<Point3D> neighborlist = latticeNeighbors(point, lattice);
	if(lattice == "bcc"){
	  list<Point3D> temp;
	  list<Point3D>::iterator niter;
	  for(niter = neighborlist.begin(); niter != neighborlist.end(); niter++) {
	    list<Point3D> secondhop = latticeNeighbors(*niter, lattice);
	    temp.splice(temp.end(), secondhop);
	  }
	  neighborlist = temp;
	}

	list<Point3D>::iterator niter;
	for(niter = neighborlist.begin(); niter != neighborlist.end(); niter++) {
	  point = *niter;
	  if((painted.find(point) != painted.end()) || (failed.find(point) != failed.end())) {
	    continue;
	  }		
	  
	  if(!floatEq(point.getZ(), drawHeight)) {
	    continue;
	  }

	  //figure out if point is in or out
	  if(pointInEdgeset(point, edgeset, tolerance)) {
	    points2paint.push_back(point);
	    if(point.getX()<0) {
	      display_gtk_message("Error: model leak during solid fill. Please report this bug!");
	      return catoms; 
	      }
	  }
	  else {
	    failed.insert(point);
	  }
	}//end for
      }//end while
      
    }

    //move to next slice
    drawHeight = nextSlice(lattice, drawHeight);

    //progress bar
    while(int((drawHeight/(size/100+1)))>bar+1) {
      bar++;
      update_widget("catomizeBar", "GtkProgressBar", ((float)drawHeight)/((float)size));
    }	
  }
  while(nextSlice(lattice, drawHeight) <= size); //end slicing
  
  update_widget("catomizeBar", "GtkProgressBar", 0);
  update_widget("catomizeBar", "GtkProgressBar", "");
  
  return catoms;
}


static bool checkMesh(modelData model)
{
  //run through all edges and make sure that each is a member of 2 faces
  //any edge which is a member of only one face is a part of a hole
 
  multiset<Pt3DPair, comparePt3DPair> counter;
  list<vector<Point3D*> >::iterator fi;
  vector<Point3D*>::iterator pi;
  
  //progress bar
  int bar = 0;
  int modelmeshsize = model.mesh.size();
  int progress=0;
  update_widget("catomizeBar", "GtkProgressBar", "Checking mesh for holes...");

  //run through faces
  for(fi = model.mesh.begin(); fi != model.mesh.end(); fi++){
  
    //for each face, iterate through points
    pi = fi->begin();

    Point3D* first = *pi;
    Point3D* prev = first;
      
    for(; pi != fi->end(); pi++){

      //for each pair of points, add an edge
      Point3D* current = *pi;

      if(prev != current){
	struct Pt3DPair ce;
	if(comparePoint3Dptr()(current, prev)){
	  ce.start.setData(*current);
	  ce.end.setData(*prev);
	}
	else{
	  ce.start.setData(*prev);
	  ce.end.setData(*current);
	}
	
	counter.insert(ce);
	prev = current;
      }

    }
  
    //finally, connect the last point to the first
    if(prev != first){
      struct Pt3DPair ce;
      if(comparePoint3Dptr()(first, prev)){
	ce.start.setData(*first);
	ce.end.setData(*prev);
      }
      else{
	ce.start.setData(*prev);
	ce.end.setData(*first);
      }
      
      counter.insert(ce);
    }

    progress++;
    //progress bar
    while(int(progress/(modelmeshsize/100+1))>bar+1) {
      bar++;
      update_widget("catomizeBar", "GtkProgressBar", ((float)progress-1)/((float)modelmeshsize));
    }
  }

  update_widget("catomizeBar", "GtkProgressBar", 0);
  update_widget("catomizeBar", "GtkProgressBar", "");

  //run through multiset, making sure each element is included more than once
  multiset<Pt3DPair, comparePt3DPair>::iterator ei;
  bool flag = true;
  for(ei = counter.begin(); ei != counter.end(); ei++) {
    //cout<<"-----------"<<endl;
    //cout<<ei->start.getX()<<", "<<ei->start.getY()<<", "<<ei->start.getZ()<<endl;
    //cout<<ei->end.getX()<<", "<<ei->end.getY()<<", "<<ei->end.getZ()<<endl;
    if(counter.count(*ei) < 2) {
      //cout<<"*"<<endl;
      flag = false;
    }
  }
  
  return flag;
}


static bool pointInEdgeset(Point3D point, set<Pt3DPair, comparePt3DPair> edgeset, double tolerance) {
  
  bool result = false;
  
  //for each line in the edgeSet, we send a ray out to the right of the point
  //if it intersects, we negate result
  //crossing an even number of times means the point is out
  //crossing an odd number of times means the point is in

  set<Pt3DPair, comparePt3DPair>::iterator ei;
  point.setZ(0);
  //cout<<"point: "<<point.getX()<<", "<<point.getY()<<", "<<point.getZ()<<endl;
  for(ei = edgeset.begin(); ei != edgeset.end(); ei++){
    Point3D lead = ei->end;
    lead.setZ(0);
    Point3D follow = ei->start;
    follow.setZ(0);
    //cout<<"line: "<<lead.getX()<<","<<lead.getY()<<" to "<<follow.getX()<<","<<follow.getY()<<endl;
    
    
    if(floatEq(lead.getY(),point.getY())) {lead.setY(point.getY());}
    if(floatEq(follow.getY(),point.getY())) {follow.setY(point.getY());}
    bool betweenY =	
      (	
       ((lead.getY() >= point.getY()) && (follow.getY() < point.getY()))
       ||
       ((follow.getY() >= point.getY()) && (lead.getY() < point.getY()))
       );
    
    
    if (betweenY) {
      bool crossX = (point.getX() < (follow.getX() - lead.getX()) * (point.getY() - lead.getY()) / (follow.getY() - lead.getY()) + lead.getX());
      if(crossX){
	result = !result;
      }
    } 
    //cout<<(result?"true":"false")<<endl;
  }
  
  return result;
}

static bool pointOnPolygon(Point3D point, double eqn[6], vector<Point3D> vertices, double tolerance) {
  //the first step is to project the face onto a plane
  //then we will use a 2d point in polygon test
  
  //we will project it along the axis with the largest normal;
  //this guarantees the largest projection
  double a = fabs(eqn[0]);
  double b = fabs(eqn[1]);
  double c = fabs(eqn[2]);
  char discard;
  if(a >= b && a >= c) {
    discard = 'x'; 
    point.setX(point.getY()); 
    point.setY(point.getZ());
  }
  else if(b >= a && b >= c) {
    discard = 'y'; 
    point.setY(point.getZ());
  }
  point.setZ(0);
  
  vector<Point3D>::iterator iter;
  for(iter = vertices.begin(); iter != vertices.end(); iter++) {
    if(discard == 'x') {
      iter->setX(iter->getY()); 
      iter->setY(iter->getZ());
    }
    else if(discard == 'y') {
      iter->setY(iter->getZ());
    }
    iter->setZ(0);
  }

  //having projected the face and the point, we apply a point-in-polygon test
  bool result = false;
  
  //for each line in the face, we send a ray out to the right of the point
  //if it intersects, we negate result
  //crossing an even number of times means the point is out
  //crossing an odd number of times means the point is in

  vector<Point3D>::iterator lead, follow;
  follow = vertices.end(); 
  follow--;
  
  for (lead = vertices.begin(); lead != vertices.end(); follow = lead++) {
    bool betweenY =	
      (	
       ((lead->getY() <= point.getY()) && (point.getY() < follow->getY()))
       ||
       ((follow->getY() <= point.getY()) && (point.getY() < lead->getY()))
      );
    if (betweenY) {
      bool crossX = (point.getX() < (follow->getX() - lead->getX()) * (point.getY() - lead->getY()) / (follow->getY() - lead->getY()) + lead->getX());
      if(crossX) 
	{
	  result = !result;
	}
    }
  }
  
  //if we are not strictly inside the polygon, we may be within tolerance of it.
  //loop thru all lines in face and check if we are close enough.
  //TODO: maybe change this, instead making each line a face?
  if(!result) {
    follow = vertices.end(); 
    follow--;   
    for (lead = vertices.begin(); lead != vertices.end(); follow = lead++) {
      double line[6];
      line[0] = lead->getX() - follow->getX();
      line[1] = follow->getX();
      line[2] = lead->getY() - follow->getY();
      line[3] = follow->getY();
      line[4] = 0;
      line[5] = 0; 			
      if (distancePointToLine(point, line)<MINWIDTH) {
	return true;
      }
    } 
  }
  
  return result;
}


static double distancePointToLine(Point3D point, double line[6]) {
  //we will project the point onto the line segment
  //then we will measure the distance from the point to the projection
  double param;
  if(floatEq(line[0]+line[2]+line[4], 0)) {
    //if the line is really a point, the point projects right onto it
    param = 0;
  }
  else {
    //find the projected length by dividing the dotproduct by the length
    //find the param by dividing by length again
    
    Point3D linevector(line[0], line[2], line[4]);
    Point3D initial(line[1], line[3], line[5]);
    param = linevector.dot(point-initial)/pow(linevector.getR(),2);
  }
 	
  //if we are off the edge one way or another, project to the endpoints
  if (param < 0) param = 0;
  if (param > 1) param = 1;	 

  Point3D projection(line[0]*param + line[1], line[2]*param + line[3], line[4]*param + line[5]);
  
  return point.distanceFrom(projection);
}

