/contrib/Netgen/libsrc/stlgeom/stlgeom.cpp
C++ | 3506 lines | 2425 code | 564 blank | 517 comment | 608 complexity | c9efdd8c1eb92e185901b6f00e324a99 MD5 | raw file
Possible License(s): GPL-2.0, LGPL-2.1, BSD-3-Clause
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);…
Large files files are truncated, but you can click here to view the full file