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 #ifndef TRIANGULATOR_HXX 17 #define TRIANGULATOR_HXX 18 19 #include <algorithm> 20 #include <limits> 21 #include <list> 22 #include <vector> 23 24 #include <iostream> 25 26 /* 27 * A structure representing a point. 28 */ 29 struct Vector3d 30 { 31 double x; 32 double y; 33 double z; 34 Vector3dVector3d35 Vector3d() : x(0.), y(0.), z(0.) { } Vector3dVector3d36 Vector3d(const double _x, const double _y, const double _z) : x(_x), y(_y), z(_z) { } 37 }; 38 39 /** 40 * Triangulator class 41 * An implementation of the ear-clipping triangulation algorithm, 42 * an O(n^2) complexity algorithm, when n is the triangulated polygon's 43 * number of points. The polygon must be simple and non-intersecting. 44 * Triangulation occurs as if the polygon were plane, which is why only 45 * its two largest dimensions are considered and the third is ignored. 46 * 47 * To do: 48 * -extend to take into account self-intersecting polygons 49 * -use a more efficient and robust algorithm (such as the decomposition into monotone pieces O(n log n) algorithm) 50 * -use a dedicated library (more efficient and/or robust) 51 */ 52 class Triangulator 53 { 54 private: 55 /** The array of input points. */ 56 std::vector<Vector3d> inputPoints; 57 58 /** 59 * The array of points, filled from the array of input points depending on the polygon's dimensions. 60 * It contains exactly the same number of points. 61 */ 62 std::vector<Vector3d> points; 63 64 /** The polygon's number of points. */ 65 int numPoints; 66 67 /** The polygons's initial number of points, including colinear vertices. */ 68 int numInitPoints; 69 70 /** 71 * Specifies whether the list of vertex indices must be flipped or not 72 * and indicates the vertices' order. If false, vertices are ordered counter-clockwise 73 * whereas they are ordered clockwise if it is true. 74 */ 75 bool flipped; 76 77 /** The list of vertex indices. */ 78 std::list<int> vertexIndices; 79 80 /** The list of actual vertex indices. */ 81 std::vector<int> actualVertexIndices; 82 83 /** The list of ear vertex indices. */ 84 std::list<int> earList; 85 86 /** The list of convex vertices. */ 87 std::list<int> convexList; 88 89 /** The list of reflex vertices. */ 90 std::list<int> reflexList; 91 92 /** The convexity flag array. */ 93 std::vector<bool> flagList; 94 95 /** The list of output triangle indices. */ 96 std::vector<int> triangleIndices; 97 98 /** The number of insertions into the ear vertex list. */ 99 int numAddEars; 100 101 /** The number of deletions from the ear vertex list. */ 102 int numDelEars; 103 104 /** The number of steps taken by the triangulation's execution. */ 105 int numSteps; 106 107 /** The number of ear tests performed. */ 108 int numEarTests; 109 110 /** The number of colinear vertices. */ 111 int numColinearVertices; 112 113 double xmin, xmax, ymin, ymax, zmin, zmax; 114 115 private: 116 117 /** 118 * Fills the array of points with a projection of the points in the array of input points. 119 **/ 120 void fillPoints(void); 121 122 /** 123 * Computes and returns the polygon's signed area. 124 * Its sign depends on the order of the polygon's vertices. 125 * If positive, the vertices are ordered counter-clockwise whereas they are ordered 126 * clockwise if it is negative. 127 * @return the polygon's area. 128 */ 129 double computeArea(void); 130 131 /** 132 * Fills the list of vertex indices, depending on their order. 133 */ 134 void fillVertexIndices(void); 135 136 /** 137 * Removes colinear vertices. 138 */ 139 void removeColinearVertices(void); 140 141 /** 142 * Removes duplicate vertices. 143 */ 144 void removeDuplicateVertices(void); 145 146 /** 147 * Fills the list of convex vertices, determining whether each vertex vi is convex or not. 148 * It also fills the convexity flag array. 149 */ 150 void fillConvexVerticesList(void); 151 152 /** 153 * Fills the list of ear vertices. 154 */ 155 void fillEarList(void); 156 157 /** 158 * Gets the vertices adjacent to vertex i. 159 * @param[in] vertex i's iterator. 160 * @param[out] a reference to vertex i-1's iterator. 161 * @param[out] a reference to vertex i+1's iterator. 162 */ 163 void getAdjacentVertices(std::list<int>::iterator vi, std::list<int>::iterator& vim1, std::list<int>::iterator& vip1); 164 165 /** 166 * Determines whether a vertex is convex or not. 167 * @param[in] the vertex's iterator. 168 * @return true if the vertex is convex, false if it is not. 169 */ 170 bool isConvex(std::list<int>::iterator vertex); 171 172 /** 173 * Determines whether a vertex is an ear vertex. 174 * @param[in] the vertex's iterator. 175 * @return true if it is an ear, false if not. 176 */ 177 bool isAnEar(std::list<int>::iterator vertex); 178 179 /** 180 * Updates a vertex's state due to the next/previous vertex 181 * having been removed. It performs an ear test and accordingly 182 * updates the ear vertex list. If the vertex becomes convex, its 183 * flag is updated as well as the reflex vertex list. 184 * @param[in] the updated vertex's iterator. 185 */ 186 void updateVertex(std::list<int>::iterator vertex); 187 188 /** 189 * Computes the dot product between edge e(i) (from vertex i to i+1) and the vector 190 * orthogonal to edge e(i-1) (from vertex i-1 to i) such that the angle from e(i-1) 191 * to this vector vector is +90°. 192 * @param[in] the index of vertex i-1. 193 * @param[in] the index of vertex i. 194 * @param[in] the index of vertex i+1. 195 * @return the dot product. 196 */ 197 double computeDotProduct(int im1, int i, int ip1); 198 199 /** 200 * Determines whether a point P is located within the triangle ABC. 201 * It considers that P and ABC both lie in the xy plane. 202 * @param[in] the triangle's first point. 203 * @param[in] the triangle's second point. 204 * @param[in] the triangle's third point. 205 * @return true if P is located within ABC, false if not. 206 */ 207 bool pointInTriangle(Vector3d A, Vector3d B, Vector3d C, Vector3d P); 208 209 /** 210 * Subtracts the second input vector from the first one and returns the result. 211 * It should be moved to a Vector3d class. 212 * @param[in] the first vector. 213 * @param[in] the second vector. 214 * @return the resulting vector. 215 */ 216 static Vector3d minus(Vector3d v0, Vector3d v1); 217 218 /** 219 * Add the first and second vectors 220 * It should be moved to a Vector3d class. 221 * @param[in] the first vector. 222 * @param[in] the second vector. 223 * @return the resulting vector. 224 */ 225 static Vector3d plus(Vector3d v0, Vector3d v1); 226 227 /** 228 * Computes and returns the dot product of two vectors. 229 * It should be moved to a Vector3d class. 230 * @param[in] the first vector. 231 * @param[in] the second vector. 232 * @return the dot product. 233 */ 234 static double dot(Vector3d v0, Vector3d v1); 235 236 /** 237 * Normalizes a 2D vector, the z coordinate is ignored. 238 * It should be moved to a Vector3d class. 239 * @param[in] the vector to normalize. 240 * @return the normalized vector. 241 */ 242 static Vector3d normalize(Vector3d v); 243 244 /** 245 * Compares whether two vertices are identical or not. 246 * Two vertices are considered identical if they x and y coordinates 247 * are equal, the z coordinate being ignored. 248 * It should be moved to a Vector3d class. 249 * @param[in] the first vector. 250 * @param[in] the second vector. 251 * @return true if the two vertices are identical, false if not. 252 */ 253 static bool compareVertices(Vector3d v0, Vector3d v1); 254 255 /** 256 * Determines whether two floating-point values are equal. 257 * @param[in] the first value. 258 * @param[in] the second value. 259 * @return true if the values are equal, false if not. 260 */ 261 static bool areEqual(double x0, double x1); 262 263 /** 264 * Computes and returns a vector p orthogonal to vector v, 265 * such that the angle from v to p is +90°. v and p are considered 266 * to lie in the xy plane. 267 * It should be moved to a Vector3d class. 268 * @param[in] the vector. 269 * @return the vector perpendicular to vector v. 270 */ 271 static Vector3d perpendicularVector(Vector3d v); 272 273 /** 274 * Computes and returns the cross product of two vectors. 275 * It should be moved to a Vector3d class. 276 * @param[in] the vector a. 277 * @param[in] the vector b. 278 * @return the cross product vector perpendicular to vectors a and b. 279 */ 280 static Vector3d cross(Vector3d a, Vector3d b); 281 282 /** 283 * Computes the multimplication of 3x3 matrix Out = A * B; 284 * @param[in] the left matrix a. 285 * @param[in] the right matrix b. 286 * @param[out] the multiplied matrix. 287 **/ 288 void matrixMatrixMul(double (&a)[3][3], double (&b)[3][3], double (&out)[3][3]); 289 290 /** 291 * Computes the multimplication of 3x3 matrix and a vector Pout = M * Pin; 292 * @param[in] the matrix m. 293 * @param[in] the vector pin. 294 * @param[out] the multiplied vector pout. 295 **/ 296 void matrixVectorMul(double (&m)[3][3], Vector3d & pin, Vector3d & pout); 297 298 /** 299 * Computes the multimplication of 3x3 matrix and a vector PinOut = M * PinOut; 300 * @param[in] the matrix m. 301 * @param[in] the vector pinOut. 302 * @param[out] the multiplied vector pinOut. 303 **/ 304 void matrixVectorMul(double (&m)[3][3], Vector3d & pinOut); 305 306 /** 307 * Transform all points from the array of points, to the cube [0,1]^3 308 **/ 309 void normalizePoints(void); 310 311 /** 312 * 3x3 Matriz diagonalizer using jacobi method 313 * Original idea from: http://www.melax.com/diag.html 314 * @param[in] the matrix A to be diagonalized. 315 * @param[out] the matrix Q which diagonalizer A, D = Q' * D * Q. 316 * @param[out] the matrix D diagonal. 317 **/ 318 void diagonalize(const double (&A)[3][3], double (&Q)[3][3], double (&D)[3][3]); 319 320 public: 321 /** Default constructor. */ 322 Triangulator(void); 323 324 /** 325 * Initializes and fills the lists used by the algorithm from the 326 * list of input points, which must have been filled beforehand. 327 */ 328 void initialize(void); 329 330 /** 331 * Adds a point to the array of input points. 332 * @param[in] the point's x-coordinate. 333 * @param[in] the point's y-coordinate. 334 * @param[in] the point's z-coordinate. 335 */ 336 void addPoint(double x, double y, double z); 337 338 /** 339 * Triangulates the input data. 340 */ 341 void triangulate(void); 342 343 /** 344 * Returns the number of triangles resulting from the algorithm's execution. 345 * @return the number of triangles. 346 */ 347 int getNumberTriangles(void); 348 349 /** 350 * Returns the array of triangle indices resulting from the algorithm's execution. 351 * @return a pointer to the array of triangle indices. 352 */ 353 int* getIndices(void); 354 355 /** 356 * Returns the number steps performed by the algorithm. 357 * @return the number of steps performed. 358 */ 359 int getNumberSteps(); 360 361 /** 362 * Returns the number of ear tests performed by the algorithm. 363 * Relevant only to debugging. 364 * @return the number of ear tests performed. 365 */ 366 int getNumberEarTests(); 367 368 /** 369 * Clears the lists of points, vertex indices, and the lists storing 370 * the algorithm's internal state (vertex indices, ear vertices and others). 371 */ 372 void clear(void); 373 }; 374 375 /** 376 * An arbitrary tolerance value used to determine colinear edges. 377 */ 378 #define TOLERANCE 0.0000001 379 380 /** 381 * An epsilon value used when comparing vertices. 382 */ 383 #define EPSILON 0.00000001 384 385 #endif 386