/contrib/Netgen/libsrc/stlgeom/stlgeom.cpp

Large files files are truncated, but you can click here to view the full file

#include <meshing.hpp>

#include "stlgeom.hpp"

namespace netgen
{

//globalen searchtree fuer gesamte geometry aktivieren
int geomsearchtreeon = 0;

int usechartnormal = 1;  

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

void STLMeshing (STLGeometry & geom,
		 Mesh & mesh)
{
  geom.Clear();
  geom.BuildEdges();
  geom.MakeAtlas(mesh);
  geom.CalcFaceNums();
  geom.AddFaceEdges();
  geom.LinkEdges();

  mesh.ClearFaceDescriptors();
  for (int i = 1; i <= geom.GetNOFaces(); i++)
    mesh.AddFaceDescriptor (FaceDescriptor (i, 1, 0, 0));
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//+++++++++++++++++++   STL GEOMETRY   ++++++++++++++++++++++++++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


  STLGeometry :: STLGeometry()
  /*
  : edges(), edgesperpoint(),
    normals(),  externaledges(),
    atlas(), chartmark(), 
    lines(), outerchartspertrig(), vicinity(), markedtrigs(), markedsegs(),
    lineendpoints(), spiralpoints(), selectedmultiedge()
  */
{
	edgedata = new STLEdgeDataList(*this);
  externaledges.SetSize(0);
  Clear();
  meshchart = 0; // initialize all ?? JS

  if (geomsearchtreeon)
    searchtree = new Box3dTree (GetBoundingBox().PMin() - Vec3d(1,1,1),
				GetBoundingBox().PMax() + Vec3d(1,1,1));
  else
    searchtree = NULL;

  status = STL_GOOD;
  statustext = "Good Geometry";
  smoothedges = NULL;
}

STLGeometry :: ~STLGeometry()
{
  delete edgedata;
}

void STLGeometry :: Save (string filename) const
{
  const char  * cfilename = filename.c_str();
  if (strlen(cfilename) < 4) 
    throw NgException ("illegal filename");

  if (strlen(cfilename) > 3 &&
      strcmp (&cfilename[strlen(cfilename)-3], "stl") == 0)
    {
      STLTopology::Save (cfilename);
    }
  else if (strlen(cfilename) > 4 &&
	   strcmp (&cfilename[strlen(cfilename)-4], "stlb") == 0)
    {
      SaveBinary (cfilename,"Binary STL Geometry");
      
    }
  else if (strlen(cfilename) > 4 &&
	   strcmp (&cfilename[strlen(cfilename)-4], "stle") == 0)
    {
      SaveSTLE (cfilename);
    }
}



int STLGeometry :: GenerateMesh (Mesh*& mesh, MeshingParameters & mparam,
				 int perfstepsstart, int perfstepsend)
{
  return STLMeshingDummy (this, mesh, mparam, perfstepsstart, perfstepsend);
}


const Refinement & STLGeometry :: GetRefinement () const
{
  return RefinementSTLGeometry (*this);
}



void STLGeometry :: STLInfo(double* data)
{
  data[0] = GetNT();

  Box<3> b = GetBoundingBox();
  data[1] = b.PMin()(0);
  data[2] = b.PMax()(0);
  data[3] = b.PMin()(1);
  data[4] = b.PMax()(1);
  data[5] = b.PMin()(2);
  data[6] = b.PMax()(2);

  int i;
 
  int cons = 1;
  for (i = 1; i <= GetNT(); i++)
    {
      if (NONeighbourTrigs(i) != 3) {cons = 0;}
    }
  data[7] = cons;
}

void STLGeometry :: MarkNonSmoothNormals()
{

  PrintFnStart("Mark Non-Smooth Normals");

  int i,j;

  markedtrigs.SetSize(GetNT());

  for (i = 1; i <= GetNT(); i++)
    {
      SetMarkedTrig(i, 0);
    }

  double dirtyangle = stlparam.yangle/180.*M_PI;

  int cnt = 0;
  int lp1,lp2;
  for (i = 1; i <= GetNT(); i++)
    {
      for (j = 1; j <= NONeighbourTrigs(i); j++)
	{
	  if (GetAngle(i, NeighbourTrig(i,j)) > dirtyangle)
	    {
	      GetTriangle(i).GetNeighbourPoints(GetTriangle(NeighbourTrig(i,j)), lp1, lp2);
	      if (!IsEdge(lp1,lp2))
		{
		  if (!IsMarkedTrig(i)) {SetMarkedTrig(i,1); cnt++;}
		}
	    }
	}
    }

  PrintMessage(5,"marked ",cnt," non-smooth trig-normals");

}

void STLGeometry :: SmoothNormals()
{
  multithread.terminate = 0;

  //  UseExternalEdges();

  BuildEdges();


  DenseMatrix m(3), hm(3);
  Vector rhs(3), sol(3), hv(3), hv2(3);

  Vec<3> ri;

  double wnb = stldoctor.smoothnormalsweight;   // neigbour normal weight
  double wgeom = 1-wnb;   // geometry normal weight


  // minimize 
  //  wgeom sum_T  \sum ri  \| ri^T (n - n_geom) \|^2  
  //  + wnb sum_SE  \| ri x (n - n_nb) \|^2
  
  int i, j, k, l;
  int nt = GetNT();
  
  PushStatusF("Smooth Normals");
    
  //int testmode;

  for (i = 1; i <= nt; i++)
    {

      SetThreadPercent( 100.0 * (double)i / (double)nt);

      const STLTriangle & trig = GetTriangle (i);
      
      m = 0;
      rhs = 0;

      // normal of geometry:
      Vec<3> ngeom = trig.GeomNormal(points);
      ngeom.Normalize();

      for (j = 1; j <= 3; j++)
	{ 
	  int pi1 = trig.PNumMod (j);
	  int pi2 = trig.PNumMod (j+1);

	  // edge vector
	  ri = GetPoint (pi2) - GetPoint (pi1);
	  
	  for (k = 0; k < 3; k++)
	    for (l = 0; l < 3; l++)
	      hm.Elem(k+1, l+1) = wgeom * ri(k) * ri(l);
	  
	  
	  for (k = 0; k < 3; k++)
	    hv(k) = ngeom(k);
	  
	  hm.Mult (hv, hv2);
	  /*
	  if (testmode)
	    (*testout) << "add vec " << hv2 << endl 
		       << " add m " << hm << endl;
	  */
	  rhs.Add (1, hv2);
	  m += hm;


	  int nbt = 0;
	  int fp1,fp2;
	  for (k = 1; k <= NONeighbourTrigs(i); k++)
	    {
	      trig.GetNeighbourPoints(GetTriangle(NeighbourTrig(i, k)),fp1,fp2);
	      if (fp1 == pi1 && fp2 == pi2)
		{
		  nbt = NeighbourTrig(i, k);
		}
	    }

	  if (!nbt)
	    {
	      cerr << "ERROR: stlgeom::Smoothnormals, nbt = 0" << endl;
	    }

	  // smoothed normal
	  Vec<3> nnb = GetTriangle(nbt).Normal();   // neighbour normal
	  nnb.Normalize();

	  if (!IsEdge(pi1,pi2)) 
	    {
	      double lr2 = ri * ri;
	      for (k = 0; k < 3; k++)
		{
		  for (l = 0; l < k; l++)
		    {
		      hm.Elem(k+1, l+1) = -wnb * ri(k) * ri(l);
		      hm.Elem(l+1, k+1) = -wnb * ri(k) * ri(l);
		    }
		  
		  hm.Elem(k+1, k+1) = wnb * (lr2 - ri(k) * ri(k));
		}
	      
	      for (k = 0; k < 3; k++)
		hv(k) = nnb(k);
	      
	      hm.Mult (hv, hv2);
	      /*
	      if (testmode)
		(*testout) << "add nb vec " << hv2 << endl 
			   << " add nb m " << hm << endl;
	      */

	      rhs.Add (1, hv2);
	      m += hm;
	    }
	}

      m.Solve (rhs, sol);
      Vec3d newn(sol(0), sol(1), sol(2));
      newn /= (newn.Length() + 1e-24);      

      GetTriangle(i).SetNormal(newn);
      // setnormal (sol);
    }

  /*
  for (i = 1; i <= nt; i++)
    SetMarkedTrig(i, 0);      		



  int crloop;
  for (crloop = 1; crloop <= 3; crloop++)
    {

  // find critical:

  Array<INDEX_2> critpairs;
  for (i = 1; i <= nt; i++)
    {
      const STLTriangle & trig = GetTriangle (i);
      
      Vec3d ngeom = GetTriangleNormal (i); // trig.Normal(points);
      ngeom /= (ngeom.Length() + 1e-24);

      for (j = 1; j <= 3; j++)
	{ 
	  int pi1 = trig.PNumMod (j);
	  int pi2 = trig.PNumMod (j+1);

	  int nbt = 0;
	  int fp1,fp2;
	  for (k = 1; k <= NONeighbourTrigs(i); k++)
	    {
	      trig.GetNeighbourPoints(GetTriangle(NeighbourTrig(i, k)),fp1,fp2);
	      if (fp1 == pi1 && fp2 == pi2)
		{
		  nbt = NeighbourTrig(i, k);
		}
	    }
	  
	  if (!nbt)
	    {
	      cerr << "ERROR: stlgeom::Smoothnormals, nbt = 0" << endl;
	    }

	  Vec3d nnb = GetTriangleNormal(nbt);   // neighbour normal
	  nnb /= (nnb.Length() + 1e-24);

	  if (!IsEdge(pi1,pi2)) 
	    {
	      if (Angle (nnb, ngeom) > 150 * M_PI/180)
		{
		  SetMarkedTrig(i, 1);      		
		  SetMarkedTrig(nbt, 1);      		
		  critpairs.Append (INDEX_2 (i, nbt));
		}
	    }

	}
    }

  if (!critpairs.Size())
    {
      break;
    }

  if (critpairs.Size())
    {

      Array<int> friends;
      double area1 = 0, area2 = 0;

      for (i = 1; i <= critpairs.Size(); i++)
	{
	  int tnr1 = critpairs.Get(i).I1();
	  int tnr2 = critpairs.Get(i).I2();
	  (*testout) << "t1 = " << tnr1 << ", t2 = " << tnr2
		     << " angle = " << Angle (GetTriangleNormal (tnr1),
					      GetTriangleNormal (tnr2))
		     << endl;

	  // who has more friends ?
	  int side;
	  area1 = 0;
	  area2 = 0;
	  for (side = 1; side <= 2; side++)
	    {
	      friends.SetSize (0);
	      friends.Append ( (side == 1) ? tnr1 : tnr2);

	      for (j = 1; j <= 3; j++)
		{
		  int fsize = friends.Size();
		  for (k = 1; k <= fsize; k++)
		    {
		      int testtnr = friends.Get(k);
		      Vec3d ntt = GetTriangleNormal(testtnr);
		      ntt /= (ntt.Length() + 1e-24);
		      
		      for (l = 1; l <= NONeighbourTrigs(testtnr); l++)
			{
			  int testnbnr = NeighbourTrig(testtnr, l);
			  Vec3d nbt = GetTriangleNormal(testnbnr);
			  nbt /= (nbt.Length() + 1e-24);

			  if (Angle (nbt, ntt) < 15 * M_PI/180)
			    {
			      int ii;
			      int found = 0;
			      for (ii = 1; ii <= friends.Size(); ii++)
				{
				  if (friends.Get(ii) == testnbnr)
				    {
				      found = 1;
				      break;
				    }
				}
			      if (!found)
				friends.Append (testnbnr);
			    }
			}
		    }
		}

	      // compute area:
	      for (k = 1; k <= friends.Size(); k++)
		{
		  double area = 
		    GetTriangle (friends.Get(k)).Area(points);

		  if (side == 1)
		    area1 += area;
		  else
		    area2 += area;
		}
	      
	    }

	  (*testout) << "area1 = " << area1 << " area2 = " << area2 << endl;
	  if (area1 < 0.1 * area2)
	    {
	      Vec3d n = GetTriangleNormal (tnr1);
	      n *= -1;
	      SetTriangleNormal(tnr1, n);
	    }
	  if (area2 < 0.1 * area1)
	    {
	      Vec3d n = GetTriangleNormal (tnr2);
	      n *= -1;
	      SetTriangleNormal(tnr2, n);
	    }
	}
    }
    }
  */

  calcedgedataanglesnew = 1;
  PopStatus();
}


int STLGeometry :: AddEdge(int ap1, int ap2)
{
  STLEdge e(ap1,ap2);
  e.SetLeftTrig(GetLeftTrig(ap1,ap2));
  e.SetRightTrig(GetRightTrig(ap1,ap2));
  return edges.Append(e);
}

void STLGeometry :: STLDoctorConfirmEdge()
{
  StoreEdgeData();
  if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT() && GetNodeOfSelTrig())
    {
      if (stldoctor.selectmode == 1)
	{
	  int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
	  int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
	  edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus (ED_CONFIRMED);
	}
      else if (stldoctor.selectmode == 3 || stldoctor.selectmode == 4)
	{
	  int i;
	  for (i = 1; i <= selectedmultiedge.Size(); i++)
	    {
	      int ap1 = selectedmultiedge.Get(i).i1;
	      int ap2 = selectedmultiedge.Get(i).i2;
	      edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus (ED_CONFIRMED);
	    }
	}
    }
}

void STLGeometry :: STLDoctorCandidateEdge()
{
  StoreEdgeData();
  if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT() && GetNodeOfSelTrig())
    {
      if (stldoctor.selectmode == 1)
	{
	  int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
	  int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
	  edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus (ED_CANDIDATE);
	}
      else if (stldoctor.selectmode == 3 || stldoctor.selectmode == 4)
	{
	  int i;
	  for (i = 1; i <= selectedmultiedge.Size(); i++)
	    {
	      int ap1 = selectedmultiedge.Get(i).i1;
	      int ap2 = selectedmultiedge.Get(i).i2;
	      edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus (ED_CANDIDATE);
	    }
	}
    }
}

void STLGeometry :: STLDoctorExcludeEdge()
{
  StoreEdgeData();
  if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT() && GetNodeOfSelTrig())
    {
      if (stldoctor.selectmode == 1)
	{
	  int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
	  int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
	  edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus(ED_EXCLUDED);
	}
      else if (stldoctor.selectmode == 3 || stldoctor.selectmode == 4)
	{
	  int i;
	  for (i = 1; i <= selectedmultiedge.Size(); i++)
	    {
	      int ap1 = selectedmultiedge.Get(i).i1;
	      int ap2 = selectedmultiedge.Get(i).i2;
	      edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus(ED_EXCLUDED);
	    }
	}
    }
}

void STLGeometry :: STLDoctorUndefinedEdge()
{
  StoreEdgeData();
  if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT() && GetNodeOfSelTrig())
    {
      if (stldoctor.selectmode == 1)
	{
	  int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
	  int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
	  edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus(ED_UNDEFINED);
	}
      else if (stldoctor.selectmode == 3 || stldoctor.selectmode == 4)
	{
	  int i;
	  for (i = 1; i <= selectedmultiedge.Size(); i++)
	    {
	      int ap1 = selectedmultiedge.Get(i).i1;
	      int ap2 = selectedmultiedge.Get(i).i2;
	      edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus(ED_UNDEFINED);
	    }
	}
    }
}

void STLGeometry :: STLDoctorSetAllUndefinedEdges()
{
  edgedata->ResetAll();
}

void STLGeometry :: STLDoctorEraseCandidateEdges()
{
  StoreEdgeData();
  edgedata->ChangeStatus(ED_CANDIDATE, ED_UNDEFINED);
}

void STLGeometry :: STLDoctorConfirmCandidateEdges()
{
  StoreEdgeData();
  edgedata->ChangeStatus(ED_CANDIDATE, ED_CONFIRMED);
}

void STLGeometry :: STLDoctorConfirmedToCandidateEdges()
{
  StoreEdgeData();
  edgedata->ChangeStatus(ED_CONFIRMED, ED_CANDIDATE);
}

void STLGeometry :: STLDoctorDirtyEdgesToCandidates()
{
  StoreEdgeData();
}

void STLGeometry :: STLDoctorLongLinesToCandidates()
{
  StoreEdgeData();
}

twoint STLGeometry :: GetNearestSelectedDefinedEdge()
{
  Point<3> pestimate = Center(GetTriangle(GetSelectTrig()).center,
  			     GetPoint(GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig())));
    //Point3d pestimate = GetTriangle(GetSelectTrig()).center;

  int i, j, en;
  Array<int> vic;
  GetVicinity(GetSelectTrig(),4,vic);
  

  twoint fedg;
  fedg.i1 = 0;
  fedg.i2 = 0;
  double mindist = 1E50;
  double dist;
  Point<3> p;

  for (i = 1; i <= vic.Size(); i++)
  {
    const STLTriangle& t = GetTriangle(vic.Get(i));
    for (j = 1; j <= 3; j++)
      {
	en = edgedata->GetEdgeNum(t.PNum(j),t.PNumMod(j+1));
	if (edgedata->Get(en).GetStatus() != ED_UNDEFINED)
	  {
	    p = pestimate;
	    dist = GetDistFromLine(GetPoint(t.PNum(j)),GetPoint(t.PNumMod(j+1)),p);
	    if (dist < mindist)
	      {
		mindist = dist;
		fedg.i1 = t.PNum(j);
		fedg.i2 = t.PNumMod(j+1);
	      }
	  }
      }
  }
  return fedg;
}
 
void STLGeometry :: BuildSelectedMultiEdge(twoint ep)
{
  if (edgedata->Size() == 0 || 
      !GetEPPSize()) 
    {
      return; 
    }

  selectedmultiedge.SetSize(0);
  int tenum = GetTopEdgeNum (ep.i1, ep.i2);

  if (edgedata->Get(tenum).GetStatus() == ED_UNDEFINED)
    {
      twoint epnew = GetNearestSelectedDefinedEdge();
      if (epnew.i1) 
	{
	  ep = epnew;
	  tenum = GetTopEdgeNum (ep.i1, ep.i2);
	}
    }

  selectedmultiedge.Append(twoint(ep));

  if (edgedata->Get(tenum).GetStatus() == ED_UNDEFINED)
    {
      return;
    }

  edgedata->BuildLineWithEdge(ep.i1,ep.i2,selectedmultiedge);
}

void STLGeometry :: BuildSelectedEdge(twoint ep)
{
  if (edgedata->Size() == 0 || 
      !GetEPPSize()) 
    {
      return; 
    }

  selectedmultiedge.SetSize(0);

  selectedmultiedge.Append(twoint(ep));
}

void STLGeometry :: BuildSelectedCluster(twoint ep)
{
  if (edgedata->Size() == 0 || 
      !GetEPPSize()) 
    {
      return; 
    }

  selectedmultiedge.SetSize(0);

  int tenum = GetTopEdgeNum (ep.i1, ep.i2);

  if (edgedata->Get(tenum).GetStatus() == ED_UNDEFINED)
    {
      twoint epnew = GetNearestSelectedDefinedEdge();
      if (epnew.i1) 
	{
	  ep = epnew;
	  tenum = GetTopEdgeNum (ep.i1, ep.i2);
	}
    }

  selectedmultiedge.Append(twoint(ep));

  if (edgedata->Get(tenum).GetStatus() == ED_UNDEFINED)
    {
      return;
    }

  edgedata->BuildClusterWithEdge(ep.i1,ep.i2,selectedmultiedge);
}

void STLGeometry :: ImportEdges()
{
  StoreEdgeData();

  PrintMessage(5, "import edges from file 'edges.ng'");
  ifstream fin("edges.ng");

  int ne;
  fin >> ne;

  Array<Point<3> > eps;

  int i;
  Point<3> p;
  for (i = 1; i <= 2*ne; i++)
    {
      fin >> p(0); 
      fin >> p(1); 
      fin >> p(2);
      eps.Append(p);
    }
  AddEdges(eps);
}

void STLGeometry :: AddEdges(const Array<Point<3> >& eps)
{
  int i;
  int ne = eps.Size()/2;
  
  Array<int> epsi;
  Box<3> bb = GetBoundingBox();
  bb.Increase(1);

  Point3dTree ptree (bb.PMin(), 
			 bb.PMax());
  Array<int> pintersect;

  double gtol = GetBoundingBox().Diam()/1.E10;
  Point<3> p;

  for (i = 1; i <= GetNP(); i++)
    {
      p = GetPoint(i);
      ptree.Insert (p, i);
    }
  
  int error = 0;
  for (i = 1; i <= 2*ne; i++)
    {
      p = eps.Get(i);
      Point3d pmin = p - Vec3d (gtol, gtol, gtol);
      Point3d pmax = p + Vec3d (gtol, gtol, gtol);
	  
      ptree.GetIntersecting (pmin, pmax, pintersect);
      if (pintersect.Size() > 1)
	{
	  PrintError("Found too much points in epsilon-dist");
	  error = 1;
	}
      else if (pintersect.Size() == 0)
	{
	  error = 1;
	  PrintError("edgepoint does not exist!");
	  PrintMessage(5,"p=",Point3d(eps.Get(i)));
	}
      else
	{
	  epsi.Append(pintersect.Get(1));
	}
    }

  if (error) return;

  int en;
  for (i = 1; i <= ne; i++)
    {
      if (epsi.Get(2*i-1) == epsi.Get(2*i)) {PrintError("Edge with zero length!");}
      else 
	{
	  en = edgedata->GetEdgeNum(epsi.Get(2*i-1),epsi.Get(2*i));
	  edgedata->Elem(en).SetStatus (ED_CONFIRMED);
	}
    }

}



void STLGeometry :: ImportExternalEdges(const char * filename)
{
  //AVL edges!!!!!!

  ifstream inf (filename);
  char ch;
  //int cnt = 0;
  int records, units, i, j;
  PrintFnStart("Import edges from ",filename);
  
  const int flen=30;
  char filter[flen+1];
  filter[flen] = 0;
  char buf[20];

  Array<Point3d> importpoints;
  Array<int> importlines;
  Array<int> importpnums;

  while (inf.good())
    {
      inf.get(ch);
      //      (*testout) << cnt << ": " << ch << endl;
      
      for (i = 0; i < flen; i++)
	filter[i] = filter[i+1];
      filter[flen-1] = ch;
      //      (*testout) << filter << endl;

      if (strcmp (filter+flen-7, "RECORDS") == 0)
	{
	  inf.get(ch);  // '='
	  inf >> records;
	}
      if (strcmp (filter+flen-5, "UNITS") == 0)
	{
	  inf.get(ch);  // '='
	  inf >> units;
	}

      if (strcmp (filter+flen-17, "EDGE NODE NUMBERS") == 0)
	{
	  int nodenr;
	  importlines.SetSize (units);
	  for (i = 1; i <= units; i++)
	    {
	      inf >> nodenr;
	      importlines.Elem(i) = nodenr;
	      //	      (*testout) << nodenr << endl;
	    }
	}

      if (strcmp (filter+flen-23, "EDGE POINT COORD IN DIR") == 0)
	{
	  int coord;

	  inf >> coord;
	  
	  importpoints.SetSize (units);

	  inf >> ch;
	  inf.putback (ch);

	  for (i = 1; i <= units; i++)
	    {
	      for (j = 0; j < 12; j++)
		inf.get (buf[j]);
	      buf[12] = 0;

	      importpoints.Elem(i).X(coord) = 1000 * atof (buf);
	    }
	}
    }

  /*
  (*testout) << "lines: " << endl;
  for (i = 1; i <= importlines.Size(); i++)
    (*testout) << importlines.Get(i) << endl;
  (*testout) << "points: " << endl;
  for (i = 1; i <= importpoints.Size(); i++)
    (*testout) << importpoints.Get(i) << endl;
  */



  importpnums.SetSize (importpoints.Size());
  

  Box3d bb (GetBoundingBox().PMin() + Vec3d (-1,-1,-1),
	    GetBoundingBox().PMax() + Vec3d (1, 1, 1));

  Point3dTree ptree (bb.PMin(), 
			 bb.PMax());


  PrintMessage(7,"stl - bb: ",bb.PMin(), " - ", bb.PMax());
  
  Box3d ebb;
  ebb.SetPoint (importpoints.Get(1));
  for (i = 1; i <= importpoints.Size(); i++)
    ebb.AddPoint (importpoints.Get(i));
  PrintMessage(7,"edgep - bb: ", ebb.PMin(), " - ", ebb.PMax());

  Array<int> pintersect;

  double gtol = GetBoundingBox().Diam()/1.E6;

  for (i = 1; i <= GetNP(); i++)
    {
      Point3d p = GetPoint(i);
      //      (*testout) << "stlpt: " << p << endl;
      ptree.Insert (p, i);
    }
  

  for (i = 1; i <= importpoints.Size(); i++)
    {
      Point3d p = importpoints.Get(i);
      Point3d pmin = p - Vec3d (gtol, gtol, gtol);
      Point3d pmax = p + Vec3d (gtol, gtol, gtol);
	  
      ptree.GetIntersecting (pmin, pmax, pintersect);
      if (pintersect.Size() > 1)
	{
	  importpnums.Elem(i) = 0;
	  PrintError("Found too many points in epsilon-dist");
	}
      else if (pintersect.Size() == 0)
	{
	  importpnums.Elem(i) = 0;
	  PrintError("Edgepoint does not exist!");
	}
      else
	{
	  importpnums.Elem(i) = pintersect.Get(1);
	}
    }

  //  if (!error) 
    {
      PrintMessage(7,"found all edge points in stl file");


      StoreEdgeData();

      int oldp = 0;

      for (i = 1; i <= importlines.Size(); i++)
	{
	  int newp = importlines.Get(i);
	  if (!importpnums.Get(abs(newp)))
	    newp = 0;

	  if (oldp && newp)
	    {
	      int en = edgedata->GetEdgeNum(importpnums.Get(oldp), 
					   importpnums.Get(abs(newp)));
	      edgedata->Elem(en).SetStatus (ED_CONFIRMED);
	    }
	  
	  if (newp < 0)
	    oldp = 0;
	  else
	    oldp = newp;
	}
    }


}



void STLGeometry :: ExportEdges()
{
  PrintFnStart("Save edges to file 'edges.ng'");

  ofstream fout("edges.ng");
  fout.precision(16);

  int n = edgedata->GetNConfEdges();
  
  fout << n << endl;

  int i;
  for (i = 1; i <= edgedata->Size(); i++)
    {
      if (edgedata->Get(i).GetStatus() == ED_CONFIRMED)
	{
	  const STLTopEdge & e = edgedata->Get(i);
	  fout << GetPoint(e.PNum(1))(0) << " " << GetPoint(e.PNum(1))(1) << " " << GetPoint(e.PNum(1))(2) << endl;
	  fout << GetPoint(e.PNum(2))(0) << " " << GetPoint(e.PNum(2))(1) << " " << GetPoint(e.PNum(2))(2) << endl;
	}
    }

}

void STLGeometry :: LoadEdgeData(const char* file)
{
  StoreEdgeData();

  PrintFnStart("Load edges from file '", file, "'");
  ifstream fin(file);

  edgedata->Read(fin);

  //  calcedgedataanglesnew = 1;
}

void STLGeometry :: SaveEdgeData(const char* file)
{
  PrintFnStart("save edges to file '", file, "'");
  ofstream fout(file);

  edgedata->Write(fout);
}







/*
void STLGeometry :: SaveExternalEdges()
{
  ofstream fout("externaledgesp3.ng");
  fout.precision(16);

  int n = NOExternalEdges();
  fout << n << endl;

  int i;
  for (i = 1; i <= n; i++)
    {
      twoint e = GetExternalEdge(i);
      fout << GetPoint(e.i1)(0) << " " << GetPoint(e.i1)(1) << " " << GetPoint(e.i1)(2) << endl;
      fout << GetPoint(e.i2)(0) << " " << GetPoint(e.i2)(1) << " " << GetPoint(e.i2)(2) << endl;
    }

}
*/
void STLGeometry :: StoreExternalEdges()
{
  storedexternaledges.SetSize(0);
  undoexternaledges = 1;
  int i;
  for (i = 1; i <= externaledges.Size(); i++)
    {
      storedexternaledges.Append(externaledges.Get(i));      
    }

}

void STLGeometry :: UndoExternalEdges()
{
  if (!undoexternaledges) 
    {
      PrintMessage(1, "undo not further possible!");
      return;
    }
  RestoreExternalEdges();
  undoexternaledges = 0;
}

void STLGeometry :: RestoreExternalEdges()
{
  externaledges.SetSize(0);
  int i;
  for (i = 1; i <= storedexternaledges.Size(); i++)
    {
      externaledges.Append(storedexternaledges.Get(i));      
    }

}


void STLGeometry :: AddExternalEdgeAtSelected()
{
  StoreExternalEdges();
  if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
    {
      int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
      int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
      if (!IsExternalEdge(ap1,ap2)) {AddExternalEdge(ap1,ap2);}
    }
}

void STLGeometry :: AddClosedLinesToExternalEdges()
{
  StoreExternalEdges();

  int i, j;
  for (i = 1; i <= GetNLines(); i++)
    {
      STLLine* l = GetLine(i);
      if (l->StartP() == l->EndP()) 
	{
	  for (j = 1; j < l->NP(); j++)
	    {
	      int ap1 = l->PNum(j);
	      int ap2 = l->PNum(j+1);

	      if (!IsExternalEdge(ap1,ap2)) {AddExternalEdge(ap1,ap2);}	      
	    }
	}
    }
}

void STLGeometry :: AddLongLinesToExternalEdges()
{
  StoreExternalEdges();

  double diamfact = stldoctor.dirtytrigfact;
  double diam = GetBoundingBox().Diam();

  int i, j;
  for (i = 1; i <= GetNLines(); i++)
    {
      STLLine* l = GetLine(i);
      if (l->GetLength(points) >= diamfact*diam) 
	{
	  for (j = 1; j < l->NP(); j++)
	    {
	      int ap1 = l->PNum(j);
	      int ap2 = l->PNum(j+1);

	      if (!IsExternalEdge(ap1,ap2)) {AddExternalEdge(ap1,ap2);}	      
	    }
	}
    }
}

void STLGeometry :: AddAllNotSingleLinesToExternalEdges()
{
  StoreExternalEdges();

  int i, j;
  for (i = 1; i <= GetNLines(); i++)
    {
      STLLine* l = GetLine(i);
      if (GetNEPP(l->StartP()) > 1 || GetNEPP(l->EndP()) > 1) 
	{
	  for (j = 1; j < l->NP(); j++)
	    {
	      int ap1 = l->PNum(j);
	      int ap2 = l->PNum(j+1);

	      if (!IsExternalEdge(ap1,ap2)) {AddExternalEdge(ap1,ap2);}	      
	    }
	}
    }
}

void STLGeometry :: DeleteDirtyExternalEdges()
{
  //delete single triangle edges and single edge-lines in clusters"
  StoreExternalEdges();

  int i, j;
  for (i = 1; i <= GetNLines(); i++)
    {
      STLLine* l = GetLine(i);
      if (l->NP() <= 3 || (l->StartP() == l->EndP() && l->NP() == 4))
	{
	  for (j = 1; j < l->NP(); j++)
	    {
	      int ap1 = l->PNum(j);
	      int ap2 = l->PNum(j+1);

	      if (IsExternalEdge(ap1,ap2)) {DeleteExternalEdge(ap1,ap2);}	      
	    }
	}
    }
}

void STLGeometry :: AddExternalEdgesFromGeomLine()
{
  StoreExternalEdges();
  if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
    {
      int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
      int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);

      if (IsEdge(ap1,ap2))
	{
	  int edgenum = IsEdgeNum(ap1,ap2);
	  if (!IsExternalEdge(ap1,ap2)) {AddExternalEdge(ap1,ap2);}
	  
	  int noend = 1;
	  int startp = ap1;
	  int laste = edgenum;
	  int np1, np2;
	  while (noend)
	    {
	      if (GetNEPP(startp) == 2)
		{
		  if (GetEdgePP(startp,1) != laste) {laste = GetEdgePP(startp,1);}
		  else {laste = GetEdgePP(startp,2);}
		  np1 = GetEdge(laste).PNum(1);
		  np2 = GetEdge(laste).PNum(2);
		  
		  if (!IsExternalEdge(np1, np2)) {AddExternalEdge(np1, np2);}
		  else {noend = 0;}
		  if (np1 != startp) {startp = np1;}
		  else {startp = np2;}
		}
	      else {noend = 0;}
	    }

	  startp = ap2;
	  laste = edgenum;
	  noend = 1;
	  while (noend)
	    {
	      if (GetNEPP(startp) == 2)
		{
		  if (GetEdgePP(startp,1) != laste) {laste = GetEdgePP(startp,1);}
		  else {laste = GetEdgePP(startp,2);}
		  np1 = GetEdge(laste).PNum(1);
		  np2 = GetEdge(laste).PNum(2);
		  
		  if (!IsExternalEdge(np1, np2)) {AddExternalEdge(np1, np2);}
		  else {noend = 0;}
		  if (np1 != startp) {startp = np1;}
		  else {startp = np2;}
		}
	      else {noend = 0;}
	    }
	  
	}

    }
  
}

void STLGeometry :: ClearEdges()
{
  edgesfound = 0;
  edges.SetSize(0);
  //edgedata->SetSize(0);
  // externaledges.SetSize(0);
  edgesperpoint.SetSize(0);
  undoexternaledges = 0;

}

void STLGeometry :: STLDoctorBuildEdges()
{
  //  if (!trigsconverted) {return;}
  ClearEdges();

  meshlines.SetSize(0);
  FindEdgesFromAngles();
}

void STLGeometry :: DeleteExternalEdgeAtSelected()
{
  StoreExternalEdges();
  if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
    {
      int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
      int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
      if (IsExternalEdge(ap1,ap2)) {DeleteExternalEdge(ap1,ap2);}
    }
}

void STLGeometry :: DeleteExternalEdgeInVicinity()
{
  StoreExternalEdges();
  if (!stldoctor.showvicinity || vicinity.Size() != GetNT()) {return;}

  int i, j, ap1, ap2;
  
  for (i = 1; i <= GetNT(); i++)
    {
      if (vicinity.Elem(i))
	{
	  for (j = 1; j <= 3; j++)
	    {
	      ap1 = GetTriangle(i).PNum(j);
	      ap2 = GetTriangle(i).PNumMod(j+1);

	      if (IsExternalEdge(ap1,ap2))
		{
		  DeleteExternalEdge(ap1,ap2);
		}
	    }
	}
    }
}

void STLGeometry :: BuildExternalEdgesFromEdges()
{
  StoreExternalEdges();

  if (GetNE() == 0) {PrintWarning("Edges possibly not generated!");}

  int i;
  externaledges.SetSize(0);

  for (i = 1; i <= GetNE(); i++)
    {
      STLEdge e = GetEdge(i);
      AddExternalEdge(e.PNum(1), e.PNum(2));
    }

}


void STLGeometry :: AddExternalEdge(int ap1, int ap2)
{
  externaledges.Append(twoint(ap1,ap2));
}

void STLGeometry :: DeleteExternalEdge(int ap1, int ap2)
{

  int i;
  int found = 0;
  for (i = 1; i <= NOExternalEdges(); i++)
    {
      if ((GetExternalEdge(i).i1 == ap1 && GetExternalEdge(i).i2 == ap2) ||
	  (GetExternalEdge(i).i1 == ap2 && GetExternalEdge(i).i2 == ap1)) {found = 1;};
      if (found && i < NOExternalEdges())
	{
	  externaledges.Elem(i) = externaledges.Get(i+1);
	}
    }
  if (!found) {PrintWarning("edge not found");}
  else
    {
      externaledges.SetSize(externaledges.Size()-1);
    }

}

int STLGeometry :: IsExternalEdge(int ap1, int ap2)
{
  int i;
  for (i = 1; i <= NOExternalEdges(); i++)
    {
      if ((GetExternalEdge(i).i1 == ap1 && GetExternalEdge(i).i2 == ap2) ||
	  (GetExternalEdge(i).i1 == ap2 && GetExternalEdge(i).i2 == ap1)) {return 1;};
    }
  return 0;
}

void STLGeometry :: DestroyDirtyTrigs()
{

  PrintFnStart("Destroy dirty triangles");
  PrintMessage(5,"original number of triangles=", GetNT());

  //destroy every triangle with other than 3 neighbours;
  int changed = 1;
  int i, j, k;
  while (changed)
    {
      changed = 0;
      Clear();

      for (i = 1; i <= GetNT(); i++)
	{
	  int dirty = NONeighbourTrigs(i) < 3;

	  for (j = 1; j <= 3; j++)
	    {
	      int pnum = GetTriangle(i).PNum(j);
	      /*
	      if (pnum == 1546)
		{
		// for (k = 1; k <=  NOTrigsPerPoint(pnum); k++)
		}
	      */
	      if (NOTrigsPerPoint(pnum) <= 2) 
		dirty = 1;
	    }
	  
	  int pi1 = GetTriangle(i).PNum(1);
	  int pi2 = GetTriangle(i).PNum(2);
	  int pi3 = GetTriangle(i).PNum(3);
	  if (pi1 == pi2 || pi1 == pi3 || pi2 == pi3)
	    {
	      PrintMessage(5,"triangle with Volume 0: ", i, "  nodes: ", pi1, ", ", pi2, ", ", pi3);
	      dirty = 1;
	    }

	  if (dirty)
	    {
	      for (k = i+1; k <= GetNT(); k++)
		{
		  trias.Elem(k-1) = trias.Get(k);
		  // readtrias: not longer permanent, JS
		  //		  readtrias.Elem(k-1) = readtrias.Get(k); 
		}
	      int size = GetNT();
	      trias.SetSize(size-1);
	      //	      readtrias.SetSize(size-1);
	      changed = 1;
	      break;
	    }
	}
    }  

  FindNeighbourTrigs();
  PrintMessage(5,"final number of triangles=", GetNT());
}

void STLGeometry :: CalcNormalsFromGeometry()
{
  int i;
  for (i = 1; i <= GetNT(); i++)
    {
      const STLTriangle & tr = GetTriangle(i);
      const Point3d& ap1 = GetPoint(tr.PNum(1));
      const Point3d& ap2 = GetPoint(tr.PNum(2));
      const Point3d& ap3 = GetPoint(tr.PNum(3));

      Vec3d normal = Cross (ap2-ap1, ap3-ap1);
      
      if (normal.Length() != 0)
	{
	  normal /= (normal.Length());		  
	}
      GetTriangle(i).SetNormal(normal);
    }
  PrintMessage(5,"Normals calculated from geometry!!!");

  calcedgedataanglesnew = 1;
}

void STLGeometry :: SetSelectTrig(int trig)
{
  stldoctor.selecttrig = trig;
}

int STLGeometry :: GetSelectTrig() const
{
  return stldoctor.selecttrig;
}

void STLGeometry :: SetNodeOfSelTrig(int n)
{
  stldoctor.nodeofseltrig = n;
}

int STLGeometry :: GetNodeOfSelTrig() const
{
  return stldoctor.nodeofseltrig;
}

void STLGeometry :: MoveSelectedPointToMiddle()
{
  if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
    {
      int p = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
      Point<3> pm(0.,0.,0.); //Middlevector;
      Point<3> p0(0.,0.,0.);
      PrintMessage(5,"original point=", Point3d(GetPoint(p)));

      int i;
      int cnt = 0;
      for (i = 1; i <= trigsperpoint.EntrySize(p); i++)
	{
	  const STLTriangle& tr = GetTriangle(trigsperpoint.Get(p,i));
	  int j;
	  for (j = 1; j <= 3; j++)
	    {
	      if (tr.PNum(j) != p)
		{
		  cnt++;
		  pm(0) += GetPoint(tr.PNum(j))(0);
		  pm(1) += GetPoint(tr.PNum(j))(1);
		  pm(2) += GetPoint(tr.PNum(j))(2);
		}
	    }
	}

      Point<3> origp = GetPoint(p);
      double fact = 0.2;

      SetPoint(p, p0 + fact*(1./(double)cnt)*(pm-p0)+(1.-fact)*(origp-p0));

      PrintMessage(5,"middle point=", Point3d (GetPoint(p)));
      
      PrintMessage(5,"moved point ", Point3d (p));

    }
}

void STLGeometry :: PrintSelectInfo()
{

  //int trig = GetSelectTrig();
  //int p = GetTriangle(trig).PNum(GetNodeOfSelTrig());
  
  PrintMessage(1,"touch triangle ", GetSelectTrig()
       , ", local node ", GetNodeOfSelTrig()
       , " (=", GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig()), ")");
  if (AtlasMade() && GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
    {
      PrintMessage(1,"           chartnum=",GetChartNr(GetSelectTrig()));
      /*      
      PointBetween(Center(Center(GetPoint(GetTriangle(270).PNum(1)),
				 GetPoint(GetTriangle(270).PNum(2))),
			  GetPoint(GetTriangle(270).PNum(3))),270,
		   Center(Center(GetPoint(GetTriangle(trig).PNum(1)),
				 GetPoint(GetTriangle(trig).PNum(2))),
			  GetPoint(GetTriangle(trig).PNum(3))),trig);
      */
      //PointBetween(Point3d(5.7818, 7.52768, 4.14879),260,Point3d(6.80292, 6.55392, 4.70184),233);
    }
}

void STLGeometry :: ShowSelectedTrigChartnum()
{
  int st = GetSelectTrig();

  if (st >= 1 && st <= GetNT() && AtlasMade())
    PrintMessage(1,"selected trig ", st, " has chartnumber ", GetChartNr(st));
}

void STLGeometry :: ShowSelectedTrigCoords()
{
  int st = GetSelectTrig();

  /*
  //testing!!!!
  Array<int> trigs;
  GetSortedTrianglesAroundPoint(GetTriangle(st).PNum(GetNodeOfSelTrig()),st,trigs);
  */

  if (st >= 1 && st <= GetNT())
    {
      PrintMessage(1, "coordinates of selected trig ", st, ":");
      PrintMessage(1, "   p1 = ", GetTriangle(st).PNum(1), " = ", 
		   Point3d (GetPoint(GetTriangle(st).PNum(1))));
      PrintMessage(1, "   p2 = ", GetTriangle(st).PNum(2), " = ", 
		   Point3d (GetPoint(GetTriangle(st).PNum(2))));
      PrintMessage(1, "   p3 = ", GetTriangle(st).PNum(3), " = ", 
		   Point3d (GetPoint(GetTriangle(st).PNum(3))));
    }
}

void STLGeometry :: LoadMarkedTrigs()
{
  PrintFnStart("load marked trigs from file 'markedtrigs.ng'");
  ifstream fin("markedtrigs.ng");

  int n;
  fin >> n;
  if (n != GetNT() || n == 0) {PrintError("Not a suitable marked-trig-file!"); return;}

  int i, m;
  for (i = 1; i <= n; i++)
    {
      fin >> m;
      SetMarkedTrig(i, m);      
    }

  fin >> n;
  if (n != 0) 
    {
 
      Point<3> ap1, ap2;
      for (i = 1; i <= n; i++)
	{
	  fin >> ap1(0); fin >> ap1(1); fin >> ap1(2);
	  fin >> ap2(0); fin >> ap2(1); fin >> ap2(2);
	  AddMarkedSeg(ap1,ap2);      
	}
    }
}

void STLGeometry :: SaveMarkedTrigs()
{
  PrintFnStart("save marked trigs to file 'markedtrigs.ng'");
  ofstream fout("markedtrigs.ng");

  int n = GetNT();
  fout << n << endl;

  int i;
  for (i = 1; i <= n; i++)
    {
      fout << IsMarkedTrig(i) << "\n";
    }

  n = GetNMarkedSegs();
  fout << n << endl;

  Point<3> ap1,ap2;
  for (i = 1; i <= n; i++)
    {
      GetMarkedSeg(i,ap1,ap2);
      fout << ap1(0) << " " << ap1(1) << " " << ap1(2) << "  ";
      fout << ap2(0) << " " << ap2(1) << " " << ap2(2) << " " << "\n";
    }

}

void STLGeometry :: NeighbourAnglesOfSelectedTrig()
{
  int st = GetSelectTrig();

  if (st >= 1 && st <= GetNT())
    {
      int i;
      PrintMessage(1,"Angle to triangle ", st, ":");
      for (i = 1; i <= NONeighbourTrigs(st); i++)
	{
	  PrintMessage(1,"   triangle ", NeighbourTrig(st,i), ": angle = ", 
		       180./M_PI*GetAngle(st, NeighbourTrig(st,i)), "°",
		       ", calculated = ", 180./M_PI*Angle(GetTriangle(st).GeomNormal(points), 
							  GetTriangle(NeighbourTrig(st,i)).GeomNormal(points)), "°");
	}
    }
}

void STLGeometry :: GetVicinity(int starttrig, int size, Array<int>& vic)
{
  if (starttrig == 0 || starttrig > GetNT()) {return;} 

  Array<int> vicarray;
  vicarray.SetSize(GetNT());

  int i;
  for (i = 1; i <= vicarray.Size(); i++)
    {
      vicarray.Elem(i) = 0;
    }
 
  vicarray.Elem(starttrig) = 1;
  
  int j = 0,k;

  Array <int> list1;
  list1.SetSize(0);
  Array <int> list2;
  list2.SetSize(0);
  list1.Append(starttrig);

  while (j < size)
    {
      j++;
      for (i = 1; i <= list1.Size(); i++)
	{
	  for (k = 1; k <= NONeighbourTrigs(i); k++)
	    {
	      int nbtrig = NeighbourTrig(list1.Get(i),k);
	      if (nbtrig && vicarray.Get(nbtrig) == 0)
		{
		  list2.Append(nbtrig);
		  vicarray.Elem(nbtrig) = 1;
		}
	    }
	}
      list1.SetSize(0);
      for (i = 1; i <= list2.Size(); i++)
	{
	  list1.Append(list2.Get(i));
	}
      list2.SetSize(0);
    }

  vic.SetSize(0);
  for (i = 1; i <= vicarray.Size(); i++)
    {
      if (vicarray.Get(i)) {vic.Append(i);}
    }
}

void STLGeometry :: CalcVicinity(int starttrig)
{
  if (starttrig == 0 || starttrig > GetNT()) {return;} 

  vicinity.SetSize(GetNT());

  if (!stldoctor.showvicinity) {return;}

  int i;
  for (i = 1; i <= vicinity.Size(); i++)
    {
      vicinity.Elem(i) = 0;
    }
 
  vicinity.Elem(starttrig) = 1;
  
  int j = 0,k;

  Array <int> list1;
  list1.SetSize(0);
  Array <int> list2;
  list2.SetSize(0);
  list1.Append(starttrig);

  //  int cnt = 1;
  while (j < stldoctor.vicinity)
    {
      j++;
      for (i = 1; i <= list1.Size(); i++)
	{
	  for (k = 1; k <= NONeighbourTrigs(i); k++)
	    {
	      int nbtrig = NeighbourTrig(list1.Get(i),k);
	      if (nbtrig && vicinity.Get(nbtrig) == 0)
		{
		  list2.Append(nbtrig);
		  vicinity.Elem(nbtrig) = 1;
		  //cnt++;
		}
	    }
	}
      list1.SetSize(0);
      for (i = 1; i <= list2.Size(); i++)
	{
	  list1.Append(list2.Get(i));
	}
      list2.SetSize(0);
    }

}

int STLGeometry :: Vicinity(int trig) const 
{
  if (trig <= vicinity.Size() && trig >=1)
    {
      return vicinity.Get(trig);
    }
  else {PrintSysError("In STLGeometry::Vicinity");}
  return 0;
}

void STLGeometry :: InitMarkedTrigs()
{
  markedtrigs.SetSize(GetNT());
  int i;
  for (i = 1; i <= GetNT(); i++)
    {
      SetMarkedTrig(i, 0);
    }
}

void STLGeometry :: MarkDirtyTrigs()
{
  PrintFnStart("mark dirty trigs");
  int i,j;

  markedtrigs.SetSize(GetNT());

  for (i = 1; i <= GetNT(); i++)
    {
      SetMarkedTrig(i, 0);
    }

  int found;
  double dirtyangle = stlparam.yangle/2./180.*M_PI;
  int cnt = 0;
  for (i = 1; i <= GetNT(); i++)
    {
      found = 0;
      for (j = 1; j <= NONeighbourTrigs(i); j++)
	{
	  if (GetAngle(i, NeighbourTrig(i,j)) > dirtyangle)
	    {
	      found++;
	    }
	}
      if (found && GetTriangle(i).MinHeight(points) < 
	  stldoctor.dirtytrigfact*GetTriangle(i).MaxLength(points))
	{
	  SetMarkedTrig(i, 1); cnt++;
	}
      /*
      else if (found == 3)
	{
	  SetMarkedTrig(i, 1); cnt++;	  
	}
      */
    }

  PrintMessage(1, "marked ", cnt, " dirty trigs");
}


void STLGeometry :: MarkTopErrorTrigs()
{
  int cnt = 0;
  markedtrigs.SetSize(GetNT());
  for (int i = 1; i <= GetNT(); i++)
    {
      const STLTriangle & trig = GetTriangle(i);

      SetMarkedTrig(i, trig.flags.toperror);
      if (trig.flags.toperror) cnt++;
    }
  PrintMessage(1,"marked ", cnt, " inconsistent triangles");
}



double STLGeometry :: CalcTrigBadness(int i)
{
  int j;
  double maxbadness = 0;
  int ap1, ap2;
  for (j = 1; j <= NONeighbourTrigs(i); j++)
    {
      GetTriangle(i).GetNeighbourPoints(GetTriangle(NeighbourTrig(i,j)), ap1, ap2);
      
      if (!IsEdge(ap1,ap2) && GetGeomAngle(i, NeighbourTrig(i,j)) > maxbadness)
	{
	  maxbadness = GetGeomAngle(i, NeighbourTrig(i,j));
	}
    }
  return maxbadness;

}

void STLGeometry :: GeomSmoothRevertedTrigs()
{
  //double revertedangle = stldoctor.smoothangle/180.*M_PI;
  double fact = stldoctor.dirtytrigfact;

  MarkRevertedTrigs();

  int i, j, k, l, p;

  for (i = 1; i <= GetNT(); i++)
    {
      if (IsMarkedTrig(i)) 
	{
	  for (j = 1; j <= 3; j++)
	    {
	      double origbadness = CalcTrigBadness(i);

	      p = GetTriangle(i).PNum(j);
	      Point<3> pm(0.,0.,0.); //Middlevector;
	      Point<3> p0(0.,0.,0.);

	      int cnt = 0;

	      for (k = 1; k <= trigsperpoint.EntrySize(p); k++)
		{
		  const STLTriangle& tr = GetTriangle(trigsperpoint.Get(p,k));
		  for (l = 1; l <= 3; l++)
		    {
		      if (tr.PNum(l) != p)
			{
			  cnt++;
			  pm(0) += GetPoint(tr.PNum(l))(0);
			  pm(1) += GetPoint(tr.PNum(l))(1);
			  pm(2) += GetPoint(tr.PNum(l))(2);
			}
		    }
		}
	      Point3d origp = GetPoint(p);
	      Point3d newp = p0 + fact*(1./(double)cnt)*(pm-p0)+(1.-fact)*(origp-p0);

	      SetPoint(p, newp);

	      if (CalcTrigBadness(i) > 0.9*origbadness) {SetPoint(p,origp); PrintDot('f');}
	      else {PrintDot('s');}
	    }
	}
    }
  MarkRevertedTrigs();
}

void STLGeometry :: MarkRevertedTrigs()
{
  int i,j;
  if (edgesperpoint.Size() != GetNP()) {BuildEdges();}

  PrintFnStart("mark reverted trigs");

  InitMarkedTrigs();

  int found;
  double revertedangle = stldoctor.smoothangle/180.*M_PI;

  int cnt = 0;
  int ap1, ap2;
  for (i = 1; i <= GetNT(); i++)
    {
      found = 0;
      for (j = 1; j <= NONeighbourTrigs(i); j++)
	{
	  GetTriangle(i).GetNeighbourPoints(GetTriangle(NeighbourTrig(i,j)), ap1, ap2);

	  if (!IsEdge(ap1,ap2))
	    {
              if (GetGeomAngle(i, NeighbourTrig(i,j)) > revertedangle)
		{
		  found = 1;
		  break;
		}
	    }
	}
      
      if (found)
	{
	  SetMarkedTrig(i, 1); cnt++;
	}
      
    }

  PrintMessage(5, "found ", cnt, " reverted trigs");


}

void STLGeometry :: SmoothDirtyTrigs()
{
  PrintFnStart("smooth dirty trigs");

  MarkDirtyTrigs();

  int i,j;
  int changed = 1;
  int ap1, ap2;
  
  while (changed)
    {
      changed = 0;
      for (i = 1; i <= GetNT(); i++)
	{
	  if (IsMarkedTrig(i))
	    {
	      int foundtrig = 0;
	      double maxlen = 0;
	      // JS: darf normalvector nicht ueber kurze Seite erben
	      maxlen = GetTriangle(i).MaxLength(GetPoints()) / 2.1; //JG: bei flachem dreieck auch kurze Seite

	      for (j = 1; j <= NONeighbourTrigs(i); j++)
		{
		  if (!IsMarkedTrig(NeighbourTrig(i,j)))
		    {
		      GetTriangle(i).GetNeighbourPoints(GetTriangle(NeighbourTrig(i,j)),ap1,ap2);
		      if (Dist(GetPoint(ap1),GetPoint(ap2)) >= maxlen)
			{
			  foundtrig = NeighbourTrig(i,j);
			  maxlen = Dist(GetPoint(ap1),GetPoint(ap2));
			}
		    }
		}
	      if (foundtrig)
		{
		  GetTriangle(i).SetNormal(GetTriangle(foundtrig).Normal());
		  changed = 1;
		  SetMarkedTrig(i,0);
		}
	    }
	}
    }

  calcedgedataanglesnew = 1;


  MarkDirtyTrigs();

  int cnt = 0;
  for (i = 1; i <= GetNT(); i++)
    {
      if (IsMarkedTrig(i)) {cnt++;}
    }

  PrintMessage(5,"NO marked dirty trigs=", cnt);

}

int STLGeometry :: IsMarkedTrig(int trig) const 
{
  if (trig <= markedtrigs.Size() && trig >=1)
    {
      return markedtrigs.Get(trig);
    }
  else {PrintSysError("In STLGeometry::IsMarkedTrig");}

  return 0;  
}

void STLGeometry :: SetMarkedTrig(int trig, int num)
{
  if (trig <= markedtrigs.Size() && trig >=1)
    {
      markedtrigs.Elem(trig) = num;
    }
  else {PrintSysError("In STLGeometry::SetMarkedTrig");}
}

void STLGeometry :: Clear()
{
  PrintFnStart("Clear");

  surfacemeshed = 0;
  surfaceoptimized = 0;
  volumemeshed = 0;

  selectedmultiedge.SetSize(0);
  meshlines.SetSize(0);
  // neighbourtrigs.SetSize(0);
  outerchartspertrig.SetSize(0);
  atlas.SetSize(0);
  ClearMarkedSegs();
  ClearSpiralPoints();
  ClearLineEndPoints();

  SetSelectTrig(0);
  SetNodeOfSelTrig(1);
  facecnt = 0;

  SetThreadPercent(100.);

  ClearEdges();
}

double STLGeometry :: Area()
{
  double ar = 0;
  int i;
  for (i = 1; i <= GetNT(); i++)
    {
      ar += GetTriangle(i).Area(points);
    }
  return ar;
}

double STLGeometry :: GetAngle(int t1, int t2)
{
  return Angle(GetTriangle(t1).Normal(),GetTriangle(t2).Normal());
}

double STLGeometry :: GetGeomAngle(int t1, int t2)
{
  Vec3d n1 = GetTriangle(t1).GeomNormal(points);
  Vec3d n2 = GetTriangle(t2).GeomNormal(points);
  return Angle(n1,n2);
}


void STLGeometry :: InitSTLGeometry(const Array<STLReadTriangle> & readtrias)
{
  PrintFnStart("Init STL Geometry");
  STLTopology::InitSTLGeometry(readtrias);

  int i, k;

  //const double geometry_tol_fact = 1E8; //distances lower than max_box_size/tol are ignored

  int np = GetNP();
  PrintMessage(5,"NO points= ", GetNP());
  normals.SetSize(GetNP());
  Array<int> normal_cnt(GetNP()); // counts number of added normals in a point

  for (i = 1; i <= np; i++)
    {
      normal_cnt.Elem(i) = 0;
      normals.Elem(i) = Vec3d (0,0,0);
    }

  for(i = 1; i <= GetNT(); i++)
    {
      //      STLReadTriangle t = GetReadTriangle(i);
      //      STLTriangle st;

      Vec<3> n = GetTriangle(i).Normal ();

      for (k = 1; k <= 3; k++)
	{
	  int pi = GetTriangle(i).PNum(k);
	  
	  normal_cnt.Elem(pi)++;
	  SetNormal(pi, GetNormal(pi) + n);
	}
    } 

  //normalize the normals
  for (i = 1; i <= GetNP(); i++)
    {
      SetNormal(i,1./(double)normal_cnt.Get(i)*GetNormal(i));
    }

  trigsconverted = 1;

  vicinity.SetSize(GetNT());
  markedtrigs.SetSize(GetNT());
  for (i = 1; i <= GetNT(); i++)
    {
      markedtrigs.Elem(i) = 0;
      vicinity.Elem(i) = 1;
    }

  ha_points.SetSize(GetNP());
  for (i = 1; i <= GetNP(); i++)
    ha_points.Elem(i) = 0;

  calcedgedataanglesnew = 0;
  edgedatastored = 0;
  edgedata->Clear();


  if (GetStatus() == STL_ERROR) return;

  CalcEdgeData();
  CalcEdgeDataAngles();

  ClearLineEndPoints();

  CheckGeometryOverlapping();
}

void STLGeometry :: TopologyChanged()
{
  calcedgedataanglesnew = 1;
}

int STLGeometry :: CheckGeometryOverlapping()
{
  int i, j, k;

  Box<3> geombox = GetBoundingBox();
  Point<3> pmin = geombox.PMin();
  Point<3> pmax = geombox.PMax();

  Box3dTree setree(pmin, pmax);
  Array<int> inters;

  int oltrigs = 0;
  markedtrigs.SetSize(GetNT());

  for (i = 1; i <= GetNT(); i++)
    SetMarkedTrig(i, 0);

  for (i = 1; i <= GetNT(); i++)
    {
      const STLTriangle & tri = GetTriangle(i);
      
      Point<3> tpmin = tri.box.PMin();
      Point<3> tpmax = tri.box.PMax();
      Vec<3> diag = tpmax - tpmin;

      tpmax = tpmax + 0.001 * diag;
      tpmin = tpmin - 0.001 * diag;

      setree.Insert (tpmin, tpmax, i);
    }

  for (i = 1; i <= GetNT(); i++)
    {
      const STLTriangle & tri = GetTriangle(i);
      
      Point<3> tpmin = tri.box.PMin();
      Point<3> tpmax = tri.box.PMax();

      setree.GetIntersecting (tpmin, tpmax, inters);

      for (j = 1; j <= inters.Size(); j++)
	{
	  const STLTriangle & tri2 = GetTriangle(inters.Get(j));

	  const Point<3> *trip1[3], *trip2[3];	
	  Point<3> hptri1[3], hptri2[3];
	  /*
	  for (k = 1; k <= 3; k++)
	    {
	      trip1[k-1] = &GetPoint (tri.PNum(k));
	      trip2[k-1] = &GetPoint (tri2.PNum(k));
	    }
	  */

	  for (k = 0; k < 3; k++)
	    {
	      hptri1[k] = GetPoint (tri[k]);
	      hptri2[k] = GetPoint (tri2[k]);
	      trip1[k] = &hptri1[k];
	      trip2[k] = &hptri2[k];
	    }

	  if (IntersectTriangleTriangle (&trip1[0], &trip2[0]))
	    {
	      oltrigs++;
	      PrintMessage(5,"Intersecting Triangles: trig ",i," with ",inters.Get(j),"!");
	      SetMarkedTrig(i, 1);
	      SetMarkedTrig(inters.Get(j), 1);
	    }
	}
    }

  PrintMessage(3,"Check Geometry Overlapping: overlapping triangles = ",oltrigs);
  return oltrigs;
}

/*
void STLGeometry :: InitSTLGeometry()
{
  STLTopology::InitSTLGeometry();

  int i, j, k;

  const double geometry_tol_fact = 1E8; //distances lower than max_box_size/tol are ignored


  trias.SetSize(0);
  points.SetSize(0);
  normals.SetSize(0);

  Array<int> normal_cnt; // counts number of added normals in a point

  Box3d bb (GetBoundingBox().PMin() + Vec3d (-1,-1,-1),
	    GetBoundingBox().PMax() + Vec3d (1, 1, 1));

  Point3dTree pointtree (bb.PMin(), 
			 bb.PMax());
  Array<int> pintersect;

  double gtol = GetBoundingBox().CalcDiam()/geometry_tol_fact;

  for(i = 1; i <= GetReadNT(); i++)
    {
      //if (i%500==499) {(*mycout) << (double)i/(double)GetReadNT()*100. << "%" << endl;}

      STLReadTriangle t = GetReadTriangle(i);
      STLTriangle st;
      Vec3d n = t.normal;

      for (k = 0; k < 3; k++)
	{
	  Point3d p = t.pts[k];

	  Point3d pmin = p - Vec3d (gtol, gtol, gtol);
	  Point3d pmax = p + Vec3d (gtol, gtol, gtol);
	  
	  pointtree.GetIntersecting (pmin, pmax, pintersect);
	  
	  if (pintersect.Size() > 1)
	    (*mycout) << "found too much  " << char(7) << endl;
	  int foundpos = 0;
	  if (pintersect.Size())
	    foundpos = pintersect.Get(1);

	  if (foundpos) 
	    {
	      normal_cnt[foundpos]++;
	      SetNormal(foundpos,GetNormal(foundpos)+n);
	      //	      (*testout) << "found p " << p << endl;
	    }
	  else
	    {
	      foundpos = AddPoint(p);
	      AddNormal(n);
	      normal_cnt.Append(1);

	      pointtree.Insert (p, foundpos);
	    }
	  //(*mycout) << "foundpos=" << foundpos << endl;
	  st.pts[k] = foundpos;
	}

      if ( (st.pts[0] == st.pts[1]) || 
	   (st.pts[0] == st.pts[2]) || 
	   (st.pts[1] == st.pts[2]) )
	{
	  (*mycout) << "ERROR: STL Triangle degenerated" << endl;
	}
      else
	{
	  // do not add ? js
	  AddTriangle(st);
	}
      //(*mycout) << "TRIG" << i << " = " << st << endl;
      
    } 
  //normal the normals
  for (i = 1; i <= GetNP(); i++)
    {
      SetNormal(i,1./(double)normal_cnt[i]*GetNormal(i));
    }

  trigsconverted = 1;

  vicinity.SetSize(GetNT());
  markedtrigs.SetSize(GetNT());
  for (i = 1; i <= GetNT(); i++)
    {
      markedtrigs.Elem(i) = 0;
      vicinity.Elem(i) = 1;
    }

  ha_points.SetSize(GetNP());
  for (i = 1; i <= GetNP(); i++)
    ha_points.Elem(i) = 0;

  calcedgedataanglesnew = 0;
  edgedatastored = 0;
  edgedata->Clear();

  CalcEdgeData();
  CalcEdgeDataAngles();

  ClearLineEndPoints();

  (*mycout) << "done" << endl;
}
*/



void STLGeometry :: SetLineEndPoint(int pn) 
{
  if (pn <1 || pn > lineendpoints.Size()) {PrintSysError("Illegal pnum in SetLineEndPoint!!!"); return; }
  lineendpoints.Elem(pn) = 1;
}

int STLGeometry :: IsLineEndPoint(int pn) 
{
  //  return 0;
  if (pn <1 || pn > lineendpoints.Size()) 
    {PrintSysError("Illegal pnum in IsLineEndPoint!!!"); return 0;}
  return lineendpoints.Get(pn);
}

void STLGeometry :: ClearLineEndPoints()
{
  lineendpoints.SetSize(GetNP());
  int i;
  for (i = 1; i <= GetNP(); i++)
    {
      lineendpoints.Elem(i) = 0;
    }
}

int STLGeometry :: IsEdge(int ap1, int ap2)
{
  int i,j;
  for (i = 1; i <= GetNEPP(ap1); i++)
    {
      for (j = 1; j <= GetNEPP(ap2); j++)
	{
	  if (GetEdgePP(ap1,i) == GetEdgePP(ap2,j)) {return 1;}
	}
    }
  return 0;
}

int STLGeometry :: IsEdgeNum(int ap1, int ap2)
{
  int i,j;
  for (i = 1; i <= GetNEPP(ap1); i++)
    {
      for (j = 1; j <= GetNEPP(ap2); j++)
	{
	  if (GetEdgePP(ap1,i) == GetEdgePP(ap2,j)) {return GetEdgePP(ap1,i);…