1 #include <meshing.hpp>
2 #include <core/register_archive.hpp>
3
4 #include "stlgeom.hpp"
5
6 namespace netgen
7 {
8
9 //globalen searchtree fuer gesamte geometry aktivieren
10 int geomsearchtreeon = 0;
11
12 int usechartnormal = 1;
13
14 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
STLMeshing(STLGeometry & geom,Mesh & mesh,const MeshingParameters & mparam,const STLParameters & stlpar)16 void STLMeshing (STLGeometry & geom,
17 Mesh & mesh,
18 const MeshingParameters& mparam,
19 const STLParameters& stlpar)
20 {
21 geom.Clear();
22 geom.BuildEdges(stlpar);
23 geom.MakeAtlas(mesh, mparam, stlpar);
24 if (multithread.terminate) { return; }
25 geom.CalcFaceNums();
26 geom.AddFaceEdges();
27 geom.LinkEdges(stlpar);
28
29 mesh.ClearFaceDescriptors();
30 for (int i = 1; i <= geom.GetNOFaces(); i++)
31 mesh.AddFaceDescriptor (FaceDescriptor (i, 1, 0, 0));
32 }
33
34 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
35 //+++++++++++++++++++ STL GEOMETRY ++++++++++++++++++++++++++++
36 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
37
38
STLGeometry()39 STLGeometry :: STLGeometry()
40 /*
41 : edges(), edgesperpoint(),
42 normals(), externaledges(),
43 atlas(), chartmark(),
44 lines(), outerchartspertrig(), vicinity(), markedtrigs(), markedsegs(),
45 lineendpoints(), spiralpoints(), selectedmultiedge()
46 */
47 {
48 edgedata = make_unique<STLEdgeDataList>(*this);
49 externaledges.SetSize(0);
50 Clear();
51 meshchart = 0; // initialize all ?? JS
52
53 if (geomsearchtreeon)
54 searchtree = new BoxTree<3> (GetBoundingBox().PMin() - Vec3d(1,1,1),
55 GetBoundingBox().PMax() + Vec3d(1,1,1));
56 else
57 searchtree = NULL;
58
59 status = STL_GOOD;
60 statustext = "Good Geometry";
61 smoothedges = NULL;
62 area = -1;
63 }
64
~STLGeometry()65 STLGeometry :: ~STLGeometry()
66 {
67 // for (auto p : atlas) delete p;
68 // delete edgedata;
69 }
70
Save(string filename) const71 void STLGeometry :: Save (string filename) const
72 {
73 const char * cfilename = filename.c_str();
74 if (strlen(cfilename) < 4)
75 throw NgException ("illegal filename");
76
77 if (strlen(cfilename) > 3 &&
78 strcmp (&cfilename[strlen(cfilename)-3], "stl") == 0)
79 {
80 STLTopology::Save (cfilename);
81 }
82 else if (strlen(cfilename) > 4 &&
83 strcmp (&cfilename[strlen(cfilename)-4], "stlb") == 0)
84 {
85 SaveBinary (cfilename,"Binary STL Geometry");
86
87 }
88 else if (strlen(cfilename) > 4 &&
89 strcmp (&cfilename[strlen(cfilename)-4], "stle") == 0)
90 {
91 SaveSTLE (cfilename);
92 }
93 }
94
95
96
97 DLL_HEADER extern STLParameters stlparam;
GenerateMesh(shared_ptr<Mesh> & mesh,MeshingParameters & mparam)98 int STLGeometry :: GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
99 {
100 STLParameters stlpar = stlparam;
101 return STLMeshingDummy (this, mesh, mparam, stlpar);
102 }
103
GetNormal(int surfind,const Point<3> & p,const PointGeomInfo * gi) const104 Vec<3> STLGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const
105 {
106 if(!gi)
107 throw Exception("STLGeometry::GetNormal without PointGeomInfo called");
108 return GetChart(GetChartNr(gi->trignum)).GetNormal();
109 }
110
CalcPointGeomInfo(int,PointGeomInfo & gi,const Point<3> & p3) const111 bool STLGeometry :: CalcPointGeomInfo(int /*surfind*/, PointGeomInfo& gi, const Point<3> & p3) const
112 {
113 Point<3> hp = p3;
114 SelectChartOfTriangle(gi.trignum);
115
116 gi.trignum = Project (hp);
117
118 if (gi.trignum) return true;
119
120 return false;
121 }
122
ProjectPointGI(int surfind,Point<3> & p,PointGeomInfo & gi) const123 bool STLGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
124 {
125 static std::mutex mutex_project_whole_surface;
126 int meshchart = GetChartNr(gi.trignum);
127 const STLChart& chart = GetChart(meshchart);
128 int trignum = chart.ProjectNormal(p);
129 if(trignum==0)
130 {
131 // non-thread-safe implementation
132 std::lock_guard<std::mutex> guard(mutex_project_whole_surface);
133 PrintMessage(7,"project failed");
134 SelectChartOfTriangle (gi.trignum); // needed because ProjectOnWholeSurface uses meshchartnv (the normal vector of selected chart)
135 trignum = ProjectOnWholeSurface(p);
136 if(trignum==0)
137 {
138 PrintMessage(7, "project on whole surface failed");
139 return false;
140 }
141 }
142 return true;
143 }
144
ProjectPoint(INDEX surfind,Point<3> & p) const145 PointGeomInfo STLGeometry :: ProjectPoint (INDEX surfind, Point<3> & p) const
146 {
147 throw Exception("ProjectPoint without PointGeomInfo not implemented");
148 }
149
150 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) const151 PointBetween (const Point<3> & p1, const Point<3> & p2, double secpoint,
152 int surfi,
153 const PointGeomInfo & gi1,
154 const PointGeomInfo & gi2,
155 Point<3> & newp, PointGeomInfo & newgi) const
156 {
157 newp = p1+secpoint*(p2-p1);
158
159 /*
160 (*testout) << "surf-between: p1 = " << p1 << ", p2 = " << p2
161 << ", gi = " << gi1 << " - " << gi2 << endl;
162 */
163
164 if (gi1.trignum > 0)
165 {
166 // ((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
167
168 Point<3> np1 = newp;
169 Point<3> np2 = newp;
170 auto ngi1 = gi1;
171 auto ngi2 = gi2;
172 // SelectChartOfTriangle (gi1.trignum);
173 int tn1 = ProjectPointGI (surfi, np1, ngi1);
174
175 // SelectChartOfTriangle (gi2.trignum);
176 int tn2 = ProjectPointGI (surfi, np2, ngi2);
177
178 newgi.trignum = tn1; //urspruengliche version
179 newp = np1; //urspruengliche version
180
181 if (!newgi.trignum)
182 { newgi.trignum = tn2; newp = np2; }
183 if (!newgi.trignum) newgi.trignum = gi1.trignum;
184 }
185 else
186 {
187 // (*testout) << "WARNING: PointBetween got geominfo = 0" << endl;
188 newp = p1+secpoint*(p2-p1);
189 newgi.trignum = 0;
190 }
191 }
192
193 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) const194 PointBetweenEdge (const Point<3> & p1, const Point<3> & p2, double secpoint,
195 int surfi1, int surfi2,
196 const EdgePointGeomInfo & gi1,
197 const EdgePointGeomInfo & gi2,
198 Point<3> & newp, EdgePointGeomInfo & newgi) const
199 {
200 /*
201 (*testout) << "edge-between: p1 = " << p1 << ", p2 = " << p2
202 << ", gi1,2 = " << gi1 << ", " << gi2 << endl;
203 */
204 /*
205 newp = Center (p1, p2);
206 ((STLGeometry&)geom).SelectChartOfTriangle (gi1.trignum);
207 newgi.trignum = geom.Project (newp);
208 */
209 int hi;
210 newgi.dist = (1.0-secpoint) * gi1.dist + secpoint*gi2.dist;
211 newgi.edgenr = gi1.edgenr;
212
213 /*
214 (*testout) << "p1 = " << p1 << ", p2 = " << p2 << endl;
215 (*testout) << "refedge: " << gi1.edgenr
216 << " d1 = " << gi1.dist << ", d2 = " << gi2.dist << endl;
217 */
218 newp = GetLine (gi1.edgenr)->GetPointInDist (GetPoints(), newgi.dist, hi);
219
220 // (*testout) << "newp = " << newp << endl;
221 }
222
STLInfo(double * data)223 void STLGeometry :: STLInfo(double* data)
224 {
225 data[0] = GetNT();
226
227 Box<3> b = GetBoundingBox();
228 data[1] = b.PMin()(0);
229 data[2] = b.PMax()(0);
230 data[3] = b.PMin()(1);
231 data[4] = b.PMax()(1);
232 data[5] = b.PMin()(2);
233 data[6] = b.PMax()(2);
234
235 int i;
236
237 int cons = 1;
238 for (i = 1; i <= GetNT(); i++)
239 {
240 if (NONeighbourTrigs(i) != 3) {cons = 0;}
241 }
242 data[7] = cons;
243 }
244
MarkNonSmoothNormals(const STLParameters & stlparam)245 void STLGeometry :: MarkNonSmoothNormals(const STLParameters& stlparam)
246 {
247
248 PrintFnStart("Mark Non-Smooth Normals");
249
250 int i,j;
251
252 markedtrigs.SetSize(GetNT());
253
254 for (i = 1; i <= GetNT(); i++)
255 {
256 SetMarkedTrig(i, 0);
257 }
258
259 double dirtyangle = stlparam.yangle/180.*M_PI;
260
261 int cnt = 0;
262 STLPointId lp1,lp2;
263 for (i = 1; i <= GetNT(); i++)
264 {
265 for (j = 1; j <= NONeighbourTrigs(i); j++)
266 {
267 if (GetAngle(i, NeighbourTrig(i,j)) > dirtyangle)
268 {
269 GetTriangle(i).GetNeighbourPoints(GetTriangle(NeighbourTrig(i,j)), lp1, lp2);
270 if (!IsEdge(lp1,lp2))
271 {
272 if (!IsMarkedTrig(i)) {SetMarkedTrig(i,1); cnt++;}
273 }
274 }
275 }
276 }
277
278 PrintMessage(5,"marked ",cnt," non-smooth trig-normals");
279
280 }
281
SmoothNormals(const STLParameters & stlparam)282 void STLGeometry :: SmoothNormals(const STLParameters& stlparam)
283 {
284 multithread.terminate = 0;
285
286 // UseExternalEdges();
287
288 BuildEdges(stlparam);
289
290
291 DenseMatrix m(3), hm(3);
292 Vector rhs(3), sol(3), hv(3), hv2(3);
293
294 Vec<3> ri;
295
296 double wnb = stldoctor.smoothnormalsweight; // neighbour normal weight
297 double wgeom = 1-wnb; // geometry normal weight
298
299
300 // minimize
301 // wgeom sum_T \sum ri \| ri^T (n - n_geom) \|^2
302 // + wnb sum_SE \| ri x (n - n_nb) \|^2
303
304 int i, j, k, l;
305 int nt = GetNT();
306
307 PushStatusF("Smooth Normals");
308
309 //int testmode;
310
311 for (i = 1; i <= nt; i++)
312 {
313
314 SetThreadPercent( 100.0 * (double)i / (double)nt);
315
316 const STLTriangle & trig = GetTriangle (i);
317
318 m = 0;
319 rhs = 0;
320
321 // normal of geometry:
322 Vec<3> ngeom = trig.GeomNormal(points);
323 ngeom.Normalize();
324
325 for (j = 1; j <= 3; j++)
326 {
327 int pi1 = trig.PNumMod (j);
328 int pi2 = trig.PNumMod (j+1);
329
330 // edge vector
331 ri = GetPoint (pi2) - GetPoint (pi1);
332
333 for (k = 0; k < 3; k++)
334 for (l = 0; l < 3; l++)
335 hm.Elem(k+1, l+1) = wgeom * ri(k) * ri(l);
336
337
338 for (k = 0; k < 3; k++)
339 hv(k) = ngeom(k);
340
341 hm.Mult (hv, hv2);
342 /*
343 if (testmode)
344 (*testout) << "add vec " << hv2 << endl
345 << " add m " << hm << endl;
346 */
347 rhs.Add (1, hv2);
348 m += hm;
349
350
351 int nbt = 0;
352 STLPointId fp1,fp2;
353 for (k = 1; k <= NONeighbourTrigs(i); k++)
354 {
355 trig.GetNeighbourPoints(GetTriangle(NeighbourTrig(i, k)),fp1,fp2);
356 if (fp1 == pi1 && fp2 == pi2)
357 {
358 nbt = NeighbourTrig(i, k);
359 }
360 }
361
362 if (!nbt)
363 {
364 cerr << "ERROR: stlgeom::Smoothnormals, nbt = 0" << endl;
365 }
366
367 // smoothed normal
368 Vec<3> nnb = GetTriangle(nbt).Normal(); // neighbour normal
369 nnb.Normalize();
370
371 if (!IsEdge(pi1,pi2))
372 {
373 double lr2 = ri * ri;
374 for (k = 0; k < 3; k++)
375 {
376 for (l = 0; l < k; l++)
377 {
378 hm.Elem(k+1, l+1) = -wnb * ri(k) * ri(l);
379 hm.Elem(l+1, k+1) = -wnb * ri(k) * ri(l);
380 }
381
382 hm.Elem(k+1, k+1) = wnb * (lr2 - ri(k) * ri(k));
383 }
384
385 for (k = 0; k < 3; k++)
386 hv(k) = nnb(k);
387
388 hm.Mult (hv, hv2);
389 /*
390 if (testmode)
391 (*testout) << "add nb vec " << hv2 << endl
392 << " add nb m " << hm << endl;
393 */
394
395 rhs.Add (1, hv2);
396 m += hm;
397 }
398 }
399
400 m.Solve (rhs, sol);
401 Vec3d newn(sol(0), sol(1), sol(2));
402 newn /= (newn.Length() + 1e-24);
403
404 GetTriangle(i).SetNormal(newn);
405 // setnormal (sol);
406 }
407
408 /*
409 for (i = 1; i <= nt; i++)
410 SetMarkedTrig(i, 0);
411
412
413
414 int crloop;
415 for (crloop = 1; crloop <= 3; crloop++)
416 {
417
418 // find critical:
419
420 NgArray<INDEX_2> critpairs;
421 for (i = 1; i <= nt; i++)
422 {
423 const STLTriangle & trig = GetTriangle (i);
424
425 Vec3d ngeom = GetTriangleNormal (i); // trig.Normal(points);
426 ngeom /= (ngeom.Length() + 1e-24);
427
428 for (j = 1; j <= 3; j++)
429 {
430 int pi1 = trig.PNumMod (j);
431 int pi2 = trig.PNumMod (j+1);
432
433 int nbt = 0;
434 int fp1,fp2;
435 for (k = 1; k <= NONeighbourTrigs(i); k++)
436 {
437 trig.GetNeighbourPoints(GetTriangle(NeighbourTrig(i, k)),fp1,fp2);
438 if (fp1 == pi1 && fp2 == pi2)
439 {
440 nbt = NeighbourTrig(i, k);
441 }
442 }
443
444 if (!nbt)
445 {
446 cerr << "ERROR: stlgeom::Smoothnormals, nbt = 0" << endl;
447 }
448
449 Vec3d nnb = GetTriangleNormal(nbt); // neighbour normal
450 nnb /= (nnb.Length() + 1e-24);
451
452 if (!IsEdge(pi1,pi2))
453 {
454 if (Angle (nnb, ngeom) > 150 * M_PI/180)
455 {
456 SetMarkedTrig(i, 1);
457 SetMarkedTrig(nbt, 1);
458 critpairs.Append (INDEX_2 (i, nbt));
459 }
460 }
461
462 }
463 }
464
465 if (!critpairs.Size())
466 {
467 break;
468 }
469
470 if (critpairs.Size())
471 {
472
473 NgArray<int> friends;
474 double area1 = 0, area2 = 0;
475
476 for (i = 1; i <= critpairs.Size(); i++)
477 {
478 int tnr1 = critpairs.Get(i).I1();
479 int tnr2 = critpairs.Get(i).I2();
480 (*testout) << "t1 = " << tnr1 << ", t2 = " << tnr2
481 << " angle = " << Angle (GetTriangleNormal (tnr1),
482 GetTriangleNormal (tnr2))
483 << endl;
484
485 // who has more friends ?
486 int side;
487 area1 = 0;
488 area2 = 0;
489 for (side = 1; side <= 2; side++)
490 {
491 friends.SetSize (0);
492 friends.Append ( (side == 1) ? tnr1 : tnr2);
493
494 for (j = 1; j <= 3; j++)
495 {
496 int fsize = friends.Size();
497 for (k = 1; k <= fsize; k++)
498 {
499 int testtnr = friends.Get(k);
500 Vec3d ntt = GetTriangleNormal(testtnr);
501 ntt /= (ntt.Length() + 1e-24);
502
503 for (l = 1; l <= NONeighbourTrigs(testtnr); l++)
504 {
505 int testnbnr = NeighbourTrig(testtnr, l);
506 Vec3d nbt = GetTriangleNormal(testnbnr);
507 nbt /= (nbt.Length() + 1e-24);
508
509 if (Angle (nbt, ntt) < 15 * M_PI/180)
510 {
511 int ii;
512 int found = 0;
513 for (ii = 1; ii <= friends.Size(); ii++)
514 {
515 if (friends.Get(ii) == testnbnr)
516 {
517 found = 1;
518 break;
519 }
520 }
521 if (!found)
522 friends.Append (testnbnr);
523 }
524 }
525 }
526 }
527
528 // compute area:
529 for (k = 1; k <= friends.Size(); k++)
530 {
531 double area =
532 GetTriangle (friends.Get(k)).Area(points);
533
534 if (side == 1)
535 area1 += area;
536 else
537 area2 += area;
538 }
539
540 }
541
542 (*testout) << "area1 = " << area1 << " area2 = " << area2 << endl;
543 if (area1 < 0.1 * area2)
544 {
545 Vec3d n = GetTriangleNormal (tnr1);
546 n *= -1;
547 SetTriangleNormal(tnr1, n);
548 }
549 if (area2 < 0.1 * area1)
550 {
551 Vec3d n = GetTriangleNormal (tnr2);
552 n *= -1;
553 SetTriangleNormal(tnr2, n);
554 }
555 }
556 }
557 }
558 */
559
560 calcedgedataanglesnew = 1;
561 PopStatus();
562 }
563
564
AddEdge(int ap1,int ap2)565 int STLGeometry :: AddEdge(int ap1, int ap2)
566 {
567 STLEdge e(ap1,ap2);
568 e.SetLeftTrig(GetLeftTrig(ap1,ap2));
569 e.SetRightTrig(GetRightTrig(ap1,ap2));
570 edges.Append(e);
571 return edges.Size();
572 }
573
STLDoctorConfirmEdge()574 void STLGeometry :: STLDoctorConfirmEdge()
575 {
576 StoreEdgeData();
577 if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT() && GetNodeOfSelTrig())
578 {
579 if (stldoctor.selectmode == 1)
580 {
581 int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
582 int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
583 edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus (ED_CONFIRMED);
584 }
585 else if (stldoctor.selectmode == 3 || stldoctor.selectmode == 4)
586 {
587 int i;
588 for (i = 1; i <= selectedmultiedge.Size(); i++)
589 {
590 int ap1 = selectedmultiedge.Get(i).i1;
591 int ap2 = selectedmultiedge.Get(i).i2;
592 edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus (ED_CONFIRMED);
593 }
594 }
595 }
596 }
597
STLDoctorCandidateEdge()598 void STLGeometry :: STLDoctorCandidateEdge()
599 {
600 StoreEdgeData();
601 if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT() && GetNodeOfSelTrig())
602 {
603 if (stldoctor.selectmode == 1)
604 {
605 int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
606 int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
607 edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus (ED_CANDIDATE);
608 }
609 else if (stldoctor.selectmode == 3 || stldoctor.selectmode == 4)
610 {
611 int i;
612 for (i = 1; i <= selectedmultiedge.Size(); i++)
613 {
614 int ap1 = selectedmultiedge.Get(i).i1;
615 int ap2 = selectedmultiedge.Get(i).i2;
616 edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus (ED_CANDIDATE);
617 }
618 }
619 }
620 }
621
STLDoctorExcludeEdge()622 void STLGeometry :: STLDoctorExcludeEdge()
623 {
624 StoreEdgeData();
625 if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT() && GetNodeOfSelTrig())
626 {
627 if (stldoctor.selectmode == 1)
628 {
629 int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
630 int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
631 edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus(ED_EXCLUDED);
632 }
633 else if (stldoctor.selectmode == 3 || stldoctor.selectmode == 4)
634 {
635 int i;
636 for (i = 1; i <= selectedmultiedge.Size(); i++)
637 {
638 int ap1 = selectedmultiedge.Get(i).i1;
639 int ap2 = selectedmultiedge.Get(i).i2;
640 edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus(ED_EXCLUDED);
641 }
642 }
643 }
644 }
645
STLDoctorUndefinedEdge()646 void STLGeometry :: STLDoctorUndefinedEdge()
647 {
648 StoreEdgeData();
649 if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT() && GetNodeOfSelTrig())
650 {
651 if (stldoctor.selectmode == 1)
652 {
653 int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
654 int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
655 edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus(ED_UNDEFINED);
656 }
657 else if (stldoctor.selectmode == 3 || stldoctor.selectmode == 4)
658 {
659 int i;
660 for (i = 1; i <= selectedmultiedge.Size(); i++)
661 {
662 int ap1 = selectedmultiedge.Get(i).i1;
663 int ap2 = selectedmultiedge.Get(i).i2;
664 edgedata->Elem(edgedata->GetEdgeNum(ap1,ap2)).SetStatus(ED_UNDEFINED);
665 }
666 }
667 }
668 }
669
STLDoctorSetAllUndefinedEdges()670 void STLGeometry :: STLDoctorSetAllUndefinedEdges()
671 {
672 edgedata->ResetAll();
673 }
674
STLDoctorEraseCandidateEdges()675 void STLGeometry :: STLDoctorEraseCandidateEdges()
676 {
677 StoreEdgeData();
678 edgedata->ChangeStatus(ED_CANDIDATE, ED_UNDEFINED);
679 }
680
STLDoctorConfirmCandidateEdges()681 void STLGeometry :: STLDoctorConfirmCandidateEdges()
682 {
683 StoreEdgeData();
684 edgedata->ChangeStatus(ED_CANDIDATE, ED_CONFIRMED);
685 }
686
STLDoctorConfirmedToCandidateEdges()687 void STLGeometry :: STLDoctorConfirmedToCandidateEdges()
688 {
689 StoreEdgeData();
690 edgedata->ChangeStatus(ED_CONFIRMED, ED_CANDIDATE);
691 }
692
STLDoctorDirtyEdgesToCandidates()693 void STLGeometry :: STLDoctorDirtyEdgesToCandidates()
694 {
695 StoreEdgeData();
696 }
697
STLDoctorLongLinesToCandidates()698 void STLGeometry :: STLDoctorLongLinesToCandidates()
699 {
700 StoreEdgeData();
701 }
702
GetNearestSelectedDefinedEdge()703 twoint STLGeometry :: GetNearestSelectedDefinedEdge()
704 {
705 Point<3> pestimate = Center(GetTriangle(GetSelectTrig()).center,
706 GetPoint(GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig())));
707 //Point3d pestimate = GetTriangle(GetSelectTrig()).center;
708
709 int i, j, en;
710 NgArray<int> vic;
711 GetVicinity(GetSelectTrig(),4,vic);
712
713
714 twoint fedg;
715 fedg.i1 = 0;
716 fedg.i2 = 0;
717 double mindist = 1E50;
718 double dist;
719 Point<3> p;
720
721 for (i = 1; i <= vic.Size(); i++)
722 {
723 const STLTriangle& t = GetTriangle(vic.Get(i));
724 for (j = 1; j <= 3; j++)
725 {
726 en = edgedata->GetEdgeNum(t.PNum(j),t.PNumMod(j+1));
727 if (edgedata->Get(en).GetStatus() != ED_UNDEFINED)
728 {
729 p = pestimate;
730 dist = GetDistFromLine(GetPoint(t.PNum(j)),GetPoint(t.PNumMod(j+1)),p);
731 if (dist < mindist)
732 {
733 mindist = dist;
734 fedg.i1 = t.PNum(j);
735 fedg.i2 = t.PNumMod(j+1);
736 }
737 }
738 }
739 }
740 return fedg;
741 }
742
BuildSelectedMultiEdge(twoint ep)743 void STLGeometry :: BuildSelectedMultiEdge(twoint ep)
744 {
745 if (edgedata->Size() == 0 ||
746 !GetEPPSize())
747 {
748 return;
749 }
750
751 selectedmultiedge.SetSize(0);
752 int tenum = GetTopEdgeNum (ep.i1, ep.i2);
753
754 if (edgedata->Get(tenum).GetStatus() == ED_UNDEFINED)
755 {
756 twoint epnew = GetNearestSelectedDefinedEdge();
757 if (epnew.i1)
758 {
759 ep = epnew;
760 tenum = GetTopEdgeNum (ep.i1, ep.i2);
761 }
762 }
763
764 selectedmultiedge.Append(twoint(ep));
765
766 if (edgedata->Get(tenum).GetStatus() == ED_UNDEFINED)
767 {
768 return;
769 }
770
771 edgedata->BuildLineWithEdge(ep.i1,ep.i2,selectedmultiedge);
772 }
773
BuildSelectedEdge(twoint ep)774 void STLGeometry :: BuildSelectedEdge(twoint ep)
775 {
776 if (edgedata->Size() == 0 ||
777 !GetEPPSize())
778 {
779 return;
780 }
781
782 selectedmultiedge.SetSize(0);
783
784 selectedmultiedge.Append(twoint(ep));
785 }
786
BuildSelectedCluster(twoint ep)787 void STLGeometry :: BuildSelectedCluster(twoint ep)
788 {
789 if (edgedata->Size() == 0 ||
790 !GetEPPSize())
791 {
792 return;
793 }
794
795 selectedmultiedge.SetSize(0);
796
797 int tenum = GetTopEdgeNum (ep.i1, ep.i2);
798
799 if (edgedata->Get(tenum).GetStatus() == ED_UNDEFINED)
800 {
801 twoint epnew = GetNearestSelectedDefinedEdge();
802 if (epnew.i1)
803 {
804 ep = epnew;
805 tenum = GetTopEdgeNum (ep.i1, ep.i2);
806 }
807 }
808
809 selectedmultiedge.Append(twoint(ep));
810
811 if (edgedata->Get(tenum).GetStatus() == ED_UNDEFINED)
812 {
813 return;
814 }
815
816 edgedata->BuildClusterWithEdge(ep.i1,ep.i2,selectedmultiedge);
817 }
818
ImportEdges()819 void STLGeometry :: ImportEdges()
820 {
821 StoreEdgeData();
822
823 PrintMessage(5, "import edges from file 'edges.ng'");
824 ifstream fin("edges.ng");
825
826 int ne;
827 fin >> ne;
828
829 NgArray<Point<3> > eps;
830
831 int i;
832 Point<3> p;
833 for (i = 1; i <= 2*ne; i++)
834 {
835 fin >> p(0);
836 fin >> p(1);
837 fin >> p(2);
838 eps.Append(p);
839 }
840 AddEdges(eps);
841 }
842
AddEdges(const NgArray<Point<3>> & eps)843 void STLGeometry :: AddEdges(const NgArray<Point<3> >& eps)
844 {
845 int i;
846 int ne = eps.Size()/2;
847
848 NgArray<int> epsi;
849 Box<3> bb = GetBoundingBox();
850 bb.Increase(1);
851
852 Point3dTree ptree (bb.PMin(),
853 bb.PMax());
854 NgArray<int> pintersect;
855
856 double gtol = GetBoundingBox().Diam()/1.E10;
857 Point<3> p;
858
859 for (i = 1; i <= GetNP(); i++)
860 {
861 p = GetPoint(i);
862 ptree.Insert (p, i);
863 }
864
865 int error = 0;
866 for (i = 1; i <= 2*ne; i++)
867 {
868 p = eps.Get(i);
869 Point3d pmin = p - Vec3d (gtol, gtol, gtol);
870 Point3d pmax = p + Vec3d (gtol, gtol, gtol);
871
872 ptree.GetIntersecting (pmin, pmax, pintersect);
873 if (pintersect.Size() > 1)
874 {
875 PrintError("Found too much points in epsilon-dist");
876 error = 1;
877 }
878 else if (pintersect.Size() == 0)
879 {
880 error = 1;
881 PrintError("edgepoint does not exist!");
882 PrintMessage(5,"p=",Point3d(eps.Get(i)));
883 }
884 else
885 {
886 epsi.Append(pintersect.Get(1));
887 }
888 }
889
890 if (error) return;
891
892 int en;
893 for (i = 1; i <= ne; i++)
894 {
895 if (epsi.Get(2*i-1) == epsi.Get(2*i)) {PrintError("Edge with zero length!");}
896 else
897 {
898 en = edgedata->GetEdgeNum(epsi.Get(2*i-1),epsi.Get(2*i));
899 edgedata->Elem(en).SetStatus (ED_CONFIRMED);
900 }
901 }
902
903 }
904
905
906
ImportExternalEdges(const char * filename)907 void STLGeometry :: ImportExternalEdges(const char * filename)
908 {
909 //AVL edges!!!!!!
910
911 ifstream inf (filename);
912 char ch;
913 //int cnt = 0;
914 int records, units, i, j;
915 PrintFnStart("Import edges from ",filename);
916
917 const int flen=30;
918 char filter[flen+1];
919 filter[flen] = 0;
920 char buf[20];
921
922 NgArray<Point3d> importpoints;
923 NgArray<int> importlines;
924 NgArray<int> importpnums;
925
926 while (inf.good())
927 {
928 inf.get(ch);
929 // (*testout) << cnt << ": " << ch << endl;
930
931 for (i = 0; i < flen; i++)
932 filter[i] = filter[i+1];
933 filter[flen-1] = ch;
934 // (*testout) << filter << endl;
935
936 if (strcmp (filter+flen-7, "RECORDS") == 0)
937 {
938 inf.get(ch); // '='
939 inf >> records;
940 }
941 if (strcmp (filter+flen-5, "UNITS") == 0)
942 {
943 inf.get(ch); // '='
944 inf >> units;
945 }
946
947 if (strcmp (filter+flen-17, "EDGE NODE NUMBERS") == 0)
948 {
949 int nodenr;
950 importlines.SetSize (units);
951 for (i = 1; i <= units; i++)
952 {
953 inf >> nodenr;
954 importlines.Elem(i) = nodenr;
955 // (*testout) << nodenr << endl;
956 }
957 }
958
959 if (strcmp (filter+flen-23, "EDGE POINT COORD IN DIR") == 0)
960 {
961 int coord;
962
963 inf >> coord;
964
965 importpoints.SetSize (units);
966
967 inf >> ch;
968 inf.putback (ch);
969
970 for (i = 1; i <= units; i++)
971 {
972 for (j = 0; j < 12; j++)
973 inf.get (buf[j]);
974 buf[12] = 0;
975
976 importpoints.Elem(i).X(coord) = 1000 * atof (buf);
977 }
978 }
979 }
980
981 /*
982 (*testout) << "lines: " << endl;
983 for (i = 1; i <= importlines.Size(); i++)
984 (*testout) << importlines.Get(i) << endl;
985 (*testout) << "points: " << endl;
986 for (i = 1; i <= importpoints.Size(); i++)
987 (*testout) << importpoints.Get(i) << endl;
988 */
989
990
991
992 importpnums.SetSize (importpoints.Size());
993
994
995 Box3d bb (GetBoundingBox().PMin() + Vec3d (-1,-1,-1),
996 GetBoundingBox().PMax() + Vec3d (1, 1, 1));
997
998 Point3dTree ptree (bb.PMin(),
999 bb.PMax());
1000
1001
1002 PrintMessage(7,"stl - bb: ",bb.PMin(), " - ", bb.PMax());
1003
1004 Box3d ebb;
1005 ebb.SetPoint (importpoints.Get(1));
1006 for (i = 1; i <= importpoints.Size(); i++)
1007 ebb.AddPoint (importpoints.Get(i));
1008 PrintMessage(7,"edgep - bb: ", ebb.PMin(), " - ", ebb.PMax());
1009
1010 NgArray<int> pintersect;
1011
1012 double gtol = GetBoundingBox().Diam()/1.E6;
1013
1014 for (i = 1; i <= GetNP(); i++)
1015 {
1016 Point3d p = GetPoint(i);
1017 // (*testout) << "stlpt: " << p << endl;
1018 ptree.Insert (p, i);
1019 }
1020
1021
1022 for (i = 1; i <= importpoints.Size(); i++)
1023 {
1024 Point3d p = importpoints.Get(i);
1025 Point3d pmin = p - Vec3d (gtol, gtol, gtol);
1026 Point3d pmax = p + Vec3d (gtol, gtol, gtol);
1027
1028 ptree.GetIntersecting (pmin, pmax, pintersect);
1029 if (pintersect.Size() > 1)
1030 {
1031 importpnums.Elem(i) = 0;
1032 PrintError("Found too many points in epsilon-dist");
1033 }
1034 else if (pintersect.Size() == 0)
1035 {
1036 importpnums.Elem(i) = 0;
1037 PrintError("Edgepoint does not exist!");
1038 }
1039 else
1040 {
1041 importpnums.Elem(i) = pintersect.Get(1);
1042 }
1043 }
1044
1045 // if (!error)
1046 {
1047 PrintMessage(7,"found all edge points in stl file");
1048
1049
1050 StoreEdgeData();
1051
1052 int oldp = 0;
1053
1054 for (i = 1; i <= importlines.Size(); i++)
1055 {
1056 int newp = importlines.Get(i);
1057 if (!importpnums.Get(abs(newp)))
1058 newp = 0;
1059
1060 if (oldp && newp)
1061 {
1062 int en = edgedata->GetEdgeNum(importpnums.Get(oldp),
1063 importpnums.Get(abs(newp)));
1064 edgedata->Elem(en).SetStatus (ED_CONFIRMED);
1065 }
1066
1067 if (newp < 0)
1068 oldp = 0;
1069 else
1070 oldp = newp;
1071 }
1072 }
1073
1074
1075 }
1076
1077
1078
ExportEdges()1079 void STLGeometry :: ExportEdges()
1080 {
1081 PrintFnStart("Save edges to file 'edges.ng'");
1082
1083 ofstream fout("edges.ng");
1084 fout.precision(16);
1085
1086 int n = edgedata->GetNConfEdges();
1087
1088 fout << n << endl;
1089
1090 int i;
1091 for (i = 1; i <= edgedata->Size(); i++)
1092 {
1093 if (edgedata->Get(i).GetStatus() == ED_CONFIRMED)
1094 {
1095 const STLTopEdge & e = edgedata->Get(i);
1096 fout << GetPoint(e.PNum(1))(0) << " " << GetPoint(e.PNum(1))(1) << " " << GetPoint(e.PNum(1))(2) << endl;
1097 fout << GetPoint(e.PNum(2))(0) << " " << GetPoint(e.PNum(2))(1) << " " << GetPoint(e.PNum(2))(2) << endl;
1098 }
1099 }
1100
1101 }
1102
LoadEdgeData(const char * file)1103 void STLGeometry :: LoadEdgeData(const char* file)
1104 {
1105 StoreEdgeData();
1106
1107 PrintFnStart("Load edges from file '", file, "'");
1108 ifstream fin(file);
1109
1110 edgedata->Read(fin);
1111
1112 // calcedgedataanglesnew = 1;
1113 }
1114
SaveEdgeData(const char * file)1115 void STLGeometry :: SaveEdgeData(const char* file)
1116 {
1117 PrintFnStart("save edges to file '", file, "'");
1118 ofstream fout(file);
1119
1120 edgedata->Write(fout);
1121 }
1122
1123
1124
1125
1126
1127
1128
1129 /*
1130 void STLGeometry :: SaveExternalEdges()
1131 {
1132 ofstream fout("externaledgesp3.ng");
1133 fout.precision(16);
1134
1135 int n = NOExternalEdges();
1136 fout << n << endl;
1137
1138 int i;
1139 for (i = 1; i <= n; i++)
1140 {
1141 twoint e = GetExternalEdge(i);
1142 fout << GetPoint(e.i1)(0) << " " << GetPoint(e.i1)(1) << " " << GetPoint(e.i1)(2) << endl;
1143 fout << GetPoint(e.i2)(0) << " " << GetPoint(e.i2)(1) << " " << GetPoint(e.i2)(2) << endl;
1144 }
1145
1146 }
1147 */
StoreExternalEdges()1148 void STLGeometry :: StoreExternalEdges()
1149 {
1150 storedexternaledges.SetSize(0);
1151 undoexternaledges = 1;
1152 int i;
1153 for (i = 1; i <= externaledges.Size(); i++)
1154 {
1155 storedexternaledges.Append(externaledges.Get(i));
1156 }
1157
1158 }
1159
UndoExternalEdges()1160 void STLGeometry :: UndoExternalEdges()
1161 {
1162 if (!undoexternaledges)
1163 {
1164 PrintMessage(1, "undo not further possible!");
1165 return;
1166 }
1167 RestoreExternalEdges();
1168 undoexternaledges = 0;
1169 }
1170
RestoreExternalEdges()1171 void STLGeometry :: RestoreExternalEdges()
1172 {
1173 externaledges.SetSize(0);
1174 int i;
1175 for (i = 1; i <= storedexternaledges.Size(); i++)
1176 {
1177 externaledges.Append(storedexternaledges.Get(i));
1178 }
1179
1180 }
1181
1182
AddExternalEdgeAtSelected()1183 void STLGeometry :: AddExternalEdgeAtSelected()
1184 {
1185 StoreExternalEdges();
1186 if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
1187 {
1188 int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
1189 int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
1190 if (!IsExternalEdge(ap1,ap2)) {AddExternalEdge(ap1,ap2);}
1191 }
1192 }
1193
AddClosedLinesToExternalEdges()1194 void STLGeometry :: AddClosedLinesToExternalEdges()
1195 {
1196 StoreExternalEdges();
1197
1198 int i, j;
1199 for (i = 1; i <= GetNLines(); i++)
1200 {
1201 STLLine* l = GetLine(i);
1202 if (l->StartP() == l->EndP())
1203 {
1204 for (j = 1; j < l->NP(); j++)
1205 {
1206 int ap1 = l->PNum(j);
1207 int ap2 = l->PNum(j+1);
1208
1209 if (!IsExternalEdge(ap1,ap2)) {AddExternalEdge(ap1,ap2);}
1210 }
1211 }
1212 }
1213 }
1214
AddLongLinesToExternalEdges()1215 void STLGeometry :: AddLongLinesToExternalEdges()
1216 {
1217 StoreExternalEdges();
1218
1219 double diamfact = stldoctor.dirtytrigfact;
1220 double diam = GetBoundingBox().Diam();
1221
1222 int i, j;
1223 for (i = 1; i <= GetNLines(); i++)
1224 {
1225 STLLine* l = GetLine(i);
1226 if (l->GetLength(points) >= diamfact*diam)
1227 {
1228 for (j = 1; j < l->NP(); j++)
1229 {
1230 int ap1 = l->PNum(j);
1231 int ap2 = l->PNum(j+1);
1232
1233 if (!IsExternalEdge(ap1,ap2)) {AddExternalEdge(ap1,ap2);}
1234 }
1235 }
1236 }
1237 }
1238
AddAllNotSingleLinesToExternalEdges()1239 void STLGeometry :: AddAllNotSingleLinesToExternalEdges()
1240 {
1241 StoreExternalEdges();
1242
1243 int i, j;
1244 for (i = 1; i <= GetNLines(); i++)
1245 {
1246 STLLine* l = GetLine(i);
1247 if (GetNEPP(l->StartP()) > 1 || GetNEPP(l->EndP()) > 1)
1248 {
1249 for (j = 1; j < l->NP(); j++)
1250 {
1251 int ap1 = l->PNum(j);
1252 int ap2 = l->PNum(j+1);
1253
1254 if (!IsExternalEdge(ap1,ap2)) {AddExternalEdge(ap1,ap2);}
1255 }
1256 }
1257 }
1258 }
1259
DeleteDirtyExternalEdges()1260 void STLGeometry :: DeleteDirtyExternalEdges()
1261 {
1262 //delete single triangle edges and single edge-lines in clusters"
1263 StoreExternalEdges();
1264
1265 int i, j;
1266 for (i = 1; i <= GetNLines(); i++)
1267 {
1268 STLLine* l = GetLine(i);
1269 if (l->NP() <= 3 || (l->StartP() == l->EndP() && l->NP() == 4))
1270 {
1271 for (j = 1; j < l->NP(); j++)
1272 {
1273 int ap1 = l->PNum(j);
1274 int ap2 = l->PNum(j+1);
1275
1276 if (IsExternalEdge(ap1,ap2)) {DeleteExternalEdge(ap1,ap2);}
1277 }
1278 }
1279 }
1280 }
1281
AddExternalEdgesFromGeomLine()1282 void STLGeometry :: AddExternalEdgesFromGeomLine()
1283 {
1284 StoreExternalEdges();
1285 if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
1286 {
1287 int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
1288 int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
1289
1290 if (IsEdge(ap1,ap2))
1291 {
1292 int edgenum = IsEdgeNum(ap1,ap2);
1293 if (!IsExternalEdge(ap1,ap2)) {AddExternalEdge(ap1,ap2);}
1294
1295 int noend = 1;
1296 int startp = ap1;
1297 int laste = edgenum;
1298 int np1, np2;
1299 while (noend)
1300 {
1301 if (GetNEPP(startp) == 2)
1302 {
1303 if (GetEdgePP(startp,1) != laste) {laste = GetEdgePP(startp,1);}
1304 else {laste = GetEdgePP(startp,2);}
1305 np1 = GetEdge(laste).PNum(1);
1306 np2 = GetEdge(laste).PNum(2);
1307
1308 if (!IsExternalEdge(np1, np2)) {AddExternalEdge(np1, np2);}
1309 else {noend = 0;}
1310 if (np1 != startp) {startp = np1;}
1311 else {startp = np2;}
1312 }
1313 else {noend = 0;}
1314 }
1315
1316 startp = ap2;
1317 laste = edgenum;
1318 noend = 1;
1319 while (noend)
1320 {
1321 if (GetNEPP(startp) == 2)
1322 {
1323 if (GetEdgePP(startp,1) != laste) {laste = GetEdgePP(startp,1);}
1324 else {laste = GetEdgePP(startp,2);}
1325 np1 = GetEdge(laste).PNum(1);
1326 np2 = GetEdge(laste).PNum(2);
1327
1328 if (!IsExternalEdge(np1, np2)) {AddExternalEdge(np1, np2);}
1329 else {noend = 0;}
1330 if (np1 != startp) {startp = np1;}
1331 else {startp = np2;}
1332 }
1333 else {noend = 0;}
1334 }
1335
1336 }
1337
1338 }
1339
1340 }
1341
ClearEdges()1342 void STLGeometry :: ClearEdges()
1343 {
1344 edgesfound = 0;
1345 edges.SetSize(0);
1346 //edgedata->SetSize(0);
1347 // externaledges.SetSize(0);
1348 edgesperpoint.SetSize(0);
1349 undoexternaledges = 0;
1350
1351 }
1352
STLDoctorBuildEdges(const STLParameters & stlparam)1353 void STLGeometry :: STLDoctorBuildEdges(const STLParameters& stlparam)
1354 {
1355 // if (!trigsconverted) {return;}
1356 ClearEdges();
1357
1358 meshlines.SetSize(0);
1359 FindEdgesFromAngles(stlparam);
1360 }
1361
DeleteExternalEdgeAtSelected()1362 void STLGeometry :: DeleteExternalEdgeAtSelected()
1363 {
1364 StoreExternalEdges();
1365 if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
1366 {
1367 int ap1 = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
1368 int ap2 = GetTriangle(GetSelectTrig()).PNumMod(GetNodeOfSelTrig()+1);
1369 if (IsExternalEdge(ap1,ap2)) {DeleteExternalEdge(ap1,ap2);}
1370 }
1371 }
1372
DeleteExternalEdgeInVicinity()1373 void STLGeometry :: DeleteExternalEdgeInVicinity()
1374 {
1375 StoreExternalEdges();
1376 if (!stldoctor.showvicinity || vicinity.Size() != GetNT()) {return;}
1377
1378 int i, j, ap1, ap2;
1379
1380 for (i = 1; i <= GetNT(); i++)
1381 {
1382 if (vicinity.Elem(i))
1383 {
1384 for (j = 1; j <= 3; j++)
1385 {
1386 ap1 = GetTriangle(i).PNum(j);
1387 ap2 = GetTriangle(i).PNumMod(j+1);
1388
1389 if (IsExternalEdge(ap1,ap2))
1390 {
1391 DeleteExternalEdge(ap1,ap2);
1392 }
1393 }
1394 }
1395 }
1396 }
1397
BuildExternalEdgesFromEdges()1398 void STLGeometry :: BuildExternalEdgesFromEdges()
1399 {
1400 StoreExternalEdges();
1401
1402 if (GetNE() == 0) {PrintWarning("Edges possibly not generated!");}
1403
1404 int i;
1405 externaledges.SetSize(0);
1406
1407 for (i = 1; i <= GetNE(); i++)
1408 {
1409 STLEdge e = GetEdge(i);
1410 AddExternalEdge(e.PNum(1), e.PNum(2));
1411 }
1412
1413 }
1414
1415
AddExternalEdge(int ap1,int ap2)1416 void STLGeometry :: AddExternalEdge(int ap1, int ap2)
1417 {
1418 externaledges.Append(twoint(ap1,ap2));
1419 }
1420
DeleteExternalEdge(int ap1,int ap2)1421 void STLGeometry :: DeleteExternalEdge(int ap1, int ap2)
1422 {
1423
1424 int i;
1425 int found = 0;
1426 for (i = 1; i <= NOExternalEdges(); i++)
1427 {
1428 if ((GetExternalEdge(i).i1 == ap1 && GetExternalEdge(i).i2 == ap2) ||
1429 (GetExternalEdge(i).i1 == ap2 && GetExternalEdge(i).i2 == ap1)) {found = 1;};
1430 if (found && i < NOExternalEdges())
1431 {
1432 externaledges.Elem(i) = externaledges.Get(i+1);
1433 }
1434 }
1435 if (!found) {PrintWarning("edge not found");}
1436 else
1437 {
1438 externaledges.SetSize(externaledges.Size()-1);
1439 }
1440
1441 }
1442
IsExternalEdge(int ap1,int ap2)1443 int STLGeometry :: IsExternalEdge(int ap1, int ap2)
1444 {
1445 int i;
1446 for (i = 1; i <= NOExternalEdges(); i++)
1447 {
1448 if ((GetExternalEdge(i).i1 == ap1 && GetExternalEdge(i).i2 == ap2) ||
1449 (GetExternalEdge(i).i1 == ap2 && GetExternalEdge(i).i2 == ap1)) {return 1;};
1450 }
1451 return 0;
1452 }
1453
DestroyDirtyTrigs()1454 void STLGeometry :: DestroyDirtyTrigs()
1455 {
1456
1457 PrintFnStart("Destroy dirty triangles");
1458 PrintMessage(5,"original number of triangles=", GetNT());
1459
1460 //destroy every triangle with other than 3 neighbours;
1461 int changed = 1;
1462 int i, j, k;
1463 while (changed)
1464 {
1465 changed = 0;
1466 Clear();
1467
1468 for (i = 1; i <= GetNT(); i++)
1469 {
1470 int dirty = NONeighbourTrigs(i) < 3;
1471
1472 for (j = 1; j <= 3; j++)
1473 {
1474 int pnum = GetTriangle(i).PNum(j);
1475 /*
1476 if (pnum == 1546)
1477 {
1478 // for (k = 1; k <= NOTrigsPerPoint(pnum); k++)
1479 }
1480 */
1481 if (NOTrigsPerPoint(pnum) <= 2)
1482 dirty = 1;
1483 }
1484
1485 int pi1 = GetTriangle(i).PNum(1);
1486 int pi2 = GetTriangle(i).PNum(2);
1487 int pi3 = GetTriangle(i).PNum(3);
1488 if (pi1 == pi2 || pi1 == pi3 || pi2 == pi3)
1489 {
1490 PrintMessage(5,"triangle with Volume 0: ", i, " nodes: ", pi1, ", ", pi2, ", ", pi3);
1491 dirty = 1;
1492 }
1493
1494 if (dirty)
1495 {
1496 for (k = i+1; k <= GetNT(); k++)
1497 {
1498 trias[k-1] = trias[k];
1499 // readtrias: not longer permanent, JS
1500 // readtrias.Elem(k-1) = readtrias.Get(k);
1501 }
1502 int size = GetNT();
1503 trias.SetSize(size-1);
1504 // readtrias.SetSize(size-1);
1505 changed = 1;
1506 break;
1507 }
1508 }
1509 }
1510
1511 FindNeighbourTrigs();
1512 PrintMessage(5,"final number of triangles=", GetNT());
1513 }
1514
CalcNormalsFromGeometry()1515 void STLGeometry :: CalcNormalsFromGeometry()
1516 {
1517 int i;
1518 for (i = 1; i <= GetNT(); i++)
1519 {
1520 const STLTriangle & tr = GetTriangle(i);
1521 const Point3d& ap1 = GetPoint(tr.PNum(1));
1522 const Point3d& ap2 = GetPoint(tr.PNum(2));
1523 const Point3d& ap3 = GetPoint(tr.PNum(3));
1524
1525 Vec3d normal = Cross (ap2-ap1, ap3-ap1);
1526
1527 if (normal.Length() != 0)
1528 {
1529 normal /= (normal.Length());
1530 }
1531 GetTriangle(i).SetNormal(normal);
1532 }
1533 PrintMessage(5,"Normals calculated from geometry!!!");
1534
1535 calcedgedataanglesnew = 1;
1536 }
1537
SetSelectTrig(int trig)1538 void STLGeometry :: SetSelectTrig(int trig)
1539 {
1540 stldoctor.selecttrig = trig;
1541 }
1542
GetSelectTrig() const1543 int STLGeometry :: GetSelectTrig() const
1544 {
1545 return stldoctor.selecttrig;
1546 }
1547
SetNodeOfSelTrig(int n)1548 void STLGeometry :: SetNodeOfSelTrig(int n)
1549 {
1550 stldoctor.nodeofseltrig = n;
1551 }
1552
GetNodeOfSelTrig() const1553 int STLGeometry :: GetNodeOfSelTrig() const
1554 {
1555 return stldoctor.nodeofseltrig;
1556 }
1557
MoveSelectedPointToMiddle()1558 void STLGeometry :: MoveSelectedPointToMiddle()
1559 {
1560 if (GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
1561 {
1562 int p = GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig());
1563 Point<3> pm(0.,0.,0.); //Middlevector;
1564 Point<3> p0(0.,0.,0.);
1565 PrintMessage(5,"original point=", Point3d(GetPoint(p)));
1566
1567 int i;
1568 int cnt = 0;
1569 for (i = 1; i <= trigsperpoint.EntrySize(p); i++)
1570 {
1571 const STLTriangle& tr = GetTriangle(trigsperpoint.Get(p,i));
1572 int j;
1573 for (j = 1; j <= 3; j++)
1574 {
1575 if (tr.PNum(j) != p)
1576 {
1577 cnt++;
1578 pm(0) += GetPoint(tr.PNum(j))(0);
1579 pm(1) += GetPoint(tr.PNum(j))(1);
1580 pm(2) += GetPoint(tr.PNum(j))(2);
1581 }
1582 }
1583 }
1584
1585 Point<3> origp = GetPoint(p);
1586 double fact = 0.2;
1587
1588 SetPoint(p, p0 + fact*(1./(double)cnt)*(pm-p0)+(1.-fact)*(origp-p0));
1589
1590 PrintMessage(5,"middle point=", Point3d (GetPoint(p)));
1591
1592 PrintMessage(5,"moved point ", Point3d (p));
1593
1594 }
1595 }
1596
PrintSelectInfo()1597 void STLGeometry :: PrintSelectInfo()
1598 {
1599
1600 //int trig = GetSelectTrig();
1601 //int p = GetTriangle(trig).PNum(GetNodeOfSelTrig());
1602
1603 PrintMessage(1,"touch triangle ", GetSelectTrig()
1604 , ", local node ", GetNodeOfSelTrig()
1605 , " (=", int(GetTriangle(GetSelectTrig()).PNum(GetNodeOfSelTrig())), ")");
1606 if (AtlasMade() && GetSelectTrig() >= 1 && GetSelectTrig() <= GetNT())
1607 {
1608 PrintMessage(1," chartnum=", int(GetChartNr(GetSelectTrig())));
1609 /*
1610 PointBetween(Center(Center(GetPoint(GetTriangle(270).PNum(1)),
1611 GetPoint(GetTriangle(270).PNum(2))),
1612 GetPoint(GetTriangle(270).PNum(3))),270,
1613 Center(Center(GetPoint(GetTriangle(trig).PNum(1)),
1614 GetPoint(GetTriangle(trig).PNum(2))),
1615 GetPoint(GetTriangle(trig).PNum(3))),trig);
1616 */
1617 //PointBetween(Point3d(5.7818, 7.52768, 4.14879),260,Point3d(6.80292, 6.55392, 4.70184),233);
1618 }
1619 }
1620
ShowSelectedTrigChartnum()1621 void STLGeometry :: ShowSelectedTrigChartnum()
1622 {
1623 int st = GetSelectTrig();
1624
1625 if (st >= 1 && st <= GetNT() && AtlasMade())
1626 PrintMessage(1,"selected trig ", st, " has chartnumber ", int(GetChartNr(st)));
1627 }
1628
ShowSelectedTrigCoords()1629 void STLGeometry :: ShowSelectedTrigCoords()
1630 {
1631 int st = GetSelectTrig();
1632
1633 /*
1634 //testing!!!!
1635 NgArray<int> trigs;
1636 GetSortedTrianglesAroundPoint(GetTriangle(st).PNum(GetNodeOfSelTrig()),st,trigs);
1637 */
1638
1639 if (st >= 1 && st <= GetNT())
1640 {
1641 PrintMessage(1, "coordinates of selected trig ", st, ":");
1642 PrintMessage(1, " p1 = ", int(GetTriangle(st).PNum(1)), " = ",
1643 Point3d (GetPoint(GetTriangle(st).PNum(1))));
1644 PrintMessage(1, " p2 = ", int(GetTriangle(st).PNum(2)), " = ",
1645 Point3d (GetPoint(GetTriangle(st).PNum(2))));
1646 PrintMessage(1, " p3 = ", int(GetTriangle(st).PNum(3)), " = ",
1647 Point3d (GetPoint(GetTriangle(st).PNum(3))));
1648 }
1649 }
1650
LoadMarkedTrigs()1651 void STLGeometry :: LoadMarkedTrigs()
1652 {
1653 PrintFnStart("load marked trigs from file 'markedtrigs.ng'");
1654 ifstream fin("markedtrigs.ng");
1655
1656 int n;
1657 fin >> n;
1658 if (n != GetNT() || n == 0) {PrintError("Not a suitable marked-trig-file!"); return;}
1659
1660 int i, m;
1661 for (i = 1; i <= n; i++)
1662 {
1663 fin >> m;
1664 SetMarkedTrig(i, m);
1665 }
1666
1667 fin >> n;
1668 if (n != 0)
1669 {
1670
1671 Point<3> ap1, ap2;
1672 for (i = 1; i <= n; i++)
1673 {
1674 fin >> ap1(0); fin >> ap1(1); fin >> ap1(2);
1675 fin >> ap2(0); fin >> ap2(1); fin >> ap2(2);
1676 AddMarkedSeg(ap1,ap2);
1677 }
1678 }
1679 }
1680
SaveMarkedTrigs()1681 void STLGeometry :: SaveMarkedTrigs()
1682 {
1683 PrintFnStart("save marked trigs to file 'markedtrigs.ng'");
1684 ofstream fout("markedtrigs.ng");
1685
1686 int n = GetNT();
1687 fout << n << endl;
1688
1689 int i;
1690 for (i = 1; i <= n; i++)
1691 {
1692 fout << IsMarkedTrig(i) << "\n";
1693 }
1694
1695 n = GetNMarkedSegs();
1696 fout << n << endl;
1697
1698 Point<3> ap1,ap2;
1699 for (i = 1; i <= n; i++)
1700 {
1701 GetMarkedSeg(i,ap1,ap2);
1702 fout << ap1(0) << " " << ap1(1) << " " << ap1(2) << " ";
1703 fout << ap2(0) << " " << ap2(1) << " " << ap2(2) << " " << "\n";
1704 }
1705
1706 }
1707
NeighbourAnglesOfSelectedTrig()1708 void STLGeometry :: NeighbourAnglesOfSelectedTrig()
1709 {
1710 int st = GetSelectTrig();
1711
1712 if (st >= 1 && st <= GetNT())
1713 {
1714 int i;
1715 PrintMessage(1,"Angle to triangle ", st, ":");
1716 for (i = 1; i <= NONeighbourTrigs(st); i++)
1717 {
1718 PrintMessage(1," triangle ", int(NeighbourTrig(st,i)), ": angle = ",
1719 180./M_PI*GetAngle(st, NeighbourTrig(st,i)), "°",
1720 ", calculated = ", 180./M_PI*Angle(GetTriangle(st).GeomNormal(points),
1721 GetTriangle(NeighbourTrig(st,i)).GeomNormal(points)), "°");
1722 }
1723 }
1724 }
1725
GetVicinity(int starttrig,int size,NgArray<int> & vic)1726 void STLGeometry :: GetVicinity(int starttrig, int size, NgArray<int>& vic)
1727 {
1728 if (starttrig == 0 || starttrig > GetNT()) {return;}
1729
1730 NgArray<int> vicarray;
1731 vicarray.SetSize(GetNT());
1732
1733 int i;
1734 for (i = 1; i <= vicarray.Size(); i++)
1735 {
1736 vicarray.Elem(i) = 0;
1737 }
1738
1739 vicarray.Elem(starttrig) = 1;
1740
1741 int j = 0,k;
1742
1743 NgArray <int> list1;
1744 list1.SetSize(0);
1745 NgArray <int> list2;
1746 list2.SetSize(0);
1747 list1.Append(starttrig);
1748
1749 while (j < size)
1750 {
1751 j++;
1752 for (i = 1; i <= list1.Size(); i++)
1753 {
1754 for (k = 1; k <= NONeighbourTrigs(i); k++)
1755 {
1756 int nbtrig = NeighbourTrig(list1.Get(i),k);
1757 if (nbtrig && vicarray.Get(nbtrig) == 0)
1758 {
1759 list2.Append(nbtrig);
1760 vicarray.Elem(nbtrig) = 1;
1761 }
1762 }
1763 }
1764 list1.SetSize(0);
1765 for (i = 1; i <= list2.Size(); i++)
1766 {
1767 list1.Append(list2.Get(i));
1768 }
1769 list2.SetSize(0);
1770 }
1771
1772 vic.SetSize(0);
1773 for (i = 1; i <= vicarray.Size(); i++)
1774 {
1775 if (vicarray.Get(i)) {vic.Append(i);}
1776 }
1777 }
1778
CalcVicinity(int starttrig)1779 void STLGeometry :: CalcVicinity(int starttrig)
1780 {
1781 if (starttrig == 0 || starttrig > GetNT()) {return;}
1782
1783 vicinity.SetSize(GetNT());
1784
1785 if (!stldoctor.showvicinity) {return;}
1786
1787 int i;
1788 for (i = 1; i <= vicinity.Size(); i++)
1789 {
1790 vicinity.Elem(i) = 0;
1791 }
1792
1793 vicinity.Elem(starttrig) = 1;
1794
1795 int j = 0,k;
1796
1797 NgArray <int> list1;
1798 list1.SetSize(0);
1799 NgArray <int> list2;
1800 list2.SetSize(0);
1801 list1.Append(starttrig);
1802
1803 // int cnt = 1;
1804 while (j < stldoctor.vicinity)
1805 {
1806 j++;
1807 for (i = 1; i <= list1.Size(); i++)
1808 {
1809 for (k = 1; k <= NONeighbourTrigs(i); k++)
1810 {
1811 int nbtrig = NeighbourTrig(list1.Get(i),k);
1812 if (nbtrig && vicinity.Get(nbtrig) == 0)
1813 {
1814 list2.Append(nbtrig);
1815 vicinity.Elem(nbtrig) = 1;
1816 //cnt++;
1817 }
1818 }
1819 }
1820 list1.SetSize(0);
1821 for (i = 1; i <= list2.Size(); i++)
1822 {
1823 list1.Append(list2.Get(i));
1824 }
1825 list2.SetSize(0);
1826 }
1827
1828 }
1829
Vicinity(int trig) const1830 int STLGeometry :: Vicinity(int trig) const
1831 {
1832 if (trig <= vicinity.Size() && trig >=1)
1833 {
1834 return vicinity.Get(trig);
1835 }
1836 else {PrintSysError("In STLGeometry::Vicinity");}
1837 return 0;
1838 }
1839
InitMarkedTrigs()1840 void STLGeometry :: InitMarkedTrigs()
1841 {
1842 markedtrigs.SetSize(GetNT());
1843 int i;
1844 for (i = 1; i <= GetNT(); i++)
1845 {
1846 SetMarkedTrig(i, 0);
1847 }
1848 }
1849
MarkDirtyTrigs(const STLParameters & stlparam)1850 void STLGeometry :: MarkDirtyTrigs(const STLParameters& stlparam)
1851 {
1852 PrintFnStart("mark dirty trigs");
1853 int i,j;
1854
1855 markedtrigs.SetSize(GetNT());
1856
1857 for (i = 1; i <= GetNT(); i++)
1858 {
1859 SetMarkedTrig(i, 0);
1860 }
1861
1862 int found;
1863 double dirtyangle = stlparam.yangle/2./180.*M_PI;
1864 int cnt = 0;
1865 for (i = 1; i <= GetNT(); i++)
1866 {
1867 found = 0;
1868 for (j = 1; j <= NONeighbourTrigs(i); j++)
1869 {
1870 if (GetAngle(i, NeighbourTrig(i,j)) > dirtyangle)
1871 {
1872 found++;
1873 }
1874 }
1875 if (found && GetTriangle(i).MinHeight(points) <
1876 stldoctor.dirtytrigfact*GetTriangle(i).MaxLength(points))
1877 {
1878 SetMarkedTrig(i, 1); cnt++;
1879 }
1880 /*
1881 else if (found == 3)
1882 {
1883 SetMarkedTrig(i, 1); cnt++;
1884 }
1885 */
1886 }
1887
1888 PrintMessage(1, "marked ", cnt, " dirty trigs");
1889 }
1890
1891
MarkTopErrorTrigs()1892 void STLGeometry :: MarkTopErrorTrigs()
1893 {
1894 int cnt = 0;
1895 markedtrigs.SetSize(GetNT());
1896 for (int i = 1; i <= GetNT(); i++)
1897 {
1898 const STLTriangle & trig = GetTriangle(i);
1899
1900 SetMarkedTrig(i, trig.flags.toperror);
1901 if (trig.flags.toperror) cnt++;
1902 }
1903 PrintMessage(1,"marked ", cnt, " inconsistent triangles");
1904 }
1905
1906
1907
CalcTrigBadness(int i)1908 double STLGeometry :: CalcTrigBadness(int i)
1909 {
1910 int j;
1911 double maxbadness = 0;
1912 STLPointId ap1, ap2;
1913 for (j = 1; j <= NONeighbourTrigs(i); j++)
1914 {
1915 GetTriangle(i).GetNeighbourPoints(GetTriangle(NeighbourTrig(i,j)), ap1, ap2);
1916
1917 if (!IsEdge(ap1,ap2) && GetGeomAngle(i, NeighbourTrig(i,j)) > maxbadness)
1918 {
1919 maxbadness = GetGeomAngle(i, NeighbourTrig(i,j));
1920 }
1921 }
1922 return maxbadness;
1923
1924 }
1925
GeomSmoothRevertedTrigs(const STLParameters & stlparam)1926 void STLGeometry :: GeomSmoothRevertedTrigs(const STLParameters& stlparam)
1927 {
1928 //double revertedangle = stldoctor.smoothangle/180.*M_PI;
1929 double fact = stldoctor.dirtytrigfact;
1930
1931 MarkRevertedTrigs(stlparam);
1932
1933 int i, j, k, l, p;
1934
1935 for (i = 1; i <= GetNT(); i++)
1936 {
1937 if (IsMarkedTrig(i))
1938 {
1939 for (j = 1; j <= 3; j++)
1940 {
1941 double origbadness = CalcTrigBadness(i);
1942
1943 p = GetTriangle(i).PNum(j);
1944 Point<3> pm(0.,0.,0.); //Middlevector;
1945 Point<3> p0(0.,0.,0.);
1946
1947 int cnt = 0;
1948
1949 for (k = 1; k <= trigsperpoint.EntrySize(p); k++)
1950 {
1951 const STLTriangle& tr = GetTriangle(trigsperpoint.Get(p,k));
1952 for (l = 1; l <= 3; l++)
1953 {
1954 if (tr.PNum(l) != p)
1955 {
1956 cnt++;
1957 pm(0) += GetPoint(tr.PNum(l))(0);
1958 pm(1) += GetPoint(tr.PNum(l))(1);
1959 pm(2) += GetPoint(tr.PNum(l))(2);
1960 }
1961 }
1962 }
1963 Point3d origp = GetPoint(p);
1964 Point3d newp = p0 + fact*(1./(double)cnt)*(pm-p0)+(1.-fact)*(origp-p0);
1965
1966 SetPoint(p, newp);
1967
1968 if (CalcTrigBadness(i) > 0.9*origbadness) {SetPoint(p,origp); PrintDot('f');}
1969 else {PrintDot('s');}
1970 }
1971 }
1972 }
1973 MarkRevertedTrigs(stlparam);
1974 }
1975
MarkRevertedTrigs(const STLParameters & stlparam)1976 void STLGeometry :: MarkRevertedTrigs(const STLParameters& stlparam)
1977 {
1978 int i,j;
1979 if (edgesperpoint.Size() != GetNP()) {BuildEdges(stlparam);}
1980
1981 PrintFnStart("mark reverted trigs");
1982
1983 InitMarkedTrigs();
1984
1985 int found;
1986 double revertedangle = stldoctor.smoothangle/180.*M_PI;
1987
1988 int cnt = 0;
1989 STLPointId ap1, ap2;
1990 for (i = 1; i <= GetNT(); i++)
1991 {
1992 found = 0;
1993 for (j = 1; j <= NONeighbourTrigs(i); j++)
1994 {
1995 GetTriangle(i).GetNeighbourPoints(GetTriangle(NeighbourTrig(i,j)), ap1, ap2);
1996
1997 if (!IsEdge(ap1,ap2))
1998 {
1999 if (GetGeomAngle(i, NeighbourTrig(i,j)) > revertedangle)
2000 {
2001 found = 1;
2002 break;
2003 }
2004 }
2005 }
2006
2007 if (found)
2008 {
2009 SetMarkedTrig(i, 1); cnt++;
2010 }
2011
2012 }
2013
2014 PrintMessage(5, "found ", cnt, " reverted trigs");
2015
2016
2017 }
2018
SmoothDirtyTrigs(const STLParameters & stlparam)2019 void STLGeometry :: SmoothDirtyTrigs(const STLParameters& stlparam)
2020 {
2021 PrintFnStart("smooth dirty trigs");
2022
2023 MarkDirtyTrigs(stlparam);
2024
2025 int i,j;
2026 int changed = 1;
2027 STLPointId ap1, ap2;
2028
2029 while (changed)
2030 {
2031 changed = 0;
2032 for (i = 1; i <= GetNT(); i++)
2033 {
2034 if (IsMarkedTrig(i))
2035 {
2036 int foundtrig = 0;
2037 double maxlen = 0;
2038 // JS: darf normalvector nicht ueber kurze Seite erben
2039 maxlen = GetTriangle(i).MaxLength(GetPoints()) / 2.1; //JG: bei flachem dreieck auch kurze Seite
2040
2041 for (j = 1; j <= NONeighbourTrigs(i); j++)
2042 {
2043 if (!IsMarkedTrig(NeighbourTrig(i,j)))
2044 {
2045 GetTriangle(i).GetNeighbourPoints(GetTriangle(NeighbourTrig(i,j)),ap1,ap2);
2046 if (Dist(GetPoint(ap1),GetPoint(ap2)) >= maxlen)
2047 {
2048 foundtrig = NeighbourTrig(i,j);
2049 maxlen = Dist(GetPoint(ap1),GetPoint(ap2));
2050 }
2051 }
2052 }
2053 if (foundtrig)
2054 {
2055 GetTriangle(i).SetNormal(GetTriangle(foundtrig).Normal());
2056 changed = 1;
2057 SetMarkedTrig(i,0);
2058 }
2059 }
2060 }
2061 }
2062
2063 calcedgedataanglesnew = 1;
2064
2065
2066 MarkDirtyTrigs(stlparam);
2067
2068 int cnt = 0;
2069 for (i = 1; i <= GetNT(); i++)
2070 {
2071 if (IsMarkedTrig(i)) {cnt++;}
2072 }
2073
2074 PrintMessage(5,"NO marked dirty trigs=", cnt);
2075
2076 }
2077
IsMarkedTrig(int trig) const2078 int STLGeometry :: IsMarkedTrig(int trig) const
2079 {
2080 if (trig <= markedtrigs.Size() && trig >=1)
2081 {
2082 return markedtrigs.Get(trig);
2083 }
2084 else {PrintSysError("In STLGeometry::IsMarkedTrig");}
2085
2086 return 0;
2087 }
2088
SetMarkedTrig(int trig,int num)2089 void STLGeometry :: SetMarkedTrig(int trig, int num)
2090 {
2091 if (trig <= markedtrigs.Size() && trig >=1)
2092 {
2093 markedtrigs.Elem(trig) = num;
2094 }
2095 else {PrintSysError("In STLGeometry::SetMarkedTrig");}
2096 }
2097
Clear()2098 void STLGeometry :: Clear()
2099 {
2100 PrintFnStart("Clear");
2101
2102 surfacemeshed = 0;
2103 surfaceoptimized = 0;
2104 volumemeshed = 0;
2105
2106 selectedmultiedge.SetSize(0);
2107 meshlines.SetSize(0);
2108 // neighbourtrigs.SetSize(0);
2109 outerchartspertrig.SetSize(0);
2110 atlas.SetSize(0);
2111 ClearMarkedSegs();
2112 ClearSpiralPoints();
2113 ClearLineEndPoints();
2114
2115 SetSelectTrig(0);
2116 SetNodeOfSelTrig(1);
2117 facecnt = 0;
2118
2119 SetThreadPercent(100.);
2120
2121 ClearEdges();
2122 }
2123
Area()2124 double STLGeometry :: Area()
2125 {
2126 if (area >= 0) return area;
2127 area = 0;
2128 for (int i = 1; i <= GetNT(); i++)
2129 area += GetTriangle(i).Area(points);
2130 return area;
2131 }
2132
GetAngle(int t1,int t2)2133 double STLGeometry :: GetAngle(int t1, int t2)
2134 {
2135 return Angle(GetTriangle(t1).Normal(),GetTriangle(t2).Normal());
2136 }
2137
GetGeomAngle(int t1,int t2)2138 double STLGeometry :: GetGeomAngle(int t1, int t2)
2139 {
2140 Vec3d n1 = GetTriangle(t1).GeomNormal(points);
2141 Vec3d n2 = GetTriangle(t2).GeomNormal(points);
2142 return Angle(n1,n2);
2143 }
2144
2145
InitSTLGeometry(const NgArray<STLReadTriangle> & readtrias)2146 void STLGeometry :: InitSTLGeometry(const NgArray<STLReadTriangle> & readtrias)
2147 {
2148 PrintFnStart("Init STL Geometry");
2149 STLTopology::InitSTLGeometry(readtrias);
2150
2151 int i, k;
2152
2153 //const double geometry_tol_fact = 1E8; //distances lower than max_box_size/tol are ignored
2154
2155 int np = GetNP();
2156 PrintMessage(5,"NO points= ", GetNP());
2157 normals.SetSize(GetNP());
2158 NgArray<int> normal_cnt(GetNP()); // counts number of added normals in a point
2159
2160 for (i = 1; i <= np; i++)
2161 {
2162 normal_cnt.Elem(i) = 0;
2163 normals.Elem(i) = Vec3d (0,0,0);
2164 }
2165
2166 for(i = 1; i <= GetNT(); i++)
2167 {
2168 // STLReadTriangle t = GetReadTriangle(i);
2169 // STLTriangle st;
2170
2171 Vec<3> n = GetTriangle(i).Normal ();
2172
2173 for (k = 1; k <= 3; k++)
2174 {
2175 int pi = GetTriangle(i).PNum(k);
2176
2177 normal_cnt.Elem(pi)++;
2178 SetNormal(pi, GetNormal(pi) + n);
2179 }
2180 }
2181
2182 //normalize the normals
2183 for (i = 1; i <= GetNP(); i++)
2184 {
2185 SetNormal(i,1./(double)normal_cnt.Get(i)*GetNormal(i));
2186 }
2187
2188 trigsconverted = 1;
2189
2190 vicinity.SetSize(GetNT());
2191 markedtrigs.SetSize(GetNT());
2192 for (i = 1; i <= GetNT(); i++)
2193 {
2194 markedtrigs.Elem(i) = 0;
2195 vicinity.Elem(i) = 1;
2196 }
2197
2198 ha_points.SetSize(GetNP());
2199 for (i = 1; i <= GetNP(); i++)
2200 ha_points.Elem(i) = 0;
2201
2202 calcedgedataanglesnew = 0;
2203 edgedatastored = 0;
2204 edgedata->Clear();
2205
2206
2207 if (GetStatus() == STL_ERROR) return;
2208
2209 CalcEdgeData();
2210 CalcEdgeDataAngles();
2211
2212 ClearLineEndPoints();
2213
2214 CheckGeometryOverlapping();
2215 }
2216
TopologyChanged()2217 void STLGeometry :: TopologyChanged()
2218 {
2219 calcedgedataanglesnew = 1;
2220 }
2221
CheckGeometryOverlapping()2222 int STLGeometry :: CheckGeometryOverlapping()
2223 {
2224 PrintMessageCR(3,"Check overlapping geometry ...");
2225
2226 Box<3> geombox = GetBoundingBox();
2227 Point<3> pmin = geombox.PMin();
2228 Point<3> pmax = geombox.PMax();
2229
2230 BoxTree<3> setree(pmin, pmax);
2231
2232 int oltrigs = 0;
2233 markedtrigs.SetSize(GetNT());
2234
2235 for (int i = 1; i <= GetNT(); i++)
2236 SetMarkedTrig(i, 0);
2237
2238 for (int i = 1; i <= GetNT(); i++)
2239 {
2240 const STLTriangle & tri = GetTriangle(i);
2241
2242 Point<3> tpmin = tri.box.PMin();
2243 Point<3> tpmax = tri.box.PMax();
2244 Vec<3> diag = tpmax - tpmin;
2245
2246 tpmax = tpmax + 0.001 * diag;
2247 tpmin = tpmin - 0.001 * diag;
2248
2249 setree.Insert (tpmin, tpmax, i);
2250 }
2251
2252
2253 {
2254 mutex inters_mutex;
2255
2256 ParallelFor( 1, GetNT()+1, [&] (int first, int next)
2257 {
2258 NgArray<int> inters;
2259 for (int i=first; i<next; i++) {
2260 const STLTriangle & tri = GetTriangle(i);
2261
2262 Point<3> tpmin = tri.box.PMin();
2263 Point<3> tpmax = tri.box.PMax();
2264
2265 setree.GetIntersecting (tpmin, tpmax, inters);
2266
2267 for (int j = 1; j <= inters.Size(); j++)
2268 {
2269 const STLTriangle & tri2 = GetTriangle(inters.Get(j));
2270
2271 const Point<3> *trip1[3], *trip2[3];
2272 Point<3> hptri1[3], hptri2[3];
2273 /*
2274 for (k = 1; k <= 3; k++)
2275 {
2276 trip1[k-1] = &GetPoint (tri.PNum(k));
2277 trip2[k-1] = &GetPoint (tri2.PNum(k));
2278 }
2279 */
2280
2281 for (int k = 0; k < 3; k++)
2282 {
2283 hptri1[k] = GetPoint (tri[k]);
2284 hptri2[k] = GetPoint (tri2[k]);
2285 trip1[k] = &hptri1[k];
2286 trip2[k] = &hptri2[k];
2287 }
2288
2289 if (IntersectTriangleTriangle (&trip1[0], &trip2[0]))
2290 {
2291 lock_guard<mutex> guard(inters_mutex);
2292 {
2293 oltrigs++;
2294 PrintMessage(5,"Intersecting Triangles: trig ",i," with ",inters.Get(j),"!");
2295 SetMarkedTrig(i, 1);
2296 SetMarkedTrig(inters.Get(j), 1);
2297 }
2298 }
2299 }
2300 }
2301 });
2302 }
2303 PrintMessage(3,"Check overlapping geometry ... ", oltrigs, " triangles overlap");
2304 return oltrigs;
2305 }
2306
2307 /*
2308 void STLGeometry :: InitSTLGeometry()
2309 {
2310 STLTopology::InitSTLGeometry();
2311
2312 int i, j, k;
2313
2314 const double geometry_tol_fact = 1E8; //distances lower than max_box_size/tol are ignored
2315
2316
2317 trias.SetSize(0);
2318 points.SetSize(0);
2319 normals.SetSize(0);
2320
2321 NgArray<int> normal_cnt; // counts number of added normals in a point
2322
2323 Box3d bb (GetBoundingBox().PMin() + Vec3d (-1,-1,-1),
2324 GetBoundingBox().PMax() + Vec3d (1, 1, 1));
2325
2326 Point3dTree pointtree (bb.PMin(),
2327 bb.PMax());
2328 NgArray<int> pintersect;
2329
2330 double gtol = GetBoundingBox().CalcDiam()/geometry_tol_fact;
2331
2332 for(i = 1; i <= GetReadNT(); i++)
2333 {
2334 //if (i%500==499) {(*mycout) << (double)i/(double)GetReadNT()*100. << "%" << endl;}
2335
2336 STLReadTriangle t = GetReadTriangle(i);
2337 STLTriangle st;
2338 Vec3d n = t.normal;
2339
2340 for (k = 0; k < 3; k++)
2341 {
2342 Point3d p = t.pts[k];
2343
2344 Point3d pmin = p - Vec3d (gtol, gtol, gtol);
2345 Point3d pmax = p + Vec3d (gtol, gtol, gtol);
2346
2347 pointtree.GetIntersecting (pmin, pmax, pintersect);
2348
2349 if (pintersect.Size() > 1)
2350 (*mycout) << "found too much " << char(7) << endl;
2351 int foundpos = 0;
2352 if (pintersect.Size())
2353 foundpos = pintersect.Get(1);
2354
2355 if (foundpos)
2356 {
2357 normal_cnt[foundpos]++;
2358 SetNormal(foundpos,GetNormal(foundpos)+n);
2359 // (*testout) << "found p " << p << endl;
2360 }
2361 else
2362 {
2363 foundpos = AddPoint(p);
2364 AddNormal(n);
2365 normal_cnt.Append(1);
2366
2367 pointtree.Insert (p, foundpos);
2368 }
2369 //(*mycout) << "foundpos=" << foundpos << endl;
2370 st.pts[k] = foundpos;
2371 }
2372
2373 if ( (st.pts[0] == st.pts[1]) ||
2374 (st.pts[0] == st.pts[2]) ||
2375 (st.pts[1] == st.pts[2]) )
2376 {
2377 (*mycout) << "ERROR: STL Triangle degenerated" << endl;
2378 }
2379 else
2380 {
2381 // do not add ? js
2382 AddTriangle(st);
2383 }
2384 //(*mycout) << "TRIG" << i << " = " << st << endl;
2385
2386 }
2387 //normal the normals
2388 for (i = 1; i <= GetNP(); i++)
2389 {
2390 SetNormal(i,1./(double)normal_cnt[i]*GetNormal(i));
2391 }
2392
2393 trigsconverted = 1;
2394
2395 vicinity.SetSize(GetNT());
2396 markedtrigs.SetSize(GetNT());
2397 for (i = 1; i <= GetNT(); i++)
2398 {
2399 markedtrigs.Elem(i) = 0;
2400 vicinity.Elem(i) = 1;
2401 }
2402
2403 ha_points.SetSize(GetNP());
2404 for (i = 1; i <= GetNP(); i++)
2405 ha_points.Elem(i) = 0;
2406
2407 calcedgedataanglesnew = 0;
2408 edgedatastored = 0;
2409 edgedata->Clear();
2410
2411 CalcEdgeData();
2412 CalcEdgeDataAngles();
2413
2414 ClearLineEndPoints();
2415
2416 (*mycout) << "done" << endl;
2417 }
2418 */
2419
2420
2421
SetLineEndPoint(int pn)2422 void STLGeometry :: SetLineEndPoint(int pn)
2423 {
2424 if (pn <1 || pn > lineendpoints.Size()) {PrintSysError("Illegal pnum in SetLineEndPoint!!!"); return; }
2425 lineendpoints.Elem(pn) = 1;
2426 }
2427
IsLineEndPoint(int pn)2428 int STLGeometry :: IsLineEndPoint(int pn)
2429 {
2430 // return 0;
2431 if (pn <1 || pn > lineendpoints.Size())
2432 {PrintSysError("Illegal pnum in IsLineEndPoint!!!"); return 0;}
2433 return lineendpoints.Get(pn);
2434 }
2435
ClearLineEndPoints()2436 void STLGeometry :: ClearLineEndPoints()
2437 {
2438 lineendpoints.SetSize(GetNP());
2439 int i;
2440 for (i = 1; i <= GetNP(); i++)
2441 {
2442 lineendpoints.Elem(i) = 0;
2443 }
2444 }
2445
IsEdge(int ap1,int ap2)2446 int STLGeometry :: IsEdge(int ap1, int ap2)
2447 {
2448 int i,j;
2449 for (i = 1; i <= GetNEPP(ap1); i++)
2450 {
2451 for (j = 1; j <= GetNEPP(ap2); j++)
2452 {
2453 if (GetEdgePP(ap1,i) == GetEdgePP(ap2,j)) {return 1;}
2454 }
2455 }
2456 return 0;
2457 }
2458
IsEdgeNum(int ap1,int ap2)2459 int STLGeometry :: IsEdgeNum(int ap1, int ap2)
2460 {
2461 int i,j;
2462 for (i = 1; i <= GetNEPP(ap1); i++)
2463 {
2464 for (j = 1; j <= GetNEPP(ap2); j++)
2465 {
2466 if (GetEdgePP(ap1,i) == GetEdgePP(ap2,j)) {return GetEdgePP(ap1,i);}
2467 }
2468 }
2469 return 0;
2470 }
2471
2472
BuildEdges(const STLParameters & stlparam)2473 void STLGeometry :: BuildEdges(const STLParameters& stlparam)
2474 {
2475 //PrintFnStart("build edges");
2476 edges.SetSize(0);
2477 meshlines.SetSize(0);
2478 FindEdgesFromAngles(stlparam);
2479 }
2480
UseExternalEdges()2481 void STLGeometry :: UseExternalEdges()
2482 {
2483 for (int i = 1; i <= NOExternalEdges(); i++)
2484 AddEdge(GetExternalEdge(i).i1,GetExternalEdge(i).i2);
2485 //BuildEdgesPerPointy();
2486 }
2487
UndoEdgeChange()2488 void STLGeometry :: UndoEdgeChange()
2489 {
2490 if (edgedatastored)
2491 {
2492 RestoreEdgeData();
2493 }
2494 else
2495 {
2496 PrintWarning("no edge undo possible");
2497 }
2498 }
2499
2500
StoreEdgeData()2501 void STLGeometry :: StoreEdgeData()
2502 {
2503 // edgedata_store = *edgedata;
2504
2505 edgedata->Store();
2506 edgedatastored = 1;
2507
2508 // put stlgeom-edgedata to stltopology edgedata
2509 /*
2510 int i;
2511 for (i = 1; i <= GetNTE(); i++)
2512 {
2513 const STLTopEdge & topedge = GetTopEdge (i);
2514 int ednum = edgedata->GetEdgeNum (topedge.PNum(1),
2515 topedge.PNum(2));
2516 topedges.Elem(i).SetStatus (edgedata->Get (ednum).status);
2517 }
2518 */
2519 }
2520
RestoreEdgeData()2521 void STLGeometry :: RestoreEdgeData()
2522 {
2523 // *edgedata = edgedata_store;
2524 edgedata->Restore();
2525 edgedatastored=0;
2526 }
2527
2528
CalcEdgeData()2529 void STLGeometry :: CalcEdgeData()
2530 {
2531 PushStatus("Calc Edge Data");
2532
2533 STLPointId np1, np2;
2534
2535 int ecnt = 0;
2536 edgedata->SetSize(GetNT()/2*3);
2537
2538 for (int i = 1; i <= GetNT(); i++)
2539 {
2540 SetThreadPercent((double)i/(double)GetNT()*100.);
2541
2542 const STLTriangle & t1 = GetTriangle(i);
2543
2544 for (int j = 1; j <= NONeighbourTrigs(i); j++)
2545 {
2546 int nbti = NeighbourTrig(i,j);
2547 if (nbti > i)
2548 {
2549 const STLTriangle & t2 = GetTriangle(nbti);
2550
2551 if (t1.IsNeighbourFrom(t2))
2552 {
2553 ecnt++; if (ecnt > edgedata->Size()) {PrintError("In Calc edge data, illegal geometry");}
2554
2555 t1.GetNeighbourPoints(t2,np1,np2);
2556
2557 /* ang = GetAngle(i,nbti);
2558 if (ang < -M_PI) {ang += 2*M_PI;}*/
2559
2560
2561 // edgedata->Add(STLEdgeData(0, np1, np2, i, nbti),ecnt);
2562 edgedata->Elem(ecnt).SetStatus(ED_UNDEFINED);
2563
2564 // edgedata->Elem(ecnt).top = this;
2565 // edgedata->Elem(ecnt).topedgenr = GetTopEdgeNum (np1, np2);
2566 }
2567 }
2568 }
2569 }
2570
2571 //BuildEdgesPerPoint();
2572 PopStatus();
2573 }
2574
CalcEdgeDataAngles()2575 void STLGeometry :: CalcEdgeDataAngles()
2576 {
2577 PrintMessageCR (5,"calc edge data angles ... ");
2578
2579 for (int i = 1; i <= GetNTE(); i++)
2580 {
2581 STLTopEdge & edge = GetTopEdge (i);
2582 double cosang =
2583 GetTriangle(edge.TrigNum(1)).Normal() *
2584 GetTriangle(edge.TrigNum(2)).Normal();
2585 edge.SetCosAngle (cosang);
2586 }
2587
2588 for (int i = 1; i <= edgedata->Size(); i++)
2589 {
2590 /*
2591 const STLEdgeData& e = edgedata->Get(i);
2592 ang = GetAngle(e.lt,e.rt);
2593 if (ang < -M_PI) {ang += 2*M_PI;}
2594 edgedata->Elem(i).angle = fabs(ang);
2595 */
2596 }
2597 PrintMessage (5,"calc edge data angles ... done");
2598 }
2599
FindEdgesFromAngles(const STLParameters & stlparam)2600 void STLGeometry :: FindEdgesFromAngles(const STLParameters& stlparam)
2601 {
2602 // PrintFnStart("find edges from angles");
2603
2604 double min_edge_angle = stlparam.yangle/180.*M_PI;
2605 double cont_min_edge_angle = stlparam.contyangle/180.*M_PI;
2606
2607 double cos_min_edge_angle = cos (min_edge_angle);
2608 double cos_cont_min_edge_angle = cos (cont_min_edge_angle);
2609
2610 if (calcedgedataanglesnew) {CalcEdgeDataAngles(); calcedgedataanglesnew = 0;}
2611
2612 for (int i = 1; i <= edgedata->Size(); i++)
2613 {
2614 STLTopEdge & sed = edgedata->Elem(i);
2615 if (sed.GetStatus() == ED_CANDIDATE ||
2616 sed.GetStatus() == ED_UNDEFINED)
2617 {
2618 if (sed.CosAngle() <= cos_min_edge_angle)
2619 {
2620 sed.SetStatus (ED_CANDIDATE);
2621 }
2622 else
2623 {
2624 sed.SetStatus(ED_UNDEFINED);
2625 }
2626 }
2627 }
2628
2629 if (stlparam.contyangle < stlparam.yangle)
2630 {
2631 int changed = 1;
2632 int its = 0;
2633 while (changed && stlparam.contyangle < stlparam.yangle)
2634 {
2635 its++;
2636 //(*mycout) << "." << flush;
2637 changed = 0;
2638 for (int i = 1; i <= edgedata->Size(); i++)
2639 {
2640 STLTopEdge & sed = edgedata->Elem(i);
2641 if (sed.CosAngle() <= cos_cont_min_edge_angle
2642 && sed.GetStatus() == ED_UNDEFINED &&
2643 (edgedata->GetNConfCandEPP(sed.PNum(1)) == 1 ||
2644 edgedata->GetNConfCandEPP(sed.PNum(2)) == 1))
2645 {
2646 changed = 1;
2647 sed.SetStatus (ED_CANDIDATE);
2648 }
2649 }
2650 }
2651 }
2652
2653 int confcand = 0;
2654 if (edgedata->GetNConfEdges() == 0)
2655 {
2656 confcand = 1;
2657 }
2658
2659 for (int i = 1; i <= edgedata->Size(); i++)
2660 {
2661 STLTopEdge & sed = edgedata->Elem(i);
2662 if (sed.GetStatus() == ED_CONFIRMED ||
2663 (sed.GetStatus() == ED_CANDIDATE && confcand))
2664 {
2665 STLEdge se(sed.PNum(1),sed.PNum(2));
2666 se.SetLeftTrig(sed.TrigNum(1));
2667 se.SetRightTrig(sed.TrigNum(2));
2668 AddEdge(se);
2669 }
2670 }
2671 BuildEdgesPerPoint();
2672
2673
2674
2675 //(*mycout) << "its for continued angle = " << its << endl;
2676 PrintMessage(5,"built ", GetNE(), " edges with yellow angle = ", stlparam.yangle, " degree");
2677
2678 }
2679
2680 /*
2681 void STLGeometry :: FindEdgesFromAngles()
2682 {
2683 double yangle = stlparam.yangle;
2684 char * savetask = multithread.task;
2685 multithread.task = "find edges";
2686
2687 const double min_edge_angle = yangle/180.*M_PI;
2688
2689 int np1, np2;
2690 double ang;
2691 int i;
2692
2693 //(*mycout) << "area=" << Area() << endl;
2694
2695 for (i = 1; i <= GetNT(); i++)
2696 {
2697 multithread.percent = (double)i/(double)GetReadNT()*100.;
2698
2699 const STLTriangle & t1 = GetTriangle(i);
2700 //NeighbourTrigs(nt,i);
2701
2702 for (int j = 1; j <= NONeighbourTrigs(i); j++)
2703 {
2704 int nbti = NeighbourTrig(i,j);
2705 if (nbti > i)
2706 {
2707 const STLTriangle & t2 = GetTriangle(nbti);
2708
2709 if (t1.IsNeighbourFrom(t2))
2710 {
2711 ang = GetAngle(i,nbti);
2712 if (ang < -M_PI*0.5) {ang += 2*M_PI;}
2713
2714 t1.GetNeighbourPoints(t2,np1,np2);
2715
2716 if (fabs(ang) >= min_edge_angle)
2717 {
2718 STLEdge se(np1,np2);
2719 se.SetLeftTrig(i);
2720 se.SetRightTrig(nbti);
2721 AddEdge(se);
2722 }
2723 }
2724 }
2725 }
2726 }
2727
2728 (*mycout) << "added " << GetNE() << " edges" << endl;
2729
2730 //BuildEdgesPerPoint();
2731
2732 multithread.percent = 100.;
2733 multithread.task = savetask;
2734
2735 }
2736 */
BuildEdgesPerPoint()2737 void STLGeometry :: BuildEdgesPerPoint()
2738 {
2739 //cout << "*** build edges per point" << endl;
2740 edgesperpoint.SetSize(GetNP());
2741
2742 //add edges to points
2743 for (int i = 1; i <= GetNE(); i++)
2744 {
2745 //(*mycout) << "EDGE " << GetEdge(i).PNum(1) << " - " << GetEdge(i).PNum(2) << endl;
2746 for (int j = 1; j <= 2; j++)
2747 {
2748 AddEdgePP(GetEdge(i).PNum(j),i);
2749 }
2750 }
2751 }
2752
AddFaceEdges()2753 void STLGeometry :: AddFaceEdges()
2754 {
2755 PrintFnStart("Add starting edges for faces");
2756
2757 //für Kugel eine STLLine hinzufügen (Vorteil: verfeinerbar, unabhängig von Auflösung der Geometrie!!!):
2758 //Grenze von 1. gefundener chart
2759
2760 NgArray<int> edgecnt;
2761 NgArray<int> chartindex;
2762 edgecnt.SetSize(GetNOFaces());
2763 chartindex.SetSize(GetNOFaces());
2764
2765 for (int i = 1; i <= GetNOFaces(); i++)
2766 {
2767 edgecnt.Elem(i) = 0;
2768 chartindex.Elem(i) = 0;
2769 }
2770
2771 for (int i = 1; i <= GetNT(); i++)
2772 {
2773 int fn = GetTriangle(i).GetFaceNum();
2774 if (!chartindex.Get(fn)) {chartindex.Elem(fn) = GetChartNr(i);}
2775 for (int j = 1; j <= 3; j++)
2776 {
2777 edgecnt.Elem(fn) += GetNEPP(GetTriangle(i).PNum(j));
2778 }
2779 }
2780
2781 for (int i = 1; i <= GetNOFaces(); i++)
2782 {
2783 if (!edgecnt.Get(i)) {PrintMessage(5,"Face", i, " has no edge!");}
2784 }
2785
2786 int changed = 0;
2787 STLPointId ap1, ap2;
2788 for (int i = 1; i <= GetNOFaces(); i++)
2789 {
2790 if (!edgecnt.Get(i))
2791 {
2792 const STLChart& c = GetChart(chartindex.Get(i));
2793 // bool foundone = false;
2794 int longest_ap1 = -1, longest_ap2 = -1;
2795 double maxlen = -1;
2796 for (int j = 1; j <= c.GetNChartT(); j++)
2797 {
2798 const STLTriangle& t1 = GetTriangle(c.GetChartTrig1(j));
2799 for (int k = 1; k <= 3; k++)
2800 {
2801 int nt = NeighbourTrig(c.GetChartTrig1(j),k);
2802 if (GetChartNr(nt) != chartindex.Get(i))
2803 {
2804 t1.GetNeighbourPoints(GetTriangle(nt),ap1,ap2);
2805 // AddEdge(ap1,ap2);
2806 double len = Dist(GetPoint(ap1), GetPoint(ap2));
2807 if (len > maxlen)
2808 {
2809 maxlen = len;
2810 longest_ap1 = ap1;
2811 longest_ap2 = ap2;
2812 }
2813 changed = 1;
2814 }
2815 }
2816 }
2817 if (maxlen > 0)
2818 AddEdge(longest_ap1,longest_ap2);
2819 }
2820
2821 }
2822
2823 if (changed) BuildEdgesPerPoint();
2824
2825 }
2826
LinkEdges(const STLParameters & stlparam)2827 void STLGeometry :: LinkEdges(const STLParameters& stlparam)
2828 {
2829 PushStatusF("Link Edges");
2830 PrintMessage(5,"have now ", GetNE(), " edges with yellow angle = ", stlparam.yangle, " degree");
2831
2832 int i;
2833
2834 lines.SetSize(0);
2835 int starte(0);
2836 int edgecnt = 0;
2837 int found;
2838 int rev(0); //indicates, that edge is inserted reverse
2839
2840 //worked edges
2841 NgArray<int> we(GetNE());
2842
2843 //setlineendpoints; wenn 180°, dann keine endpunkte
2844 //nur punkte mit 2 edges kommen in frage, da bei mehr oder weniger punkten ohnehin ein meshpoint hinkommt
2845
2846 Vec3d v1,v2;
2847 double cos_eca = cos(stlparam.edgecornerangle/180.*M_PI);
2848 int ecnt = 0;
2849 int lp1, lp2;
2850 if (stlparam.edgecornerangle < 180)
2851 {
2852 for (i = 1; i <= GetNP(); i++)
2853 {
2854 if (GetNEPP(i) == 2)
2855 {
2856 if (GetEdge(GetEdgePP(i,1)).PNum(2) == GetEdge(GetEdgePP(i,2)).PNum(1) ||
2857 GetEdge(GetEdgePP(i,1)).PNum(1) == GetEdge(GetEdgePP(i,2)).PNum(2))
2858 {
2859 lp1 = 1; lp2 = 2;
2860 }
2861 else
2862 {
2863 lp1 = 2; lp2 = 1;
2864 }
2865
2866 v1 = Vec3d(GetPoint(GetEdge(GetEdgePP(i,1)).PNum(1)),
2867 GetPoint(GetEdge(GetEdgePP(i,1)).PNum(2)));
2868 v2 = Vec3d(GetPoint(GetEdge(GetEdgePP(i,2)).PNum(lp1)),
2869 GetPoint(GetEdge(GetEdgePP(i,2)).PNum(lp2)));
2870 if ((v1*v2)/sqrt(v1.Length2()*v2.Length2()) < cos_eca)
2871 {
2872 //(*testout) << "add edgepoint " << i << endl;
2873 SetLineEndPoint(i);
2874 ecnt++;
2875 }
2876 }
2877 }
2878 }
2879 PrintMessage(5, "added ", ecnt, " mesh_points due to edge corner angle (",
2880 stlparam.edgecornerangle, " degree)");
2881
2882 for (i = 1; i <= GetNE(); i++) {we.Elem(i) = 0;}
2883
2884 while(edgecnt < GetNE())
2885 {
2886 SetThreadPercent((double)edgecnt/(double)GetNE()*100.);
2887
2888 STLLine* line = new STLLine(this);
2889
2890 //find start edge
2891 int j = 1;
2892 found = 0;
2893 //try second time, if only rings are left!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2894 int second = 0;
2895
2896 //find a starting edge at point with 1 or more than 2 edges or at lineendpoint
2897 while (!found && j<=GetNE())
2898 {
2899 if (!we.Get(j))
2900 {
2901 if (GetNEPP(GetEdge(j).PNum(1)) != 2 || IsLineEndPoint(GetEdge(j).PNum(1)))
2902 {
2903 starte = j;
2904 found = 1;
2905 rev = 0;
2906 }
2907 else
2908 if (GetNEPP(GetEdge(j).PNum(2)) != 2 || IsLineEndPoint(GetEdge(j).PNum(2)))
2909 {
2910 starte = j;
2911 found = 1;
2912 rev = 1;
2913 }
2914 else if (second)
2915 {
2916 starte = j;
2917 found = 1;
2918 rev = 0; //0 or 1 are possible
2919 }
2920 }
2921 j++;
2922 if (!second && j == GetNE()) {second = 1; j = 1;}
2923 }
2924
2925 if (!found) {PrintSysError("No starting edge found, edgecnt=", edgecnt, ", GETNE=", GetNE());}
2926
2927 line->AddPoint(GetEdge(starte).PNum(1+rev));
2928 line->AddPoint(GetEdge(starte).PNum(2-rev));
2929 if (!rev)
2930 {
2931 line->AddLeftTrig(GetEdge(starte).LeftTrig());
2932 line->AddRightTrig(GetEdge(starte).RightTrig());
2933 }
2934 else
2935 {
2936 line->AddLeftTrig(GetEdge(starte).RightTrig());
2937 line->AddRightTrig(GetEdge(starte).LeftTrig());
2938 }
2939 edgecnt++; we.Elem(starte) = 1;
2940
2941 //add segments to line as long as segments other than starting edge are found or lineendpoint is reached
2942 found = 1;
2943 int other;
2944 while(found)
2945 {
2946 found = 0;
2947 int fp = GetEdge(starte).PNum(2-rev);
2948 if (GetNEPP(fp) == 2 && !IsLineEndPoint(fp))
2949 {
2950 //find the "other" edge of point fp
2951 other = 0;
2952 if (GetEdgePP(fp,1) == starte) {other = 1;}
2953
2954 starte = GetEdgePP(fp,1+other);
2955
2956 //falls ring -> aufhoeren !!!!!!!!!!!
2957 if (!we.Elem(starte))
2958 {
2959 found = 1;
2960 rev = 0;
2961 if (GetEdge(starte).PNum(2) == fp) {rev = 1;}
2962 else if (GetEdge(starte).PNum(1) != fp) {PrintSysError("In Link Edges!");}
2963
2964 line->AddPoint(GetEdge(starte).PNum(2-rev));
2965 if (!rev)
2966 {
2967 line->AddLeftTrig(GetEdge(starte).LeftTrig());
2968 line->AddRightTrig(GetEdge(starte).RightTrig());
2969 }
2970 else
2971 {
2972 line->AddLeftTrig(GetEdge(starte).RightTrig());
2973 line->AddRightTrig(GetEdge(starte).LeftTrig());
2974 }
2975 edgecnt++; we.Elem(starte) = 1;
2976 }
2977 }
2978 }
2979 AddLine(line);
2980 }
2981 PrintMessage(5,"number of lines generated = ", GetNLines());
2982
2983 //check, which lines must have at least one midpoint
2984 INDEX_2_HASHTABLE<int> lineht(GetNLines()+1);
2985
2986 for (i = 1; i <= GetNLines(); i++)
2987 {
2988 if (GetLine(i)->StartP() == GetLine(i)->EndP())
2989 {
2990 GetLine(i)->DoSplit();
2991 }
2992 }
2993
2994 for (i = 1; i <= GetNLines(); i++)
2995 {
2996 INDEX_2 lineep (GetLine(i)->StartP(),GetLine(i)->EndP());
2997 lineep.Sort();
2998
2999 if (lineht.Used (lineep))
3000 {
3001 GetLine(i)->DoSplit();
3002 int other = lineht.Get(lineep);
3003 GetLine(other)->DoSplit();
3004 }
3005 else
3006 {
3007 lineht.Set (lineep, i);
3008 }
3009 }
3010
3011 for (i = 1; i <= GetNLines(); i++)
3012 {
3013 STLLine* line = GetLine(i);
3014 for (int ii = 1; ii <= line->GetNS(); ii++)
3015 {
3016 int ap1, ap2;
3017 line->GetSeg(ii,ap1,ap2);
3018 // (*mycout) << "SEG " << p1 << " - " << p2 << endl;
3019 }
3020 }
3021
3022 PopStatus();
3023 }
3024
GetNOBodys()3025 int STLGeometry :: GetNOBodys()
3026 {
3027 int markedtrigs1 = 0;
3028 int starttrig = 1;
3029 int i, k, nnt;
3030 int bodycnt = 0;
3031
3032 NgArray<int> bodynum(GetNT());
3033
3034 for (i = 1; i <= GetNT(); i++)
3035 bodynum.Elem(i)=0;
3036
3037
3038 while (markedtrigs1 < GetNT())
3039 {
3040 for (i = starttrig; i <= GetNT(); i++)
3041 {
3042 if (!bodynum.Get(i))
3043 {
3044 starttrig = i;
3045 break;
3046 }
3047 }
3048 //add all triangles around starttriangle, which is reachable without going over an edge
3049 NgArray<int> todolist;
3050 NgArray<int> nextlist;
3051 bodycnt++;
3052 markedtrigs1++;
3053 bodynum.Elem(starttrig) = bodycnt;
3054 todolist.Append(starttrig);
3055
3056 while(todolist.Size())
3057 {
3058 for (i = 1; i <= todolist.Size(); i++)
3059 {
3060 //const STLTriangle& tt = GetTriangle(todolist.Get(i));
3061 for (k = 1; k <= NONeighbourTrigs(todolist.Get(i)); k++)
3062 {
3063 nnt = NeighbourTrig(todolist.Get(i),k);
3064 if (!bodynum.Get(nnt))
3065 {
3066 nextlist.Append(nnt);
3067 bodynum.Elem(nnt) = bodycnt;
3068 markedtrigs1++;
3069 }
3070 }
3071 }
3072
3073 todolist.SetSize(0);
3074 for (i = 1; i <= nextlist.Size(); i++)
3075 {
3076 todolist.Append(nextlist.Get(i));
3077 }
3078 nextlist.SetSize(0);
3079 }
3080 }
3081 PrintMessage(3, "Geometry has ", bodycnt, " separated bodys");
3082
3083 return bodycnt;
3084 }
3085
CalcFaceNums()3086 void STLGeometry :: CalcFaceNums()
3087 {
3088 int markedtrigs1 = 0;
3089 int starttrig(0);
3090 int laststarttrig = 1;
3091 int i, k, nnt;
3092 facecnt = 0;
3093
3094
3095 for (i = 1; i <= GetNT(); i++)
3096 GetTriangle(i).SetFaceNum(0);
3097
3098
3099 while (markedtrigs1 < GetNT())
3100 {
3101 for (i = laststarttrig; i <= GetNT(); i++)
3102 {
3103 if (!GetTriangle(i).GetFaceNum())
3104 {
3105 starttrig = i;
3106 laststarttrig = i;
3107 break;
3108 }
3109 }
3110 //add all triangles around starttriangle, which is reachable without going over an edge
3111 NgArray<int> todolist;
3112 NgArray<int> nextlist;
3113 facecnt++;
3114 markedtrigs1++;
3115 GetTriangle(starttrig).SetFaceNum(facecnt);
3116 todolist.Append(starttrig);
3117 STLPointId ap1, ap2;
3118
3119 while(todolist.Size())
3120 {
3121 for (i = 1; i <= todolist.Size(); i++)
3122 {
3123 const STLTriangle& tt = GetTriangle(todolist.Get(i));
3124 for (k = 1; k <= NONeighbourTrigs(todolist.Get(i)); k++)
3125 {
3126 nnt = NeighbourTrig(todolist.Get(i),k);
3127 STLTriangle& nt = GetTriangle(nnt);
3128 if (!nt.GetFaceNum())
3129 {
3130 tt.GetNeighbourPoints(nt,ap1,ap2);
3131 if (!IsEdge(ap1,ap2))
3132 {
3133 nextlist.Append(nnt);
3134 nt.SetFaceNum(facecnt);
3135 markedtrigs1++;
3136 }
3137 }
3138 }
3139 }
3140
3141 todolist.SetSize(0);
3142 for (i = 1; i <= nextlist.Size(); i++)
3143 {
3144 todolist.Append(nextlist.Get(i));
3145 }
3146 nextlist.SetSize(0);
3147 }
3148 }
3149 GetNOBodys();
3150 PrintMessage(3,"generated ", facecnt, " faces");
3151 }
3152
ClearSpiralPoints()3153 void STLGeometry :: ClearSpiralPoints()
3154 {
3155 spiralpoints.SetSize(GetNP());
3156 int i;
3157 for (i = 1; i <= spiralpoints.Size(); i++)
3158 {
3159 spiralpoints.Elem(i) = 0;
3160 }
3161 }
3162
3163
BuildSmoothEdges()3164 void STLGeometry :: BuildSmoothEdges ()
3165 {
3166 if (smoothedges) delete smoothedges;
3167
3168 smoothedges = new INDEX_2_HASHTABLE<int> (GetNE()/10 + 1);
3169
3170
3171 // Jack: Ok ?
3172 // UseExternalEdges();
3173
3174 PushStatusF("Build Smooth Edges");
3175
3176 int nt = GetNT();
3177 Vec3d ng1, ng2;
3178
3179 for (int i = 1; i <= nt; i++)
3180 {
3181 if (multithread.terminate)
3182 {PopStatus();return;}
3183
3184 SetThreadPercent(100.0 * (double)i / (double)nt);
3185
3186 const STLTriangle & trig = GetTriangle (i);
3187
3188 ng1 = trig.GeomNormal(points);
3189 ng1 /= (ng1.Length() + 1e-24);
3190
3191 for (int j = 1; j <= 3; j++)
3192 {
3193 int nbt = NeighbourTrig (i, j);
3194
3195 ng2 = GetTriangle(nbt).GeomNormal(points);
3196 ng2 /= (ng2.Length() + 1e-24);
3197
3198 STLPointId pi1, pi2;
3199 trig.GetNeighbourPoints(GetTriangle(nbt), pi1, pi2);
3200
3201 if (!IsEdge(pi1,pi2))
3202 {
3203 if (ng1 * ng2 < 0)
3204 {
3205 PrintMessage(7,"smoothedge found");
3206 INDEX_2 i2(pi1, pi2);
3207 i2.Sort();
3208 smoothedges->Set (i2, 1);
3209 }
3210 }
3211 }
3212 }
3213 PopStatus();
3214 }
3215
3216
3217
3218
3219
IsSmoothEdge(int pi1,int pi2) const3220 bool STLGeometry :: IsSmoothEdge (int pi1, int pi2) const
3221 {
3222 if (!smoothedges)
3223 return false;
3224 INDEX_2 i2(pi1, pi2);
3225 i2.Sort();
3226 return smoothedges->Used (i2);
3227 }
3228
3229
3230
3231 /*
3232 //function is not used now
3233 int IsInArray(int n, const NgArray<int>& ia)
3234 {
3235 int i;
3236 for (i = 1; i <= ia.Size(); i++)
3237 {
3238 if (ia.Get(i) == n) {return 1;}
3239 }
3240 return 0;
3241 }
3242 */
3243
AddConeAndSpiralEdges(const STLParameters & stlparam)3244 void STLGeometry :: AddConeAndSpiralEdges(const STLParameters& stlparam)
3245 {
3246 PrintMessage(5,"have now ", GetNE(), " edges with yellow angle = ", stlparam.yangle, " degree");
3247
3248 PrintFnStart("AddConeAndSpiralEdges");
3249
3250 // int i,j,k,n;
3251 // int changed = 0;
3252
3253 //check edges, where inner chart and no outer chart come together without an edge
3254 STLPointId np1, np2;
3255 int cnt = 0;
3256
3257 for (ChartId i = 1; i <= GetNOCharts(); i++)
3258 {
3259 STLChart& chart = GetChart(i);
3260 for (int j = 1; j <= chart.GetNChartT(); j++)
3261 {
3262 STLTrigId t = chart.GetChartTrig1(j);
3263 const STLTriangle& tt = GetTriangle(t);
3264
3265 for (int k = 1; k <= 3; k++)
3266 {
3267 STLTrigId nt = NeighbourTrig(t,k);
3268 if (GetChartNr(nt) != i && !TrigIsInOC(nt,i))
3269 {
3270 tt.GetNeighbourPoints(GetTriangle(nt),np1,np2);
3271 if (!IsEdge(np1,np2))
3272 {
3273 STLEdge se(np1,np2);
3274 se.SetLeftTrig(t);
3275 se.SetRightTrig(nt);
3276 int edgenum = AddEdge(se);
3277 AddEdgePP(np1,edgenum);
3278 AddEdgePP(np2,edgenum);
3279 //changed = 1;
3280 PrintWarning("Found a spiral like structure: chart=", int(i),
3281 ", trig=", int(t), ", p1=", int(np1), ", p2=", int(np2));
3282 cnt++;
3283 }
3284 }
3285 }
3286 }
3287
3288 }
3289
3290 PrintMessage(5, "found ", cnt, " spiral like structures");
3291 PrintMessage(5, "added ", cnt, " edges due to spiral like structures");
3292
3293 cnt = 0;
3294 int edgecnt = 0;
3295
3296 Array<STLTrigId> trigsaroundp;
3297 NgArray<int> chartpointchecked(GetNP()); //gets number of chart, if in this chart already checked
3298 chartpointchecked = 0;
3299
3300
3301 int onoc, notonoc, tpp, pn;
3302 STLPointId ap1, ap2;
3303 int tn1, tn2, l, problem;
3304
3305 if (!stldoctor.conecheck) {PrintWarning("++++++++++++ \ncone checking deactivated by user!!!!!\n+++++++++++++++"); return ;}
3306
3307 PushStatus("Find Critical Points");
3308
3309 int addedges = 0;
3310
3311 for (ChartId i = 1; i <= GetNOCharts(); i++)
3312 {
3313 SetThreadPercent((double)i/(double)GetNOCharts()*100.);
3314 if (multithread.terminate)
3315 {PopStatus();return;}
3316
3317 STLChart& chart = GetChart(i);
3318 for (int j = 1; j <= chart.GetNChartT(); j++)
3319 {
3320 STLTrigId t = chart.GetChartTrig1(j);
3321 const STLTriangle& tt = GetTriangle(t);
3322
3323 for (int k = 1; k <= 3; k++)
3324 {
3325 pn = tt.PNum(k);
3326 if (chartpointchecked.Get(pn) == i)
3327 {continue;}
3328
3329 int checkpoint = 0;
3330 for (int n = 1; n <= trigsperpoint.EntrySize(pn); n++)
3331 {
3332 if (trigsperpoint.Get(pn,n) != t &&
3333 GetChartNr(trigsperpoint.Get(pn,n)) != i &&
3334 !TrigIsInOC(trigsperpoint.Get(pn,n),i)) {checkpoint = 1;};
3335 }
3336 if (checkpoint)
3337 {
3338 chartpointchecked.Elem(pn) = i;
3339
3340 int worked = 0;
3341 int spworked = 0;
3342 GetSortedTrianglesAroundPoint(pn,t,trigsaroundp);
3343 trigsaroundp.Append(t);
3344
3345 problem = 0;
3346 for (int l = 2; l <= trigsaroundp.Size()-1; l++)
3347 {
3348 tn1 = trigsaroundp[l-2];
3349 tn2 = trigsaroundp[l-1];
3350 const STLTriangle& t1 = GetTriangle(tn1);
3351 const STLTriangle& t2 = GetTriangle(tn2);
3352 t1.GetNeighbourPoints(t2, ap1, ap2);
3353 if (IsEdge(ap1,ap2)) break;
3354
3355 if (GetChartNr(tn2) != i && !TrigIsInOC(tn2,i)) {problem = 1;}
3356 }
3357
3358 if (problem)
3359 {
3360 for (int l = 2; l <= trigsaroundp.Size()-1; l++)
3361 {
3362 tn1 = trigsaroundp[l-2];
3363 tn2 = trigsaroundp[l-1];
3364 const STLTriangle& t1 = GetTriangle(tn1);
3365 const STLTriangle& t2 = GetTriangle(tn2);
3366 t1.GetNeighbourPoints(t2, ap1, ap2);
3367 if (IsEdge(ap1,ap2)) break;
3368
3369 if ((GetChartNr(tn1) == i && GetChartNr(tn2) != i && TrigIsInOC(tn2,i)) ||
3370 (GetChartNr(tn2) == i && GetChartNr(tn1) != i && TrigIsInOC(tn1,i)))
3371 {
3372 if (addedges || !GetNEPP(pn))
3373 {
3374 STLEdge se(ap1,ap2);
3375 se.SetLeftTrig(tn1);
3376 se.SetRightTrig(tn2);
3377 int edgenum = AddEdge(se);
3378 AddEdgePP(ap1,edgenum);
3379 AddEdgePP(ap2,edgenum);
3380 edgecnt++;
3381 }
3382 if (!addedges && !GetSpiralPoint(pn))
3383 {
3384 SetSpiralPoint(pn);
3385 spworked = 1;
3386 }
3387 worked = 1;
3388 }
3389 }
3390 }
3391 //backwards:
3392 problem = 0;
3393 for (int l = trigsaroundp.Size()-1; l >= 2; l--)
3394 {
3395 tn1 = trigsaroundp[l];
3396 tn2 = trigsaroundp[l-1];
3397 const STLTriangle& t1 = GetTriangle(tn1);
3398 const STLTriangle& t2 = GetTriangle(tn2);
3399 t1.GetNeighbourPoints(t2, ap1, ap2);
3400 if (IsEdge(ap1,ap2)) break;
3401
3402 if (GetChartNr(tn2) != i && !TrigIsInOC(tn2,i)) {problem = 1;}
3403 }
3404 if (problem)
3405 for (int l = trigsaroundp.Size()-1; l >= 2; l--)
3406 {
3407 tn1 = trigsaroundp[l];
3408 tn2 = trigsaroundp[l-1];
3409 const STLTriangle& t1 = GetTriangle(tn1);
3410 const STLTriangle& t2 = GetTriangle(tn2);
3411 t1.GetNeighbourPoints(t2, ap1, ap2);
3412 if (IsEdge(ap1,ap2)) break;
3413
3414 if ((GetChartNr(tn1) == i && GetChartNr(tn2) != i && TrigIsInOC(tn2,i)) ||
3415 (GetChartNr(tn2) == i && GetChartNr(tn1) != i && TrigIsInOC(tn1,i)))
3416 {
3417 if (addedges || !GetNEPP(pn))
3418 {
3419 STLEdge se(ap1,ap2);
3420 se.SetLeftTrig(tn1);
3421 se.SetRightTrig(tn2);
3422 int edgenum = AddEdge(se);
3423 AddEdgePP(ap1,edgenum);
3424 AddEdgePP(ap2,edgenum);
3425 edgecnt++;
3426 }
3427 if (!addedges && !GetSpiralPoint(pn))
3428 {
3429 SetSpiralPoint(pn);
3430 spworked = 1;
3431 //if (GetNEPP(pn) == 0) {(*mycout) << "ERROR: spiralpoint with no edge found!" << endl;}
3432 }
3433 worked = 1;
3434 }
3435 }
3436
3437 if (worked)
3438 {
3439 //(*testout) << "set edgepoint due to spirals: pn=" << i << endl;
3440 SetLineEndPoint(pn);
3441 }
3442 if (spworked)
3443 {
3444 /*
3445 (*mycout) << "Warning: Critical Point " << tt.PNum(k)
3446 << "( chart " << i << ", trig " << t
3447 << ") has been neutralized!!!" << endl;
3448 */
3449 cnt++;
3450 }
3451 // markedpoints.Elem(tt.PNum(k)) = 1;
3452 }
3453 }
3454 }
3455 }
3456 PrintMessage(5, "found ", cnt, " critical points!");
3457 PrintMessage(5, "added ", edgecnt, " edges due to critical points!");
3458
3459 PopStatus();
3460
3461 //search points where inner chart and outer chart and "no chart" trig come together at edge-point
3462
3463 PrintMessage(7,"search for special chart points");
3464 for (ChartId i = 1; i <= GetNOCharts(); i++)
3465 {
3466 STLChart& chart = GetChart(i);
3467 for (int j = 1; j <= chart.GetNChartT(); j++)
3468 {
3469 STLTrigId t = chart.GetChartTrig1(j);
3470 const STLTriangle& tt = GetTriangle(t);
3471
3472 for (int k = 1; k <= 3; k++)
3473 {
3474 pn = tt.PNum(k);
3475 if (GetNEPP(pn) == 2)
3476 {
3477 onoc = 0;
3478 notonoc = 0;
3479 for (int n = 1; n <= trigsperpoint.EntrySize(pn); n++)
3480 {
3481 tpp = trigsperpoint.Get(pn,n);
3482 if (tpp != t && GetChartNr(tpp) != i)
3483 {
3484 if (TrigIsInOC(tpp,i)) {onoc = 1;}
3485 if (!TrigIsInOC(tpp,i)) {notonoc = 1;}
3486 }
3487 }
3488 if (onoc && notonoc && !IsLineEndPoint(pn))
3489 {
3490 GetSortedTrianglesAroundPoint(pn,t,trigsaroundp);
3491 int here = 1; //we start on this side of edge, !here = there
3492 int thereOC = 0;
3493 int thereNotOC = 0;
3494 for (l = 2; l <= trigsaroundp.Size(); l++)
3495 {
3496 GetTriangle(trigsaroundp[l-2]).
3497 GetNeighbourPoints(GetTriangle(trigsaroundp[l-1]), ap1, ap2);
3498 if (IsEdge(ap1,ap2)) {here = (here+1)%2;}
3499 if (!here && TrigIsInOC(trigsaroundp[l-1],i)) {thereOC = 1;}
3500 if (!here && !TrigIsInOC(trigsaroundp[l-1],i)) {thereNotOC = 1;}
3501 }
3502 if (thereOC && thereNotOC)
3503 {
3504 //(*mycout) << "Special OCICnotC - point " << pn << " found!" << endl;
3505 //(*testout) << "set edgepoint due to spirals: pn=" << i << endl;
3506 SetLineEndPoint(pn);
3507 }
3508 }
3509 }
3510 }
3511 }
3512 }
3513 PrintMessage(5,"have now ", GetNE(), " edges with yellow angle = ", stlparam.yangle, " degree");
3514 }
3515
3516 //get trigs at a point, started with starttrig, then every left
GetSortedTrianglesAroundPoint(STLPointId p,STLTrigId starttrig,Array<STLTrigId> & trigs)3517 void STLGeometry :: GetSortedTrianglesAroundPoint(STLPointId p, STLTrigId starttrig, Array<STLTrigId>& trigs)
3518 {
3519 STLTrigId acttrig = starttrig;
3520 trigs.SetAllocSize(trigsperpoint.EntrySize(p));
3521 trigs.SetSize(0);
3522 trigs.Append(acttrig);
3523 int locindex1(0), locindex2(0);
3524 //(*mycout) << "trigs around point " << p << endl;
3525
3526 int end = 0;
3527 while (!end)
3528 {
3529 const STLTriangle& at = GetTriangle(acttrig);
3530 for (int i = 1; i <= trigsperpoint.EntrySize(p); i++)
3531 {
3532 STLTrigId t = trigsperpoint.Get(p,i);
3533 const STLTriangle& nt = GetTriangle(t);
3534 if (at.IsNeighbourFrom(nt))
3535 {
3536 STLPointId ap1, ap2;
3537 at.GetNeighbourPoints(nt, ap1, ap2);
3538 if (ap2 == p) {Swap(ap1,ap2);}
3539 if (ap1 != p) {PrintSysError("In GetSortedTrianglesAroundPoint!!!");}
3540
3541 for (int j = 1; j <= 3; j++)
3542 {
3543 if (at.PNum(j) == ap1) {locindex1 = j;};
3544 if (at.PNum(j) == ap2) {locindex2 = j;};
3545 }
3546 if ((locindex2+1)%3+1 == locindex1)
3547 {
3548 if (t != starttrig)
3549 {
3550 trigs.Append(t);
3551 // (*mycout) << "trig " << t << endl;
3552 acttrig = t;
3553 }
3554 else
3555 {
3556 end = 1;
3557 }
3558 break;
3559 }
3560 }
3561 }
3562 }
3563
3564 }
3565
3566 /*
3567 int STLGeometry :: NeighbourTrig(int trig, int nr) const
3568 {
3569 return neighbourtrigs.Get(trig,nr);
3570 }
3571 */
3572
3573
3574
SmoothGeometry()3575 void STLGeometry :: SmoothGeometry ()
3576 {
3577 int i, j, k;
3578
3579 double maxerr0, maxerr;
3580
3581 for (i = 1; i <= GetNP(); i++)
3582 {
3583 if (GetNEPP(i)) continue;
3584
3585 maxerr0 = 0;
3586 for (j = 1; j <= NOTrigsPerPoint(i); j++)
3587 {
3588 int tnum = TrigPerPoint(i, j);
3589 double err = Angle (GetTriangle(tnum).Normal (),
3590 GetTriangle(tnum).GeomNormal(GetPoints()));
3591 if (err > maxerr0)
3592 maxerr0 = err;
3593 }
3594
3595 Point3d pi = GetPoint (i);
3596 if (maxerr0 < 1.1) continue; // about 60 degree
3597
3598 maxerr0 /= 2; // should be at least halfen
3599
3600 for (k = 1; k <= NOTrigsPerPoint(i); k++)
3601 {
3602 const STLTriangle & trig = GetTriangle (TrigPerPoint (i, k));
3603 Point3d c = Center(GetPoint (trig.PNum(1)),
3604 GetPoint (trig.PNum(2)),
3605 GetPoint (trig.PNum(3)));
3606
3607 Point3d np = pi + 0.1 * (c - pi);
3608 SetPoint (i, np);
3609
3610 maxerr = 0;
3611 for (j = 1; j <= NOTrigsPerPoint(i); j++)
3612 {
3613 int tnum = TrigPerPoint(i, j);
3614 double err = Angle (GetTriangle(tnum).Normal (),
3615 GetTriangle(tnum).GeomNormal(GetPoints()));
3616 if (err > maxerr)
3617 maxerr = err;
3618 }
3619
3620 if (maxerr < maxerr0)
3621 {
3622 pi = np;
3623 }
3624 }
3625
3626 SetPoint (i, pi);
3627 }
3628 }
3629
WriteChartToFile(ChartId chartnumber,string filename)3630 void STLGeometry :: WriteChartToFile( ChartId chartnumber, string filename )
3631 {
3632 PrintMessage(1,"write chart ", int(chartnumber), " to ", filename);
3633 Array<int> trignums;
3634
3635 if (chartnumber >= 1 && chartnumber <= GetNOCharts())
3636 {
3637 const STLChart& chart = GetChart(chartnumber);
3638
3639 for (int j = 1; j <= chart.GetNChartT(); j++)
3640 trignums.Append(chart.GetChartTrig1(j));
3641
3642 for (int j = 1; j <= chart.GetNOuterT(); j++)
3643 trignums.Append(chart.GetOuterTrig1(j));
3644
3645 QuickSort(trignums);
3646 STLGeometry geo;
3647 NgArray<STLReadTriangle> readtrigs;
3648 const auto & first_trig = GetTriangle(chart.GetChartTrig1(1));
3649 auto normal = first_trig.Normal();
3650 Box<3> box{Box<3>::EMPTY_BOX};
3651
3652 for(auto j : trignums)
3653 {
3654 const auto& trig = GetTriangle(j);
3655 Point<3> pts[3];
3656 for(auto k : Range(3))
3657 {
3658 pts[k] = GetPoint(trig[k]);
3659 box.Add(pts[k]);
3660 }
3661 Vec3d normal = Cross( pts[1]-pts[0], pts[2]-pts[0] );
3662 readtrigs.Append(STLReadTriangle(pts, trig.Normal()));
3663 }
3664 auto dist = box.PMax() - box.PMin();
3665 auto extra_point = GetPoint(first_trig[0]) - dist.Length()*normal;
3666
3667 NgArray<int> acttrigs(GetNT());
3668 acttrigs = -1;
3669 for (int j = 1; j <= chart.GetNT(); j++)
3670 acttrigs.Elem(chart.GetTrig1(j)) = chartnumber;
3671
3672 for (int j = 1; j <= chart.GetNT(); j++)
3673 {
3674 auto t = chart.GetTrig1(j);
3675 const auto & tt = GetTriangle(t);
3676 for (int k = 1; k <= 3; k++)
3677 {
3678 int nt = NeighbourTrig(t,k);
3679 if (acttrigs.Get(nt) != chartnumber)
3680 {
3681 STLPointId np1, np2;
3682 tt.GetNeighbourPoints(GetTriangle(nt),np1,np2);
3683
3684 Point<3> pts[3];
3685 pts[0] = GetPoint(np2);
3686 pts[1] = GetPoint(np1);
3687 pts[2] = extra_point;
3688 Vec3d normal = -Cross( pts[2]-pts[0], pts[1]-pts[0] );
3689 readtrigs.Append(STLReadTriangle(pts, normal));
3690 }
3691 }
3692 }
3693
3694 geo.InitSTLGeometry(readtrigs);
3695 geo.Save(filename);
3696 }
3697 }
3698
3699
3700
3701 class STLGeometryRegister : public GeometryRegister
3702 {
3703 public:
3704 virtual NetgenGeometry * Load (string filename) const;
3705 };
3706
Load(string filename) const3707 NetgenGeometry * STLGeometryRegister :: Load (string filename) const
3708 {
3709 const char * cfilename = filename.c_str();
3710
3711 if (strcmp (&cfilename[strlen(cfilename)-3], "stl") == 0)
3712 {
3713 PrintMessage (1, "Load STL geometry file ", cfilename);
3714
3715 ifstream infile(cfilename);
3716
3717 STLGeometry * hgeom = STLGeometry :: Load (infile);
3718 hgeom -> edgesfound = 0;
3719 return hgeom;
3720 }
3721 else if (strcmp (&cfilename[strlen(cfilename)-4], "stlb") == 0)
3722 {
3723 PrintMessage (1, "Load STL binary geometry file ", cfilename);
3724
3725 ifstream infile(cfilename);
3726
3727 STLGeometry * hgeom = STLGeometry :: LoadBinary (infile);
3728 hgeom -> edgesfound = 0;
3729 return hgeom;
3730 }
3731 else if (strcmp (&cfilename[strlen(cfilename)-3], "nao") == 0)
3732 {
3733 PrintMessage (1, "Load naomi (F. Kickinger) geometry file ", cfilename);
3734
3735 ifstream infile(cfilename);
3736
3737 STLGeometry * hgeom = STLGeometry :: LoadNaomi (infile);
3738 hgeom -> edgesfound = 0;
3739 return hgeom;
3740 }
3741
3742
3743 return NULL;
3744 }
3745
3746
3747 class STLInit
3748 {
3749 public:
STLInit()3750 STLInit()
3751 {
3752 geometryregister.Append (new STLGeometryRegister);
3753 }
3754 };
3755
3756 STLInit stlinit;
3757
3758 static RegisterClassForArchive<STLGeometry, NetgenGeometry, STLTopology> stlgeo;
3759 }
3760