netgen/libsrc/stlgeom/stltool.cpp
2019-10-01 12:34:23 +02:00

1495 lines
33 KiB
C++

#include <mystdlib.h>
#include <myadt.hpp>
#include <linalg.hpp>
#include <gprim.hpp>
#include <meshing.hpp>
#include "stlgeom.hpp"
namespace netgen
{
//add a point into a pointlist, return pointnumber
int AddPointIfNotExists(NgArray<Point3d>& ap, const Point3d& p, double eps)
{
double eps2 = sqr(eps);
for (int i = 1; i <= ap.Size(); i++)
if (Dist2(ap.Get(i),p) <= eps2 )
return i;
ap.Append(p);
return ap.Size();
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
double GetDistFromLine(const Point<3> & lp1, const Point<3> & lp2,
Point<3> & p)
{
Vec3d vn = lp2 - lp1;
Vec3d v1 = p - lp1;
Vec3d v2 = lp2 - p;
Point3d pold = p;
if (v2 * vn <= 0) {p = lp2; return (pold - p).Length();}
if (v1 * vn <= 0) {p = lp1; return (pold - p).Length();}
double vnl = vn.Length();
if (vnl == 0) {return Dist(lp1,p);}
vn /= vnl;
p = lp1 + (v1 * vn) * vn;
return (pold - p).Length();
};
double GetDistFromInfiniteLine(const Point<3>& lp1, const Point<3>& lp2, const Point<3>& p)
{
Vec3d vn(lp1, lp2);
Vec3d v1(lp1, p);
double vnl = vn.Length();
if (vnl == 0)
{
return Dist (lp1, p);
}
else
{
return Cross (vn, v1).Length() / vnl;
}
};
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Binary IO-Manipulation
void FIOReadInt(istream& ios, int& i)
{
const int ilen = sizeof(int);
char buf[ilen];
for (int j = 0; j < ilen; j++)
ios.get(buf[j]);
memcpy(&i, &buf, ilen);
}
void FIOWriteInt(ostream& ios, const int& i)
{
const int ilen = sizeof(int);
char buf[ilen];
memcpy(&buf, &i, ilen);
for (int j = 0; j < ilen; j++)
ios << buf[j];
}
void FIOReadDouble(istream& ios, double& i)
{
const int ilen = sizeof(double);
char buf[ilen];
for (int j = 0; j < ilen; j++)
ios.get(buf[j]);
memcpy(&i, &buf, ilen);
}
void FIOWriteDouble(ostream& ios, const double& i)
{
const int ilen = sizeof(double);
char buf[ilen];
memcpy(&buf, &i, ilen);
for (int j = 0; j < ilen; j++)
ios << buf[j];
}
void FIOReadFloat(istream& ios, float& i)
{
const int ilen = sizeof(float);
char buf[ilen];
int j;
for (j = 0; j < ilen; j++)
{
ios.get(buf[j]);
}
memcpy(&i, &buf, ilen);
}
void FIOWriteFloat(ostream& ios, const float& i)
{
const int ilen = sizeof(float);
char buf[ilen];
memcpy(&buf, &i, ilen);
for (int j = 0; j < ilen; j++)
ios << buf[j];
}
void FIOReadString(istream& ios, char* str, int len)
{
for (int j = 0; j < len; j++)
ios.get(str[j]);
}
//read string and add terminating 0
void FIOReadStringE(istream& ios, char* str, int len)
{
for (int j = 0; j < len; j++)
ios.get(str[j]);
str[len] = 0;
}
void FIOWriteString(ostream& ios, char* str, int len)
{
for (int j = 0; j < len; j++)
ios << str[j];
}
/*
void FIOReadInt(istream& ios, int& i)
{
const int ilen = sizeof(int);
char buf[ilen];
int j;
for (j = 0; j < ilen; j++)
{
ios.get(buf[ilen-j-1]);
}
memcpy(&i, &buf, ilen);
}
void FIOWriteInt(ostream& ios, const int& i)
{
const int ilen = sizeof(int);
char buf[ilen];
memcpy(&buf, &i, ilen);
int j;
for (j = 0; j < ilen; j++)
{
ios << buf[ilen-j-1];
}
}
void FIOReadDouble(istream& ios, double& i)
{
const int ilen = sizeof(double);
char buf[ilen];
int j;
for (j = 0; j < ilen; j++)
{
ios.get(buf[ilen-j-1]);
}
memcpy(&i, &buf, ilen);
}
void FIOWriteDouble(ostream& ios, const double& i)
{
const int ilen = sizeof(double);
char buf[ilen];
memcpy(&buf, &i, ilen);
int j;
for (j = 0; j < ilen; j++)
{
ios << buf[ilen-j-1];
}
}
void FIOReadFloat(istream& ios, float& i)
{
const int ilen = sizeof(float);
char buf[ilen];
int j;
for (j = 0; j < ilen; j++)
{
ios.get(buf[ilen-j-1]);
}
memcpy(&i, &buf, ilen);
}
void FIOWriteFloat(ostream& ios, const float& i)
{
const int ilen = sizeof(float);
char buf[ilen];
memcpy(&buf, &i, ilen);
int j;
for (j = 0; j < ilen; j++)
{
ios << buf[ilen-j-1];
}
}
void FIOReadString(istream& ios, char* str, int len)
{
int j;
for (j = 0; j < len; j++)
{
ios.get(str[j]);
}
}
//read string and add terminating 0
void FIOReadStringE(istream& ios, char* str, int len)
{
int j;
for (j = 0; j < len; j++)
{
ios.get(str[j]);
}
str[len] = 0;
}
void FIOWriteString(ostream& ios, char* str, int len)
{
int j;
for (j = 0; j < len; j++)
{
ios << str[j];
}
}
*/
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
STLReadTriangle :: STLReadTriangle (const Point<3> * apts,
const Vec<3> & anormal)
{
pts[0] = apts[0];
pts[1] = apts[1];
pts[2] = apts[2];
normal = anormal;
}
STLTriangle :: STLTriangle(const STLPointId * apts)
{
pts[0] = apts[0];
pts[1] = apts[1];
pts[2] = apts[2];
facenum = 0;
}
int STLTriangle :: IsNeighbourFrom(const STLTriangle& t) const
{
//triangles must have same orientation!!!
for(int i = 0; i <= 2; i++)
for(int j = 0; j <= 2; j++)
if (t.pts[(i+1)%3] == pts[j] &&
t.pts[i] == pts[(j+1)%3])
return 1;
return 0;
}
int STLTriangle :: IsWrongNeighbourFrom(const STLTriangle& t) const
{
//triangles have not same orientation!!!
for(int i = 0; i <= 2; i++)
for(int j = 0; j <= 2; j++)
if (t.pts[(i+1)%3] == pts[(j+1)%3] &&
t.pts[i] == pts[j])
return 1;
return 0;
}
void STLTriangle :: GetNeighbourPoints(const STLTriangle& t, STLPointId & p1, STLPointId & p2) const
{
for(int i = 1; i <= 3; i++)
for(int j = 1; j <= 3; j++)
if (t.PNumMod(i+1) == PNumMod(j) &&
t.PNumMod(i) == PNumMod(j+1))
{
p1 = PNumMod(j);
p2 = PNumMod(j+1);
return;
}
PrintSysError("Get neighbourpoints failed!");
}
int STLTriangle :: GetNeighbourPointsAndOpposite(const STLTriangle& t, STLPointId & p1,
STLPointId & p2, STLPointId & po) const
{
for(int i = 1; i <= 3; i++)
for(int j = 1; j <= 3; j++)
if (t.PNumMod(i+1) == PNumMod(j) &&
t.PNumMod(i) == PNumMod(j+1))
{
p1 = PNumMod(j);
p2 = PNumMod(j+1);
po = PNumMod(j+2);
return 1;
}
return 0;
}
Vec<3> STLTriangle :: GeomNormal(const Array<Point<3>,STLPointId>& ap) const
{
const Point<3> & p1 = ap[PNum(1)];
const Point<3> & p2 = ap[PNum(2)];
const Point<3> & p3 = ap[PNum(3)];
return Cross(p2-p1, p3-p1);
}
void STLTriangle :: SetNormal (const Vec<3> & n)
{
double len = n.Length();
if (len > 0)
{
normal = n;
normal.Normalize();
}
else
{
normal = Vec<3> (1, 0, 0);
}
}
void STLTriangle :: ChangeOrientation()
{
normal *= -1;
Swap(pts[0],pts[1]);
}
double STLTriangle :: Area(const Array<Point<3>,STLPointId>& ap) const
{
return 0.5 * Cross(ap[PNum(2)]-ap[PNum(1)],
ap[PNum(3)]-ap[PNum(1)]).Length();
}
double STLTriangle :: MinHeight(const Array<Point<3>,STLPointId>& ap) const
{
double ml = MaxLength(ap);
if (ml != 0) {return 2.*Area(ap)/ml;}
PrintWarning("max Side Length of a triangle = 0!!!");
return 0;
}
double STLTriangle :: MaxLength(const Array<Point<3>,STLPointId>& ap) const
{
return max3(Dist(ap[PNum(1)],ap[PNum(2)]),
Dist(ap[PNum(2)],ap[PNum(3)]),
Dist(ap[PNum(3)],ap[PNum(1)]));
}
void STLTriangle :: ProjectInPlain(const Array<Point<3>,STLPointId>& ap,
const Vec<3> & n, Point<3> & pp) const
{
const Point<3> & p1 = ap[PNum(1)];
const Point<3> & p2 = ap[PNum(2)];
const Point<3> & p3 = ap[PNum(3)];
Vec<3> v1 = p2 - p1;
Vec<3> v2 = p3 - p1;
Vec<3> nt = Cross(v1, v2);
double c = - (p1(0)*nt(0) + p1(1)*nt(1) + p1(2)*nt(2));
double prod = n * nt;
if (fabs(prod) == 0)
{
pp = Point<3>(1.E20,1.E20,1.E20);
return;
}
double nfact = -(pp(0)*nt(0) + pp(1)*nt(1) + pp(2)*nt(2) + c) / (prod);
pp = pp + (nfact) * n;
}
int STLTriangle :: ProjectInPlain (const Array<Point<3>,STLPointId>& ap,
const Vec<3> & nproj,
Point<3> & pp, Vec<3> & lam) const
{
const Point<3> & p1 = ap[PNum(1)];
const Point<3> & p2 = ap[PNum(2)];
const Point<3> & p3 = ap[PNum(3)];
Vec<3> v1 = p2-p1;
Vec<3> v2 = p3-p1;
Mat<3> mat;
for (int i = 0; i < 3; i++)
{
mat(i,0) = v1(i);
mat(i,1) = v2(i);
mat(i,2) = nproj(i);
}
int err = 0;
mat.Solve (pp-p1, lam);
// int err = SolveLinearSystem (v1, v2, nproj, pp-p1, lam);
if (!err)
{
// pp = p1 + lam(0) * v1 + lam(1) * v2;
pp(0) = p1(0) + lam(0) * v1(0) + lam(1) * v2(0);
pp(1) = p1(1) + lam(0) * v1(1) + lam(1) * v2(1);
pp(2) = p1(2) + lam(0) * v1(2) + lam(1) * v2(2);
}
return err;
}
void STLTriangle :: ProjectInPlain(const Array<Point<3>,STLPointId>& ap,
Point<3> & pp) const
{
const Point<3> & p1 = ap[PNum(1)];
const Point<3> & p2 = ap[PNum(2)];
const Point<3> & p3 = ap[PNum(3)];
Vec<3> v1 = p2 - p1;
Vec<3> v2 = p3 - p1;
Vec<3> nt = Cross(v1, v2);
double c = - (p1(0)*nt(0) + p1(1)*nt(1) + p1(2)*nt(2));
double prod = nt * nt;
double nfact = -(pp(0)*nt(0) + pp(1)*nt(1) + pp(2)*nt(2) + c) / (prod);
pp = pp + (nfact) * nt;
}
bool STLTriangle :: PointInside(const Array<Point<3>,STLPointId> & ap,
const Point<3> & pp) const
{
const Point<3> & p1 = ap[PNum(1)];
const Point<3> & p2 = ap[PNum(2)];
const Point<3> & p3 = ap[PNum(3)];
Vec<3> v1 = p2 - p1;
Vec<3> v2 = p3 - p1;
Vec<3> v = pp - p1;
double det, l1, l2;
Vec<3> ex, ey, ez;
ez = GeomNormal(ap);
ez /= ez.Length();
ex = v1;
ex /= ex.Length();
ey = Cross (ez, ex);
Vec<2> v1p(v1*ex, v1*ey);
Vec<2> v2p(v2*ex, v2*ey);
Vec<2> vp(v*ex, v*ey);
det = v2p(1) * v1p(0) - v2p(0) * v1p(1);
if (fabs(det) == 0) {return 0;}
l2 = (vp(1) * v1p(0) - vp(0) * v1p(1)) / det;
if (v1p(0) != 0.)
{
l1 = (vp(0) - l2 * v2p(0)) / v1p(0);
}
else if (v1p(1) != 0.)
{
l1 = (vp(1) - l2 * v2p(1)) / v1p(1);
}
else {return 0;}
if (l1 >= -1E-10 && l2 >= -1E-10 && l1 + l2 <= 1.+1E-10) {return 1;}
return 0;
}
double STLTriangle :: GetNearestPoint(const Array<Point<3>,STLPointId>& ap,
Point<3> & p3d) const
{
Point<3> p = p3d;
ProjectInPlain(ap, p);
double dist = (p - p3d).Length();
if (PointInside(ap, p)) {p3d = p; return dist;}
else
{
Point<3> pf = 0.0;
double nearest = 1E50;
//int fi = 0;
for (int j = 1; j <= 3; j++)
{
p = p3d;
dist = GetDistFromLine(ap[PNum(j)], ap[PNumMod(j+1)], p);
if (dist < nearest)
{
nearest = dist;
pf = p;
}
}
p3d = pf;
return nearest;
}
}
bool STLTriangle :: HasEdge(STLPointId p1, STLPointId p2) const
{
for (int i = 1; i <= 3; i++)
if (p1 == PNum(i) && p2 == PNumMod(i+1))
return true;
return false;
}
ostream& operator<<(ostream& os, const STLTriangle& t)
{
os << "[";
os << t[0] << ",";
os << t[1] << ",";
os << t[2] << "]";
return os;
};
STLTopEdge :: STLTopEdge ()
{
pts[0] = pts[1] = 0;
trigs[0] = trigs[1] = 0;
cosangle = 1;
status = ED_UNDEFINED;
}
STLTopEdge :: STLTopEdge (STLPointId p1, STLPointId p2, int trig1, int trig2)
{
pts[0] = p1;
pts[1] = p2;
trigs[0] = trig1;
trigs[1] = trig2;
cosangle = 1;
status = ED_UNDEFINED;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//+++++++++++++++++++ STL CHART +++++++++++++++++++++++++++++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
STLChart :: STLChart(STLGeometry * ageometry, const STLParameters& astlparam)
: geometry(ageometry), stlparam(astlparam)
{
// charttrigs = new NgArray<int> (0,0);
// outertrigs = new NgArray<int> (0,0);
// ilimit = new NgArray<twoint> (0,0);
// olimit = new NgArray<twoint> (0,0);
geometry = ageometry;
if ( stlparam.usesearchtree == 1)
{
Box<3> box = geometry->GetBoundingBox();
box.Increase (0.2*box.Diam()+1e-12);
searchtree = new BoxTree<3,STLTrigId> (box);
/*
searchtree = new BoxTree<3> (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
geometry->GetBoundingBox().PMax() + Vec3d(1,1,1));
*/
}
else
searchtree = NULL;
}
STLChart :: ~STLChart()
{
delete searchtree;
}
void STLChart :: AddChartTrig(STLTrigId i)
{
// static int timer = NgProfiler::CreateTimer ("STLChart::AddChartTrig");
// NgProfiler::RegionTimer reg(timer);
charttrigs.Append(i);
const STLTriangle & trig = geometry->GetTriangle(i);
const Point<3> & p1 = geometry->GetPoint (trig.PNum(1));
const Point<3> & p2 = geometry->GetPoint (trig.PNum(2));
const Point<3> & p3 = geometry->GetPoint (trig.PNum(3));
/*
Point3d pmin(p1), pmax(p1);
pmin.SetToMin (p2);
pmin.SetToMin (p3);
pmax.SetToMax (p2);
pmax.SetToMax (p3);
*/
/*
Box<3> box(p1);
box.Add(p2);
box.Add(p3);
*/
Box<3> box(p1,p2,p3);
if (!geomsearchtreeon && (stlparam.usesearchtree == 1))
// {searchtree->Insert (pmin, pmax, i);}
{
searchtree->Insert (box, i);
}
}
void STLChart :: AddOuterTrig(STLTrigId i)
{
// static int timer = NgProfiler::CreateTimer ("STLChart::AddOuterTrig");
// NgProfiler::RegionTimer reg(timer);
outertrigs.Append(i);
const STLTriangle & trig = geometry->GetTriangle(i);
const Point3d & p1 = geometry->GetPoint (trig.PNum(1));
const Point3d & p2 = geometry->GetPoint (trig.PNum(2));
const Point3d & p3 = geometry->GetPoint (trig.PNum(3));
Point3d pmin(p1), pmax(p1);
pmin.SetToMin (p2);
pmin.SetToMin (p3);
pmax.SetToMax (p2);
pmax.SetToMax (p3);
if (!geomsearchtreeon && (stlparam.usesearchtree==1))
{searchtree->Insert (pmin, pmax, i);}
}
bool STLChart :: IsInWholeChart(int nr) const
{
return charttrigs.Contains(nr) || outertrigs.Contains(nr);
}
void STLChart :: GetTrianglesInBox (const Point3d & pmin,
const Point3d & pmax,
NgArray<STLTrigId> & trias) const
{
if (geomsearchtreeon) {PrintMessage(5,"geomsearchtreeon is set!!!");}
if (searchtree)
searchtree -> GetIntersecting (pmin, pmax, trias);
else
{
Box<3> box1(pmin, pmax);
box1.Increase (1e-2*box1.Diam());
trias.SetSize(0);
int nt = GetNT();
for (int i = 1; i <= nt; i++)
{
STLTrigId trignum = GetTrig1(i);
const STLTriangle & trig = geometry->GetTriangle(trignum);
Box<3> box2(geometry->GetPoint (trig.PNum(1)),
geometry->GetPoint (trig.PNum(2)),
geometry->GetPoint (trig.PNum(3)));
if (box1.Intersect (box2))
trias.Append (trignum);
}
}
}
//trigs may contain the same triangle double
void STLChart :: MoveToOuterChart(const NgArray<int>& trigs)
{
if (!trigs.Size()) return;
for (int i = 1; i <= trigs.Size(); i++)
{
if (charttrigs[trigs.Get(i)-1] != -1)
AddOuterTrig(charttrigs[trigs.Get(i)-1]);
charttrigs[trigs.Get(i)-1] = -1;
}
DelChartTrigs(trigs);
}
//trigs may contain the same triangle double
void STLChart :: DelChartTrigs(const NgArray<int>& trigs)
{
if (!trigs.Size()) return;
for (int i = 1; i <= trigs.Size(); i++)
charttrigs[trigs.Get(i)-1] = -1;
int cnt = 0;
for (int i = 1; i <= charttrigs.Size(); i++)
{
if (charttrigs[i-1] == -1)
cnt++;
if (cnt != 0 && i < charttrigs.Size())
charttrigs[i-cnt] = charttrigs[i];
}
int i = charttrigs.Size() - trigs.Size();
charttrigs.SetSize(i);
if (!geomsearchtreeon && stlparam.usesearchtree == 1)
{
PrintMessage(7, "Warning: unsecure routine due to first use of searchtrees!!!");
//bould new searchtree!!!
searchtree = new BoxTree<3,STLTrigId> (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
geometry->GetBoundingBox().PMax() + Vec3d(1,1,1));
for (int i = 1; i <= charttrigs.Size(); i++)
{
const STLTriangle & trig = geometry->GetTriangle(i);
const Point3d & p1 = geometry->GetPoint (trig.PNum(1));
const Point3d & p2 = geometry->GetPoint (trig.PNum(2));
const Point3d & p3 = geometry->GetPoint (trig.PNum(3));
Point3d pmin(p1), pmax(p1);
pmin.SetToMin (p2);
pmin.SetToMin (p3);
pmax.SetToMax (p2);
pmax.SetToMax (p3);
searchtree->Insert (pmin, pmax, i);
}
}
}
void STLChart :: SetNormal (const Point<3> & apref, const Vec<3> & anormal)
{
pref = apref;
normal = anormal;
double len = normal.Length();
if (len) normal /= len;
else normal = Vec<3> (1, 0, 0);
t1 = normal.GetNormal ();
t2 = Cross (normal, t1);
}
void STLChart :: BuildInnerSearchTree()
{
Box<2> chart_bbox(Box<2>::EMPTY_BOX);
for (STLTrigId trigid : charttrigs)
{
for (STLPointId pi : (*geometry)[trigid].PNums())
{
Point<3> p = (*geometry)[pi];
Point<2> p2d = Project2d(p);
chart_bbox.Add(p2d);
}
}
chart_bbox.Increase (1e-2*chart_bbox.Diam());
inner_searchtree = make_unique<BoxTree<2,STLTrigId>> (chart_bbox);
for (STLTrigId trigid : charttrigs)
{
Box<2> bbox(Box<2>::EMPTY_BOX);
for (STLPointId pi : (*geometry)[trigid].PNums())
{
Point<3> p = (*geometry)[pi];
Point<2> p2d = Project2d(p);
bbox.Add(p2d);
}
inner_searchtree->Insert (bbox, trigid);
}
}
STLTrigId STLChart :: ProjectNormal (Point<3> & p3d) const
{
int nt = GetNT();
double lamtol = 1e-6;
QuadraticFunction3d quadfun(p3d, GetNormal());
int starttrig = 1;
if (inner_searchtree)
{
NgArray<STLTrigId> trigs;
starttrig = GetNChartT()+1;
Point<2> p2d = Project2d (p3d);
inner_searchtree->GetIntersecting(p2d, p2d, trigs);
for (STLTrigId i : trigs)
{
auto & trig = geometry->GetTriangle(i);
const Point<3> & c = trig.center;
if (quadfun.Eval(c) > sqr (trig.rad))
continue;
Point<3> p = p3d;
Vec<3> lam;
int err = trig.ProjectInPlain(geometry->GetPoints(), GetNormal(), p, lam);
bool inside = (err == 0 && lam(0) > -lamtol &&
lam(1) > -lamtol && (1-lam(0)-lam(1)) > -lamtol);
if (inside)
{
p3d = p;
return i;
}
}
}
for (int j = starttrig; j <= nt; j++)
{
STLTrigId i = GetTrig1(j);
auto & trig = geometry->GetTriangle(i);
const Point<3> & c = trig.center;
if (quadfun.Eval(c) > sqr (trig.rad))
continue;
Point<3> p = p3d;
Vec<3> lam;
int err = trig.ProjectInPlain(geometry->GetPoints(), GetNormal(), p, lam);
bool inside = (err == 0 && lam(0) > -lamtol &&
lam(1) > -lamtol && (1-lam(0)-lam(1)) > -lamtol);
if (inside)
{
p3d = p;
return i;
}
}
return 0;
}
/*
Point<2> STLChart :: Project2d (const Point<3> & p3d) const
{
Vec<3> v = p3d-pref;
return Point<2> (t1 * v, t2 * v);
}
*/
/*
Point3d p1, p2, center;
double rad;
int i1, i2;
public:
*/
/*
STLBoundarySeg ::
STLBoundarySeg (int ai1, int ai2, const NgArray<Point<3> > & points,
const STLChart * chart)
{
i1 = ai1;
i2 = ai2;
p1 = points.Get(i1);
p2 = points.Get(i2);
center = ::netgen::Center (p1, p2);
rad = Dist (p1, center);
p2d1 = chart->Project2d (p1);
p2d2 = chart->Project2d (p2);
boundingbox.Set (p2d1);
boundingbox.Add (p2d2);
}
*/
void STLBoundarySeg :: Swap ()
{
::netgen::Swap (i1, i2);
::netgen::Swap (p1, p2);
}
STLBoundary :: STLBoundary (STLGeometry * ageometry)
: geometry(ageometry)
{ ; }
/*
void STLBoundary :: AddOrDelSegment(const STLBoundarySeg & seg)
{
bool found = false;
for (int i = 1; i <= boundary.Size(); i++)
{
if (found) { boundary.Elem(i-1) = boundary.Get(i); }
if (boundary.Get(i) == seg) { found = true; }
}
if (!found)
{
boundary.Append(seg);
}
else
{
boundary.SetSize(boundary.Size()-1);
}
}
*/
void STLBoundary ::AddTriangle(const STLTriangle & t)
{
// static Timer timer("STLBoundary::AddTriangle"); RegionTimer reg(timer);
// static int timer_old = NgProfiler::CreateTimer ("STLChart::AddTriangle_old");
// static int timer_new = NgProfiler::CreateTimer ("STLChart::AddTriangle_new");
// NgProfiler::StartTimer (timer_old);
#ifdef ADDTRIGOLD
int i;
int found1 = 0;
int found2 = 0;
int found3 = 0;
//int offset = 0;
STLBoundarySeg seg1(t[0],t[1], geometry->GetPoints(), chart);
STLBoundarySeg seg2(t[1],t[2], geometry->GetPoints(), chart);
STLBoundarySeg seg3(t[2],t[0], geometry->GetPoints(), chart);
seg1.SetSmoothEdge (geometry->IsSmoothEdge (seg1.I1(), seg1.I2()));
seg2.SetSmoothEdge (geometry->IsSmoothEdge (seg2.I1(), seg2.I2()));
seg3.SetSmoothEdge (geometry->IsSmoothEdge (seg3.I1(), seg3.I2()));
/*
for (i = 1; i <= boundary.Size(); i++)
{
if (offset) {boundary.Elem(i-offset) = boundary.Get(i);}
if (boundary.Get(i) == seg1) {found1 = 1; offset++;}
if (boundary.Get(i) == seg2) {found2 = 1; offset++;}
if (boundary.Get(i) == seg3) {found3 = 1; offset++;}
}
if (offset)
{
boundary.SetSize(boundary.Size()-offset);
}
*/
for (i = boundary.Size(); i >= 1; i--)
{
if (boundary.Get(i) == seg1)
{ boundary.DeleteElement (i); found1 = 1; }
else if (boundary.Get(i) == seg2)
{ boundary.DeleteElement (i); found2 = 1; }
else if (boundary.Get(i) == seg3)
{ boundary.DeleteElement (i); found3 = 1; }
}
if (!found1)
{
seg1.Swap();
boundary.Append(seg1);
/*
int newnr;
if (freelist.Size())
{
newnr = freelist.Last();
freelist.DeleteLast();
boundary[newnr] = seg1;
}
else
{
boundary.Append(seg1);
newnr = boundary.Size();
}
// cout << "tree add el " << boundary.Size() << endl;
if (searchtree)
{
// cout << "add " << boundary.Size() << endl;
searchtree->Insert (seg1.BoundingBox(), newnr);
}
*/
}
if (!found2)
{
seg2.Swap();
boundary.Append(seg2);
/*
int newnr;
if (freelist.Size())
{
newnr = freelist.Last();
freelist.DeleteLast();
boundary[newnr] = seg2;
}
else
{
boundary.Append(seg2);
newnr = boundary.Size();
}
// boundary.Append(seg2);
// cout << "tree add el " << boundary.Size() << endl;
if (searchtree)
{
// cout << "add " << boundary.Size() << endl;
searchtree->Insert (seg2.BoundingBox(), newnr);
}
*/
}
if (!found3)
{
seg3.Swap();
boundary.Append(seg3);
/*
int newnr;
if (freelist.Size())
{
newnr = freelist.Last();
freelist.DeleteLast();
boundary[newnr] = seg3;
}
else
{
boundary.Append(seg3);
newnr = boundary.Size();
}
// cout << "tree add el " << boundary.Size() << endl;
if (searchtree)
{
// cout << "add " << boundary.Size() << endl;
searchtree->Insert (seg3.BoundingBox(), newnr);
}
*/
}
#endif
// NgProfiler::StopTimer (timer_old);
// NgProfiler::StartTimer (timer_new);
INDEX_2 segs[3];
segs[0] = INDEX_2(t[0], t[1]);
segs[1] = INDEX_2(t[1], t[2]);
segs[2] = INDEX_2(t[2], t[0]);
if(!searchtree)
BuildSearchTree();
for (auto seg : segs)
{
STLBoundarySeg bseg(seg[0], seg[1], geometry->GetPoints(), chart);
bseg.SetSmoothEdge (geometry->IsSmoothEdge (seg[0],seg[1]));
INDEX_2 op(seg[1], seg[0]);
if (boundary_ht.Used(op))
{
boundary_ht.Delete(op);
if (searchtree)
searchtree->DeleteElement(op);
}
else
{
boundary_ht[seg] = bseg;
if (searchtree)
searchtree->Insert (bseg.BoundingBox(), seg);
}
}
}
bool STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3> & sn,
double sinchartangle, int divisions, Array<Point<3>,STLPointId>& points, double eps)
{
if (usechartnormal)
return TestSegChartNV (p1, p2, sn);
#ifdef NONE
// for statistics
{
int i;
static NgArray<int> cntclass;
static int cnt = 0;
static int cnti = 0, cnto = 0;
static long int cntsegs = 0;
if (cntclass.Size() == 0)
{
cntclass.SetSize (20);
for (i = 1; i <= cntclass.Size(); i++)
cntclass.Elem(i) = 0;
}
cntsegs += NOSegments();
int cla = int (log (double(NOSegments()+1)) / log(2.0));
if (cla < 1) cla = 1;
if (cla > cntclass.Size()) cla = cntclass.Size();
cntclass.Elem(cla)++;
cnt++;
if (divisions)
cnti++;
else
cnto++;
if (cnt > 100000)
{
cnt = 0;
/*
(*testout) << "TestSeg-calls for classes:" << endl;
(*testout) << cnti << " inner calls, " << cnto << " outercalls" << endl;
(*testout) << "total testes segments: " << cntsegs << endl;
for (i = 1; i <= cntclass.Size(); i++)
{
(*testout) << int (exp (i * log(2.0))) << " bnd segs: " << cntclass.Get(i) << endl;
}
*/
}
}
#endif
int i,j,k;
Point<3> seg1p/*, seg2p*/;
Point<3> sp1,sp2;
double lambda1, lambda2, vlen2;
Vec<3> vptpl;
double sinchartangle2 = sqr(sinchartangle);
double scal;
int possible;
//double maxval = -1;
//double maxvalnew = -1;
double scalp1 = p1(0) * sn(0) + p1(1) * sn(1) + p1(2) * sn(2);
double scalp2 = p2(0) * sn(0) + p2(1) * sn(1) + p2(2) * sn(2);
double minl = min2(scalp1, scalp2);
double maxl = max2(scalp1, scalp2);
Point<3> c = Center (p1, p2);
double dist1 = Dist (c, p1);
/*
int nseg = NOSegments();
for (j = 1; j <= nseg; j++)
{
const STLBoundarySeg & seg = GetSegment(j);
*/
for(auto [i2, seg] : boundary_ht)
{
if (seg.IsSmoothEdge())
continue;
sp1 = seg.P1();
sp2 = seg.P2();
// Test, ob Spiral Konfikt moeglich
possible = 1;
double scalsp1 = sp1(0) * sn(0) + sp1(1) * sn(1) + sp1(2) * sn(2);
double scalsp2 = sp2(0) * sn(0) + sp2(1) * sn(1) + sp2(2) * sn(2);
double minsl = min2(scalsp1, scalsp2);
double maxsl = max2(scalsp1, scalsp2);
double maxdiff = max2 (maxsl - minl, maxl - minsl);
/*
Point3d sc = Center (sp1, sp2);
double mindist = Dist(c, sc) - dist1 - GetSegment(j).Radius();
if (maxdiff < sinchartangle * mindist)
{
possible = 0;
}
*/
double hscal = maxdiff + sinchartangle * (dist1 + seg.Radius());
if (hscal * hscal < sinchartangle * Dist2(c, seg.center ))
possible = 0;
/*
if (possible)
{
double mindist2ex = MinDistLL2 (p1, p2, sp1, sp2);
if (maxdiff * maxdiff < sinchartangle2 * mindist2ex)
possible = 0;
}
*/
if (possible)
{
LinearPolynomial2V lp (scalp1 - scalsp1,
scalp2 - scalp1,
-(scalsp2 - scalsp1));
QuadraticPolynomial2V slp;
slp.Square (lp);
Vec3d v (p1, sp1);
Vec3d vl (p1, p2);
Vec3d vsl (sp1, sp2);
QuadraticPolynomial2V qp (v.Length2(),
-2 * (v * vl),
2 * (v * vsl),
vl.Length2(),
-2 * (vl * vsl),
vsl.Length2());
slp.Add (-sinchartangle2, qp);
double hv = slp.MaxUnitSquare();
if (hv > eps) return 0;
/*
if (hv > maxvalnew)
maxvalnew = hv;
*/
}
// if (possible && 0)
if (false)
for (i = 0; i <= divisions; i++)
{
lambda1 = (double)i/(double)divisions;
seg1p = Point3d(p1(0)*lambda1+p2(0)*(1.-lambda1),
p1(1)*lambda1+p2(1)*(1.-lambda1),
p1(2)*lambda1+p2(2)*(1.-lambda1));
for (k = 0; k <= divisions; k++)
{
lambda2 = (double)k/(double)divisions;
vptpl = Vec3d(sp1(0)*lambda2+sp2(0)*(1.-lambda2)-seg1p(0),
sp1(1)*lambda2+sp2(1)*(1.-lambda2)-seg1p(1),
sp1(2)*lambda2+sp2(2)*(1.-lambda2)-seg1p(2));
vlen2 = vptpl.Length2();
// if (vlen2 > 0)
{
scal = vptpl * sn;
double hv = scal*scal - sinchartangle2*vlen2;
/*
if (hv > maxval)
maxval = hv;
*/
if (hv > eps) return 0;
}
}
}
}
return 1;
// return (maxvalnew < eps);
}
void STLBoundary :: BuildSearchTree()
{
Box<2> box2d(Box<2>::EMPTY_BOX);
Box<3> box3d = geometry->GetBoundingBox();
for (size_t i = 0; i < 8; i++)
box2d.Add ( chart->Project2d (box3d.GetPointNr(i)));
searchtree = make_unique<BoxTree<2,INDEX_2>> (box2d);
// searchtree = nullptr;
}
void STLBoundary :: DeleteSearchTree()
{
searchtree = nullptr;
}
// checks, whether 2d projection intersects
bool STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
const Vec3d& sn)
{
// static int timerquick = NgProfiler::CreateTimer ("TestSegChartNV-searchtree");
// static Timer timer("TestSegChartNV"); RegionTimer reg(timer);
Point<2> p2d1 = chart->Project2d (p1);
Point<2> p2d2 = chart->Project2d (p2);
Box<2> box2d;
box2d.Set (p2d1);
box2d.Add (p2d2);
Line2d l1 (p2d1, p2d2);
double eps = 1e-6;
auto hasIntersection = [&] (auto i2) NETGEN_LAMBDA_INLINE
{
const STLBoundarySeg & seg = boundary_ht[i2];
if (seg.IsSmoothEdge()) return false;
if (!box2d.Intersect (seg.BoundingBox())) return false;
const Point<2> & sp1 = seg.P2D1();
const Point<2> & sp2 = seg.P2D2();
Line2d l2 (sp1, sp2);
double lam1, lam2;
int err = CrossPointBarycentric (l1, l2, lam1, lam2);
bool in1 = (lam1 > eps) && (lam1 < 1-eps);
bool on1 = (lam1 > -eps) && (lam1 < 1 + eps);
bool in2 = (lam2 > eps) && (lam2 < 1-eps);
bool on2 = (lam2 > -eps) && (lam2 < 1 + eps);
if(!err && ((on1 && in2) || (on2 && in1)))
return true;
return false;
};
if (searchtree)
{
bool has_intersection = false;
searchtree -> GetFirstIntersecting (box2d.PMin(), box2d.PMax(),
[&] (auto i2) NETGEN_LAMBDA_INLINE
{
has_intersection = hasIntersection(i2);
return has_intersection;
});
return !has_intersection;
}
else
{
for(auto [i2, seg] : boundary_ht)
if(hasIntersection(i2))
return false;
return true;
}
}
STLDoctorParams :: STLDoctorParams()
{
drawmeshededges = 1;
geom_tol_fact = 1E-6;
longlinefact = 0;
showexcluded = 1;
selectmode = 0;
edgeselectmode = 0;
useexternaledges = 0;
showfaces = 0;
showtouchedtrigchart = 1;
showedgecornerpoints = 1;
conecheck = 1;
spiralcheck = 1;
selecttrig = 0;
nodeofseltrig = 1;
selectwithmouse = 1;
showmarkedtrigs = 1;
dirtytrigfact = 0.001;
smoothangle = 90;
smoothnormalsweight = 0.2;
vicinity = 0;
showvicinity = 0;
}
STLDoctorParams stldoctor;
void STLDoctorParams :: Print (ostream & ost) const
{
ost << "STL doctor parameters:" << endl
<< "selecttrig = " << selecttrig << endl
<< "selectlocalpoint = " << nodeofseltrig << endl
<< "selectwithmouse = " << selectwithmouse << endl
<< "showmarkedtrigs = " << showmarkedtrigs << endl
<< "dirtytrigfact = " << dirtytrigfact << endl
<< "smoothangle = " << smoothangle << endl;
}
STLParameters :: STLParameters()
{
yangle = 30;
contyangle = 20;
edgecornerangle = 60;
chartangle = 15;
outerchartangle = 70;
usesearchtree = 0;
atlasminh = 1E-4;
resthsurfcurvfac = 2;
resthsurfcurvenable = 0;
resthatlasfac = 2;
resthatlasenable = 1;
resthchartdistfac = 1.2;
resthchartdistenable = 1;
resthlinelengthfac = 0.5;
resthlinelengthenable = 1;
resthcloseedgefac = 1;
resthcloseedgeenable = 1;
resthedgeanglefac = 1;
resthedgeangleenable = 0;
resthsurfmeshcurvfac = 1;
resthsurfmeshcurvenable = 0;
recalc_h_opt = 1;
}
void STLParameters :: Print (ostream & ost) const
{
ost << "STL parameters:" << endl
<< "yellow angle = " << yangle << endl
<< "continued yellow angle = " << contyangle << endl
<< "edgecornerangle = " << edgecornerangle << endl
<< "chartangle = " << chartangle << endl
<< "outerchartangle = " << outerchartangle << endl
<< "restrict h due to ..., enable and safety factor: " << endl
<< "surface curvature: " << resthsurfcurvenable
<< ", fac = " << resthsurfcurvfac << endl
<< "atlas surface curvature: " << resthatlasenable
<< ", fac = " << resthatlasfac << endl
<< "chart distance: " << resthchartdistenable
<< ", fac = " << resthchartdistfac << endl
<< "line length: " << resthlinelengthenable
<< ", fac = " << resthlinelengthfac << endl
<< "close edges: " << resthcloseedgeenable
<< ", fac = " << resthcloseedgefac << endl
<< "edge angle: " << resthedgeangleenable
<< ", fac = " << resthedgeanglefac << endl;
}
DLL_HEADER extern STLParameters stlparam;
STLParameters stlparam;
}