1 /*
2 This file is a part of the RepSnapper project.
3 Copyright (C) 2010 Kulitorum
4 Copyright (C) 2011-12 martin.dieringer@gmx.de
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License along
17 with this program; if not, write to the Free Software Foundation, Inc.,
18 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
19 */
20
21 #include "shape.h"
22 #include "files.h"
23 #include "ui/progress.h"
24 #include "settings.h"
25 #include "clipping.h"
26 #include "render.h"
27
28 #ifdef _OPENMP
29 #include <omp.h>
30 #endif
31
32 // Constructor
Shape()33 Shape::Shape()
34 : slow_drawing(false), gl_List(-1)
35 {
36 Min.set(0,0,0);
37 Max.set(200,200,200);
38 CalcBBox();
39 }
40
clear()41 void Shape::clear() {
42 triangles.clear();
43 if (gl_List>=0)
44 glDeleteLists(gl_List,1);
45 gl_List = -1;
46 };
47
setTriangles(const vector<Triangle> & triangles_)48 void Shape::setTriangles(const vector<Triangle> &triangles_)
49 {
50 triangles = triangles_;
51
52 CalcBBox();
53 double vol = volume();
54 if (vol < 0) {
55 invertNormals();
56 vol = -vol;
57 }
58
59 //PlaceOnPlatform();
60 cerr << _("Shape has volume ") << volume() << _(" mm^3 and ")
61 << triangles.size() << _(" triangles") << endl;
62 }
63
64
saveBinarySTL(Glib::ustring filename) const65 int Shape::saveBinarySTL(Glib::ustring filename) const
66 {
67 if (!File::saveBinarySTL(filename, triangles, transform3D.transform))
68 return -1;
69 return 0;
70
71 }
72
hasAdjacentTriangleTo(const Triangle & triangle,double sqdistance) const73 bool Shape::hasAdjacentTriangleTo(const Triangle &triangle, double sqdistance) const
74 {
75 bool haveadj = false;
76 int count = (int)triangles.size();
77 #ifdef _OPENMP
78 #pragma omp parallel for schedule(dynamic)
79 #endif
80 for (int i = 0; i < count; i++)
81 if (!haveadj)
82 if (triangle.isConnectedTo(triangles[i],sqdistance)) {
83 haveadj = true;
84 }
85 return haveadj;
86 }
87
88
89 // recursively build a list of triangles for a shape
addtoshape(uint i,const vector<vector<uint>> & adj,vector<uint> & tr,vector<bool> & done)90 void addtoshape(uint i, const vector< vector<uint> > &adj,
91 vector<uint> &tr, vector<bool> &done)
92 {
93 if (!done[i]) {
94 // add this index to tr
95 tr.push_back(i);
96 done[i] = true;
97 for (uint j = 0; j < adj[i].size(); j++) {
98 if (adj[i][j]!=i)
99 addtoshape(adj[i][j], adj, tr, done);
100 }
101 }
102 // we have a complete list of adjacent triangles indices
103 }
104
splitshapes(vector<Shape * > & shapes,ViewProgress * progress)105 void Shape::splitshapes(vector<Shape*> &shapes, ViewProgress *progress)
106 {
107 int n_tr = (int)triangles.size();
108 if (progress) progress->start(_("Split Shapes"), n_tr);
109 int progress_steps = max(1,(int)(n_tr/100));
110 vector<bool> done(n_tr);
111 bool cont = true;
112 // make list of adjacent triangles for each triangle
113 vector< vector<uint> > adj(n_tr);
114 if (progress) progress->set_label(_("Split: Sorting Triangles ..."));
115 #ifdef _OPENMP
116 omp_lock_t progress_lock;
117 omp_init_lock(&progress_lock);
118 #pragma omp parallel for schedule(dynamic)
119 #endif
120 for (int i = 0; i < n_tr; i++) {
121 if (progress && i%progress_steps==0) {
122 #ifdef _OPENMP
123 omp_set_lock(&progress_lock);
124 #endif
125 cont = progress->update(i);
126 #ifdef _OPENMP
127 omp_unset_lock(&progress_lock);
128 #endif
129 }
130 vector<uint> trv;
131 for (int j = 0; j < n_tr; j++) {
132 if (i!=j) {
133 bool add = false;
134 if (j<i) // maybe(!) we have it already
135 for (uint k = 0; k<adj[j].size(); k++) {
136 if ((int)adj[j][k] == i) {
137 add = true; break;
138 }
139 }
140 add |= (triangles[i].isConnectedTo(triangles[j], 0.01));
141 if (add) trv.push_back(j);
142 }
143 }
144 adj[i] = trv;
145 if (!cont) i=n_tr;
146 }
147
148 if (progress) progress->set_label(_("Split: Building shapes ..."));
149
150
151 // triangle indices of shapes
152 vector< vector<uint> > shape_tri;
153
154 for (int i = 0; i < n_tr; i++) done[i] = false;
155 for (int i = 0; i < n_tr; i++) {
156 if (progress && i%progress_steps==0)
157 cont = progress->update(i);
158 if (!done[i]){
159 cerr << _("Shape ") << shapes.size()+1 << endl;
160 vector<uint> current;
161 addtoshape(i, adj, current, done);
162 Shape *shape = new Shape();
163 shapes.push_back(shape);
164 shapes.back()->triangles.resize(current.size());
165 for (uint i = 0; i < current.size(); i++)
166 shapes.back()->triangles[i] = triangles[current[i]];
167 shapes.back()->CalcBBox();
168 }
169 if (!cont) i=n_tr;
170 }
171
172 if (progress) progress->stop("_(Done)");
173 }
174
cube(Vector3d Min,Vector3d Max)175 vector<Triangle> cube(Vector3d Min, Vector3d Max)
176 {
177 vector<Triangle> c;
178 const Vector3d diag = Max-Min;
179 const Vector3d dX(diag.x(),0,0);
180 const Vector3d dY(0,diag.y(),0);
181 const Vector3d dZ(0,0,diag.z());
182 // front
183 c.push_back(Triangle(Min, Min+dX, Min+dX+dZ));
184 c.push_back(Triangle(Min, Min+dX+dZ, Min+dZ));
185 // back
186 c.push_back(Triangle(Min+dY, Min+dY+dX+dZ, Min+dY+dX));
187 c.push_back(Triangle(Min+dY, Min+dY+dZ, Min+dY+dX+dZ));
188 // left
189 c.push_back(Triangle(Min, Min+dZ, Min+dY+dZ));
190 c.push_back(Triangle(Min, Min+dY+dZ, Min+dY));
191 // right
192 c.push_back(Triangle(Min+dX, Min+dX+dY+dZ, Min+dX+dZ));
193 c.push_back(Triangle(Min+dX, Min+dX+dY, Min+dX+dY+dZ));
194 // bottom
195 c.push_back(Triangle(Min, Min+dX+dY, Min+dX));
196 c.push_back(Triangle(Min, Min+dY, Min+dX+dY));
197 // top
198 c.push_back(Triangle(Min+dZ, Min+dZ+dX, Min+dZ+dX+dY));
199 c.push_back(Triangle(Min+dZ, Min+dZ+dX+dY, Min+dZ+dY));
200 return c;
201 }
202
makeHollow(double wallthickness)203 void Shape::makeHollow(double wallthickness)
204 {
205 invertNormals();
206 const Vector3d wall(wallthickness,wallthickness,wallthickness);
207 Matrix4d invT = transform3D.getInverse();
208 vector<Triangle> cubet = cube(invT*Min-wall, invT*Max+wall);
209 triangles.insert(triangles.end(),cubet.begin(),cubet.end());
210 CalcBBox();
211 }
212
invertNormals()213 void Shape::invertNormals()
214 {
215 for (uint i = 0; i < triangles.size(); i++)
216 triangles[i].invertNormal();
217 }
218
219 // doesn't work
repairNormals(double sqdistance)220 void Shape::repairNormals(double sqdistance)
221 {
222 for (uint i = 0; i < triangles.size(); i++) {
223 vector<uint> adjacent;
224 uint numadj=0, numwrong=0;
225 for (uint j = i+1; j < triangles.size(); j++) {
226 if (i!=j) {
227 if (triangles[i].isConnectedTo(triangles[j], sqdistance)) {
228 numadj++;
229 if (triangles[i].wrongOrientationWith(triangles[j], sqdistance)) {
230 numwrong++;
231 triangles[j].invertNormal();
232 }
233 }
234 }
235 }
236 //cerr << i<< ": " << numadj << " - " << numwrong << endl;
237 //if (numwrong > numadj/2) triangles[i].invertNormal();
238 }
239 }
240
mirror()241 void Shape::mirror()
242 {
243 const Vector3d mCenter = transform3D.getInverse() * Center;
244 for (uint i = 0; i < triangles.size(); i++)
245 triangles[i].mirrorX(mCenter);
246 CalcBBox();
247 }
248
volume() const249 double Shape::volume() const
250 {
251 double vol=0;
252 for (uint i = 0; i < triangles.size(); i++)
253 vol+=triangles[i].projectedvolume(transform3D.transform);
254 return vol;
255 }
256
getSTLsolid() const257 string Shape::getSTLsolid() const
258 {
259 stringstream sstr;
260 sstr << "solid " << filename <<endl;
261 for (uint i = 0; i < triangles.size(); i++)
262 sstr << triangles[i].getSTLfacet(transform3D.transform);
263 sstr << "endsolid " << filename <<endl;
264 return sstr.str();
265 }
266
addTriangles(const vector<Triangle> & tr)267 void Shape::addTriangles(const vector<Triangle> &tr)
268 {
269 triangles.insert(triangles.end(), tr.begin(), tr.end());
270 CalcBBox();
271 }
272
getTriangles(const Matrix4d & T) const273 vector<Triangle> Shape::getTriangles(const Matrix4d &T) const
274 {
275 vector<Triangle> tr(triangles.size());
276 for (uint i = 0; i < triangles.size(); i++) {
277 tr[i] = triangles[i].transformed(T*transform3D.transform);
278 }
279 return tr;
280 }
281
282
trianglesSteeperThan(double angle) const283 vector<Triangle> Shape::trianglesSteeperThan(double angle) const
284 {
285 vector<Triangle> tr;
286 for (uint i = 0; i < triangles.size(); i++) {
287 // negative angles are triangles facing downwards
288 const double tangle = -triangles[i].slopeAngle(transform3D.transform);
289 if (tangle >= angle)
290 tr.push_back(triangles[i]);
291 }
292 return tr;
293 }
294
295
FitToVolume(const Vector3d & vol)296 void Shape::FitToVolume(const Vector3d &vol)
297 {
298 Vector3d diag = Max-Min;
299 const double sc_x = diag.x() / vol.x();
300 const double sc_y = diag.y() / vol.y();
301 const double sc_z = diag.z() / vol.z();
302 double max_sc = max(max(sc_x, sc_y),sc_z);
303 if (max_sc > 1.)
304 Scale(1./max_sc, true);
305 }
306
Scale(double in_scale_factor,bool calcbbox)307 void Shape::Scale(double in_scale_factor, bool calcbbox)
308 {
309 transform3D.move(-Center);
310 transform3D.scale(in_scale_factor);
311 transform3D.move(Center);
312 if (calcbbox)
313 CalcBBox();
314 }
315
ScaleX(double x)316 void Shape::ScaleX(double x)
317 {
318 transform3D.move(-Center);
319 transform3D.scale_x(x);
320 transform3D.move(Center);
321 }
ScaleY(double x)322 void Shape::ScaleY(double x)
323 {
324 transform3D.move(-Center);
325 transform3D.scale_y(x);
326 transform3D.move(Center);
327 }
ScaleZ(double x)328 void Shape::ScaleZ(double x)
329 {
330 transform3D.move(-Center);
331 transform3D.scale_z(x);
332 transform3D.move(Center);
333 }
334
CalcBBox()335 void Shape::CalcBBox()
336 {
337 Min.set(INFTY,INFTY,INFTY);
338 Max.set(-INFTY,-INFTY,-INFTY);
339 for(size_t i = 0; i < triangles.size(); i++) {
340 triangles[i].AccumulateMinMax (Min, Max, transform3D.transform);
341 }
342 Center = (Max + Min) / 2;
343 if (gl_List>=0)
344 glDeleteLists(gl_List,1);
345 gl_List = -1;
346 }
347
scaledCenter() const348 Vector3d Shape::scaledCenter() const
349 {
350 return Center * transform3D.get_scale();
351 }
352
353 struct SNorm {
354 Vector3d normal;
355 double area;
operator <SNorm356 bool operator<(const SNorm &other) const {return (area<other.area);};
357 } ;
358
getMostUsedNormals() const359 vector<Vector3d> Shape::getMostUsedNormals() const
360 {
361 vector<struct SNorm> normals;
362 // vector<Vector3d> normals;
363 // vector<double> area;
364 uint ntr = triangles.size();
365 vector<bool> done(ntr);
366 normals.reserve(ntr);
367 for(size_t i=0;i<ntr;i++)
368 {
369 bool havenormal = false;
370 int numnorm = normals.size();
371 #ifdef _OPENMP
372 #pragma omp parallel for schedule(dynamic)
373 #endif
374 for (int n = 0; n < numnorm; n++) {
375 if ( (normals[n].normal -
376 triangles[i].transformed(transform3D.transform).Normal)
377 .squared_length() < 0.000001) {
378 havenormal = true;
379 normals[n].area += triangles[i].area();
380 done[i] = true;
381 }
382 }
383 if (!havenormal){
384 SNorm n;
385 n.normal = triangles[i].transformed(transform3D.transform).Normal;
386 n.area = triangles[i].area();
387 normals.push_back(n);
388 done[i] = true;
389 }
390 }
391 std::sort(normals.rbegin(),normals.rend());
392 //cerr << normals.size() << endl;
393 vector<Vector3d> nv(normals.size());
394 for (uint n=0; n < normals.size(); n++) nv[n] = normals[n].normal;
395 return nv;
396 }
397
398
OptimizeRotation()399 void Shape::OptimizeRotation()
400 {
401 // CenterAroundXY();
402 vector<Vector3d> normals = getMostUsedNormals();
403 // cycle through most-used normals?
404
405 Vector3d N;
406 Vector3d Z(0,0,-1);
407 double angle=0;
408 int count = (int)normals.size();
409 for (int n=0; n < count; n++) {
410 //cerr << n << normals[n] << endl;
411 N = normals[n];
412 angle = acos(N.dot(Z));
413 if (angle>0) {
414 Vector3d axis = N.cross(Z);
415 if (axis.squared_length()>0.1) {
416 Rotate(axis,angle);
417 break;
418 }
419 }
420 }
421 CalcBBox();
422 PlaceOnPlatform();
423 }
424
divideAtZ(double z,Shape * upper,Shape * lower,const Matrix4d & T) const425 int Shape::divideAtZ(double z, Shape *upper, Shape *lower, const Matrix4d &T) const
426 {
427 vector<Poly> polys;
428 vector<Poly> supportpolys;
429 double max_grad;
430 bool ok = getPolygonsAtZ(T, z, polys, max_grad, supportpolys, -1);
431 if (!ok) return 0;
432 vector< vector<Triangle> > surfs;
433 triangulate(polys, surfs);
434
435 vector<Triangle> surf;
436 for (uint i=0; i<surfs.size(); i++)
437 surf.insert(surf.end(), surfs[i].begin(), surfs[i].end());
438
439 lower->triangles.insert(lower->triangles.end(),surf.begin(),surf.end());
440 for (guint i=0; i<surf.size(); i++) surf[i].invertNormal();
441 upper->triangles.insert(upper->triangles.end(),surf.begin(),surf.end());
442 vector<Triangle> toboth;
443 for (guint i=0; i< triangles.size(); i++) {
444 Triangle tt = triangles[i].transformed(T*transform3D.transform);
445 if (tt.A.z() < z && tt.B.z() < z && tt.C.z() < z )
446 lower->triangles.push_back(tt);
447 else if (tt.A.z() > z && tt.B.z() > z && tt.C.z() > z )
448 upper->triangles.push_back(tt);
449 else
450 toboth.push_back(tt);
451 }
452 vector<Triangle> uppersplit,lowersplit;
453 for (guint i=0; i< toboth.size(); i++) {
454 toboth[i].SplitAtPlane(z, uppersplit, lowersplit);
455 }
456 upper->triangles.insert(upper->triangles.end(),
457 uppersplit.begin(),uppersplit.end());
458 lower->triangles.insert(lower->triangles.end(),
459 lowersplit.begin(),lowersplit.end());
460 upper->CalcBBox();
461 lower->CalcBBox();
462 lower->Rotate(Vector3d(0,1,0),M_PI);
463 upper->move(Vector3d(10+Max.x()-Min.x(),0,0));
464 lower->move(Vector3d(2*(10+Max.x()-Min.x()),0,0));
465 upper->PlaceOnPlatform();
466 lower->PlaceOnPlatform();
467 return 2;
468 }
469
PlaceOnPlatform()470 void Shape::PlaceOnPlatform()
471 {
472 transform3D.move(Vector3d(0,0,-Min.z()));
473 }
474
475 // Rotate and adjust for the user - not a pure rotation by any means
Rotate(const Vector3d & axis,const double & angle)476 void Shape::Rotate(const Vector3d & axis, const double & angle)
477 {
478 transform3D.rotate(Center, axis, angle);
479 return;
480 // CenterAroundXY();
481 // // do a real rotation because matrix transform gives errors when slicing
482 // int count = (int)triangles.size();
483 // #ifdef _OPENMP
484 // #pragma omp parallel for schedule(dynamic)
485 // #endif
486 // for (int i=0; i < count ; i++)
487 // {
488 // triangles[i].rotate(axis, angle);
489 // }
490 // PlaceOnPlatform();
491 }
492
493 // this is primitive, it just rotates triangle vertices
Twist(double angle)494 void Shape::Twist(double angle)
495 {
496 CalcBBox();
497 double h = Max.z()-Min.z();
498 double hangle=0;
499 Vector3d axis(0,0,1);
500 int count = (int)triangles.size();
501 #ifdef _OPENMP
502 #pragma omp parallel for schedule(dynamic)
503 #endif
504 for (int i=0; i<count; i++) {
505 for (size_t j=0; j<3; j++)
506 {
507 hangle = angle * (triangles[i][j].z() - Min.z()) / h;
508 triangles[i][j] = triangles[i][j].rotate(hangle,axis);
509 }
510 triangles[i].calcNormal();
511 }
512 CalcBBox();
513 }
514
515 // void Shape::CenterAroundXY()
516 // {
517 // CalcBBox();
518 // return;
519
520 // /* // this moves all triangles
521 // Vector3d displacement = transform3D.getTranslation() - Center;
522 // int count = (int)triangles.size();
523 // #ifdef _OPENMP
524 // #pragma omp parallel for schedule(dynamic)
525 // #endif
526 // for(int i=0; i<count ; i++)
527 // {
528 // triangles[i].Translate(displacement);
529 // }
530 // transform3D.move(-displacement);
531 // */
532
533 // //cerr << "DISPL " << displacement << endl;
534 // //CalcBBox();
535 // // Min -= displacement;
536 // // Max -= displacement;
537 // // Center -= displacement;
538 // }
539
540 /*
541 Poly Shape::getOutline(const Matrix4d &T, double maxlen) const
542 {
543 Matrix4d transform = T * transform3D.transform ;
544 vector<Vector2d> points(triangles.size()*3);
545 for (uint i = 0; i<triangles.size(); i++) {
546 for (uint j = 0; j<3; j++) {
547 points[i*3 + j] = Vector2d(triangles[i][j].x(),triangles[i][j].y());
548 }
549 }
550 Poly hull = concaveHull2D(points, maxlen);
551 return hull;
552 }
553 */
554
getLineSequences(const vector<Segment> lines,vector<vector<uint>> & connectedlines)555 bool getLineSequences(const vector<Segment> lines, vector< vector<uint> > &connectedlines)
556 {
557 uint nlines = lines.size();
558 //cerr << "lines size " << nlines << endl;
559 if (nlines==0) return true;
560 vector<bool> linedone(nlines);
561 for (uint l=0; l < nlines; l++) linedone[l] = false;
562 vector<uint> sequence;
563 uint donelines = 0;
564 vector<uint> connections;
565 while (donelines < nlines) {
566 connections.clear();
567 for (uint l=0; l < nlines; l++) { // add next connecting line
568 if ( !linedone[l] &&
569 ( (sequence.size()==0) ||
570 (lines[l].start == lines[sequence.back()].end) ) ) {
571 connections.push_back(l);
572 //cerr << "found connection" << endl;
573 }
574 }
575 if (connections.size()>0) {
576 //cerr << "found " << connections.size() << " connections" << endl;
577 sequence.push_back(connections[0]);
578 linedone[connections[0]] =true;
579 donelines++;
580 if (lines[sequence.front()].start == lines[sequence.back()].end) {
581 //cerr << "closed sequence" << endl;
582 connectedlines.push_back(sequence);
583 sequence.clear();
584 }
585 } else { // (!foundconnection) { // new sequence
586 //cerr << "sequence size " << sequence.size() << endl;
587 connectedlines.push_back(sequence);
588 sequence.clear();
589 for (uint l=0; l < nlines; l++) { // add next best undone line
590 if (!linedone[l]) {
591 sequence.push_back(l);
592 linedone[l] = true;
593 donelines++;
594 break;
595 }
596 }
597 }
598 }
599 if (sequence.size()>0)
600 connectedlines.push_back(sequence);
601 //cerr << "found "<< connectedlines.size() << " sequences" << endl;
602 return true;
603 }
604
getPolygonsAtZ(const Matrix4d & T,double z,vector<Poly> & polys,double & max_gradient,vector<Poly> & supportpolys,double max_supportangle,double thickness) const605 bool Shape::getPolygonsAtZ(const Matrix4d &T, double z,
606 vector<Poly> &polys,
607 double &max_gradient,
608 vector<Poly> &supportpolys,
609 double max_supportangle,
610 double thickness) const
611 {
612 vector<Vector2d> vertices;
613 vector<Triangle> support_triangles;
614 vector<Segment> lines = getCutlines(T, z, vertices, max_gradient,
615 support_triangles, max_supportangle, thickness);
616 //cerr << vertices.size() << " " << lines.size() << endl;
617 if (!CleanupSharedSegments(lines)) return false;
618 //cerr << vertices.size() << " " << lines.size() << endl;
619 if (!CleanupConnectSegments(vertices,lines,true)) return false;
620 //cerr << vertices.size() << " " << lines.size() << endl;
621 vector< vector<uint> > connectedlines; // sequence of connected lines indices
622 if (!getLineSequences(lines, connectedlines)) return false;
623 for (uint i=0; i<connectedlines.size();i++){
624 Poly poly(z);
625 for (uint j = 0; j < connectedlines[i].size();j++){
626 poly.addVertex(vertices[lines[connectedlines[i][j]].start]);
627 }
628 if (lines[connectedlines[i].back()].end !=
629 lines[connectedlines[i].front()].start )
630 poly.addVertex(vertices[lines[connectedlines[i].back()].end]);
631 //cerr << "poly size " << poly.size() << endl;
632 poly.calcHole();
633 polys.push_back(poly);
634 }
635
636 for (uint i = 0; i < support_triangles.size(); i++) {
637 Poly p(z);
638 // keep only part of triangle above z
639 Vector2d lineStart;
640 Vector2d lineEnd;
641 // support_triangles are already transformed
642 int num_cutpoints = support_triangles[i].CutWithPlane(z, Matrix4d::IDENTITY,
643 lineStart, lineEnd);
644 if (num_cutpoints == 0) {
645 for (uint j = 0; j < 3; j++) {
646 p.addVertex(Vector2d(support_triangles[i][j].x(),
647 support_triangles[i][j].y()));
648 }
649 }
650 else if (num_cutpoints > 1) {
651 // add points of triangle above z
652 for (uint j = 0; j < 3; j++) {
653 if (support_triangles[i][j].z() > z) {
654 p.addVertex(Vector2d(support_triangles[i][j].x(),
655 support_triangles[i][j].y()));
656 }
657 }
658 bool reverse = false;
659 // test for order if we get 4 points (2 until now)
660 if (p.size() > 1) {
661 Vector2d i0, i1;
662 const int is = intersect2D_Segments(p[1], lineStart, lineEnd, p[0],
663 i0, i1);
664 if (is > 0 && is < 3) {
665 reverse = true;
666 }
667 }
668 // add cutline points
669 if (reverse) {
670 p.addVertex(lineEnd);
671 p.addVertex(lineStart);
672 } else {
673 p.addVertex(lineStart);
674 p.addVertex(lineEnd);
675 }
676 }
677 if (p.isHole()) p.reverse();
678 supportpolys.push_back(p);
679 }
680
681 // remove polygon areas from triangles
682 // Clipping clipp;
683 // clipp.clear();
684 // clipp.addPolys(supportpolys, subject);
685 // clipp.addPolys(polys, clip);
686 // supportpolys = clipp.subtract(CL::pftPositive,CL::pftPositive);
687
688 return true;
689 }
690
691
find_vertex(const vector<Vector2d> & vertices,const Vector2d & v,double delta=0.0001)692 int find_vertex(const vector<Vector2d> &vertices,
693 const Vector2d &v, double delta = 0.0001)
694 {
695 int found = -1;
696 int count = (int)vertices.size();
697 #ifdef _OPENMP
698 #pragma omp parallel for schedule(dynamic)
699 #endif
700 for (int i=0; i<count; i++) {
701 if (found != -1) continue;
702 if ( (v-vertices[i]).squared_length() < delta ) {
703 found = i;
704 #ifndef _OPENMP
705 break;
706 #endif
707 }
708 }
709 return found;
710 }
711
getCutlines(const Matrix4d & T,double z,vector<Vector2d> & vertices,double & max_gradient,vector<Triangle> & support_triangles,double supportangle,double thickness) const712 vector<Segment> Shape::getCutlines(const Matrix4d &T, double z,
713 vector<Vector2d> &vertices,
714 double &max_gradient,
715 vector<Triangle> &support_triangles,
716 double supportangle,
717 double thickness) const
718 {
719 Vector2d lineStart;
720 Vector2d lineEnd;
721 vector<Segment> lines;
722 // we know our own tranform:
723 Matrix4d transform = T * transform3D.transform ;
724
725 int count = (int)triangles.size();
726 // #ifdef _OPENMP
727 // #pragma omp parallel for schedule(dynamic)
728 // #endif
729 for (int i = 0; i < count; i++)
730 {
731 Segment line(-1,-1);
732 int num_cutpoints = triangles[i].CutWithPlane(z, transform, lineStart, lineEnd);
733 if (num_cutpoints == 0) {
734 if (supportangle >= 0 && thickness > 0) {
735 if (triangles[i].isInZrange(z-thickness, z, transform)) {
736 const double slope = -triangles[i].slopeAngle(transform);
737 if (slope >= supportangle) {
738 support_triangles.push_back(triangles[i].transformed(transform));
739 }
740 }
741 }
742 continue;
743 }
744 if (num_cutpoints > 0) {
745 int havev = find_vertex(vertices, lineStart);
746 if (havev >= 0)
747 line.start = havev;
748 else {
749 line.start = vertices.size();
750 vertices.push_back(lineStart);
751 }
752 if (abs(triangles[i].Normal.z()) > max_gradient)
753 max_gradient = abs(triangles[i].Normal.z());
754 if (supportangle >= 0) {
755 const double slope = -triangles[i].slopeAngle(transform);
756 if (slope >= supportangle)
757 support_triangles.push_back(triangles[i].transformed(transform));
758 }
759 }
760 if (num_cutpoints > 1) {
761 int havev = find_vertex(vertices, lineEnd);
762 if (havev >= 0)
763 line.end = havev;
764 else {
765 line.end = vertices.size();
766 vertices.push_back(lineEnd);
767 }
768 }
769 // Check segment normal against triangle normal. Flip segment, as needed.
770 if (line.start != -1 && line.end != -1 && line.end != line.start)
771 { // if we found a intersecting triangle
772 Vector3d Norm = triangles[i].transformed(transform).Normal;
773 Vector2d triangleNormal = Vector2d(Norm.x(), Norm.y());
774 Vector2d segment = (lineEnd - lineStart);
775 Vector2d segmentNormal(-segment.y(),segment.x());
776 triangleNormal.normalize();
777 segmentNormal.normalize();
778 if( (triangleNormal-segmentNormal).squared_length() > 0.2){
779 // if normals do not align, flip the segment
780 int iswap=line.start;line.start=line.end;line.end=iswap;
781 }
782 // cerr << "line "<<line.start << "-"<<line.end << endl;
783 lines.push_back(line);
784 }
785 }
786 return lines;
787 }
788
789
790 // called from Model::draw
draw(const Settings & settings,bool highlight,uint max_triangles)791 void Shape::draw(const Settings &settings, bool highlight, uint max_triangles)
792 {
793 //cerr << "Shape::draw" << endl;
794 // polygons
795 glEnable(GL_LIGHTING);
796
797 Vector4f no_mat(0.0f, 0.0f, 0.0f, 1.0f);
798 Vector4f low_mat(0.2f, 0.2f, 0.2f, 1.0f);
799 // float mat_ambient[] = {0.7f, 0.7f, 0.7f, 1.0f};
800 // float mat_ambient_color[] = {0.8f, 0.8f, 0.2f, 1.0f};
801 Vector4f mat_diffuse(0.1f, 0.5f, 0.8f, 1.0f);
802 Vector4f mat_specular(1.0f, 1.0f, 1.0f, 1.0f);
803 // float no_shininess = 0.0f;
804 // float low_shininess = 5.0f;
805 // float high_shininess = 100.0f;
806 // float mat_emission[] = {0.3f, 0.2f, 0.2f, 0.0f};
807
808
809 //for (uint i = 0; i < 4; i++) {
810 mat_diffuse = settings.get_colour("Display","PolygonColour");
811 //}
812
813 if (highlight)
814 mat_diffuse.array[3] += 0.3*(1.-mat_diffuse.array[3]);
815
816 // invert colours if partial draw (preview mode)
817 if (max_triangles > 0) {
818 for (uint i = 0; i < 3; i++)
819 mat_diffuse.array[i] = 1.-mat_diffuse.array[i];
820 mat_diffuse[3] = 0.9;
821 }
822
823 mat_specular.array[0] = mat_specular.array[1] = mat_specular.array[2] = settings.get_double("Display","Highlight");;
824
825 /* draw sphere in first row, first column
826 * diffuse reflection only; no ambient or specular
827 */
828 glMaterialfv(GL_FRONT, GL_AMBIENT, low_mat);
829 glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
830 glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
831 glMaterialf(GL_FRONT, GL_SHININESS, 90); // 0..128
832 glMaterialfv(GL_FRONT, GL_EMISSION, no_mat);
833
834 // glEnable (GL_POLYGON_OFFSET_FILL);
835
836 if(settings.get_boolean("Display","DisplayPolygons"))
837 {
838 glEnable(GL_CULL_FACE);
839 glEnable(GL_DEPTH_TEST);
840 glEnable(GL_BLEND);
841 // glDepthMask(GL_TRUE);
842 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); //define blending factors
843 draw_geometry(max_triangles);
844 }
845
846 glDisable (GL_POLYGON_OFFSET_FILL);
847
848 // WireFrame
849 if(settings.get_boolean("Display","DisplayWireframe"))
850 {
851 if(!settings.get_boolean("Display","DisplayWireframeShaded"))
852 glDisable(GL_LIGHTING);
853
854
855 //for (uint i = 0; i < 4; i++)
856 mat_diffuse = settings.get_colour("Display","WireframeColour");
857 glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
858
859 glColor4fv(mat_diffuse);
860 for(size_t i=0;i<triangles.size();i++)
861 {
862 glBegin(GL_LINE_LOOP);
863 glLineWidth(1);
864 glNormal3dv((GLdouble*)&(triangles[i].Normal));
865 glVertex3dv((GLdouble*)&(triangles[i].A));
866 glVertex3dv((GLdouble*)&(triangles[i].B));
867 glVertex3dv((GLdouble*)&(triangles[i].C));
868 glEnd();
869 }
870 }
871
872 glDisable(GL_LIGHTING);
873
874 // normals
875 if(settings.get_boolean("Display","DisplayNormals"))
876 {
877 glColor4fv(settings.get_colour("Display","NormalsColour"));
878 glBegin(GL_LINES);
879 double nlength = settings.get_double("Display","NormalsLength");
880 for(size_t i=0;i<triangles.size();i++)
881 {
882 Vector3d center = (triangles[i].A+triangles[i].B+triangles[i].C)/3.0;
883 glVertex3dv((GLdouble*)¢er);
884 Vector3d N = center + (triangles[i].Normal*nlength);
885 glVertex3dv((GLdouble*)&N);
886 }
887 glEnd();
888 }
889
890 // Endpoints
891 if(settings.get_boolean("Display","DisplayEndpoints"))
892 {
893 glColor4fv(settings.get_colour("Display","EndpointsColour"));
894 glPointSize(settings.get_double("Display","EndPointSize"));
895 glBegin(GL_POINTS);
896 for(size_t i=0;i<triangles.size();i++)
897 {
898 glVertex3dv((GLdouble*)&(triangles[i].A));
899 glVertex3dv((GLdouble*)&(triangles[i].B));
900 glVertex3dv((GLdouble*)&(triangles[i].C));
901 }
902 glEnd();
903 }
904 glDisable(GL_DEPTH_TEST);
905
906 }
907
908 // the bounding box is in real coordinates (not transformed)
drawBBox() const909 void Shape::drawBBox() const
910 {
911 const double minz = max(0.,Min.z()); // draw above zero plane only
912
913 // Draw bbox
914 glColor3f(1,0.2,0.2);
915 glLineWidth(1);
916 glBegin(GL_LINE_LOOP);
917 glVertex3f(Min.x(), Min.y(), minz);
918 glVertex3f(Min.x(), Max.y(), minz);
919 glVertex3f(Max.x(), Max.y(), minz);
920 glVertex3f(Max.x(), Min.y(), minz);
921 glEnd();
922 glBegin(GL_LINE_LOOP);
923 glVertex3f(Min.x(), Min.y(), Max.z());
924 glVertex3f(Min.x(), Max.y(), Max.z());
925 glVertex3f(Max.x(), Max.y(), Max.z());
926 glVertex3f(Max.x(), Min.y(), Max.z());
927 glEnd();
928 glBegin(GL_LINES);
929 glVertex3f(Min.x(), Min.y(), minz);
930 glVertex3f(Min.x(), Min.y(), Max.z());
931 glVertex3f(Min.x(), Max.y(), minz);
932 glVertex3f(Min.x(), Max.y(), Max.z());
933 glVertex3f(Max.x(), Max.y(), minz);
934 glVertex3f(Max.x(), Max.y(), Max.z());
935 glVertex3f(Max.x(), Min.y(), minz);
936 glVertex3f(Max.x(), Min.y(), Max.z());
937 glEnd();
938 /*// show center:
939 glBegin(GL_LINES);
940 glVertex3f(Min.x(), Min.y(), minz);
941 glVertex3f(Max.x(), Max.y(), Max.z());
942 glVertex3f(Max.x(), Min.y(), minz);
943 glVertex3f(Min.x(), Max.y(), Max.z());
944 glEnd();
945 glPointSize(10);
946 glBegin(GL_POINTS);
947 glVertex3dv(Center);
948 glEnd();
949 */
950
951 glColor3f(1,0.6,0.6);
952 ostringstream val;
953 val.precision(1);
954 Vector3d pos;
955 val << fixed << (Max.x()-Min.x());
956 pos = Vector3d((Max.x()+Min.x())/2.,Min.y(),Max.z());
957 Render::draw_string(pos,val.str());
958 val.str("");
959 val << fixed << (Max.y()-Min.y());
960 pos = Vector3d(Min.x(),(Max.y()+Min.y())/2.,Max.z());
961 Render::draw_string(pos,val.str());
962 val.str("");
963 val << fixed << (Max.z()-minz);
964 pos = Vector3d(Min.x(),Min.y(),(Max.z()+minz)/2.);
965 Render::draw_string(pos,val.str());
966
967 }
968
969
draw_geometry(uint max_triangles)970 void Shape::draw_geometry(uint max_triangles)
971 {
972
973 bool listDraw = (max_triangles == 0); // not in preview mode
974 bool haveList = gl_List >= 0;
975
976 if (!listDraw && haveList) {
977 if (gl_List>=0)
978 glDeleteLists(gl_List,1);
979 gl_List = -1;
980 haveList = false;
981 }
982 if (listDraw && !haveList) {
983 gl_List = glGenLists(1);
984 glNewList(gl_List, GL_COMPILE);
985 }
986 if (!listDraw || !haveList) {
987 uint step = 1;
988 if (max_triangles>0) step = floor(triangles.size()/max_triangles);
989 step = max((uint)1,step);
990
991 glBegin(GL_TRIANGLES);
992 for(size_t i=0;i<triangles.size();i+=step)
993 {
994 glNormal3dv(triangles[i].Normal);
995 glVertex3dv(triangles[i].A);
996 glVertex3dv(triangles[i].B);
997 glVertex3dv(triangles[i].C);
998 }
999 glEnd();
1000 }
1001 if (listDraw && !haveList) {
1002 glEndList();
1003 }
1004
1005 if (listDraw && gl_List >= 0) { // have stored list
1006 Glib::TimeVal starttime, endtime;
1007 if (!slow_drawing) {
1008 starttime.assign_current_time();
1009 }
1010 glCallList(gl_List);
1011 if (!slow_drawing) {
1012 endtime.assign_current_time();
1013 Glib::TimeVal usedtime = endtime-starttime;
1014 if (usedtime.as_double() > 0.2) slow_drawing = true;
1015 }
1016 return;
1017 }
1018 }
1019
1020 /*
1021 * sometimes we find adjacent polygons with shared boundary
1022 * points and lines; these cause grief and slowness in
1023 * LinkSegments, so try to identify and join those polygons
1024 * now.
1025 */
1026 // ??? as only coincident lines are removed, couldn't this be
1027 // done easier by just running through all lines and finding them?
CleanupSharedSegments(vector<Segment> & lines)1028 bool CleanupSharedSegments(vector<Segment> &lines)
1029 {
1030 #if 1 // just remove coincident lines
1031 vector<int> lines_to_delete;
1032 int count = (int)lines.size();
1033 #ifdef _OPENMP
1034 //#pragma omp parallel for schedule(dynamic)
1035 #endif
1036 for (int j = 0; j < count; j++) {
1037 const Segment &jr = lines[j];
1038 for (uint k = j + 1; k < lines.size(); k++)
1039 {
1040 const Segment &kr = lines[k];
1041 if ((jr.start == kr.start && jr.end == kr.end) ||
1042 (jr.end == kr.start && jr.start == kr.end))
1043 {
1044 lines_to_delete.push_back (j); // remove both???
1045 //lines_to_delete.push_back (k);
1046 break;
1047 }
1048 }
1049 }
1050 // we need to remove from the end first to avoid disturbing
1051 // the order of removed items
1052 std::sort(lines_to_delete.begin(), lines_to_delete.end());
1053 for (int r = lines_to_delete.size() - 1; r >= 0; r--)
1054 {
1055 lines.erase(lines.begin() + lines_to_delete[r]);
1056 }
1057 return true;
1058
1059 #endif
1060
1061
1062 #if 0
1063 vector<int> vertex_counts; // how many lines have each point
1064 vertex_counts.resize (vertices.size());
1065
1066 for (uint i = 0; i < lines.size(); i++)
1067 {
1068 vertex_counts[lines[i].start]++;
1069 vertex_counts[lines[i].end]++;
1070 }
1071
1072 // ideally all points have an entrance and
1073 // an exit, if we have co-incident lines, then
1074 // we have more than one; do we ?
1075 std::vector<int> duplicate_points;
1076 for (uint i = 0; i < vertex_counts.size(); i++)
1077 if (vertex_counts[i] > 2) // no more than 2 lines should share a point
1078 duplicate_points.push_back (i);
1079
1080 if (duplicate_points.size() == 0)
1081 return true; // all sane
1082
1083 for (uint i = 0; i < duplicate_points.size(); i++)
1084 {
1085 std::vector<int> dup_lines;
1086
1087 // find all line segments with this point in use
1088 for (uint j = 0; j < lines.size(); j++)
1089 {
1090 if (lines[j].start == duplicate_points[i] ||
1091 lines[j].end == duplicate_points[i])
1092 dup_lines.push_back (j);
1093 }
1094
1095 // identify and eliminate identical line segments
1096 // NB. hopefully by here dup_lines.size is small.
1097 std::vector<int> lines_to_delete;
1098 for (uint j = 0; j < dup_lines.size(); j++)
1099 {
1100 const Segment &jr = lines[dup_lines[j]];
1101 for (uint k = j + 1; k < dup_lines.size(); k++)
1102 {
1103 const Segment &kr = lines[dup_lines[k]];
1104 if ((jr.start == kr.start && jr.end == kr.end) ||
1105 (jr.end == kr.start && jr.start == kr.end))
1106 {
1107 lines_to_delete.push_back (dup_lines[j]);
1108 lines_to_delete.push_back (dup_lines[k]);
1109 }
1110 }
1111 }
1112 // we need to remove from the end first to avoid disturbing
1113 // the order of removed items
1114 std::sort(lines_to_delete.begin(), lines_to_delete.end());
1115 for (int r = lines_to_delete.size() - 1; r >= 0; r--)
1116 {
1117 lines.erase(lines.begin() + lines_to_delete[r]);
1118 }
1119 }
1120 return true;
1121 #endif
1122 }
1123
1124 /*
1125 * Unfortunately, finding connections via co-incident points detected by
1126 * the PointHash is not perfect. For reasons unknown (probably rounding
1127 * errors), this is often not enough. We fall-back to finding a nearest
1128 * match from any detached points and joining them, with new synthetic
1129 * segments.
1130 */
CleanupConnectSegments(const vector<Vector2d> & vertices,vector<Segment> & lines,bool connect_all)1131 bool CleanupConnectSegments(const vector<Vector2d> &vertices, vector<Segment> &lines, bool connect_all)
1132 {
1133 vector<int> vertex_types;
1134 vertex_types.resize (vertices.size());
1135 // vector<int> vertex_counts;
1136 // vertex_counts.resize (vertices.size());
1137
1138 // which vertices are referred to, and how much:
1139 int count = lines.size();
1140 for (int i = 0; i < count; i++)
1141 {
1142 vertex_types[lines[i].start]++;
1143 vertex_types[lines[i].end]--;
1144 // vertex_counts[lines[i].start]++;
1145 // vertex_counts[lines[i].end]++;
1146 }
1147
1148 // the vertex_types count should be zero for all connected lines,
1149 // positive for those ending no-where, and negative for
1150 // those starting no-where.
1151 std::vector<int> detached_points; // points with only one line
1152 count = vertex_types.size();
1153 // #ifdef _OPENMP
1154 // #pragma omp parallel for schedule(dynamic)
1155 // #endif
1156 for (int i = 0; i < count; i++)
1157 {
1158 if (vertex_types[i])
1159 {
1160 #if CUTTING_PLANE_DEBUG
1161 cout << "detached point " << i << "\t" << vertex_types[i] << " refs at " << vertices[i].x() << "\t" << vertices[i].y() << "\n";
1162 #endif
1163 detached_points.push_back (i); // i = vertex index
1164 }
1165 }
1166
1167 // Lets hope we have an even number of detached points
1168 if (detached_points.size() % 2) {
1169 cout << "oh dear - an odd number of detached points => an un-pairable impossibility\n";
1170 return false;
1171 }
1172
1173 // pair them nicely to their matching type
1174 count = detached_points.size();
1175 // #ifdef _OPENMP
1176 // #pragma omp parallel for schedule(dynamic)
1177 // #endif
1178 for (int i = 0; i < count; i++)
1179 {
1180 double nearest_dist_sq = (std::numeric_limits<double>::max)();
1181 uint nearest = 0;
1182 int n = detached_points[i]; // vertex index of detached point i
1183 if (n < 0) // handled already
1184 continue;
1185
1186 const Vector2d &p = vertices[n]; // the real detached point
1187 // find nearest other detached point to the detached point n:
1188 for (int j = i + 1; j < count; j++)
1189 {
1190 int pt = detached_points[j];
1191 if (pt < 0)
1192 continue; // already connected
1193
1194 // don't connect a start to a start, or end to end
1195 if (vertex_types[n] == vertex_types[pt])
1196 continue;
1197
1198 const Vector2d &q = vertices[pt]; // the real other point
1199 double dist_sq = (p-q).squared_length(); //pow (p.x() - q.x(), 2) + pow (p.y() - q.y(), 2);
1200 if (dist_sq < nearest_dist_sq)
1201 {
1202 nearest_dist_sq = dist_sq;
1203 nearest = j;
1204 }
1205 }
1206 //assert (nearest != 0);
1207 if (nearest == 0) continue;
1208
1209 // allow points 10mm apart to be joined, not more
1210 if (!connect_all && nearest_dist_sq > 100.0) {
1211 cout << "oh dear - the nearest connecting point is " << sqrt (nearest_dist_sq) << "mm away - not connecting\n";
1212 continue; //return false;
1213 }
1214 // warning above 1mm apart
1215 if (!connect_all && nearest_dist_sq > 1.0) {
1216 cout << "warning - the nearest connecting point is " << sqrt (nearest_dist_sq) << "mm away - connecting anyway\n";
1217 }
1218
1219 #if CUTTING_PLANE_DEBUG
1220 cout << "add join of length " << sqrt (nearest_dist_sq) << "\n" ;
1221 #endif
1222 Segment seg(n, detached_points[nearest]);
1223 if (vertex_types[n] > 0) // already had start but no end at this point
1224 seg.Swap();
1225 lines.push_back(seg);
1226 detached_points[nearest] = -1;
1227 }
1228
1229 return true;
1230 }
1231
1232
info() const1233 string Shape::info() const
1234 {
1235 ostringstream ostr;
1236 ostr <<"Shape with "<<triangles.size() << " triangles "
1237 << "min/max/center: "<<Min<<Max <<Center ;
1238 return ostr.str();
1239 }
1240
1241
1242
1243
1244