Bug Fix #13725 triangulator max initialization
[scilab.git] / scilab / modules / graphic_objects / src / cpp / Triangulator.cpp
1 /*
2  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2011-2012 - DIGITEO - Manuel Juliachs
4  *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13  *
14  */
15
16 #include "Triangulator.hxx"
17
18 #include <math.h>
19
20 void Triangulator::determineSmallestAxis(void)
21 {
22     double minval = 0.;
23
24     Vector3d min = inputPoints[0];
25     Vector3d max = min;
26
27     for (int i = 1 ; i < numPoints; i++)
28     {
29         if (inputPoints[i].x < min.x)
30         {
31             min.x = inputPoints[i].x;
32         }
33         if (inputPoints[i].y < min.y)
34         {
35             min.y = inputPoints[i].y;
36         }
37         if (inputPoints[i].z < min.z)
38         {
39             min.z = inputPoints[i].z;
40         }
41
42         if (inputPoints[i].x > max.x)
43         {
44             max.x = inputPoints[i].x;
45         }
46         if (inputPoints[i].y > max.y)
47         {
48             max.y = inputPoints[i].y;
49         }
50         if (inputPoints[i].z > max.z)
51         {
52             max.z = inputPoints[i].z;
53         }
54     }
55
56     max = minus(max, min);
57
58     if (max.x < max.y)
59     {
60         minval = max.x;
61         smallestAxis = 0;
62         largestAxes[0] = 1;
63         largestAxes[1] = 2;
64     }
65     else
66     {
67         minval = max.y;
68         smallestAxis = 1;
69         largestAxes[0] = 0;
70         largestAxes[1] = 2;
71     }
72
73     if (max.z < minval)
74     {
75         smallestAxis = 2;
76         largestAxes[0] = 0;
77         largestAxes[1] = 1;
78     }
79 }
80
81 void Triangulator::fillPoints(void)
82 {
83     Vector3d point;
84
85     for (int i = 0; i < numPoints; i++)
86     {
87         if (smallestAxis == 0)
88         {
89             /* Smallest axis: X */
90             point.x = inputPoints[i].y;
91             point.y = inputPoints[i].z;
92         }
93         else if (smallestAxis == 1)
94         {
95             /* Smallest axis: Y */
96             point.x = inputPoints[i].x;
97             point.y = inputPoints[i].z;
98         }
99         else
100         {
101             /* Smallest axis: Z (2) */
102             point.x = inputPoints[i].x;
103             point.y = inputPoints[i].y;
104         }
105
106         point.z = 0.0;
107
108         points.push_back(point);
109     }
110 }
111
112 double Triangulator::computeArea(void)
113 {
114     double area = 0.0;
115
116     for (int i = 0; i < numPoints; i++)
117     {
118         int ip1 = (i + 1) % numPoints;
119
120         area += ((points[i].x * points[ip1].y) - (points[i].y * points[ip1].x));
121     }
122
123     area *= 0.5;
124
125     return area;
126 }
127
128 void Triangulator::fillVertexIndices(void)
129 {
130     if (flipped)
131     {
132         for (int i = numPoints - 1; i >= 0; i--)
133         {
134             vertexIndices.push_back(i);
135         }
136     }
137     else
138     {
139         for (int i = 0; i < numPoints; i++)
140         {
141             vertexIndices.push_back(i);
142         }
143     }
144 }
145
146 void Triangulator::removeColinearVertices(void)
147 {
148     double dp = 0.;
149     std::list<int>::iterator vi, vim1, vip1;
150
151     std::vector<Vector3d> sievedPoints;
152     std::list<int> tmpVertexIndices;
153
154     int numColinear = 0;
155     int index = 0;
156
157     for (vi = vertexIndices.begin(); vi != vertexIndices.end(); vi++)
158     {
159         getAdjacentVertices(vi, vim1, vip1);
160
161         dp = computeDotProduct(*vim1, *vi, *vip1);
162
163         if ((!compareVertices(points[*vim1], points[*vi]) && !compareVertices(points[*vi], points[*vip1])) &&
164                 fabs(dp) < TOLERANCE)
165         {
166             numColinear++;
167         }
168         else
169         {
170             sievedPoints.push_back(points[*vi]);
171             actualVertexIndices.push_back(*vi);
172         }
173
174         index++;
175     }
176
177     points.clear();
178
179     if (flipped)
180     {
181         /* Reverse the actual vertex indices list */
182         std::vector<int> tmpList;
183
184         for (size_t i = 0; i < actualVertexIndices.size(); i++)
185         {
186             tmpList.push_back(actualVertexIndices[actualVertexIndices.size() - i - 1]);
187         }
188
189         actualVertexIndices.clear();
190
191         for (size_t i = 0; i < tmpList.size(); i++)
192         {
193             actualVertexIndices.push_back(tmpList[i]);
194         }
195
196         tmpList.clear();
197     }
198
199     if (flipped)
200     {
201         for (size_t i = 0; i < sievedPoints.size(); i++)
202         {
203             points.push_back(sievedPoints[sievedPoints.size() - i - 1]);
204         }
205     }
206     else
207     {
208         for (size_t i = 0; i < sievedPoints.size(); i++)
209         {
210             points.push_back(sievedPoints[i]);
211         }
212     }
213
214     /* Must be updated */
215     numPoints = (int)points.size();
216
217     sievedPoints.clear();
218
219     numColinearVertices = numColinear;
220 }
221
222 void Triangulator::removeDuplicateVertices(void)
223 {
224     int numDuplicateVertices = 0;
225     std::vector<int> duplicateFlagArray;
226
227     std::vector<Vector3d> sievedPoints;
228     std::vector<int> tmpActualVertexIndices;
229
230     sievedPoints.clear();
231     tmpActualVertexIndices.clear();
232
233     duplicateFlagArray.resize(points.size());
234
235     for (size_t i = 0; i < points.size(); i++)
236     {
237         int ic = ((int)i + 1) % (int)points.size();
238         int icm1 = (int)i;
239
240         Vector3d vi = points[icm1];
241         Vector3d vip1 = points[ic];
242
243         if (compareVertices(vi, vip1))
244         {
245             numDuplicateVertices++;
246             duplicateFlagArray[ic] = 1;
247         }
248         else
249         {
250             duplicateFlagArray[ic] = 0;
251         }
252     }
253
254     for (size_t i = 0; i < points.size(); i++)
255     {
256         if (duplicateFlagArray[i] == 0)
257         {
258             /* Keep as it is not a duplicate */
259             sievedPoints.push_back(points[i]);
260             tmpActualVertexIndices.push_back(actualVertexIndices[i]);
261         }
262     }
263
264     actualVertexIndices.clear();
265     points.clear();
266
267     /* Copy the two new lists */
268     for (size_t i = 0; i < tmpActualVertexIndices.size(); i++)
269     {
270         actualVertexIndices.push_back(tmpActualVertexIndices[i]);
271     }
272
273     for (size_t i = 0; i < sievedPoints.size(); i++)
274     {
275         points.push_back(sievedPoints[i]);
276     }
277
278     /* Must be updated */
279     numPoints = (int)points.size();
280
281     duplicateFlagArray.clear();
282     sievedPoints.clear();
283     tmpActualVertexIndices.clear();
284 }
285
286 void Triangulator::fillConvexVerticesList(void)
287 {
288     std::list<int>::iterator vi;
289
290     flagList.resize(vertexIndices.size());
291
292     for (vi = vertexIndices.begin(); vi != vertexIndices.end(); vi++)
293     {
294         if (isConvex(vi))
295         {
296             convexList.push_back(*vi);
297             flagList[*vi] = true;
298         }
299         else
300         {
301             reflexList.push_back(*vi);
302             flagList[*vi] = false;
303         }
304     }
305 }
306
307 void Triangulator::fillEarList(void)
308 {
309     std::list<int>::iterator vi;
310     bool res = false;
311
312     for (vi = vertexIndices.begin(); vi != vertexIndices.end(); vi++)
313     {
314         if (flagList[*vi])
315         {
316             res = isAnEar(vi);
317
318             if (res)
319             {
320                 earList.push_back(*vi);
321             }
322         }
323     }
324 }
325
326 void Triangulator::getAdjacentVertices(std::list<int>::iterator vi, std::list<int>::iterator& vim1, std::list<int>::iterator& vip1)
327 {
328     if (*vi == vertexIndices.front())
329     {
330         vim1 = vertexIndices.end();
331         vim1--;
332     }
333     else
334     {
335         vim1 = vi;
336         vim1--;
337     }
338
339     if (*vi == vertexIndices.back())
340     {
341         vip1 = vertexIndices.begin();
342     }
343     else
344     {
345         vip1 = vi;
346         vip1++;
347     }
348 }
349
350 bool Triangulator::isConvex(std::list<int>::iterator vertex)
351 {
352     double dp = 0.0;
353     std::list<int>::iterator pred, succ;
354
355     getAdjacentVertices(vertex, pred, succ);
356
357     dp = computeDotProduct(*pred, *vertex, *succ);
358
359     if (dp >= 0.0)
360     {
361         return true;
362     }
363     else
364     {
365         return false;
366     }
367 }
368
369 bool Triangulator::isAnEar(std::list<int>::iterator vertex)
370 {
371     bool isEar = true;
372     std::list<int>::iterator pred, succ;
373     std::list<int>::iterator vi;
374
375     Vector3d v0, v1, v2;
376
377     getAdjacentVertices(vertex, pred, succ);
378
379     v0 = points[*pred];
380     v1 = points[*vertex];
381     v2 = points[*succ];
382
383     for (vi = reflexList.begin(); vi != reflexList.end(); vi++)
384     {
385         if (*vi == *pred || *vi == *vertex || *vi == *succ)
386         {
387             continue;
388         }
389         else
390         {
391             bool res;
392
393             res = pointInTriangle(v0, v1, v2, points[*vi]);
394
395             if (res)
396             {
397                 isEar = false;
398                 break;
399             }
400         }
401
402     }
403
404     numEarTests++;
405
406     return isEar;
407 }
408
409 /* To do: streamline */
410 void Triangulator::updateVertex(std::list<int>::iterator vertex)
411 {
412     bool res = false;
413
414     if (flagList[*vertex])
415     {
416         /*
417          * Convex vertex: remove it from the ear list if not an ear any more,
418          * else add it if not already present.
419          */
420         res = isAnEar(vertex);
421
422         if (!res)
423         {
424             earList.remove(*vertex);
425             numDelEars++;
426         }
427         else
428         {
429             std::list<int>::iterator foundEar;
430
431             foundEar = find(earList.begin(), earList.end(), *vertex);
432
433             if (foundEar == earList.end())
434             {
435                 earList.push_front(*vertex);
436                 numAddEars++;
437             }
438         }
439     }
440     else
441     {
442         /*
443          * Non-convex vertex: may become convex, so its flag must be updated, as well
444          * as the reflex vertex list.
445          * Also determine whether it has become an ear and update the ear list accordingly.
446          */
447         if (isConvex(vertex))
448         {
449             flagList[*vertex] = true;
450         }
451
452         if (flagList[*vertex])
453         {
454             res = isAnEar(vertex);
455
456             if (res)
457             {
458                 std::list<int>::iterator foundEar;
459
460                 foundEar = find(earList.begin(), earList.end(), *vertex);
461
462                 if (foundEar == earList.end())
463                 {
464                     earList.push_front(*vertex);
465                     numAddEars++;
466                 }
467             }
468
469             reflexList.remove(*vertex);
470         }
471     }
472 }
473
474 double Triangulator::computeDotProduct(int im1, int i, int ip1)
475 {
476     double dp = 0.;
477     Vector3d eim1p;
478
479     Vector3d eim1 = minus(points[i], points[im1]);
480     Vector3d ei = minus(points[ip1], points[i]);
481
482     /* Normalize */
483     eim1 = normalize(eim1);
484     ei = normalize(ei);
485
486     /* Ought to use cross product */
487     eim1p = perpendicularVector(eim1);
488
489     dp = dot(eim1p, ei);
490
491     return dp;
492 }
493
494 bool Triangulator::pointInTriangle(Vector3d A, Vector3d B, Vector3d C, Vector3d P)
495 {
496     double dot00 = 0., dot01 = 0., dot02 = 0., dot11 = 0., dot12 = 0.;
497     double invDenom = 0.;
498     double u = 0., v = 0.;
499
500     Vector3d v0, v1, v2;
501
502     /* Compute vectors */
503     v0 = minus(C, A);
504     v1 = minus(B, A);
505     v2 = minus(P, A);
506
507     /* Compute dot products */
508     dot00 = dot(v0, v0);
509     dot01 = dot(v0, v1);
510     dot02 = dot(v0, v2);
511     dot11 = dot(v1, v1);
512     dot12 = dot(v1, v2);
513
514     /* Compute barycentric coordinates */
515     invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
516     u = (dot11 * dot02 - dot01 * dot12) * invDenom;
517     v = (dot00 * dot12 - dot01 * dot02) * invDenom;
518
519     /* Check if point is in triangle */
520     return (u >= 0) && (v >= 0) && (u + v < 1);
521 }
522
523 Vector3d Triangulator::minus(Vector3d v0, Vector3d v1)
524 {
525     Vector3d res;
526
527     res.x = v0.x - v1.x;
528     res.y = v0.y - v1.y;
529     res.z = v0.z - v1.z;
530
531     return res;
532 }
533
534 double Triangulator::dot(Vector3d v0, Vector3d v1)
535 {
536     return v0.x * v1.x + v0.y * v1.y + v0.z * v1.z;
537 }
538
539 Vector3d Triangulator::normalize(Vector3d v)
540 {
541     double n = sqrt(v.x * v.x + v.y * v.y);
542
543     if (n < EPSILON)
544     {
545         n = 1.0;
546     }
547
548     v.x /= n;
549     v.y /= n;
550
551     return v;
552 }
553
554 Vector3d Triangulator::perpendicularVector(Vector3d v)
555 {
556     Vector3d perp;
557     perp.x = -v.y;
558     perp.y = v.x;
559     perp.z = 0.0;
560
561     return perp;
562 }
563
564 bool Triangulator::compareVertices(Vector3d v0, Vector3d v1)
565 {
566     if (areEqual(v0.x, v1.x) && areEqual(v0.y, v1.y))
567     {
568         return true;
569     }
570     else
571     {
572         return false;
573     }
574 }
575
576 bool Triangulator::areEqual(double x0, double x1)
577 {
578     double maxAbs = fabs(x0) > fabs(x1) ? fabs(x0) : fabs(x1);
579
580     if (fabs(x0 - x1) <= EPSILON)
581     {
582         return true;
583     }
584     else if (fabs(x0 - x1) <= EPSILON * maxAbs)
585     {
586         return true;
587     }
588     else
589     {
590         return false;
591     }
592 }
593
594 Triangulator::Triangulator(void)
595 {
596     numPoints = 0;
597     numInitPoints = 0;
598     flipped = false;
599     numAddEars = 0;
600     numDelEars = 0;
601     numSteps = 0;
602     numEarTests = 0;
603     numColinearVertices = 0;
604
605     xmin = ymin = zmin = std::numeric_limits<double>::max();
606     xmax = ymax = zmax = std::numeric_limits<double>::lowest();
607
608     inputPoints.clear();
609     points.clear();
610     vertexIndices.clear();
611     actualVertexIndices.clear();
612     earList.clear();
613     convexList.clear();
614     reflexList.clear();
615     flagList.clear();
616     triangleIndices.clear();
617 }
618
619 void Triangulator::initialize(void)
620 {
621     const double xscale = xmax - xmin;
622     const double yscale = ymax - ymin;
623     const double zscale = zmax - zmin;
624     // we scale-translate the point in the cube [0,1]^3 to avoid error with floating point operations
625     for (std::vector<Vector3d>::iterator i = inputPoints.begin(), e = inputPoints.end(); i != e; ++i)
626     {
627         if (xscale)
628         {
629             i->x = (i->x - xmin) / xscale;
630         }
631         else
632         {
633             i->x = 1;
634         }
635         if (yscale)
636         {
637             i->y = (i->y - ymin) / yscale;
638         }
639         else
640         {
641             i->y = 1;
642         }
643         if (zscale)
644         {
645             i->z = (i->z - zmin) / zscale;
646         }
647         else
648         {
649             i->z = 1;
650         }
651     }
652
653     double area = 0.;
654
655     numPoints = (int)inputPoints.size();
656
657     determineSmallestAxis();
658     fillPoints();
659     area = computeArea();
660
661     if (area < 0.0)
662     {
663         flipped = true;
664     }
665     else
666     {
667         flipped = false;
668     }
669
670     fillVertexIndices();
671
672     numInitPoints = numPoints;
673
674     /*
675      * An additional colinear vertices removal pass should be performed
676      * after the duplicate removal pass, as some adjacent edges may become colinear.
677      */
678     removeColinearVertices();
679     removeDuplicateVertices();
680
681     /* Vertex indices must be re-filled */
682     vertexIndices.clear();
683     fillVertexIndices();
684
685     fillConvexVerticesList();
686     fillEarList();
687 }
688
689 void Triangulator::addPoint(double x, double y, double z)
690 {
691     Vector3d point(x, y, z);
692
693     if (x < xmin)
694     {
695         xmin = x;
696     }
697     if (x > xmax)
698     {
699         xmax = x;
700     }
701
702     if (y < ymin)
703     {
704         ymin = y;
705     }
706     if (y > ymax)
707     {
708         ymax = y;
709     }
710
711     if (z < zmin)
712     {
713         zmin = z;
714     }
715     if (z > zmax)
716     {
717         zmax = z;
718     }
719
720     /*point.x = x;
721     point.y = y;
722     point.z = z;*/
723
724     inputPoints.push_back(point);
725 }
726
727 void Triangulator::triangulate(void)
728 {
729     int triIndex = 0;
730     std::list<int>::iterator it;
731     std::list<int>::iterator vertex, pred, succ;
732
733     numSteps = 0;
734
735     while (vertexIndices.size() >= 3 && earList.size() > 0)
736     {
737         int v0 = 0, v1 = 0, v2 = 0;
738         int v0actual = 0, v1actual = 0, v2actual = 0;
739         int vertexIndex = 0;
740
741         it = earList.begin();
742
743         /* If not found, we should break out of the loop. To be checked. */
744         vertex = find(vertexIndices.begin(), vertexIndices.end(), *it);
745         vertexIndex = *vertex;
746
747         getAdjacentVertices(vertex, pred, succ);
748
749         /* Remove */
750         vertexIndices.remove(*vertex);
751         earList.pop_front();
752
753         numDelEars++;
754
755
756         triIndex = *pred;
757         v0 = triIndex;
758         v1 = vertexIndex;
759         triIndex = *succ;
760         v2 = triIndex;
761
762         v0actual = actualVertexIndices[v0];
763         v1actual = actualVertexIndices[v1];
764         v2actual = actualVertexIndices[v2];
765
766         triangleIndices.push_back(v0actual);
767         triangleIndices.push_back(v1actual);
768         triangleIndices.push_back(v2actual);
769
770         /* Update the predecessor vertex */
771         updateVertex(pred);
772
773         /* Update the successor vertex */
774         updateVertex(succ);
775
776         numSteps++;
777     }
778 }
779
780 int Triangulator::getNumberTriangles(void)
781 {
782     return (int)(triangleIndices.size() / 3);
783 }
784
785 int* Triangulator::getIndices(void)
786 {
787     return &triangleIndices[0];
788 }
789
790 int Triangulator::getNumberSteps(void)
791 {
792     return numSteps;
793 }
794
795 int Triangulator::getNumberEarTests(void)
796 {
797     return numEarTests;
798 }
799
800 void Triangulator::clear(void)
801 {
802     inputPoints.clear();
803     points.clear();
804     numPoints = 0;
805     numInitPoints = 0;
806
807     vertexIndices.clear();
808     actualVertexIndices.clear();
809     earList.clear();
810     convexList.clear();
811     reflexList.clear();
812     flagList.clear();
813     triangleIndices.clear();
814
815     numAddEars = 0;
816     numDelEars = 0;
817     numSteps = 0;
818     numEarTests = 0;
819     numColinearVertices = 0;
820 }
821