#include <mystdlib.h>

#include <myadt.hpp>
#include <linalg.hpp>
#include <gprim.hpp>

#include <meshing.hpp>

#include "stlgeom.hpp"

namespace netgen
{


//add a point into a pointlist, return pointnumber
int AddPointIfNotExists(Array<Point3d>& ap, const Point3d& p, double eps)
{
  double eps2 = sqr(eps);
  for (int i = 1; i <= ap.Size(); i++)
    if (Dist2(ap.Get(i),p) <= eps2 ) 
      return i;
  return ap.Append(p);
}

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

double GetDistFromLine(const Point<3> & lp1, const Point<3> & lp2, 
		       Point<3> & p)
{
  Vec3d vn = lp2 - lp1;
  Vec3d v1 = p - lp1;
  Vec3d v2 = lp2 - p;

  Point3d pold = p;

  if (v2 * vn <= 0) {p = lp2; return (pold - p).Length();}
  if (v1 * vn <= 0) {p = lp1; return (pold - p).Length();}
    
  double vnl = vn.Length();
  if (vnl == 0) {return Dist(lp1,p);}

  vn /= vnl;
  p = lp1 + (v1 * vn) * vn;
  return (pold - p).Length();
};

double GetDistFromInfiniteLine(const Point<3>& lp1, const Point<3>& lp2, const Point<3>& p)
{
  Vec3d vn(lp1, lp2);
  Vec3d v1(lp1, p);

  double vnl = vn.Length();

  if (vnl == 0)
    {
      return Dist (lp1, p);
    }
  else
    {
      return Cross (vn, v1).Length() / vnl;
    }
};



//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Binary IO-Manipulation



void FIOReadInt(istream& ios, int& i)
{
  const int ilen = sizeof(int);
  
  char buf[ilen];
  for (int j = 0; j < ilen; j++)
    ios.get(buf[j]);
  memcpy(&i, &buf, ilen);
}

void FIOWriteInt(ostream& ios, const int& i)
{
  const int ilen = sizeof(int);
  
  char buf[ilen];
  memcpy(&buf, &i, ilen);

  for (int j = 0; j < ilen; j++)
    ios << buf[j];
}

void FIOReadDouble(istream& ios, double& i)
{
  const int ilen = sizeof(double);
  
  char buf[ilen];
  for (int j = 0; j < ilen; j++)
    ios.get(buf[j]);

  memcpy(&i, &buf, ilen);
}

void FIOWriteDouble(ostream& ios, const double& i)
{
  const int ilen = sizeof(double);
  
  char buf[ilen];
  memcpy(&buf, &i, ilen);

  for (int j = 0; j < ilen; j++)
    ios << buf[j];
}

void FIOReadFloat(istream& ios, float& i)
{
  const int ilen = sizeof(float);
  
  char buf[ilen];
  int j;
  for (j = 0; j < ilen; j++)
    {
      ios.get(buf[j]);
    }
  memcpy(&i, &buf, ilen);
}

void FIOWriteFloat(ostream& ios, const float& i)
{
  const int ilen = sizeof(float);
  
  char buf[ilen];
  memcpy(&buf, &i, ilen);

  for (int j = 0; j < ilen; j++)
    ios << buf[j];
}

void FIOReadString(istream& ios, char* str, int len)
{
  for (int j = 0; j < len; j++)
    ios.get(str[j]);
}

//read string and add terminating 0
void FIOReadStringE(istream& ios, char* str, int len)
{
  for (int j = 0; j < len; j++)
    ios.get(str[j]);
  str[len] = 0;
}

void FIOWriteString(ostream& ios, char* str, int len)
{
  for (int j = 0; j < len; j++)
    ios << str[j];
}


/*
void FIOReadInt(istream& ios, int& i)
{
  const int ilen = sizeof(int);
  
  char buf[ilen];
  int j;
  for (j = 0; j < ilen; j++)
    {
      ios.get(buf[ilen-j-1]);
    }
  memcpy(&i, &buf, ilen);
}

void FIOWriteInt(ostream& ios, const int& i)
{
  const int ilen = sizeof(int);
  
  char buf[ilen];
  memcpy(&buf, &i, ilen);

  int j;
  for (j = 0; j < ilen; j++)
    {
      ios << buf[ilen-j-1];
    }
}

void FIOReadDouble(istream& ios, double& i)
{
  const int ilen = sizeof(double);
  
  char buf[ilen];
  int j;
  for (j = 0; j < ilen; j++)
    {
      ios.get(buf[ilen-j-1]);
    }
  memcpy(&i, &buf, ilen);
}

void FIOWriteDouble(ostream& ios, const double& i)
{
  const int ilen = sizeof(double);
  
  char buf[ilen];
  memcpy(&buf, &i, ilen);

  int j;
  for (j = 0; j < ilen; j++)
    {
      ios << buf[ilen-j-1];
    }
}

void FIOReadFloat(istream& ios, float& i)
{
  const int ilen = sizeof(float);
  
  char buf[ilen];
  int j;
  for (j = 0; j < ilen; j++)
    {
      ios.get(buf[ilen-j-1]);
    }
  memcpy(&i, &buf, ilen);
}

void FIOWriteFloat(ostream& ios, const float& i)
{
  const int ilen = sizeof(float);
  
  char buf[ilen];
  memcpy(&buf, &i, ilen);

  int j;
  for (j = 0; j < ilen; j++)
    {
      ios << buf[ilen-j-1];
    }
}

void FIOReadString(istream& ios, char* str, int len)
{
  int j;
  for (j = 0; j < len; j++)
    {
      ios.get(str[j]);
    }
}

//read string and add terminating 0
void FIOReadStringE(istream& ios, char* str, int len)
{
  int j;
  for (j = 0; j < len; j++)
    {
      ios.get(str[j]);
    }
  str[len] = 0;
}

void FIOWriteString(ostream& ios, char* str, int len)
{
  int j;
  for (j = 0; j < len; j++)
    {
      ios << str[j];
    }
}
*/

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

STLReadTriangle :: STLReadTriangle (const Point<3> * apts,
				    const Vec<3> & anormal)
{
  pts[0] = apts[0];
  pts[1] = apts[1];
  pts[2] = apts[2]; 
  normal = anormal;
}



STLTriangle :: STLTriangle(const int * apts)
{
  pts[0] = apts[0];
  pts[1] = apts[1];
  pts[2] = apts[2];

  facenum = 0;
}

int STLTriangle :: IsNeighbourFrom(const STLTriangle& t) const
{
  //triangles must have same orientation!!!

  for(int i = 0; i <= 2; i++)
    for(int j = 0; j <= 2; j++)
      if (t.pts[(i+1)%3] == pts[j] && 
	  t.pts[i] == pts[(j+1)%3])

	return 1;

  return 0;      
}

int STLTriangle :: IsWrongNeighbourFrom(const STLTriangle& t) const
{
  //triangles have not same orientation!!!
  for(int i = 0; i <= 2; i++)
    for(int j = 0; j <= 2; j++)
      if (t.pts[(i+1)%3] == pts[(j+1)%3] &&
	  t.pts[i] == pts[j])
	
	return 1;

  return 0;      
}

void STLTriangle :: GetNeighbourPoints(const STLTriangle& t, int& p1, int& p2) const
{
  for(int i = 1; i <= 3; i++)
    for(int j = 1; j <= 3; j++)
      if (t.PNumMod(i+1) == PNumMod(j) &&
	  t.PNumMod(i) == PNumMod(j+1))
	{
	  p1 = PNumMod(j); 
	  p2 = PNumMod(j+1); 
	  return;
	}

  PrintSysError("Get neighbourpoints failed!");
}

int STLTriangle :: GetNeighbourPointsAndOpposite(const STLTriangle& t, int& p1, int& p2, int& po) const
{
  for(int i = 1; i <= 3; i++)
    for(int j = 1; j <= 3; j++)
      if (t.PNumMod(i+1) == PNumMod(j) &&
	  t.PNumMod(i) == PNumMod(j+1))
	{
	  p1 = PNumMod(j); 
	  p2 = PNumMod(j+1); 
	  po = PNumMod(j+2); 
	  return 1;
	}
  
  return 0;
}

Vec<3> STLTriangle :: GeomNormal(const Array<Point<3> >& ap) const
{
  const Point<3> & p1 = ap.Get(PNum(1));
  const Point<3> & p2 = ap.Get(PNum(2));
  const Point<3> & p3 = ap.Get(PNum(3));
  
  return Cross(p2-p1, p3-p1);
}


void STLTriangle :: SetNormal (const Vec<3> & n)
{
  double len = n.Length();
  if (len > 0)
    {
      normal = n;
      normal.Normalize();
    }
  else
    {
      normal = Vec<3> (1, 0, 0);
    }
}


void STLTriangle :: ChangeOrientation()
{ 
  normal *= -1;
  Swap(pts[0],pts[1]); 
}



double STLTriangle :: Area(const Array<Point<3> >& ap) const
{
  return 0.5 * Cross(ap.Get(PNum(2))-ap.Get(PNum(1)), 
		     ap.Get(PNum(3))-ap.Get(PNum(1))).Length();
}

double STLTriangle :: MinHeight(const Array<Point<3> >& ap) const
{
  double ml = MaxLength(ap);
  if (ml != 0) {return 2.*Area(ap)/ml;}
  PrintWarning("max Side Length of a triangle = 0!!!");
  return 0;
}

double STLTriangle :: MaxLength(const Array<Point<3> >& ap) const
{
  return max3(Dist(ap.Get(PNum(1)),ap.Get(PNum(2))),
	      Dist(ap.Get(PNum(2)),ap.Get(PNum(3))),
	      Dist(ap.Get(PNum(3)),ap.Get(PNum(1))));
}

void STLTriangle :: ProjectInPlain(const Array<Point<3> >& ap, 
				   const Vec<3> & n, Point<3> & pp) const
{
  const Point<3> & p1 = ap.Get(PNum(1));
  const Point<3> & p2 = ap.Get(PNum(2));
  const Point<3> & p3 = ap.Get(PNum(3));
  
  Vec<3> v1 = p2 - p1;
  Vec<3> v2 = p3 - p1;
  Vec<3> nt = Cross(v1, v2);

  double c = - (p1(0)*nt(0) + p1(1)*nt(1) + p1(2)*nt(2));

  double prod = n * nt;  

  if (fabs(prod) == 0) 
    {
      pp = Point<3>(1.E20,1.E20,1.E20); 
      return; 
    }

  double nfact = -(pp(0)*nt(0) + pp(1)*nt(1) + pp(2)*nt(2) + c) / (prod);
  pp = pp + (nfact) * n;

}


int STLTriangle :: ProjectInPlain (const Array<Point<3> >& ap, 
				   const Vec<3> & nproj, 
				   Point<3> & pp, Vec<3> & lam) const
{
  const Point<3> & p1 = ap.Get(PNum(1));
  const Point<3> & p2 = ap.Get(PNum(2));
  const Point<3> & p3 = ap.Get(PNum(3));
  
  Vec<3> v1 = p2-p1;
  Vec<3> v2 = p3-p1;

  Mat<3> mat;
  for (int i = 0; i < 3; i++)
    {
      mat(i,0) = v1(i);
      mat(i,1) = v2(i);
      mat(i,2) = nproj(i);
    }

  int err = 0;
  mat.Solve (pp-p1, lam);
  //  int err = SolveLinearSystem (v1, v2, nproj, pp-p1, lam);

  if (!err)
    {
      //      pp = p1 + lam(0) * v1 + lam(1) * v2;

      pp(0) = p1(0) + lam(0) * v1(0) + lam(1) * v2(0);
      pp(1) = p1(1) + lam(0) * v1(1) + lam(1) * v2(1);
      pp(2) = p1(2) + lam(0) * v1(2) + lam(1) * v2(2);
    }
  return err;
}





void STLTriangle :: ProjectInPlain(const Array<Point<3> >& ap, 
				   Point<3> & pp) const
{
  const Point<3> & p1 = ap.Get(PNum(1));
  const Point<3> & p2 = ap.Get(PNum(2));
  const Point<3> & p3 = ap.Get(PNum(3));
  
  Vec<3> v1 = p2 - p1;
  Vec<3> v2 = p3 - p1;
  Vec<3> nt = Cross(v1, v2);

  double c = - (p1(0)*nt(0) + p1(1)*nt(1) + p1(2)*nt(2));
  
  double prod = nt * nt;  

  double nfact = -(pp(0)*nt(0) + pp(1)*nt(1) + pp(2)*nt(2) + c) / (prod);

  pp = pp + (nfact) * nt;
}

int STLTriangle :: PointInside(const Array<Point<3> > & ap, 
			       const Point<3> & pp) const
{
  const Point<3> & p1 = ap.Get(PNum(1));
  const Point<3> & p2 = ap.Get(PNum(2));
  const Point<3> & p3 = ap.Get(PNum(3));
  
  Vec<3> v1 = p2 - p1;
  Vec<3> v2 = p3 - p1;
  Vec<3> v  = pp - p1;
  double det, l1, l2;
  Vec<3> ex, ey, ez;


  ez = GeomNormal(ap);
  ez /= ez.Length();
  ex = v1;
  ex /= ex.Length();
  ey = Cross (ez, ex);
  
  Vec<2> v1p(v1*ex, v1*ey);
  Vec<2> v2p(v2*ex, v2*ey);
  Vec<2> vp(v*ex, v*ey);

  det = v2p(1) * v1p(0) - v2p(0) * v1p(1);

  if (fabs(det) == 0) {return 0;}
  
  l2 = (vp(1) * v1p(0) - vp(0) * v1p(1)) / det;
  
  if (v1p(0) != 0.)
    {
      l1 = (vp(0) - l2 * v2p(0)) / v1p(0);
    }
  else if (v1p(1) != 0.)
    {
      l1 = (vp(1) - l2 * v2p(1)) / v1p(1);
    }
  else {return 0;}
  
  if (l1 >= -1E-10 && l2 >= -1E-10 && l1 + l2 <= 1.+1E-10) {return 1;}
  return 0; 
}

double STLTriangle :: GetNearestPoint(const Array<Point<3> >& ap, 
				      Point<3> & p3d) const
{
  Point<3> p = p3d;
  ProjectInPlain(ap, p);
  double dist = (p - p3d).Length();

  if (PointInside(ap, p)) {p3d = p; return dist;}
  else
    {
      Point<3> pf = 0.0;
      double nearest = 1E50;
      //int fi = 0;
      for (int j = 1; j <= 3; j++)
	{
	  p = p3d;
	  dist = GetDistFromLine(ap.Get(PNum(j)), ap.Get(PNumMod(j+1)), p);
	  if (dist < nearest)
	    {
	      nearest = dist; 
	      pf = p;
	    }
	}
      p3d = pf;
      return nearest;
    }
}

int STLTriangle :: HasEdge(int p1, int p2) const
{
  int i;
  for (i = 1; i <= 3; i++)
    {
      if (p1 == PNum(i) && p2 == PNumMod(i+1)) {return 1;}
    }
  return 0;
}

ostream& operator<<(ostream& os, const STLTriangle& t)
{
  os << "[";
  os << t[0] << ",";
  os << t[1] << ",";
  os << t[2] << "]";

  return os;
};



STLTopEdge :: STLTopEdge ()
{
  pts[0] = pts[1] = 0;
  trigs[0] = trigs[1] = 0;
  cosangle = 1;
  status = ED_UNDEFINED;
}

STLTopEdge :: STLTopEdge (int p1, int p2, int trig1, int trig2)
{ 
  pts[0] = p1; 
  pts[1] = p2; 
  trigs[0] = trig1; 
  trigs[1] = trig2; 
  cosangle = 1;
  status = ED_UNDEFINED;
}




//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//+++++++++++++++++++   STL CHART   +++++++++++++++++++++++++++++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

STLChart :: STLChart(STLGeometry * ageometry)
{
  charttrigs = new Array<int> (0,0);
  outertrigs = new Array<int> (0,0);
  ilimit = new Array<twoint> (0,0);
  olimit = new Array<twoint> (0,0);

  geometry = ageometry;

  if ( stlparam.usesearchtree == 1)
    searchtree = new Box3dTree (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
				geometry->GetBoundingBox().PMax() + Vec3d(1,1,1));
  else
    searchtree = NULL;
}

void STLChart :: AddChartTrig(int i)
{
  charttrigs->Append(i);
  
  const STLTriangle & trig = geometry->GetTriangle(i);
  const Point3d & p1 = geometry->GetPoint (trig.PNum(1));
  const Point3d & p2 = geometry->GetPoint (trig.PNum(2));
  const Point3d & p3 = geometry->GetPoint (trig.PNum(3));

  Point3d pmin(p1), pmax(p1);
  pmin.SetToMin (p2);
  pmin.SetToMin (p3);
  pmax.SetToMax (p2);
  pmax.SetToMax (p3);
  
  if (!geomsearchtreeon && (stlparam.usesearchtree == 1))
    {searchtree->Insert (pmin, pmax, i);}
}

void STLChart :: AddOuterTrig(int i)
{
  outertrigs->Append(i);

  const STLTriangle & trig = geometry->GetTriangle(i);
  const Point3d & p1 = geometry->GetPoint (trig.PNum(1));
  const Point3d & p2 = geometry->GetPoint (trig.PNum(2));
  const Point3d & p3 = geometry->GetPoint (trig.PNum(3));

  Point3d pmin(p1), pmax(p1);
  pmin.SetToMin (p2);
  pmin.SetToMin (p3);
  pmax.SetToMax (p2);
  pmax.SetToMax (p3);
  
  if (!geomsearchtreeon && (stlparam.usesearchtree==1))
    {searchtree->Insert (pmin, pmax, i);}
}

int STLChart :: IsInWholeChart(int nr) const
{
  for (int i = 1; i <= charttrigs->Size(); i++)
    if (charttrigs->Get(i) == nr) return 1;

  for (int i = 1; i <= outertrigs->Size(); i++)
    if (outertrigs->Get(i) == nr) return 1;

  return 0;
}

void STLChart :: GetTrianglesInBox (const Point3d & pmin,
				    const Point3d & pmax,
				    Array<int> & trias) const
{
  if (geomsearchtreeon) {PrintMessage(5,"geomsearchtreeon is set!!!");}

  if (searchtree)
    searchtree -> GetIntersecting (pmin, pmax, trias);
  else
    {
      Box3d box1(pmin, pmax);
      box1.Increase (1e-4);
      Box3d box2;

      trias.SetSize(0);
      
      int nt = GetNT();
      for (int i = 1; i <= nt; i++)
	{
	  int trignum = GetTrig(i);
	  const STLTriangle & trig = geometry->GetTriangle(trignum);
	  box2.SetPoint (geometry->GetPoint (trig.PNum(1)));
	  box2.AddPoint (geometry->GetPoint (trig.PNum(2)));
	  box2.AddPoint (geometry->GetPoint (trig.PNum(3)));
	  
	  if (box1.Intersect (box2))
	    trias.Append (trignum);
	}
    }
}

//trigs may contain the same triangle double
void STLChart :: MoveToOuterChart(const Array<int>& trigs)
{
  if (!trigs.Size()) {return;}
  for (int i = 1; i <= trigs.Size(); i++)
    {
      if (charttrigs->Get(trigs.Get(i)) != -1) 
	{AddOuterTrig(charttrigs->Get(trigs.Get(i)));}
      charttrigs->Elem(trigs.Get(i)) = -1;
    }
  DelChartTrigs(trigs);
}

//trigs may contain the same triangle double
void STLChart :: DelChartTrigs(const Array<int>& trigs)
{
  if (!trigs.Size()) {return;}

  for (int i = 1; i <= trigs.Size(); i++)
    {
      charttrigs->Elem(trigs.Get(i)) = -1;
    }

  int cnt = 0;
  for (int i = 1; i <= charttrigs->Size(); i++)
    {
      if (charttrigs->Elem(i) == -1)
	{
	  cnt++;
	}
      if (cnt != 0 && i < charttrigs->Size())
	{
	  charttrigs->Elem(i-cnt+1) = charttrigs->Get(i+1);
	}
    }
  int i = charttrigs->Size() - trigs.Size();
  charttrigs->SetSize(i);

  if (!geomsearchtreeon && stlparam.usesearchtree == 1)
    {
      PrintMessage(7, "Warning: unsecure routine due to first use of searchtrees!!!");
      //bould new searchtree!!!
      searchtree = new Box3dTree (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
				  geometry->GetBoundingBox().PMax() + Vec3d(1,1,1));

      for (int i = 1; i <= charttrigs->Size(); i++)
	{
	  const STLTriangle & trig = geometry->GetTriangle(i);
	  const Point3d & p1 = geometry->GetPoint (trig.PNum(1));
	  const Point3d & p2 = geometry->GetPoint (trig.PNum(2));
	  const Point3d & p3 = geometry->GetPoint (trig.PNum(3));
	  
	  Point3d pmin(p1), pmax(p1);
	  pmin.SetToMin (p2);
	  pmin.SetToMin (p3);
	  pmax.SetToMax (p2);
	  pmax.SetToMax (p3);
	  
	  searchtree->Insert (pmin, pmax, i);	  
	}
    }
}


void STLChart :: SetNormal (const Point<3> & apref, const Vec<3> & anormal)
{
  pref = apref;
  normal = anormal;
  double len = normal.Length();
  if (len) normal /= len;
  else normal = Vec<3> (1, 0, 0);

  t1 = normal.GetNormal ();
  t2 = Cross (normal, t1);
}

/*
Point<2> STLChart :: Project2d (const Point<3> & p3d) const
{
  Vec<3> v = p3d-pref;
  return Point<2> (t1 * v, t2 * v);
}
*/



/*
  Point3d p1, p2, center;
  double rad;
  int i1, i2;
public:
*/

/*
STLBoundarySeg :: 
STLBoundarySeg (int ai1, int ai2, const Array<Point<3> > & points,
		const STLChart * chart)
{
  i1 = ai1;
  i2 = ai2; 
  p1 = points.Get(i1);
  p2 = points.Get(i2);
  center = ::netgen::Center (p1, p2);
  rad = Dist (p1, center);

  p2d1 = chart->Project2d (p1);
  p2d2 = chart->Project2d (p2);
  
  boundingbox.Set (p2d1);
  boundingbox.Add (p2d2);
}
*/

void STLBoundarySeg :: Swap ()
{
  ::netgen::Swap (i1, i2);
  ::netgen::Swap (p1, p2);
}



STLBoundary :: STLBoundary (STLGeometry * ageometry)
  : // boundary(), 
    geometry(ageometry)
{
  ;
}


void STLBoundary :: AddOrDelSegment(const STLBoundarySeg & seg)
{
  int i;
  int found = 0;
  for (i = 1; i <= boundary.Size(); i++)
    {
      if (found) {boundary.Elem(i-1) = boundary.Get(i);}
      if (boundary.Get(i) == seg) {found = 1;}
    }
  if (!found) 
    {
      boundary.Append(seg);
    }
  else 
    {
      boundary.SetSize(boundary.Size()-1);
    }
}

void STLBoundary ::AddTriangle(const STLTriangle & t)
{
  int i;
  int found1 = 0;
  int found2 = 0;
  int found3 = 0;
  //int offset = 0;
  

  STLBoundarySeg seg1(t[0],t[1], geometry->GetPoints(), chart);
  STLBoundarySeg seg2(t[1],t[2], geometry->GetPoints(), chart);
  STLBoundarySeg seg3(t[2],t[0], geometry->GetPoints(), chart);

  seg1.SetSmoothEdge (geometry->IsSmoothEdge (seg1.I1(), seg1.I2()));
  seg2.SetSmoothEdge (geometry->IsSmoothEdge (seg2.I1(), seg2.I2()));
  seg3.SetSmoothEdge (geometry->IsSmoothEdge (seg3.I1(), seg3.I2()));

  /*
  for (i = 1; i <= boundary.Size(); i++)
    {
      if (offset) {boundary.Elem(i-offset) = boundary.Get(i);}
      if (boundary.Get(i) == seg1) {found1 = 1; offset++;}
      if (boundary.Get(i) == seg2) {found2 = 1; offset++;}
      if (boundary.Get(i) == seg3) {found3 = 1; offset++;}
    }

  if (offset)
    {
      boundary.SetSize(boundary.Size()-offset);
    }    
  */
  for (i = boundary.Size(); i >= 1; i--)
    {
      if (boundary.Get(i) == seg1) 
	{ boundary.DeleteElement (i); found1 = 1; } 
      else if (boundary.Get(i) == seg2) 
	{ boundary.DeleteElement (i); found2 = 1; } 
      else if (boundary.Get(i) == seg3) 
	{ boundary.DeleteElement (i); found3 = 1; } 
    }

  if (!found1) {seg1.Swap(); boundary.Append(seg1);}
  if (!found2) {seg2.Swap(); boundary.Append(seg2);}
  if (!found3) {seg3.Swap(); boundary.Append(seg3);}
}

int STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3> & sn, 
			   double sinchartangle, int divisions, Array<Point<3> >& points, double eps)
{
  if (usechartnormal)
    return TestSegChartNV (p1, p2, sn);

#ifdef NONE
  // for statistics
  {
    int i;
    static Array<int> cntclass;
    static int cnt = 0;
    static int cnti = 0, cnto = 0;
    static long int cntsegs = 0;
    if (cntclass.Size() == 0)
      {
	cntclass.SetSize (20);
	for (i = 1; i <= cntclass.Size(); i++)
	  cntclass.Elem(i) = 0;
      }
    
    cntsegs += NOSegments();
    int cla = int (log (double(NOSegments()+1)) / log(2.0));
    if (cla < 1) cla = 1;
    if (cla > cntclass.Size()) cla = cntclass.Size();
    cntclass.Elem(cla)++;
    cnt++;
    if (divisions)
      cnti++;
    else
      cnto++;
    if (cnt > 100000) 
      {
	cnt = 0;
	/*
	(*testout) << "TestSeg-calls for classes:" << endl;
	(*testout) << cnti << " inner calls, " << cnto << " outercalls" << endl;
	(*testout) << "total testes segments: " << cntsegs << endl;
	for (i = 1; i <= cntclass.Size(); i++)
	  {
	    (*testout) << int (exp (i * log(2.0))) << " bnd segs: " << cntclass.Get(i) << endl;
	  }
	*/
      }
  }
#endif

  int i,j,k;
  Point<3> seg1p/*, seg2p*/;
  Point<3> sp1,sp2;
  double lambda1, lambda2, vlen2;
  Vec<3> vptpl;
  double sinchartangle2 = sqr(sinchartangle);
  double scal;
  int possible;

  //double maxval = -1;
  //double maxvalnew = -1;



  double scalp1 = p1(0) * sn(0) + p1(1) * sn(1) + p1(2) * sn(2);
  double scalp2 = p2(0) * sn(0) + p2(1) * sn(1) + p2(2) * sn(2);
  double minl = min2(scalp1, scalp2);
  double maxl = max2(scalp1, scalp2);
  Point<3> c = Center (p1, p2);
  double dist1 = Dist (c, p1);
 
  int nseg = NOSegments();
  for (j = 1; j <= nseg; j++)
    {
      const STLBoundarySeg & seg = GetSegment(j);


      if (seg.IsSmoothEdge())
	continue;


      sp1 = seg.P1();
      sp2 = seg.P2();

      // Test, ob Spiral Konfikt moeglich
      
      possible = 1;

      double scalsp1 = sp1(0) * sn(0) + sp1(1) * sn(1) + sp1(2) * sn(2);
      double scalsp2 = sp2(0) * sn(0) + sp2(1) * sn(1) + sp2(2) * sn(2);

      double minsl = min2(scalsp1, scalsp2);
      double maxsl = max2(scalsp1, scalsp2);
      
      double maxdiff = max2 (maxsl - minl, maxl - minsl);
      
      /*
      Point3d sc = Center (sp1, sp2);
      double mindist = Dist(c, sc) - dist1 - GetSegment(j).Radius();
      if (maxdiff < sinchartangle * mindist)
	{
	  possible = 0;
	}
      */
       
      double hscal = maxdiff + sinchartangle * (dist1 + seg.Radius());
      if (hscal * hscal < sinchartangle * Dist2(c, seg.center ))
	possible = 0;


      /*      
      if (possible)
	{
	  double mindist2ex = MinDistLL2 (p1, p2, sp1, sp2);
	  if (maxdiff * maxdiff < sinchartangle2 * mindist2ex)
	    possible = 0;
	}
      */

      if (possible)
      	{
	  LinearPolynomial2V lp (scalp1 - scalsp1,
				 scalp2 - scalp1,
				 -(scalsp2 - scalsp1));
	  QuadraticPolynomial2V slp;
	  slp.Square (lp);
	  
      
	  Vec3d v (p1, sp1);
	  Vec3d vl (p1, p2);
	  Vec3d vsl (sp1, sp2);
      
	  QuadraticPolynomial2V qp (v.Length2(),
				    -2 * (v * vl),
				    2 * (v * vsl),
				    vl.Length2(),
				    -2 * (vl * vsl),
				    vsl.Length2());
	  
	  slp.Add (-sinchartangle2, qp);

	  double hv = slp.MaxUnitSquare();

	  if (hv > eps) return 0;
	  /*
	  if (hv > maxvalnew)
	    maxvalnew = hv;
	  */
	}
      

      // if (possible && 0)
      if (false)

	for (i = 0; i <= divisions; i++)
	  {
	    
	    lambda1 = (double)i/(double)divisions;
	    seg1p = Point3d(p1(0)*lambda1+p2(0)*(1.-lambda1),
			    p1(1)*lambda1+p2(1)*(1.-lambda1),
			    p1(2)*lambda1+p2(2)*(1.-lambda1));
	    

	    
	    for (k = 0; k <= divisions; k++)
	      {
		lambda2 = (double)k/(double)divisions;
		vptpl = Vec3d(sp1(0)*lambda2+sp2(0)*(1.-lambda2)-seg1p(0),
			      sp1(1)*lambda2+sp2(1)*(1.-lambda2)-seg1p(1),
			      sp1(2)*lambda2+sp2(2)*(1.-lambda2)-seg1p(2));
		
		vlen2 = vptpl.Length2();

		//		if (vlen2 > 0)
		  {
		    scal = vptpl * sn;
		    double hv = scal*scal - sinchartangle2*vlen2;



		    /*
		    if (hv > maxval)
		      maxval = hv;
		    */
		    if (hv > eps) return 0;
		  }
	      } 
	  }
    }
  
  return 1;
  //  return (maxvalnew < eps);
}



// checks, whether 2d projection intersects
int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2, 
				  const Vec3d& sn)
{
  int timer = NgProfiler::CreateTimer ("TestSegChartNV");
  NgProfiler::StartTimer (timer);

  int nseg = NOSegments();

  Point<2> p2d1 = chart->Project2d (p1);
  Point<2> p2d2 = chart->Project2d (p2);

  Box<2> box2d;
  box2d.Set (p2d1);
  box2d.Add (p2d2);

  Line2d l1 (p2d1, p2d2);

  double eps = 1e-3;
  bool ok = true;

  /*
  static long int cnt = 0;
  static long int totnseg = 0;
  totnseg += nseg;
  cnt++;
  if ( (cnt % 100000) == 0)
  cout << "avg nseg = " << double(totnseg)/cnt << endl;
  */

  for (int j = 1; j <= nseg; j++)
    {
      const STLBoundarySeg & seg = GetSegment(j);

      if (seg.IsSmoothEdge()) continue;
      if (!box2d.Intersect (seg.BoundingBox())) continue;

      const Point<2> & sp1 = seg.P2D1();
      const Point<2> & sp2 = seg.P2D2();
      
      Line2d l2 (sp1, sp2);
      double lam1, lam2;
      
      int err = CrossPointBarycentric (l1, l2, lam1, lam2);

      if (!err && lam1 > eps && lam1 < 1-eps &&
	  lam2 > eps && lam2 < 1-eps)
        {
          ok = false;
          break;
        }
    }

  NgProfiler::StopTimer (timer);

  return ok;
}



STLDoctorParams :: STLDoctorParams()
{
  drawmeshededges = 1;
  geom_tol_fact = 1E-6;
  longlinefact = 0;
  showexcluded = 1;

  selectmode = 0;
  edgeselectmode = 0;
  useexternaledges = 0;
  showfaces = 0;
  showtouchedtrigchart = 1;
  showedgecornerpoints = 1;
  conecheck = 1;
  spiralcheck = 1;
  selecttrig = 0;
  nodeofseltrig = 1;
  selectwithmouse = 1;
  showmarkedtrigs = 1;
  dirtytrigfact = 0.001;
  smoothangle = 90;
  smoothnormalsweight = 0.2;
  vicinity = 0;
  showvicinity = 0;
}



STLDoctorParams stldoctor;

void STLDoctorParams :: Print (ostream & ost) const
{
  ost << "STL doctor parameters:" << endl
      << "selecttrig = " << selecttrig << endl
      << "selectlocalpoint = " << nodeofseltrig << endl
      << "selectwithmouse = " << selectwithmouse << endl
      << "showmarkedtrigs = " << showmarkedtrigs << endl
      << "dirtytrigfact = " << dirtytrigfact << endl
      << "smoothangle = " << smoothangle << endl;
}


STLParameters ::   STLParameters()
{
  yangle = 30;
  contyangle = 20;
  edgecornerangle = 60;
  chartangle = 15;
  outerchartangle = 70;
     
  usesearchtree = 0;
  atlasminh = 1E-4;
  resthsurfcurvfac = 2;
  resthsurfcurvenable = 0;
  resthatlasfac = 2;
  resthatlasenable = 1;
  resthchartdistfac = 1.2;
  resthchartdistenable = 1;
  resthlinelengthfac = 0.5;
  resthlinelengthenable = 1;
  resthcloseedgefac = 1;
  resthcloseedgeenable = 1;
  resthedgeanglefac = 1;
  resthedgeangleenable = 0;
  resthsurfmeshcurvfac = 1;
  resthsurfmeshcurvenable = 0;
  recalc_h_opt = 1;
}

void STLParameters :: Print (ostream & ost) const
{
  ost << "STL parameters:" << endl
      << "yellow angle = " << yangle << endl
      << "continued yellow angle = " << contyangle << endl
      << "edgecornerangle = " << edgecornerangle << endl
      << "chartangle = " << chartangle << endl
      << "outerchartangle = " << outerchartangle << endl
      << "restrict h due to ..., enable and safety factor: " << endl
      << "surface curvature: " << resthsurfcurvenable
      << ", fac = " << resthsurfcurvfac << endl
      << "atlas surface curvature: " << resthatlasenable
      << ", fac = " << resthatlasfac << endl
      << "chart distance: " << resthchartdistenable
      << ", fac = " << resthchartdistfac << endl
      << "line length: " << resthlinelengthenable
      << ", fac = " << resthlinelengthfac << endl
      << "close edges: " << resthcloseedgeenable
      << ", fac = " << resthcloseedgefac << endl
      << "edge angle: " << resthedgeangleenable
      << ", fac = " << resthedgeanglefac << endl;
}


STLParameters stlparam;


}