mirror of
https://github.com/NGSolve/netgen.git
synced 2025-01-01 00:30:34 +05:00
de7ffc5906
Current initialization of the global geometryregister suffers from a classic 'initialization order fiasco'. Depending on the order the compilation units are loaded/linked, the initialization of the global geometryregisterarray is not guaranteed to happen (and indeed often does not happen) before it is used. This leads to entries being appended before it's initialized (usually 'suceeding, but potentially causing memory corruption if the segment at that point isn't zeroed), initialization then happening halfway through (wiping the initial entries) and then the last entries being the only ones that show up. The net effect is either a crash at startup, or several geometry types seeming to be missing. Eg, step files will oad, but STL files are just ignored. The bug is actively observed on, eg, Linux. This patch implements a simple 'initialize at first access' convention for the array, eliminating the ordering problem. I've not reviewed the rest of the source for other potential examples of the fiasco pattern; this fixes only the geometryregister, since that was actively biting.
3760 lines
84 KiB
C++
3760 lines
84 KiB
C++
#include <meshing.hpp>
|
|
#include <core/register_archive.hpp>
|
|
|
|
#include "stlgeom.hpp"
|
|
|
|
namespace netgen
|
|
{
|
|
|
|
//globalen searchtree fuer gesamte geometry aktivieren
|
|
int geomsearchtreeon = 0;
|
|
|
|
int usechartnormal = 1;
|
|
|
|
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
|
|
void STLMeshing (STLGeometry & geom,
|
|
Mesh & mesh,
|
|
const MeshingParameters& mparam,
|
|
const STLParameters& stlpar)
|
|
{
|
|
geom.Clear();
|
|
geom.BuildEdges(stlpar);
|
|
geom.MakeAtlas(mesh, mparam, stlpar);
|
|
if (multithread.terminate) { return; }
|
|
geom.CalcFaceNums();
|
|
geom.AddFaceEdges();
|
|
geom.LinkEdges(stlpar);
|
|
|
|
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 = make_unique<STLEdgeDataList>(*this);
|
|
externaledges.SetSize(0);
|
|
Clear();
|
|
meshchart = 0; // initialize all ?? JS
|
|
|
|
if (geomsearchtreeon)
|
|
searchtree = new BoxTree<3> (GetBoundingBox().PMin() - Vec3d(1,1,1),
|
|
GetBoundingBox().PMax() + Vec3d(1,1,1));
|
|
else
|
|
searchtree = NULL;
|
|
|
|
status = STL_GOOD;
|
|
statustext = "Good Geometry";
|
|
smoothedges = NULL;
|
|
area = -1;
|
|
}
|
|
|
|
STLGeometry :: ~STLGeometry()
|
|
{
|
|
// for (auto p : atlas) delete p;
|
|
// delete edgedata;
|
|
}
|
|
|
|
void STLGeometry :: Save (const filesystem::path & filename) const
|
|
{
|
|
string ext = ToLower(filename.extension());
|
|
|
|
if (ext == ".stl")
|
|
{
|
|
STLTopology::Save (filename);
|
|
return;
|
|
}
|
|
else if (ext == ".stlb")
|
|
{
|
|
SaveBinary (filename,"Binary STL Geometry");
|
|
return;
|
|
}
|
|
else if (ext == ".stle")
|
|
{
|
|
SaveSTLE (filename);
|
|
return;
|
|
}
|
|
|
|
throw Exception ("Unknown target format: " + filename.string());
|
|
}
|
|
|
|
|
|
|
|
DLL_HEADER extern STLParameters stlparam;
|
|
int STLGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
|
|
{
|
|
STLParameters stlpar = stlparam;
|
|
return STLMeshingDummy (this, mesh, mparam, stlpar);
|
|
}
|
|
|
|
Vec<3> STLGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const
|
|
{
|
|
if(!gi)
|
|
throw Exception("STLGeometry::GetNormal without PointGeomInfo called");
|
|
return GetChart(GetChartNr(gi->trignum)).GetNormal();
|
|
}
|
|
|
|
bool STLGeometry :: CalcPointGeomInfo(int /*surfind*/, PointGeomInfo& gi, const Point<3> & p3) const
|
|
{
|
|
Point<3> hp = p3;
|
|
SelectChartOfTriangle(gi.trignum);
|
|
|
|
gi.trignum = Project (hp);
|
|
|
|
if (gi.trignum) return true;
|
|
|
|
return false;
|
|
}
|
|
|
|
bool STLGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
|
|
{
|
|
static std::mutex mutex_project_whole_surface;
|
|
int meshchart = GetChartNr(gi.trignum);
|
|
const STLChart& chart = GetChart(meshchart);
|
|
int trignum = chart.ProjectNormal(p);
|
|
if(trignum==0)
|
|
{
|
|
// non-thread-safe implementation
|
|
std::lock_guard<std::mutex> guard(mutex_project_whole_surface);
|
|
PrintMessage(7,"project failed");
|
|
SelectChartOfTriangle (gi.trignum); // needed because ProjectOnWholeSurface uses meshchartnv (the normal vector of selected chart)
|
|
trignum = ProjectOnWholeSurface(p);
|
|
if(trignum==0)
|
|
{
|
|
PrintMessage(7, "project on whole surface failed");
|
|
return false;
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
PointGeomInfo STLGeometry :: ProjectPoint (INDEX surfind, Point<3> & p) const
|
|
{
|
|
throw Exception("ProjectPoint without PointGeomInfo not implemented");
|
|
}
|
|
|
|
void STLGeometry ::
|
|
PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
|
|
int surfi,
|
|
const PointGeomInfo & gi1,
|
|
const PointGeomInfo & gi2,
|
|
Point<3> & newp, PointGeomInfo & newgi) const
|
|
{
|
|
newp = p1+secpoint*(p2-p1);
|
|
|
|
/*
|
|
(*testout) << "surf-between: p1 = " << p1 << ", p2 = " << p2
|
|
<< ", gi = " << gi1 << " - " << gi2 << endl;
|
|
*/
|
|
|
|
if (gi1.trignum > 0)
|
|
{
|
|
// ((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
|
|
|
|
Point<3> np1 = newp;
|
|
Point<3> np2 = newp;
|
|
auto ngi1 = gi1;
|
|
auto ngi2 = gi2;
|
|
// SelectChartOfTriangle (gi1.trignum);
|
|
int tn1 = ProjectPointGI (surfi, np1, ngi1);
|
|
|
|
// SelectChartOfTriangle (gi2.trignum);
|
|
int tn2 = ProjectPointGI (surfi, np2, ngi2);
|
|
|
|
newgi.trignum = tn1; //urspruengliche version
|
|
newp = np1; //urspruengliche version
|
|
|
|
if (!newgi.trignum)
|
|
{ newgi.trignum = tn2; newp = np2; }
|
|
if (!newgi.trignum) newgi.trignum = gi1.trignum;
|
|
}
|
|
else
|
|
{
|
|
// (*testout) << "WARNING: PointBetween got geominfo = 0" << endl;
|
|
newp = p1+secpoint*(p2-p1);
|
|
newgi.trignum = 0;
|
|
}
|
|
}
|
|
|
|
void STLGeometry ::
|
|
PointBetweenEdge (const Point<3> & p1, const Point<3> & p2, double secpoint,
|
|
int surfi1, int surfi2,
|
|
const EdgePointGeomInfo & gi1,
|
|
const EdgePointGeomInfo & gi2,
|
|
Point<3> & newp, EdgePointGeomInfo & newgi) const
|
|
{
|
|
/*
|
|
(*testout) << "edge-between: p1 = " << p1 << ", p2 = " << p2
|
|
<< ", gi1,2 = " << gi1 << ", " << gi2 << endl;
|
|
*/
|
|
/*
|
|
newp = Center (p1, p2);
|
|
((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
|
|
newgi.trignum = geom.Project (newp);
|
|
*/
|
|
int hi;
|
|
newgi.dist = (1.0-secpoint) * gi1.dist + secpoint*gi2.dist;
|
|
newgi.edgenr = gi1.edgenr;
|
|
|
|
/*
|
|
(*testout) << "p1 = " << p1 << ", p2 = " << p2 << endl;
|
|
(*testout) << "refedge: " << gi1.edgenr
|
|
<< " d1 = " << gi1.dist << ", d2 = " << gi2.dist << endl;
|
|
*/
|
|
newp = GetLine (gi1.edgenr)->GetPointInDist (GetPoints(), newgi.dist, hi);
|
|
|
|
// (*testout) << "newp = " << newp << endl;
|
|
}
|
|
|
|
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(const STLParameters& stlparam)
|
|
{
|
|
|
|
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;
|
|
STLPointId 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(const STLParameters& stlparam)
|
|
{
|
|
multithread.terminate = 0;
|
|
|
|
// UseExternalEdges();
|
|
|
|
BuildEdges(stlparam);
|
|
|
|
|
|
DenseMatrix m(3), hm(3);
|
|
Vector rhs(3), sol(3), hv(3), hv2(3);
|
|
|
|
Vec<3> ri;
|
|
|
|
double wnb = stldoctor.smoothnormalsweight; // neighbour 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;
|
|
STLPointId 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:
|
|
|
|
NgArray<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())
|
|
{
|
|
|
|
NgArray<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));
|
|
edges.Append(e);
|
|
return edges.Size();
|
|
}
|
|
|
|
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;
|
|
NgArray<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;
|
|
|
|
NgArray<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 NgArray<Point<3> >& eps)
|
|
{
|
|
int i;
|
|
int ne = eps.Size()/2;
|
|
|
|
NgArray<int> epsi;
|
|
Box<3> bb = GetBoundingBox();
|
|
bb.Increase(1);
|
|
|
|
Point3dTree ptree (bb.PMin(),
|
|
bb.PMax());
|
|
NgArray<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];
|
|
|
|
NgArray<Point3d> importpoints;
|
|
NgArray<int> importlines;
|
|
NgArray<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());
|
|
|
|
NgArray<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 filesystem::path & filename)
|
|
{
|
|
StoreEdgeData();
|
|
|
|
PrintFnStart("Load edges from file '", filename, "'");
|
|
ifstream fin(filename);
|
|
|
|
edgedata->Read(fin);
|
|
|
|
// calcedgedataanglesnew = 1;
|
|
}
|
|
|
|
void STLGeometry :: SaveEdgeData(const filesystem::path & filename)
|
|
{
|
|
PrintFnStart("save edges to file '", filename, "'");
|
|
ofstream fout(filename);
|
|
|
|
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(const STLParameters& stlparam)
|
|
{
|
|
// if (!trigsconverted) {return;}
|
|
ClearEdges();
|
|
|
|
meshlines.SetSize(0);
|
|
FindEdgesFromAngles(stlparam);
|
|
}
|
|
|
|
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[k-1] = trias[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()
|
|
, " (=", int(GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig())), ")");
|
|
if (AtlasMade() && GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
|
|
{
|
|
PrintMessage(1," chartnum=", int(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 ", int(GetChartNr(st)));
|
|
}
|
|
|
|
void STLGeometry :: ShowSelectedTrigCoords()
|
|
{
|
|
int st = GetSelectTrig();
|
|
|
|
/*
|
|
//testing!!!!
|
|
NgArray<int> trigs;
|
|
GetSortedTrianglesAroundPoint(GetTriangle(st).PNum(GetNodeOfSelTrig()),st,trigs);
|
|
*/
|
|
|
|
if (st >= 1 && st <= GetNT())
|
|
{
|
|
PrintMessage(1, "coordinates of selected trig ", st, ":");
|
|
PrintMessage(1, " p1 = ", int(GetTriangle(st).PNum(1)), " = ",
|
|
Point3d (GetPoint(GetTriangle(st).PNum(1))));
|
|
PrintMessage(1, " p2 = ", int(GetTriangle(st).PNum(2)), " = ",
|
|
Point3d (GetPoint(GetTriangle(st).PNum(2))));
|
|
PrintMessage(1, " p3 = ", int(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 ", int(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, NgArray<int>& vic)
|
|
{
|
|
if (starttrig == 0 || starttrig > GetNT()) {return;}
|
|
|
|
NgArray<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;
|
|
|
|
NgArray <int> list1;
|
|
list1.SetSize(0);
|
|
NgArray <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;
|
|
|
|
NgArray <int> list1;
|
|
list1.SetSize(0);
|
|
NgArray <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(const STLParameters& stlparam)
|
|
{
|
|
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;
|
|
STLPointId 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(const STLParameters& stlparam)
|
|
{
|
|
//double revertedangle = stldoctor.smoothangle/180.*M_PI;
|
|
double fact = stldoctor.dirtytrigfact;
|
|
|
|
MarkRevertedTrigs(stlparam);
|
|
|
|
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(stlparam);
|
|
}
|
|
|
|
void STLGeometry :: MarkRevertedTrigs(const STLParameters& stlparam)
|
|
{
|
|
int i,j;
|
|
if (edgesperpoint.Size() != GetNP()) {BuildEdges(stlparam);}
|
|
|
|
PrintFnStart("mark reverted trigs");
|
|
|
|
InitMarkedTrigs();
|
|
|
|
int found;
|
|
double revertedangle = stldoctor.smoothangle/180.*M_PI;
|
|
|
|
int cnt = 0;
|
|
STLPointId 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(const STLParameters& stlparam)
|
|
{
|
|
PrintFnStart("smooth dirty trigs");
|
|
|
|
MarkDirtyTrigs(stlparam);
|
|
|
|
int i,j;
|
|
int changed = 1;
|
|
STLPointId 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(stlparam);
|
|
|
|
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()
|
|
{
|
|
if (area >= 0) return area;
|
|
area = 0;
|
|
for (int i = 1; i <= GetNT(); i++)
|
|
area += GetTriangle(i).Area(points);
|
|
return area;
|
|
}
|
|
|
|
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 NgArray<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());
|
|
NgArray<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()
|
|
{
|
|
PrintMessageCR(3,"Check overlapping geometry ...");
|
|
|
|
Box<3> geombox = GetBoundingBox();
|
|
Point<3> pmin = geombox.PMin();
|
|
Point<3> pmax = geombox.PMax();
|
|
|
|
BoxTree<3> setree(pmin, pmax);
|
|
|
|
int oltrigs = 0;
|
|
markedtrigs.SetSize(GetNT());
|
|
|
|
for (int i = 1; i <= GetNT(); i++)
|
|
SetMarkedTrig(i, 0);
|
|
|
|
for (int 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);
|
|
}
|
|
|
|
|
|
{
|
|
mutex inters_mutex;
|
|
|
|
ParallelFor( 1, GetNT()+1, [&] (int first, int next)
|
|
{
|
|
NgArray<int> inters;
|
|
for (int i=first; i<next; i++) {
|
|
const STLTriangle & tri = GetTriangle(i);
|
|
|
|
Point<3> tpmin = tri.box.PMin();
|
|
Point<3> tpmax = tri.box.PMax();
|
|
|
|
setree.GetIntersecting (tpmin, tpmax, inters);
|
|
|
|
for (int 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 (int 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]))
|
|
{
|
|
lock_guard<mutex> guard(inters_mutex);
|
|
{
|
|
oltrigs++;
|
|
PrintMessage(5,"Intersecting Triangles: trig ",i," with ",inters.Get(j),"!");
|
|
SetMarkedTrig(i, 1);
|
|
SetMarkedTrig(inters.Get(j), 1);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
});
|
|
}
|
|
PrintMessage(3,"Check overlapping geometry ... ", oltrigs, " triangles overlap");
|
|
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);
|
|
|
|
NgArray<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());
|
|
NgArray<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);}
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
|
|
void STLGeometry :: BuildEdges(const STLParameters& stlparam)
|
|
{
|
|
//PrintFnStart("build edges");
|
|
edges.SetSize(0);
|
|
meshlines.SetSize(0);
|
|
FindEdgesFromAngles(stlparam);
|
|
}
|
|
|
|
void STLGeometry :: UseExternalEdges()
|
|
{
|
|
for (int i = 1; i <= NOExternalEdges(); i++)
|
|
AddEdge(GetExternalEdge(i).i1,GetExternalEdge(i).i2);
|
|
//BuildEdgesPerPointy();
|
|
}
|
|
|
|
void STLGeometry :: UndoEdgeChange()
|
|
{
|
|
if (edgedatastored)
|
|
{
|
|
RestoreEdgeData();
|
|
}
|
|
else
|
|
{
|
|
PrintWarning("no edge undo possible");
|
|
}
|
|
}
|
|
|
|
|
|
void STLGeometry :: StoreEdgeData()
|
|
{
|
|
// edgedata_store = *edgedata;
|
|
|
|
edgedata->Store();
|
|
edgedatastored = 1;
|
|
|
|
// put stlgeom-edgedata to stltopology edgedata
|
|
/*
|
|
int i;
|
|
for (i = 1; i <= GetNTE(); i++)
|
|
{
|
|
const STLTopEdge & topedge = GetTopEdge (i);
|
|
int ednum = edgedata->GetEdgeNum (topedge.PNum(1),
|
|
topedge.PNum(2));
|
|
topedges.Elem(i).SetStatus (edgedata->Get (ednum).status);
|
|
}
|
|
*/
|
|
}
|
|
|
|
void STLGeometry :: RestoreEdgeData()
|
|
{
|
|
// *edgedata = edgedata_store;
|
|
edgedata->Restore();
|
|
edgedatastored=0;
|
|
}
|
|
|
|
|
|
void STLGeometry :: CalcEdgeData()
|
|
{
|
|
PushStatus("Calc Edge Data");
|
|
|
|
STLPointId np1, np2;
|
|
|
|
int ecnt = 0;
|
|
edgedata->SetSize(GetNT()/2*3);
|
|
|
|
for (int i = 1; i <= GetNT(); i++)
|
|
{
|
|
SetThreadPercent((double)i/(double)GetNT()*100.);
|
|
|
|
const STLTriangle & t1 = GetTriangle(i);
|
|
|
|
for (int j = 1; j <= NONeighbourTrigs(i); j++)
|
|
{
|
|
int nbti = NeighbourTrig(i,j);
|
|
if (nbti > i)
|
|
{
|
|
const STLTriangle & t2 = GetTriangle(nbti);
|
|
|
|
if (t1.IsNeighbourFrom(t2))
|
|
{
|
|
ecnt++; if (ecnt > edgedata->Size()) {PrintError("In Calc edge data, illegal geometry");}
|
|
|
|
t1.GetNeighbourPoints(t2,np1,np2);
|
|
|
|
/* ang = GetAngle(i,nbti);
|
|
if (ang < -M_PI) {ang += 2*M_PI;}*/
|
|
|
|
|
|
// edgedata->Add(STLEdgeData(0, np1, np2, i, nbti),ecnt);
|
|
edgedata->Elem(ecnt).SetStatus(ED_UNDEFINED);
|
|
|
|
// edgedata->Elem(ecnt).top = this;
|
|
// edgedata->Elem(ecnt).topedgenr = GetTopEdgeNum (np1, np2);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
//BuildEdgesPerPoint();
|
|
PopStatus();
|
|
}
|
|
|
|
void STLGeometry :: CalcEdgeDataAngles()
|
|
{
|
|
PrintMessageCR (5,"calc edge data angles ... ");
|
|
|
|
for (int i = 1; i <= GetNTE(); i++)
|
|
{
|
|
STLTopEdge & edge = GetTopEdge (i);
|
|
double cosang =
|
|
GetTriangle(edge.TrigNum(1)).Normal() *
|
|
GetTriangle(edge.TrigNum(2)).Normal();
|
|
edge.SetCosAngle (cosang);
|
|
}
|
|
|
|
for (int i = 1; i <= edgedata->Size(); i++)
|
|
{
|
|
/*
|
|
const STLEdgeData& e = edgedata->Get(i);
|
|
ang = GetAngle(e.lt,e.rt);
|
|
if (ang < -M_PI) {ang += 2*M_PI;}
|
|
edgedata->Elem(i).angle = fabs(ang);
|
|
*/
|
|
}
|
|
PrintMessage (5,"calc edge data angles ... done");
|
|
}
|
|
|
|
void STLGeometry :: FindEdgesFromAngles(const STLParameters& stlparam)
|
|
{
|
|
// PrintFnStart("find edges from angles");
|
|
|
|
double min_edge_angle = stlparam.yangle/180.*M_PI;
|
|
double cont_min_edge_angle = stlparam.contyangle/180.*M_PI;
|
|
|
|
double cos_min_edge_angle = cos (min_edge_angle);
|
|
double cos_cont_min_edge_angle = cos (cont_min_edge_angle);
|
|
|
|
if (calcedgedataanglesnew) {CalcEdgeDataAngles(); calcedgedataanglesnew = 0;}
|
|
|
|
for (int i = 1; i <= edgedata->Size(); i++)
|
|
{
|
|
STLTopEdge & sed = edgedata->Elem(i);
|
|
if (sed.GetStatus() == ED_CANDIDATE ||
|
|
sed.GetStatus() == ED_UNDEFINED)
|
|
{
|
|
if (sed.CosAngle() <= cos_min_edge_angle)
|
|
{
|
|
sed.SetStatus (ED_CANDIDATE);
|
|
}
|
|
else
|
|
{
|
|
sed.SetStatus(ED_UNDEFINED);
|
|
}
|
|
}
|
|
}
|
|
|
|
if (stlparam.contyangle < stlparam.yangle)
|
|
{
|
|
int changed = 1;
|
|
int its = 0;
|
|
while (changed && stlparam.contyangle < stlparam.yangle)
|
|
{
|
|
its++;
|
|
//(*mycout) << "." << flush;
|
|
changed = 0;
|
|
for (int i = 1; i <= edgedata->Size(); i++)
|
|
{
|
|
STLTopEdge & sed = edgedata->Elem(i);
|
|
if (sed.CosAngle() <= cos_cont_min_edge_angle
|
|
&& sed.GetStatus() == ED_UNDEFINED &&
|
|
(edgedata->GetNConfCandEPP(sed.PNum(1)) == 1 ||
|
|
edgedata->GetNConfCandEPP(sed.PNum(2)) == 1))
|
|
{
|
|
changed = 1;
|
|
sed.SetStatus (ED_CANDIDATE);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
int confcand = 0;
|
|
if (edgedata->GetNConfEdges() == 0)
|
|
{
|
|
confcand = 1;
|
|
}
|
|
|
|
for (int i = 1; i <= edgedata->Size(); i++)
|
|
{
|
|
STLTopEdge & sed = edgedata->Elem(i);
|
|
if (sed.GetStatus() == ED_CONFIRMED ||
|
|
(sed.GetStatus() == ED_CANDIDATE && confcand))
|
|
{
|
|
STLEdge se(sed.PNum(1),sed.PNum(2));
|
|
se.SetLeftTrig(sed.TrigNum(1));
|
|
se.SetRightTrig(sed.TrigNum(2));
|
|
AddEdge(se);
|
|
}
|
|
}
|
|
BuildEdgesPerPoint();
|
|
|
|
|
|
|
|
//(*mycout) << "its for continued angle = " << its << endl;
|
|
PrintMessage(5,"built ", GetNE(), " edges with yellow angle = ", stlparam.yangle, " degree");
|
|
|
|
}
|
|
|
|
/*
|
|
void STLGeometry :: FindEdgesFromAngles()
|
|
{
|
|
double yangle = stlparam.yangle;
|
|
char * savetask = multithread.task;
|
|
multithread.task = "find edges";
|
|
|
|
const double min_edge_angle = yangle/180.*M_PI;
|
|
|
|
int np1, np2;
|
|
double ang;
|
|
int i;
|
|
|
|
//(*mycout) << "area=" << Area() << endl;
|
|
|
|
for (i = 1; i <= GetNT(); i++)
|
|
{
|
|
multithread.percent = (double)i/(double)GetReadNT()*100.;
|
|
|
|
const STLTriangle & t1 = GetTriangle(i);
|
|
//NeighbourTrigs(nt,i);
|
|
|
|
for (int j = 1; j <= NONeighbourTrigs(i); j++)
|
|
{
|
|
int nbti = NeighbourTrig(i,j);
|
|
if (nbti > i)
|
|
{
|
|
const STLTriangle & t2 = GetTriangle(nbti);
|
|
|
|
if (t1.IsNeighbourFrom(t2))
|
|
{
|
|
ang = GetAngle(i,nbti);
|
|
if (ang < -M_PI*0.5) {ang += 2*M_PI;}
|
|
|
|
t1.GetNeighbourPoints(t2,np1,np2);
|
|
|
|
if (fabs(ang) >= min_edge_angle)
|
|
{
|
|
STLEdge se(np1,np2);
|
|
se.SetLeftTrig(i);
|
|
se.SetRightTrig(nbti);
|
|
AddEdge(se);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
(*mycout) << "added " << GetNE() << " edges" << endl;
|
|
|
|
//BuildEdgesPerPoint();
|
|
|
|
multithread.percent = 100.;
|
|
multithread.task = savetask;
|
|
|
|
}
|
|
*/
|
|
void STLGeometry :: BuildEdgesPerPoint()
|
|
{
|
|
//cout << "*** build edges per point" << endl;
|
|
edgesperpoint.SetSize(GetNP());
|
|
|
|
//add edges to points
|
|
for (int i = 1; i <= GetNE(); i++)
|
|
{
|
|
//(*mycout) << "EDGE " << GetEdge(i).PNum(1) << " - " << GetEdge(i).PNum(2) << endl;
|
|
for (int j = 1; j <= 2; j++)
|
|
{
|
|
AddEdgePP(GetEdge(i).PNum(j),i);
|
|
}
|
|
}
|
|
}
|
|
|
|
void STLGeometry :: AddFaceEdges()
|
|
{
|
|
PrintFnStart("Add starting edges for faces");
|
|
|
|
//für Kugel eine STLLine hinzufügen (Vorteil: verfeinerbar, unabhängig von Auflösung der Geometrie!!!):
|
|
//Grenze von 1. gefundener chart
|
|
|
|
NgArray<int> edgecnt;
|
|
NgArray<int> chartindex;
|
|
edgecnt.SetSize(GetNOFaces());
|
|
chartindex.SetSize(GetNOFaces());
|
|
|
|
for (int i = 1; i <= GetNOFaces(); i++)
|
|
{
|
|
edgecnt.Elem(i) = 0;
|
|
chartindex.Elem(i) = 0;
|
|
}
|
|
|
|
for (int i = 1; i <= GetNT(); i++)
|
|
{
|
|
int fn = GetTriangle(i).GetFaceNum();
|
|
if (!chartindex.Get(fn)) {chartindex.Elem(fn) = GetChartNr(i);}
|
|
for (int j = 1; j <= 3; j++)
|
|
{
|
|
edgecnt.Elem(fn) += GetNEPP(GetTriangle(i).PNum(j));
|
|
}
|
|
}
|
|
|
|
for (int i = 1; i <= GetNOFaces(); i++)
|
|
{
|
|
if (!edgecnt.Get(i)) {PrintMessage(5,"Face", i, " has no edge!");}
|
|
}
|
|
|
|
int changed = 0;
|
|
STLPointId ap1, ap2;
|
|
for (int i = 1; i <= GetNOFaces(); i++)
|
|
{
|
|
if (!edgecnt.Get(i))
|
|
{
|
|
const STLChart& c = GetChart(chartindex.Get(i));
|
|
// bool foundone = false;
|
|
int longest_ap1 = -1, longest_ap2 = -1;
|
|
double maxlen = -1;
|
|
for (int j = 1; j <= c.GetNChartT(); j++)
|
|
{
|
|
const STLTriangle& t1 = GetTriangle(c.GetChartTrig1(j));
|
|
for (int k = 1; k <= 3; k++)
|
|
{
|
|
int nt = NeighbourTrig(c.GetChartTrig1(j),k);
|
|
if (GetChartNr(nt) != chartindex.Get(i))
|
|
{
|
|
t1.GetNeighbourPoints(GetTriangle(nt),ap1,ap2);
|
|
// AddEdge(ap1,ap2);
|
|
double len = Dist(GetPoint(ap1), GetPoint(ap2));
|
|
if (len > maxlen)
|
|
{
|
|
maxlen = len;
|
|
longest_ap1 = ap1;
|
|
longest_ap2 = ap2;
|
|
}
|
|
changed = 1;
|
|
}
|
|
}
|
|
}
|
|
if (maxlen > 0)
|
|
AddEdge(longest_ap1,longest_ap2);
|
|
}
|
|
|
|
}
|
|
|
|
if (changed) BuildEdgesPerPoint();
|
|
|
|
}
|
|
|
|
void STLGeometry :: LinkEdges(const STLParameters& stlparam)
|
|
{
|
|
PushStatusF("Link Edges");
|
|
PrintMessage(5,"have now ", GetNE(), " edges with yellow angle = ", stlparam.yangle, " degree");
|
|
|
|
int i;
|
|
|
|
lines.SetSize(0);
|
|
int starte(0);
|
|
int edgecnt = 0;
|
|
int found;
|
|
int rev(0); //indicates, that edge is inserted reverse
|
|
|
|
//worked edges
|
|
NgArray<int> we(GetNE());
|
|
|
|
//setlineendpoints; wenn 180°, dann keine endpunkte
|
|
//nur punkte mit 2 edges kommen in frage, da bei mehr oder weniger punkten ohnehin ein meshpoint hinkommt
|
|
|
|
Vec3d v1,v2;
|
|
double cos_eca = cos(stlparam.edgecornerangle/180.*M_PI);
|
|
int ecnt = 0;
|
|
int lp1, lp2;
|
|
if (stlparam.edgecornerangle < 180)
|
|
{
|
|
for (i = 1; i <= GetNP(); i++)
|
|
{
|
|
if (GetNEPP(i) == 2)
|
|
{
|
|
if (GetEdge(GetEdgePP(i,1)).PNum(2) == GetEdge(GetEdgePP(i,2)).PNum(1) ||
|
|
GetEdge(GetEdgePP(i,1)).PNum(1) == GetEdge(GetEdgePP(i,2)).PNum(2))
|
|
{
|
|
lp1 = 1; lp2 = 2;
|
|
}
|
|
else
|
|
{
|
|
lp1 = 2; lp2 = 1;
|
|
}
|
|
|
|
v1 = Vec3d(GetPoint(GetEdge(GetEdgePP(i,1)).PNum(1)),
|
|
GetPoint(GetEdge(GetEdgePP(i,1)).PNum(2)));
|
|
v2 = Vec3d(GetPoint(GetEdge(GetEdgePP(i,2)).PNum(lp1)),
|
|
GetPoint(GetEdge(GetEdgePP(i,2)).PNum(lp2)));
|
|
if ((v1*v2)/sqrt(v1.Length2()*v2.Length2()) < cos_eca)
|
|
{
|
|
//(*testout) << "add edgepoint " << i << endl;
|
|
SetLineEndPoint(i);
|
|
ecnt++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
PrintMessage(5, "added ", ecnt, " mesh_points due to edge corner angle (",
|
|
stlparam.edgecornerangle, " degree)");
|
|
|
|
for (i = 1; i <= GetNE(); i++) {we.Elem(i) = 0;}
|
|
|
|
while(edgecnt < GetNE())
|
|
{
|
|
SetThreadPercent((double)edgecnt/(double)GetNE()*100.);
|
|
|
|
STLLine* line = new STLLine(this);
|
|
|
|
//find start edge
|
|
int j = 1;
|
|
found = 0;
|
|
//try second time, if only rings are left!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
int second = 0;
|
|
|
|
//find a starting edge at point with 1 or more than 2 edges or at lineendpoint
|
|
while (!found && j<=GetNE())
|
|
{
|
|
if (!we.Get(j))
|
|
{
|
|
if (GetNEPP(GetEdge(j).PNum(1)) != 2 || IsLineEndPoint(GetEdge(j).PNum(1)))
|
|
{
|
|
starte = j;
|
|
found = 1;
|
|
rev = 0;
|
|
}
|
|
else
|
|
if (GetNEPP(GetEdge(j).PNum(2)) != 2 || IsLineEndPoint(GetEdge(j).PNum(2)))
|
|
{
|
|
starte = j;
|
|
found = 1;
|
|
rev = 1;
|
|
}
|
|
else if (second)
|
|
{
|
|
starte = j;
|
|
found = 1;
|
|
rev = 0; //0 or 1 are possible
|
|
}
|
|
}
|
|
j++;
|
|
if (!second && j == GetNE()) {second = 1; j = 1;}
|
|
}
|
|
|
|
if (!found) {PrintSysError("No starting edge found, edgecnt=", edgecnt, ", GETNE=", GetNE());}
|
|
|
|
line->AddPoint(GetEdge(starte).PNum(1+rev));
|
|
line->AddPoint(GetEdge(starte).PNum(2-rev));
|
|
if (!rev)
|
|
{
|
|
line->AddLeftTrig(GetEdge(starte).LeftTrig());
|
|
line->AddRightTrig(GetEdge(starte).RightTrig());
|
|
}
|
|
else
|
|
{
|
|
line->AddLeftTrig(GetEdge(starte).RightTrig());
|
|
line->AddRightTrig(GetEdge(starte).LeftTrig());
|
|
}
|
|
edgecnt++; we.Elem(starte) = 1;
|
|
|
|
//add segments to line as long as segments other than starting edge are found or lineendpoint is reached
|
|
found = 1;
|
|
int other;
|
|
while(found)
|
|
{
|
|
found = 0;
|
|
int fp = GetEdge(starte).PNum(2-rev);
|
|
if (GetNEPP(fp) == 2 && !IsLineEndPoint(fp))
|
|
{
|
|
//find the "other" edge of point fp
|
|
other = 0;
|
|
if (GetEdgePP(fp,1) == starte) {other = 1;}
|
|
|
|
starte = GetEdgePP(fp,1+other);
|
|
|
|
//falls ring -> aufhoeren !!!!!!!!!!!
|
|
if (!we.Elem(starte))
|
|
{
|
|
found = 1;
|
|
rev = 0;
|
|
if (GetEdge(starte).PNum(2) == fp) {rev = 1;}
|
|
else if (GetEdge(starte).PNum(1) != fp) {PrintSysError("In Link Edges!");}
|
|
|
|
line->AddPoint(GetEdge(starte).PNum(2-rev));
|
|
if (!rev)
|
|
{
|
|
line->AddLeftTrig(GetEdge(starte).LeftTrig());
|
|
line->AddRightTrig(GetEdge(starte).RightTrig());
|
|
}
|
|
else
|
|
{
|
|
line->AddLeftTrig(GetEdge(starte).RightTrig());
|
|
line->AddRightTrig(GetEdge(starte).LeftTrig());
|
|
}
|
|
edgecnt++; we.Elem(starte) = 1;
|
|
}
|
|
}
|
|
}
|
|
AddLine(line);
|
|
}
|
|
PrintMessage(5,"number of lines generated = ", GetNLines());
|
|
|
|
//check, which lines must have at least one midpoint
|
|
INDEX_2_HASHTABLE<int> lineht(GetNLines()+1);
|
|
|
|
for (i = 1; i <= GetNLines(); i++)
|
|
{
|
|
if (GetLine(i)->StartP() == GetLine(i)->EndP())
|
|
{
|
|
GetLine(i)->DoSplit();
|
|
}
|
|
}
|
|
|
|
for (i = 1; i <= GetNLines(); i++)
|
|
{
|
|
INDEX_2 lineep (GetLine(i)->StartP(),GetLine(i)->EndP());
|
|
lineep.Sort();
|
|
|
|
if (lineht.Used (lineep))
|
|
{
|
|
GetLine(i)->DoSplit();
|
|
int other = lineht.Get(lineep);
|
|
GetLine(other)->DoSplit();
|
|
}
|
|
else
|
|
{
|
|
lineht.Set (lineep, i);
|
|
}
|
|
}
|
|
|
|
for (i = 1; i <= GetNLines(); i++)
|
|
{
|
|
STLLine* line = GetLine(i);
|
|
for (int ii = 1; ii <= line->GetNS(); ii++)
|
|
{
|
|
int ap1, ap2;
|
|
line->GetSeg(ii,ap1,ap2);
|
|
// (*mycout) << "SEG " << p1 << " - " << p2 << endl;
|
|
}
|
|
}
|
|
|
|
PopStatus();
|
|
}
|
|
|
|
int STLGeometry :: GetNOBodys()
|
|
{
|
|
int markedtrigs1 = 0;
|
|
int starttrig = 1;
|
|
int i, k, nnt;
|
|
int bodycnt = 0;
|
|
|
|
NgArray<int> bodynum(GetNT());
|
|
|
|
for (i = 1; i <= GetNT(); i++)
|
|
bodynum.Elem(i)=0;
|
|
|
|
|
|
while (markedtrigs1 < GetNT())
|
|
{
|
|
for (i = starttrig; i <= GetNT(); i++)
|
|
{
|
|
if (!bodynum.Get(i))
|
|
{
|
|
starttrig = i;
|
|
break;
|
|
}
|
|
}
|
|
//add all triangles around starttriangle, which is reachable without going over an edge
|
|
NgArray<int> todolist;
|
|
NgArray<int> nextlist;
|
|
bodycnt++;
|
|
markedtrigs1++;
|
|
bodynum.Elem(starttrig) = bodycnt;
|
|
todolist.Append(starttrig);
|
|
|
|
while(todolist.Size())
|
|
{
|
|
for (i = 1; i <= todolist.Size(); i++)
|
|
{
|
|
//const STLTriangle& tt = GetTriangle(todolist.Get(i));
|
|
for (k = 1; k <= NONeighbourTrigs(todolist.Get(i)); k++)
|
|
{
|
|
nnt = NeighbourTrig(todolist.Get(i),k);
|
|
if (!bodynum.Get(nnt))
|
|
{
|
|
nextlist.Append(nnt);
|
|
bodynum.Elem(nnt) = bodycnt;
|
|
markedtrigs1++;
|
|
}
|
|
}
|
|
}
|
|
|
|
todolist.SetSize(0);
|
|
for (i = 1; i <= nextlist.Size(); i++)
|
|
{
|
|
todolist.Append(nextlist.Get(i));
|
|
}
|
|
nextlist.SetSize(0);
|
|
}
|
|
}
|
|
PrintMessage(3, "Geometry has ", bodycnt, " separated bodys");
|
|
|
|
return bodycnt;
|
|
}
|
|
|
|
void STLGeometry :: CalcFaceNums()
|
|
{
|
|
int markedtrigs1 = 0;
|
|
int starttrig(0);
|
|
int laststarttrig = 1;
|
|
int i, k, nnt;
|
|
facecnt = 0;
|
|
|
|
|
|
for (i = 1; i <= GetNT(); i++)
|
|
GetTriangle(i).SetFaceNum(0);
|
|
|
|
|
|
while (markedtrigs1 < GetNT())
|
|
{
|
|
for (i = laststarttrig; i <= GetNT(); i++)
|
|
{
|
|
if (!GetTriangle(i).GetFaceNum())
|
|
{
|
|
starttrig = i;
|
|
laststarttrig = i;
|
|
break;
|
|
}
|
|
}
|
|
//add all triangles around starttriangle, which is reachable without going over an edge
|
|
NgArray<int> todolist;
|
|
NgArray<int> nextlist;
|
|
facecnt++;
|
|
markedtrigs1++;
|
|
GetTriangle(starttrig).SetFaceNum(facecnt);
|
|
todolist.Append(starttrig);
|
|
STLPointId ap1, ap2;
|
|
|
|
while(todolist.Size())
|
|
{
|
|
for (i = 1; i <= todolist.Size(); i++)
|
|
{
|
|
const STLTriangle& tt = GetTriangle(todolist.Get(i));
|
|
for (k = 1; k <= NONeighbourTrigs(todolist.Get(i)); k++)
|
|
{
|
|
nnt = NeighbourTrig(todolist.Get(i),k);
|
|
STLTriangle& nt = GetTriangle(nnt);
|
|
if (!nt.GetFaceNum())
|
|
{
|
|
tt.GetNeighbourPoints(nt,ap1,ap2);
|
|
if (!IsEdge(ap1,ap2))
|
|
{
|
|
nextlist.Append(nnt);
|
|
nt.SetFaceNum(facecnt);
|
|
markedtrigs1++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
todolist.SetSize(0);
|
|
for (i = 1; i <= nextlist.Size(); i++)
|
|
{
|
|
todolist.Append(nextlist.Get(i));
|
|
}
|
|
nextlist.SetSize(0);
|
|
}
|
|
}
|
|
GetNOBodys();
|
|
PrintMessage(3,"generated ", facecnt, " faces");
|
|
}
|
|
|
|
void STLGeometry :: ClearSpiralPoints()
|
|
{
|
|
spiralpoints.SetSize(GetNP());
|
|
int i;
|
|
for (i = 1; i <= spiralpoints.Size(); i++)
|
|
{
|
|
spiralpoints.Elem(i) = 0;
|
|
}
|
|
}
|
|
|
|
|
|
void STLGeometry :: BuildSmoothEdges ()
|
|
{
|
|
if (smoothedges) delete smoothedges;
|
|
|
|
smoothedges = new INDEX_2_HASHTABLE<int> (GetNE()/10 + 1);
|
|
|
|
|
|
// Jack: Ok ?
|
|
// UseExternalEdges();
|
|
|
|
PushStatusF("Build Smooth Edges");
|
|
|
|
int nt = GetNT();
|
|
Vec3d ng1, ng2;
|
|
|
|
for (int i = 1; i <= nt; i++)
|
|
{
|
|
if (multithread.terminate)
|
|
{PopStatus();return;}
|
|
|
|
SetThreadPercent(100.0 * (double)i / (double)nt);
|
|
|
|
const STLTriangle & trig = GetTriangle (i);
|
|
|
|
ng1 = trig.GeomNormal(points);
|
|
ng1 /= (ng1.Length() + 1e-24);
|
|
|
|
for (int j = 1; j <= 3; j++)
|
|
{
|
|
int nbt = NeighbourTrig (i, j);
|
|
|
|
ng2 = GetTriangle(nbt).GeomNormal(points);
|
|
ng2 /= (ng2.Length() + 1e-24);
|
|
|
|
STLPointId pi1, pi2;
|
|
trig.GetNeighbourPoints(GetTriangle(nbt), pi1, pi2);
|
|
|
|
if (!IsEdge(pi1,pi2))
|
|
{
|
|
if (ng1 * ng2 < 0)
|
|
{
|
|
PrintMessage(7,"smoothedge found");
|
|
INDEX_2 i2(pi1, pi2);
|
|
i2.Sort();
|
|
smoothedges->Set (i2, 1);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
PopStatus();
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
bool STLGeometry :: IsSmoothEdge (int pi1, int pi2) const
|
|
{
|
|
if (!smoothedges)
|
|
return false;
|
|
INDEX_2 i2(pi1, pi2);
|
|
i2.Sort();
|
|
return smoothedges->Used (i2);
|
|
}
|
|
|
|
|
|
|
|
/*
|
|
//function is not used now
|
|
int IsInArray(int n, const NgArray<int>& ia)
|
|
{
|
|
int i;
|
|
for (i = 1; i <= ia.Size(); i++)
|
|
{
|
|
if (ia.Get(i) == n) {return 1;}
|
|
}
|
|
return 0;
|
|
}
|
|
*/
|
|
|
|
void STLGeometry :: AddConeAndSpiralEdges(const STLParameters& stlparam)
|
|
{
|
|
PrintMessage(5,"have now ", GetNE(), " edges with yellow angle = ", stlparam.yangle, " degree");
|
|
|
|
PrintFnStart("AddConeAndSpiralEdges");
|
|
|
|
// int i,j,k,n;
|
|
// int changed = 0;
|
|
|
|
//check edges, where inner chart and no outer chart come together without an edge
|
|
STLPointId np1, np2;
|
|
int cnt = 0;
|
|
|
|
for (ChartId i = 1; i <= GetNOCharts(); i++)
|
|
{
|
|
STLChart& chart = GetChart(i);
|
|
for (int j = 1; j <= chart.GetNChartT(); j++)
|
|
{
|
|
STLTrigId t = chart.GetChartTrig1(j);
|
|
const STLTriangle& tt = GetTriangle(t);
|
|
|
|
for (int k = 1; k <= 3; k++)
|
|
{
|
|
STLTrigId nt = NeighbourTrig(t,k);
|
|
if (GetChartNr(nt) != i && !TrigIsInOC(nt,i))
|
|
{
|
|
tt.GetNeighbourPoints(GetTriangle(nt),np1,np2);
|
|
if (!IsEdge(np1,np2))
|
|
{
|
|
STLEdge se(np1,np2);
|
|
se.SetLeftTrig(t);
|
|
se.SetRightTrig(nt);
|
|
int edgenum = AddEdge(se);
|
|
AddEdgePP(np1,edgenum);
|
|
AddEdgePP(np2,edgenum);
|
|
//changed = 1;
|
|
PrintWarning("Found a spiral like structure: chart=", int(i),
|
|
", trig=", int(t), ", p1=", int(np1), ", p2=", int(np2));
|
|
cnt++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
PrintMessage(5, "found ", cnt, " spiral like structures");
|
|
PrintMessage(5, "added ", cnt, " edges due to spiral like structures");
|
|
|
|
cnt = 0;
|
|
int edgecnt = 0;
|
|
|
|
Array<STLTrigId> trigsaroundp;
|
|
NgArray<int> chartpointchecked(GetNP()); //gets number of chart, if in this chart already checked
|
|
chartpointchecked = 0;
|
|
|
|
|
|
int onoc, notonoc, tpp, pn;
|
|
STLPointId ap1, ap2;
|
|
int tn1, tn2, l, problem;
|
|
|
|
if (!stldoctor.conecheck) {PrintWarning("++++++++++++ \ncone checking deactivated by user!!!!!\n+++++++++++++++"); return ;}
|
|
|
|
PushStatus("Find Critical Points");
|
|
|
|
int addedges = 0;
|
|
|
|
for (ChartId i = 1; i <= GetNOCharts(); i++)
|
|
{
|
|
SetThreadPercent((double)i/(double)GetNOCharts()*100.);
|
|
if (multithread.terminate)
|
|
{PopStatus();return;}
|
|
|
|
STLChart& chart = GetChart(i);
|
|
for (int j = 1; j <= chart.GetNChartT(); j++)
|
|
{
|
|
STLTrigId t = chart.GetChartTrig1(j);
|
|
const STLTriangle& tt = GetTriangle(t);
|
|
|
|
for (int k = 1; k <= 3; k++)
|
|
{
|
|
pn = tt.PNum(k);
|
|
if (chartpointchecked.Get(pn) == i)
|
|
{continue;}
|
|
|
|
int checkpoint = 0;
|
|
for (int n = 1; n <= trigsperpoint.EntrySize(pn); n++)
|
|
{
|
|
if (trigsperpoint.Get(pn,n) != t &&
|
|
GetChartNr(trigsperpoint.Get(pn,n)) != i &&
|
|
!TrigIsInOC(trigsperpoint.Get(pn,n),i)) {checkpoint = 1;};
|
|
}
|
|
if (checkpoint)
|
|
{
|
|
chartpointchecked.Elem(pn) = i;
|
|
|
|
int worked = 0;
|
|
int spworked = 0;
|
|
GetSortedTrianglesAroundPoint(pn,t,trigsaroundp);
|
|
trigsaroundp.Append(t);
|
|
|
|
problem = 0;
|
|
for (int l = 2; l <= trigsaroundp.Size()-1; l++)
|
|
{
|
|
tn1 = trigsaroundp[l-2];
|
|
tn2 = trigsaroundp[l-1];
|
|
const STLTriangle& t1 = GetTriangle(tn1);
|
|
const STLTriangle& t2 = GetTriangle(tn2);
|
|
t1.GetNeighbourPoints(t2, ap1, ap2);
|
|
if (IsEdge(ap1,ap2)) break;
|
|
|
|
if (GetChartNr(tn2) != i && !TrigIsInOC(tn2,i)) {problem = 1;}
|
|
}
|
|
|
|
if (problem)
|
|
{
|
|
for (int l = 2; l <= trigsaroundp.Size()-1; l++)
|
|
{
|
|
tn1 = trigsaroundp[l-2];
|
|
tn2 = trigsaroundp[l-1];
|
|
const STLTriangle& t1 = GetTriangle(tn1);
|
|
const STLTriangle& t2 = GetTriangle(tn2);
|
|
t1.GetNeighbourPoints(t2, ap1, ap2);
|
|
if (IsEdge(ap1,ap2)) break;
|
|
|
|
if ((GetChartNr(tn1) == i && GetChartNr(tn2) != i && TrigIsInOC(tn2,i)) ||
|
|
(GetChartNr(tn2) == i && GetChartNr(tn1) != i && TrigIsInOC(tn1,i)))
|
|
{
|
|
if (addedges || !GetNEPP(pn))
|
|
{
|
|
STLEdge se(ap1,ap2);
|
|
se.SetLeftTrig(tn1);
|
|
se.SetRightTrig(tn2);
|
|
int edgenum = AddEdge(se);
|
|
AddEdgePP(ap1,edgenum);
|
|
AddEdgePP(ap2,edgenum);
|
|
edgecnt++;
|
|
}
|
|
if (!addedges && !GetSpiralPoint(pn))
|
|
{
|
|
SetSpiralPoint(pn);
|
|
spworked = 1;
|
|
}
|
|
worked = 1;
|
|
}
|
|
}
|
|
}
|
|
//backwards:
|
|
problem = 0;
|
|
for (int l = trigsaroundp.Size()-1; l >= 2; l--)
|
|
{
|
|
tn1 = trigsaroundp[l];
|
|
tn2 = trigsaroundp[l-1];
|
|
const STLTriangle& t1 = GetTriangle(tn1);
|
|
const STLTriangle& t2 = GetTriangle(tn2);
|
|
t1.GetNeighbourPoints(t2, ap1, ap2);
|
|
if (IsEdge(ap1,ap2)) break;
|
|
|
|
if (GetChartNr(tn2) != i && !TrigIsInOC(tn2,i)) {problem = 1;}
|
|
}
|
|
if (problem)
|
|
for (int l = trigsaroundp.Size()-1; l >= 2; l--)
|
|
{
|
|
tn1 = trigsaroundp[l];
|
|
tn2 = trigsaroundp[l-1];
|
|
const STLTriangle& t1 = GetTriangle(tn1);
|
|
const STLTriangle& t2 = GetTriangle(tn2);
|
|
t1.GetNeighbourPoints(t2, ap1, ap2);
|
|
if (IsEdge(ap1,ap2)) break;
|
|
|
|
if ((GetChartNr(tn1) == i && GetChartNr(tn2) != i && TrigIsInOC(tn2,i)) ||
|
|
(GetChartNr(tn2) == i && GetChartNr(tn1) != i && TrigIsInOC(tn1,i)))
|
|
{
|
|
if (addedges || !GetNEPP(pn))
|
|
{
|
|
STLEdge se(ap1,ap2);
|
|
se.SetLeftTrig(tn1);
|
|
se.SetRightTrig(tn2);
|
|
int edgenum = AddEdge(se);
|
|
AddEdgePP(ap1,edgenum);
|
|
AddEdgePP(ap2,edgenum);
|
|
edgecnt++;
|
|
}
|
|
if (!addedges && !GetSpiralPoint(pn))
|
|
{
|
|
SetSpiralPoint(pn);
|
|
spworked = 1;
|
|
//if (GetNEPP(pn) == 0) {(*mycout) << "ERROR: spiralpoint with no edge found!" << endl;}
|
|
}
|
|
worked = 1;
|
|
}
|
|
}
|
|
|
|
if (worked)
|
|
{
|
|
//(*testout) << "set edgepoint due to spirals: pn=" << i << endl;
|
|
SetLineEndPoint(pn);
|
|
}
|
|
if (spworked)
|
|
{
|
|
/*
|
|
(*mycout) << "Warning: Critical Point " << tt.PNum(k)
|
|
<< "( chart " << i << ", trig " << t
|
|
<< ") has been neutralized!!!" << endl;
|
|
*/
|
|
cnt++;
|
|
}
|
|
// markedpoints.Elem(tt.PNum(k)) = 1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
PrintMessage(5, "found ", cnt, " critical points!");
|
|
PrintMessage(5, "added ", edgecnt, " edges due to critical points!");
|
|
|
|
PopStatus();
|
|
|
|
//search points where inner chart and outer chart and "no chart" trig come together at edge-point
|
|
|
|
PrintMessage(7,"search for special chart points");
|
|
for (ChartId i = 1; i <= GetNOCharts(); i++)
|
|
{
|
|
STLChart& chart = GetChart(i);
|
|
for (int j = 1; j <= chart.GetNChartT(); j++)
|
|
{
|
|
STLTrigId t = chart.GetChartTrig1(j);
|
|
const STLTriangle& tt = GetTriangle(t);
|
|
|
|
for (int k = 1; k <= 3; k++)
|
|
{
|
|
pn = tt.PNum(k);
|
|
if (GetNEPP(pn) == 2)
|
|
{
|
|
onoc = 0;
|
|
notonoc = 0;
|
|
for (int n = 1; n <= trigsperpoint.EntrySize(pn); n++)
|
|
{
|
|
tpp = trigsperpoint.Get(pn,n);
|
|
if (tpp != t && GetChartNr(tpp) != i)
|
|
{
|
|
if (TrigIsInOC(tpp,i)) {onoc = 1;}
|
|
if (!TrigIsInOC(tpp,i)) {notonoc = 1;}
|
|
}
|
|
}
|
|
if (onoc && notonoc && !IsLineEndPoint(pn))
|
|
{
|
|
GetSortedTrianglesAroundPoint(pn,t,trigsaroundp);
|
|
int here = 1; //we start on this side of edge, !here = there
|
|
int thereOC = 0;
|
|
int thereNotOC = 0;
|
|
for (l = 2; l <= trigsaroundp.Size(); l++)
|
|
{
|
|
GetTriangle(trigsaroundp[l-2]).
|
|
GetNeighbourPoints(GetTriangle(trigsaroundp[l-1]), ap1, ap2);
|
|
if (IsEdge(ap1,ap2)) {here = (here+1)%2;}
|
|
if (!here && TrigIsInOC(trigsaroundp[l-1],i)) {thereOC = 1;}
|
|
if (!here && !TrigIsInOC(trigsaroundp[l-1],i)) {thereNotOC = 1;}
|
|
}
|
|
if (thereOC && thereNotOC)
|
|
{
|
|
//(*mycout) << "Special OCICnotC - point " << pn << " found!" << endl;
|
|
//(*testout) << "set edgepoint due to spirals: pn=" << i << endl;
|
|
SetLineEndPoint(pn);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
PrintMessage(5,"have now ", GetNE(), " edges with yellow angle = ", stlparam.yangle, " degree");
|
|
}
|
|
|
|
//get trigs at a point, started with starttrig, then every left
|
|
void STLGeometry :: GetSortedTrianglesAroundPoint(STLPointId p, STLTrigId starttrig, Array<STLTrigId>& trigs)
|
|
{
|
|
STLTrigId acttrig = starttrig;
|
|
trigs.SetAllocSize(trigsperpoint.EntrySize(p));
|
|
trigs.SetSize(0);
|
|
trigs.Append(acttrig);
|
|
int locindex1(0), locindex2(0);
|
|
//(*mycout) << "trigs around point " << p << endl;
|
|
|
|
int end = 0;
|
|
while (!end)
|
|
{
|
|
const STLTriangle& at = GetTriangle(acttrig);
|
|
for (int i = 1; i <= trigsperpoint.EntrySize(p); i++)
|
|
{
|
|
STLTrigId t = trigsperpoint.Get(p,i);
|
|
const STLTriangle& nt = GetTriangle(t);
|
|
if (at.IsNeighbourFrom(nt))
|
|
{
|
|
STLPointId ap1, ap2;
|
|
at.GetNeighbourPoints(nt, ap1, ap2);
|
|
if (ap2 == p) {Swap(ap1,ap2);}
|
|
if (ap1 != p) {PrintSysError("In GetSortedTrianglesAroundPoint!!!");}
|
|
|
|
for (int j = 1; j <= 3; j++)
|
|
{
|
|
if (at.PNum(j) == ap1) {locindex1 = j;};
|
|
if (at.PNum(j) == ap2) {locindex2 = j;};
|
|
}
|
|
if ((locindex2+1)%3+1 == locindex1)
|
|
{
|
|
if (t != starttrig)
|
|
{
|
|
trigs.Append(t);
|
|
// (*mycout) << "trig " << t << endl;
|
|
acttrig = t;
|
|
}
|
|
else
|
|
{
|
|
end = 1;
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
int STLGeometry :: NeighbourTrig(int trig, int nr) const
|
|
{
|
|
return neighbourtrigs.Get(trig,nr);
|
|
}
|
|
*/
|
|
|
|
|
|
|
|
void STLGeometry :: SmoothGeometry ()
|
|
{
|
|
int i, j, k;
|
|
|
|
double maxerr0, maxerr;
|
|
|
|
for (i = 1; i <= GetNP(); i++)
|
|
{
|
|
if (GetNEPP(i)) continue;
|
|
|
|
maxerr0 = 0;
|
|
for (j = 1; j <= NOTrigsPerPoint(i); j++)
|
|
{
|
|
int tnum = TrigPerPoint(i, j);
|
|
double err = Angle (GetTriangle(tnum).Normal (),
|
|
GetTriangle(tnum).GeomNormal(GetPoints()));
|
|
if (err > maxerr0)
|
|
maxerr0 = err;
|
|
}
|
|
|
|
Point3d pi = GetPoint (i);
|
|
if (maxerr0 < 1.1) continue; // about 60 degree
|
|
|
|
maxerr0 /= 2; // should be at least halfen
|
|
|
|
for (k = 1; k <= NOTrigsPerPoint(i); k++)
|
|
{
|
|
const STLTriangle & trig = GetTriangle (TrigPerPoint (i, k));
|
|
Point3d c = Center(GetPoint (trig.PNum(1)),
|
|
GetPoint (trig.PNum(2)),
|
|
GetPoint (trig.PNum(3)));
|
|
|
|
Point3d np = pi + 0.1 * (c - pi);
|
|
SetPoint (i, np);
|
|
|
|
maxerr = 0;
|
|
for (j = 1; j <= NOTrigsPerPoint(i); j++)
|
|
{
|
|
int tnum = TrigPerPoint(i, j);
|
|
double err = Angle (GetTriangle(tnum).Normal (),
|
|
GetTriangle(tnum).GeomNormal(GetPoints()));
|
|
if (err > maxerr)
|
|
maxerr = err;
|
|
}
|
|
|
|
if (maxerr < maxerr0)
|
|
{
|
|
pi = np;
|
|
}
|
|
}
|
|
|
|
SetPoint (i, pi);
|
|
}
|
|
}
|
|
|
|
void STLGeometry :: WriteChartToFile( ChartId chartnumber, filesystem::path filename )
|
|
{
|
|
PrintMessage(1,"write chart ", int(chartnumber), " to ", filename);
|
|
Array<int> trignums;
|
|
|
|
if (chartnumber >= 1 && chartnumber <= GetNOCharts())
|
|
{
|
|
const STLChart& chart = GetChart(chartnumber);
|
|
|
|
for (int j = 1; j <= chart.GetNChartT(); j++)
|
|
trignums.Append(chart.GetChartTrig1(j));
|
|
|
|
for (int j = 1; j <= chart.GetNOuterT(); j++)
|
|
trignums.Append(chart.GetOuterTrig1(j));
|
|
|
|
QuickSort(trignums);
|
|
STLGeometry geo;
|
|
NgArray<STLReadTriangle> readtrigs;
|
|
const auto & first_trig = GetTriangle(chart.GetChartTrig1(1));
|
|
auto normal = first_trig.Normal();
|
|
Box<3> box{Box<3>::EMPTY_BOX};
|
|
|
|
for(auto j : trignums)
|
|
{
|
|
const auto& trig = GetTriangle(j);
|
|
Point<3> pts[3];
|
|
for(auto k : Range(3))
|
|
{
|
|
pts[k] = GetPoint(trig[k]);
|
|
box.Add(pts[k]);
|
|
}
|
|
Vec3d normal = Cross( pts[1]-pts[0], pts[2]-pts[0] );
|
|
readtrigs.Append(STLReadTriangle(pts, trig.Normal()));
|
|
}
|
|
auto dist = box.PMax() - box.PMin();
|
|
auto extra_point = GetPoint(first_trig[0]) - dist.Length()*normal;
|
|
|
|
NgArray<int> acttrigs(GetNT());
|
|
acttrigs = -1;
|
|
for (int j = 1; j <= chart.GetNT(); j++)
|
|
acttrigs.Elem(chart.GetTrig1(j)) = chartnumber;
|
|
|
|
for (int j = 1; j <= chart.GetNT(); j++)
|
|
{
|
|
auto t = chart.GetTrig1(j);
|
|
const auto & tt = GetTriangle(t);
|
|
for (int k = 1; k <= 3; k++)
|
|
{
|
|
int nt = NeighbourTrig(t,k);
|
|
if (acttrigs.Get(nt) != chartnumber)
|
|
{
|
|
STLPointId np1, np2;
|
|
tt.GetNeighbourPoints(GetTriangle(nt),np1,np2);
|
|
|
|
Point<3> pts[3];
|
|
pts[0] = GetPoint(np2);
|
|
pts[1] = GetPoint(np1);
|
|
pts[2] = extra_point;
|
|
Vec3d normal = -Cross( pts[2]-pts[0], pts[1]-pts[0] );
|
|
readtrigs.Append(STLReadTriangle(pts, normal));
|
|
}
|
|
}
|
|
}
|
|
|
|
geo.InitSTLGeometry(readtrigs);
|
|
geo.Save(filename);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
class STLGeometryRegister : public GeometryRegister
|
|
{
|
|
public:
|
|
virtual NetgenGeometry * Load (const filesystem::path & filename) const;
|
|
};
|
|
|
|
NetgenGeometry * STLGeometryRegister :: Load (const filesystem::path & filename) const
|
|
{
|
|
string ext = ToLower(filename.extension());
|
|
|
|
if (ext == ".stl")
|
|
{
|
|
PrintMessage (1, "Load STL geometry file ", filename);
|
|
|
|
ifstream infile(filename);
|
|
|
|
STLGeometry * hgeom = STLGeometry :: Load (infile);
|
|
hgeom -> edgesfound = 0;
|
|
return hgeom;
|
|
}
|
|
else if (ext == ".stlb")
|
|
{
|
|
PrintMessage (1, "Load STL binary geometry file ", filename);
|
|
|
|
ifstream infile(filename);
|
|
|
|
STLGeometry * hgeom = STLGeometry :: LoadBinary (infile);
|
|
hgeom -> edgesfound = 0;
|
|
return hgeom;
|
|
}
|
|
else if (ext == ".nao")
|
|
{
|
|
PrintMessage (1, "Load naomi (F. Kickinger) geometry file ", filename);
|
|
|
|
ifstream infile(filename);
|
|
|
|
STLGeometry * hgeom = STLGeometry :: LoadNaomi (infile);
|
|
hgeom -> edgesfound = 0;
|
|
return hgeom;
|
|
}
|
|
|
|
|
|
return NULL;
|
|
}
|
|
|
|
|
|
class STLInit
|
|
{
|
|
public:
|
|
STLInit()
|
|
{
|
|
GeometryRegisterArray& gra = FetchGeometryRegisterArray();
|
|
gra.Append (new STLGeometryRegister);
|
|
}
|
|
};
|
|
|
|
STLInit stlinit;
|
|
|
|
static RegisterClassForArchive<STLGeometry, NetgenGeometry, STLTopology> stlgeo;
|
|
}
|