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