Fix T72340: Version bump for recent Userdef changes
[blender.git] / source / blender / bmesh / intern / bmesh_polygon.c
1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  */
16
17 /** \file
18  * \ingroup bmesh
19  *
20  * This file contains code for dealing
21  * with polygons (normal/area calculation,
22  * tessellation, etc)
23  */
24
25 #include "DNA_listBase.h"
26 #include "DNA_modifier_types.h"
27 #include "DNA_meshdata_types.h"
28
29 #include "MEM_guardedalloc.h"
30
31 #include "BLI_alloca.h"
32 #include "BLI_math.h"
33 #include "BLI_memarena.h"
34 #include "BLI_polyfill_2d.h"
35 #include "BLI_polyfill_2d_beautify.h"
36 #include "BLI_linklist.h"
37 #include "BLI_heap.h"
38
39 #include "bmesh.h"
40 #include "bmesh_tools.h"
41
42 #include "BKE_customdata.h"
43
44 #include "intern/bmesh_private.h"
45
46 /**
47  * \brief COMPUTE POLY NORMAL (BMFace)
48  *
49  * Same as #normal_poly_v3 but operates directly on a bmesh face.
50  */
51 static float bm_face_calc_poly_normal(const BMFace *f, float n[3])
52 {
53   BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
54   BMLoop *l_iter = l_first;
55   const float *v_prev = l_first->prev->v->co;
56   const float *v_curr = l_first->v->co;
57
58   zero_v3(n);
59
60   /* Newell's Method */
61   do {
62     add_newell_cross_v3_v3v3(n, v_prev, v_curr);
63
64     l_iter = l_iter->next;
65     v_prev = v_curr;
66     v_curr = l_iter->v->co;
67
68   } while (l_iter != l_first);
69
70   return normalize_v3(n);
71 }
72
73 /**
74  * \brief COMPUTE POLY NORMAL (BMFace)
75  *
76  * Same as #bm_face_calc_poly_normal
77  * but takes an array of vertex locations.
78  */
79 static float bm_face_calc_poly_normal_vertex_cos(const BMFace *f,
80                                                  float r_no[3],
81                                                  float const (*vertexCos)[3])
82 {
83   BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
84   BMLoop *l_iter = l_first;
85   const float *v_prev = vertexCos[BM_elem_index_get(l_first->prev->v)];
86   const float *v_curr = vertexCos[BM_elem_index_get(l_first->v)];
87
88   zero_v3(r_no);
89
90   /* Newell's Method */
91   do {
92     add_newell_cross_v3_v3v3(r_no, v_prev, v_curr);
93
94     l_iter = l_iter->next;
95     v_prev = v_curr;
96     v_curr = vertexCos[BM_elem_index_get(l_iter->v)];
97   } while (l_iter != l_first);
98
99   return normalize_v3(r_no);
100 }
101
102 /**
103  * \brief COMPUTE POLY CENTER (BMFace)
104  */
105 static void bm_face_calc_poly_center_median_vertex_cos(const BMFace *f,
106                                                        float r_cent[3],
107                                                        float const (*vertexCos)[3])
108 {
109   const BMLoop *l_first, *l_iter;
110
111   zero_v3(r_cent);
112
113   /* Newell's Method */
114   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
115   do {
116     add_v3_v3(r_cent, vertexCos[BM_elem_index_get(l_iter->v)]);
117   } while ((l_iter = l_iter->next) != l_first);
118   mul_v3_fl(r_cent, 1.0f / f->len);
119 }
120
121 /**
122  * For tools that insist on using triangles, ideally we would cache this data.
123  *
124  * \param use_fixed_quad: When true,
125  * always split quad along (0 -> 2) regardless of concave corners,
126  * (as done in #BM_mesh_calc_tessellation).
127  * \param r_loops: Store face loop pointers, (f->len)
128  * \param r_index: Store triangle triples, indices into \a r_loops,  `((f->len - 2) * 3)`
129  */
130 void BM_face_calc_tessellation(const BMFace *f,
131                                const bool use_fixed_quad,
132                                BMLoop **r_loops,
133                                uint (*r_index)[3])
134 {
135   BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
136   BMLoop *l_iter;
137
138   if (f->len == 3) {
139     *r_loops++ = (l_iter = l_first);
140     *r_loops++ = (l_iter = l_iter->next);
141     *r_loops++ = (l_iter->next);
142
143     r_index[0][0] = 0;
144     r_index[0][1] = 1;
145     r_index[0][2] = 2;
146   }
147   else if (f->len == 4 && use_fixed_quad) {
148     *r_loops++ = (l_iter = l_first);
149     *r_loops++ = (l_iter = l_iter->next);
150     *r_loops++ = (l_iter = l_iter->next);
151     *r_loops++ = (l_iter->next);
152
153     r_index[0][0] = 0;
154     r_index[0][1] = 1;
155     r_index[0][2] = 2;
156
157     r_index[1][0] = 0;
158     r_index[1][1] = 2;
159     r_index[1][2] = 3;
160   }
161   else {
162     float axis_mat[3][3];
163     float(*projverts)[2] = BLI_array_alloca(projverts, f->len);
164     int j;
165
166     axis_dominant_v3_to_m3(axis_mat, f->no);
167
168     j = 0;
169     l_iter = l_first;
170     do {
171       mul_v2_m3v3(projverts[j], axis_mat, l_iter->v->co);
172       r_loops[j] = l_iter;
173       j++;
174     } while ((l_iter = l_iter->next) != l_first);
175
176     /* complete the loop */
177     BLI_polyfill_calc(projverts, f->len, -1, r_index);
178   }
179 }
180
181 /**
182  * Return a point inside the face.
183  */
184 void BM_face_calc_point_in_face(const BMFace *f, float r_co[3])
185 {
186   const BMLoop *l_tri[3];
187
188   if (f->len == 3) {
189     const BMLoop *l = BM_FACE_FIRST_LOOP(f);
190     ARRAY_SET_ITEMS(l_tri, l, l->next, l->prev);
191   }
192   else {
193     /* tessellation here seems overkill when in many cases this will be the center,
194      * but without this we can't be sure the point is inside a concave face. */
195     const int tottri = f->len - 2;
196     BMLoop **loops = BLI_array_alloca(loops, f->len);
197     uint(*index)[3] = BLI_array_alloca(index, tottri);
198     int j;
199     int j_best = 0; /* use as fallback when unset */
200     float area_best = -1.0f;
201
202     BM_face_calc_tessellation(f, false, loops, index);
203
204     for (j = 0; j < tottri; j++) {
205       const float *p1 = loops[index[j][0]]->v->co;
206       const float *p2 = loops[index[j][1]]->v->co;
207       const float *p3 = loops[index[j][2]]->v->co;
208       const float area = area_squared_tri_v3(p1, p2, p3);
209       if (area > area_best) {
210         j_best = j;
211         area_best = area;
212       }
213     }
214
215     ARRAY_SET_ITEMS(
216         l_tri, loops[index[j_best][0]], loops[index[j_best][1]], loops[index[j_best][2]]);
217   }
218
219   mid_v3_v3v3v3(r_co, l_tri[0]->v->co, l_tri[1]->v->co, l_tri[2]->v->co);
220 }
221
222 /**
223  * get the area of the face
224  */
225 float BM_face_calc_area(const BMFace *f)
226 {
227   /* inline 'area_poly_v3' logic, avoid creating a temp array */
228   const BMLoop *l_iter, *l_first;
229   float n[3];
230
231   zero_v3(n);
232   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
233   do {
234     add_newell_cross_v3_v3v3(n, l_iter->v->co, l_iter->next->v->co);
235   } while ((l_iter = l_iter->next) != l_first);
236   return len_v3(n) * 0.5f;
237 }
238
239 /**
240  * Get the area of the face in world space.
241  */
242 float BM_face_calc_area_with_mat3(const BMFace *f, const float mat3[3][3])
243 {
244   /* inline 'area_poly_v3' logic, avoid creating a temp array */
245   const BMLoop *l_iter, *l_first;
246   float co[3];
247   float n[3];
248
249   zero_v3(n);
250   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
251   mul_v3_m3v3(co, mat3, l_iter->v->co);
252   do {
253     float co_next[3];
254     mul_v3_m3v3(co_next, mat3, l_iter->next->v->co);
255     add_newell_cross_v3_v3v3(n, co, co_next);
256     copy_v3_v3(co, co_next);
257   } while ((l_iter = l_iter->next) != l_first);
258   return len_v3(n) * 0.5f;
259 }
260
261 /**
262  * get the area of UV face
263  */
264 float BM_face_calc_area_uv(const BMFace *f, int cd_loop_uv_offset)
265 {
266   /* inline 'area_poly_v2' logic, avoid creating a temp array */
267   const BMLoop *l_iter, *l_first;
268
269   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
270   /* The Trapezium Area Rule */
271   float cross = 0.0f;
272   do {
273     const MLoopUV *luv = BM_ELEM_CD_GET_VOID_P(l_iter, cd_loop_uv_offset);
274     const MLoopUV *luv_next = BM_ELEM_CD_GET_VOID_P(l_iter->next, cd_loop_uv_offset);
275     cross += (luv_next->uv[0] - luv->uv[0]) * (luv_next->uv[1] + luv->uv[1]);
276   } while ((l_iter = l_iter->next) != l_first);
277   return fabsf(cross * 0.5f);
278 }
279
280 /**
281  * compute the perimeter of an ngon
282  */
283 float BM_face_calc_perimeter(const BMFace *f)
284 {
285   const BMLoop *l_iter, *l_first;
286   float perimeter = 0.0f;
287
288   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
289   do {
290     perimeter += len_v3v3(l_iter->v->co, l_iter->next->v->co);
291   } while ((l_iter = l_iter->next) != l_first);
292
293   return perimeter;
294 }
295
296 /**
297  * Calculate the perimeter of a ngon in world space.
298  */
299 float BM_face_calc_perimeter_with_mat3(const BMFace *f, const float mat3[3][3])
300 {
301   const BMLoop *l_iter, *l_first;
302   float co[3];
303   float perimeter = 0.0f;
304
305   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
306   mul_v3_m3v3(co, mat3, l_iter->v->co);
307   do {
308     float co_next[3];
309     mul_v3_m3v3(co_next, mat3, l_iter->next->v->co);
310     perimeter += len_v3v3(co, co_next);
311     copy_v3_v3(co, co_next);
312   } while ((l_iter = l_iter->next) != l_first);
313
314   return perimeter;
315 }
316
317 /**
318  * Utility function to calculate the edge which is most different from the other two.
319  *
320  * \return The first edge index, where the second vertex is ``(index + 1) % 3``.
321  */
322 static int bm_vert_tri_find_unique_edge(BMVert *verts[3])
323 {
324   /* find the most 'unique' loop, (greatest difference to others) */
325 #if 1
326   /* optimized version that avoids sqrt */
327   float difs[3];
328   for (int i_prev = 1, i_curr = 2, i_next = 0; i_next < 3; i_prev = i_curr, i_curr = i_next++) {
329     const float *co = verts[i_curr]->co;
330     const float *co_other[2] = {verts[i_prev]->co, verts[i_next]->co};
331     float proj_dir[3];
332     mid_v3_v3v3(proj_dir, co_other[0], co_other[1]);
333     sub_v3_v3(proj_dir, co);
334
335     float proj_pair[2][3];
336     project_v3_v3v3(proj_pair[0], co_other[0], proj_dir);
337     project_v3_v3v3(proj_pair[1], co_other[1], proj_dir);
338     difs[i_next] = len_squared_v3v3(proj_pair[0], proj_pair[1]);
339   }
340 #else
341   const float lens[3] = {
342       len_v3v3(verts[0]->co, verts[1]->co),
343       len_v3v3(verts[1]->co, verts[2]->co),
344       len_v3v3(verts[2]->co, verts[0]->co),
345   };
346   const float difs[3] = {
347       fabsf(lens[1] - lens[2]),
348       fabsf(lens[2] - lens[0]),
349       fabsf(lens[0] - lens[1]),
350   };
351 #endif
352
353   int order[3] = {0, 1, 2};
354   axis_sort_v3(difs, order);
355
356   return order[0];
357 }
358
359 /**
360  * Calculate a tangent from any 3 vertices.
361  *
362  * The tangent aligns to the most *unique* edge
363  * (the edge most unlike the other two).
364  *
365  * \param r_tangent: Calculated unit length tangent (return value).
366  */
367 void BM_vert_tri_calc_tangent_edge(BMVert *verts[3], float r_tangent[3])
368 {
369   const int index = bm_vert_tri_find_unique_edge(verts);
370
371   sub_v3_v3v3(r_tangent, verts[index]->co, verts[(index + 1) % 3]->co);
372
373   normalize_v3(r_tangent);
374 }
375
376 /**
377  * Calculate a tangent from any 3 vertices,
378  *
379  * The tangent follows the center-line formed by the most unique edges center
380  * and the opposite vertex.
381  *
382  * \param r_tangent: Calculated unit length tangent (return value).
383  */
384 void BM_vert_tri_calc_tangent_edge_pair(BMVert *verts[3], float r_tangent[3])
385 {
386   const int index = bm_vert_tri_find_unique_edge(verts);
387
388   const float *v_a = verts[index]->co;
389   const float *v_b = verts[(index + 1) % 3]->co;
390   const float *v_other = verts[(index + 2) % 3]->co;
391
392   mid_v3_v3v3(r_tangent, v_a, v_b);
393   sub_v3_v3v3(r_tangent, v_other, r_tangent);
394
395   normalize_v3(r_tangent);
396 }
397
398 /**
399  * Compute the tangent of the face, using the longest edge.
400  */
401 void BM_face_calc_tangent_edge(const BMFace *f, float r_tangent[3])
402 {
403   const BMLoop *l_long = BM_face_find_longest_loop((BMFace *)f);
404
405   sub_v3_v3v3(r_tangent, l_long->v->co, l_long->next->v->co);
406
407   normalize_v3(r_tangent);
408 }
409
410 /**
411  * Compute the tangent of the face, using the two longest disconnected edges.
412  *
413  * \param r_tangent: Calculated unit length tangent (return value).
414  */
415 void BM_face_calc_tangent_edge_pair(const BMFace *f, float r_tangent[3])
416 {
417   if (f->len == 3) {
418     BMVert *verts[3];
419
420     BM_face_as_array_vert_tri((BMFace *)f, verts);
421
422     BM_vert_tri_calc_tangent_edge_pair(verts, r_tangent);
423   }
424   else if (f->len == 4) {
425     /* Use longest edge pair */
426     BMVert *verts[4];
427     float vec[3], vec_a[3], vec_b[3];
428
429     BM_face_as_array_vert_quad((BMFace *)f, verts);
430
431     sub_v3_v3v3(vec_a, verts[3]->co, verts[2]->co);
432     sub_v3_v3v3(vec_b, verts[0]->co, verts[1]->co);
433     add_v3_v3v3(r_tangent, vec_a, vec_b);
434
435     sub_v3_v3v3(vec_a, verts[0]->co, verts[3]->co);
436     sub_v3_v3v3(vec_b, verts[1]->co, verts[2]->co);
437     add_v3_v3v3(vec, vec_a, vec_b);
438     /* use the longest edge length */
439     if (len_squared_v3(r_tangent) < len_squared_v3(vec)) {
440       copy_v3_v3(r_tangent, vec);
441     }
442   }
443   else {
444     /* For ngons use two longest disconnected edges */
445     BMLoop *l_long = BM_face_find_longest_loop((BMFace *)f);
446     BMLoop *l_long_other = NULL;
447
448     float len_max_sq = 0.0f;
449     float vec_a[3], vec_b[3];
450
451     BMLoop *l_iter = l_long->prev->prev;
452     BMLoop *l_last = l_long->next;
453
454     do {
455       const float len_sq = len_squared_v3v3(l_iter->v->co, l_iter->next->v->co);
456       if (len_sq >= len_max_sq) {
457         l_long_other = l_iter;
458         len_max_sq = len_sq;
459       }
460     } while ((l_iter = l_iter->prev) != l_last);
461
462     sub_v3_v3v3(vec_a, l_long->next->v->co, l_long->v->co);
463     sub_v3_v3v3(vec_b, l_long_other->v->co, l_long_other->next->v->co);
464     add_v3_v3v3(r_tangent, vec_a, vec_b);
465
466     /* Edges may not be opposite side of the ngon,
467      * this could cause problems for ngons with multiple-aligned edges of the same length.
468      * Fallback to longest edge. */
469     if (UNLIKELY(normalize_v3(r_tangent) == 0.0f)) {
470       normalize_v3_v3(r_tangent, vec_a);
471     }
472   }
473 }
474
475 /**
476  * Compute the tangent of the face, using the edge farthest away from any vertex in the face.
477  *
478  * \param r_tangent: Calculated unit length tangent (return value).
479  */
480 void BM_face_calc_tangent_edge_diagonal(const BMFace *f, float r_tangent[3])
481 {
482   BMLoop *l_iter, *l_first;
483
484   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
485
486   /* In case of degenerate faces. */
487   zero_v3(r_tangent);
488
489   /* warning: O(n^2) loop here, take care! */
490   float dist_max_sq = 0.0f;
491   do {
492     BMLoop *l_iter_other = l_iter->next;
493     BMLoop *l_iter_last = l_iter->prev;
494     do {
495       BLI_assert(!ELEM(l_iter->v->co, l_iter_other->v->co, l_iter_other->next->v->co));
496       float co_other[3], vec[3];
497       closest_to_line_segment_v3(
498           co_other, l_iter->v->co, l_iter_other->v->co, l_iter_other->next->v->co);
499       sub_v3_v3v3(vec, l_iter->v->co, co_other);
500
501       const float dist_sq = len_squared_v3(vec);
502       if (dist_sq > dist_max_sq) {
503         dist_max_sq = dist_sq;
504         copy_v3_v3(r_tangent, vec);
505       }
506     } while ((l_iter_other = l_iter_other->next) != l_iter_last);
507   } while ((l_iter = l_iter->next) != l_first);
508
509   normalize_v3(r_tangent);
510 }
511
512 /**
513  * Compute the tangent of the face, using longest distance between vertices on the face.
514  *
515  * \note The logic is almost identical to #BM_face_calc_tangent_edge_diagonal
516  */
517 void BM_face_calc_tangent_vert_diagonal(const BMFace *f, float r_tangent[3])
518 {
519   BMLoop *l_iter, *l_first;
520
521   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
522
523   /* In case of degenerate faces. */
524   zero_v3(r_tangent);
525
526   /* warning: O(n^2) loop here, take care! */
527   float dist_max_sq = 0.0f;
528   do {
529     BMLoop *l_iter_other = l_iter->next;
530     do {
531       float vec[3];
532       sub_v3_v3v3(vec, l_iter->v->co, l_iter_other->v->co);
533
534       const float dist_sq = len_squared_v3(vec);
535       if (dist_sq > dist_max_sq) {
536         dist_max_sq = dist_sq;
537         copy_v3_v3(r_tangent, vec);
538       }
539     } while ((l_iter_other = l_iter_other->next) != l_iter);
540   } while ((l_iter = l_iter->next) != l_first);
541
542   normalize_v3(r_tangent);
543 }
544
545 /**
546  * Compute a meaningful direction along the face (use for gizmo axis).
547  *
548  * \note Callers shouldn't depend on the *exact* method used here.
549  */
550 void BM_face_calc_tangent_auto(const BMFace *f, float r_tangent[3])
551 {
552   if (f->len == 3) {
553     /* most 'unique' edge of a triangle */
554     BMVert *verts[3];
555     BM_face_as_array_vert_tri((BMFace *)f, verts);
556     BM_vert_tri_calc_tangent_edge(verts, r_tangent);
557   }
558   else if (f->len == 4) {
559     /* longest edge pair of a quad */
560     BM_face_calc_tangent_edge_pair((BMFace *)f, r_tangent);
561   }
562   else {
563     /* longest edge of an ngon */
564     BM_face_calc_tangent_edge((BMFace *)f, r_tangent);
565   }
566 }
567
568 /**
569  * expands bounds (min/max must be initialized).
570  */
571 void BM_face_calc_bounds_expand(const BMFace *f, float min[3], float max[3])
572 {
573   const BMLoop *l_iter, *l_first;
574   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
575   do {
576     minmax_v3v3_v3(min, max, l_iter->v->co);
577   } while ((l_iter = l_iter->next) != l_first);
578 }
579
580 /**
581  * computes center of face in 3d.  uses center of bounding box.
582  */
583 void BM_face_calc_center_bounds(const BMFace *f, float r_cent[3])
584 {
585   const BMLoop *l_iter, *l_first;
586   float min[3], max[3];
587
588   INIT_MINMAX(min, max);
589
590   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
591   do {
592     minmax_v3v3_v3(min, max, l_iter->v->co);
593   } while ((l_iter = l_iter->next) != l_first);
594
595   mid_v3_v3v3(r_cent, min, max);
596 }
597
598 /**
599  * computes the center of a face, using the mean average
600  */
601 void BM_face_calc_center_median(const BMFace *f, float r_cent[3])
602 {
603   const BMLoop *l_iter, *l_first;
604
605   zero_v3(r_cent);
606
607   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
608   do {
609     add_v3_v3(r_cent, l_iter->v->co);
610   } while ((l_iter = l_iter->next) != l_first);
611   mul_v3_fl(r_cent, 1.0f / (float)f->len);
612 }
613
614 /**
615  * computes the center of a face, using the mean average
616  * weighted by edge length
617  */
618 void BM_face_calc_center_median_weighted(const BMFace *f, float r_cent[3])
619 {
620   const BMLoop *l_iter;
621   const BMLoop *l_first;
622   float totw = 0.0f;
623   float w_prev;
624
625   zero_v3(r_cent);
626
627   l_iter = l_first = BM_FACE_FIRST_LOOP(f);
628   w_prev = BM_edge_calc_length(l_iter->prev->e);
629   do {
630     const float w_curr = BM_edge_calc_length(l_iter->e);
631     const float w = (w_curr + w_prev);
632     madd_v3_v3fl(r_cent, l_iter->v->co, w);
633     totw += w;
634     w_prev = w_curr;
635   } while ((l_iter = l_iter->next) != l_first);
636
637   if (totw != 0.0f) {
638     mul_v3_fl(r_cent, 1.0f / (float)totw);
639   }
640 }
641
642 /**
643  * \brief POLY ROTATE PLANE
644  *
645  * Rotates a polygon so that it's
646  * normal is pointing towards the mesh Z axis
647  */
648 void poly_rotate_plane(const float normal[3], float (*verts)[3], const uint nverts)
649 {
650   float mat[3][3];
651   float co[3];
652   uint i;
653
654   co[2] = 0.0f;
655
656   axis_dominant_v3_to_m3(mat, normal);
657   for (i = 0; i < nverts; i++) {
658     mul_v2_m3v3(co, mat, verts[i]);
659     copy_v3_v3(verts[i], co);
660   }
661 }
662
663 /**
664  * updates face and vertex normals incident on an edge
665  */
666 void BM_edge_normals_update(BMEdge *e)
667 {
668   BMIter iter;
669   BMFace *f;
670
671   BM_ITER_ELEM (f, &iter, e, BM_FACES_OF_EDGE) {
672     BM_face_normal_update(f);
673   }
674
675   BM_vert_normal_update(e->v1);
676   BM_vert_normal_update(e->v2);
677 }
678
679 static void bm_loop_normal_accum(const BMLoop *l, float no[3])
680 {
681   float vec1[3], vec2[3], fac;
682
683   /* Same calculation used in BM_mesh_normals_update */
684   sub_v3_v3v3(vec1, l->v->co, l->prev->v->co);
685   sub_v3_v3v3(vec2, l->next->v->co, l->v->co);
686   normalize_v3(vec1);
687   normalize_v3(vec2);
688
689   fac = saacos(-dot_v3v3(vec1, vec2));
690
691   madd_v3_v3fl(no, l->f->no, fac);
692 }
693
694 bool BM_vert_calc_normal_ex(const BMVert *v, const char hflag, float r_no[3])
695 {
696   int len = 0;
697
698   zero_v3(r_no);
699
700   if (v->e) {
701     const BMEdge *e = v->e;
702     do {
703       if (e->l) {
704         const BMLoop *l = e->l;
705         do {
706           if (l->v == v) {
707             if (BM_elem_flag_test(l->f, hflag)) {
708               bm_loop_normal_accum(l, r_no);
709               len++;
710             }
711           }
712         } while ((l = l->radial_next) != e->l);
713       }
714     } while ((e = bmesh_disk_edge_next(e, v)) != v->e);
715   }
716
717   if (len) {
718     normalize_v3(r_no);
719     return true;
720   }
721   else {
722     return false;
723   }
724 }
725
726 bool BM_vert_calc_normal(const BMVert *v, float r_no[3])
727 {
728   int len = 0;
729
730   zero_v3(r_no);
731
732   if (v->e) {
733     const BMEdge *e = v->e;
734     do {
735       if (e->l) {
736         const BMLoop *l = e->l;
737         do {
738           if (l->v == v) {
739             bm_loop_normal_accum(l, r_no);
740             len++;
741           }
742         } while ((l = l->radial_next) != e->l);
743       }
744     } while ((e = bmesh_disk_edge_next(e, v)) != v->e);
745   }
746
747   if (len) {
748     normalize_v3(r_no);
749     return true;
750   }
751   else {
752     return false;
753   }
754 }
755
756 void BM_vert_normal_update_all(BMVert *v)
757 {
758   int len = 0;
759
760   zero_v3(v->no);
761
762   if (v->e) {
763     const BMEdge *e = v->e;
764     do {
765       if (e->l) {
766         const BMLoop *l = e->l;
767         do {
768           if (l->v == v) {
769             BM_face_normal_update(l->f);
770             bm_loop_normal_accum(l, v->no);
771             len++;
772           }
773         } while ((l = l->radial_next) != e->l);
774       }
775     } while ((e = bmesh_disk_edge_next(e, v)) != v->e);
776   }
777
778   if (len) {
779     normalize_v3(v->no);
780   }
781 }
782
783 /**
784  * update a vert normal (but not the faces incident on it)
785  */
786 void BM_vert_normal_update(BMVert *v)
787 {
788   BM_vert_calc_normal(v, v->no);
789 }
790
791 /**
792  * \brief BMESH UPDATE FACE NORMAL
793  *
794  * Updates the stored normal for the
795  * given face. Requires that a buffer
796  * of sufficient length to store projected
797  * coordinates for all of the face's vertices
798  * is passed in as well.
799  */
800
801 float BM_face_calc_normal(const BMFace *f, float r_no[3])
802 {
803   BMLoop *l;
804
805   /* common cases first */
806   switch (f->len) {
807     case 4: {
808       const float *co1 = (l = BM_FACE_FIRST_LOOP(f))->v->co;
809       const float *co2 = (l = l->next)->v->co;
810       const float *co3 = (l = l->next)->v->co;
811       const float *co4 = (l->next)->v->co;
812
813       return normal_quad_v3(r_no, co1, co2, co3, co4);
814     }
815     case 3: {
816       const float *co1 = (l = BM_FACE_FIRST_LOOP(f))->v->co;
817       const float *co2 = (l = l->next)->v->co;
818       const float *co3 = (l->next)->v->co;
819
820       return normal_tri_v3(r_no, co1, co2, co3);
821     }
822     default: {
823       return bm_face_calc_poly_normal(f, r_no);
824     }
825   }
826 }
827 void BM_face_normal_update(BMFace *f)
828 {
829   BM_face_calc_normal(f, f->no);
830 }
831
832 /* exact same as 'BM_face_calc_normal' but accepts vertex coords */
833 float BM_face_calc_normal_vcos(const BMesh *bm,
834                                const BMFace *f,
835                                float r_no[3],
836                                float const (*vertexCos)[3])
837 {
838   BMLoop *l;
839
840   /* must have valid index data */
841   BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
842   (void)bm;
843
844   /* common cases first */
845   switch (f->len) {
846     case 4: {
847       const float *co1 = vertexCos[BM_elem_index_get((l = BM_FACE_FIRST_LOOP(f))->v)];
848       const float *co2 = vertexCos[BM_elem_index_get((l = l->next)->v)];
849       const float *co3 = vertexCos[BM_elem_index_get((l = l->next)->v)];
850       const float *co4 = vertexCos[BM_elem_index_get((l->next)->v)];
851
852       return normal_quad_v3(r_no, co1, co2, co3, co4);
853     }
854     case 3: {
855       const float *co1 = vertexCos[BM_elem_index_get((l = BM_FACE_FIRST_LOOP(f))->v)];
856       const float *co2 = vertexCos[BM_elem_index_get((l = l->next)->v)];
857       const float *co3 = vertexCos[BM_elem_index_get((l->next)->v)];
858
859       return normal_tri_v3(r_no, co1, co2, co3);
860     }
861     default: {
862       return bm_face_calc_poly_normal_vertex_cos(f, r_no, vertexCos);
863     }
864   }
865 }
866
867 /**
868  * Calculate a normal from a vertex cloud.
869  *
870  * \note We could make a higher quality version that takes all vertices into account.
871  * Currently it finds 4 outer most points returning it's normal.
872  */
873 void BM_verts_calc_normal_from_cloud_ex(
874     BMVert **varr, int varr_len, float r_normal[3], float r_center[3], int *r_index_tangent)
875 {
876   const float varr_len_inv = 1.0f / (float)varr_len;
877
878   /* Get the center point and collect vector array since we loop over these a lot. */
879   float center[3] = {0.0f, 0.0f, 0.0f};
880   for (int i = 0; i < varr_len; i++) {
881     madd_v3_v3fl(center, varr[i]->co, varr_len_inv);
882   }
883
884   /* Find the 'co_a' point from center. */
885   int co_a_index = 0;
886   const float *co_a = NULL;
887   {
888     float dist_sq_max = -1.0f;
889     for (int i = 0; i < varr_len; i++) {
890       const float dist_sq_test = len_squared_v3v3(varr[i]->co, center);
891       if (!(dist_sq_test <= dist_sq_max)) {
892         co_a = varr[i]->co;
893         co_a_index = i;
894         dist_sq_max = dist_sq_test;
895       }
896     }
897   }
898
899   float dir_a[3];
900   sub_v3_v3v3(dir_a, co_a, center);
901   normalize_v3(dir_a);
902
903   const float *co_b = NULL;
904   float dir_b[3] = {0.0f, 0.0f, 0.0f};
905   {
906     float dist_sq_max = -1.0f;
907     for (int i = 0; i < varr_len; i++) {
908       if (varr[i]->co == co_a) {
909         continue;
910       }
911       float dir_test[3];
912       sub_v3_v3v3(dir_test, varr[i]->co, center);
913       project_plane_normalized_v3_v3v3(dir_test, dir_test, dir_a);
914       const float dist_sq_test = len_squared_v3(dir_test);
915       if (!(dist_sq_test <= dist_sq_max)) {
916         co_b = varr[i]->co;
917         dist_sq_max = dist_sq_test;
918         copy_v3_v3(dir_b, dir_test);
919       }
920     }
921   }
922
923   if (varr_len <= 3) {
924     normal_tri_v3(r_normal, center, co_a, co_b);
925     goto finally;
926   }
927
928   normalize_v3(dir_b);
929
930   const float *co_a_opposite = NULL;
931   const float *co_b_opposite = NULL;
932
933   {
934     float dot_a_min = FLT_MAX;
935     float dot_b_min = FLT_MAX;
936     for (int i = 0; i < varr_len; i++) {
937       const float *co_test = varr[i]->co;
938       float dot_test;
939
940       if (co_test != co_a) {
941         dot_test = dot_v3v3(dir_a, co_test);
942         if (dot_test < dot_a_min) {
943           dot_a_min = dot_test;
944           co_a_opposite = co_test;
945         }
946       }
947
948       if (co_test != co_b) {
949         dot_test = dot_v3v3(dir_b, co_test);
950         if (dot_test < dot_b_min) {
951           dot_b_min = dot_test;
952           co_b_opposite = co_test;
953         }
954       }
955     }
956   }
957
958   normal_quad_v3(r_normal, co_a, co_b, co_a_opposite, co_b_opposite);
959
960 finally:
961   if (r_center != NULL) {
962     copy_v3_v3(r_center, center);
963   }
964   if (r_index_tangent != NULL) {
965     *r_index_tangent = co_a_index;
966   }
967 }
968
969 void BM_verts_calc_normal_from_cloud(BMVert **varr, int varr_len, float r_normal[3])
970 {
971   BM_verts_calc_normal_from_cloud_ex(varr, varr_len, r_normal, NULL, NULL);
972 }
973
974 /**
975  * Calculates the face subset normal.
976  */
977 float BM_face_calc_normal_subset(const BMLoop *l_first, const BMLoop *l_last, float r_no[3])
978 {
979   const float *v_prev, *v_curr;
980
981   /* Newell's Method */
982   const BMLoop *l_iter = l_first;
983   const BMLoop *l_term = l_last->next;
984
985   zero_v3(r_no);
986
987   v_prev = l_last->v->co;
988   do {
989     v_curr = l_iter->v->co;
990     add_newell_cross_v3_v3v3(r_no, v_prev, v_curr);
991     v_prev = v_curr;
992   } while ((l_iter = l_iter->next) != l_term);
993
994   return normalize_v3(r_no);
995 }
996
997 /* exact same as 'BM_face_calc_normal' but accepts vertex coords */
998 void BM_face_calc_center_median_vcos(const BMesh *bm,
999                                      const BMFace *f,
1000                                      float r_cent[3],
1001                                      float const (*vertexCos)[3])
1002 {
1003   /* must have valid index data */
1004   BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
1005   (void)bm;
1006
1007   bm_face_calc_poly_center_median_vertex_cos(f, r_cent, vertexCos);
1008 }
1009
1010 /**
1011  * \brief Face Flip Normal
1012  *
1013  * Reverses the winding of a face.
1014  * \note This updates the calculated normal.
1015  */
1016 void BM_face_normal_flip_ex(BMesh *bm,
1017                             BMFace *f,
1018                             const int cd_loop_mdisp_offset,
1019                             const bool use_loop_mdisp_flip)
1020 {
1021   bmesh_kernel_loop_reverse(bm, f, cd_loop_mdisp_offset, use_loop_mdisp_flip);
1022   negate_v3(f->no);
1023 }
1024
1025 void BM_face_normal_flip(BMesh *bm, BMFace *f)
1026 {
1027   const int cd_loop_mdisp_offset = CustomData_get_offset(&bm->ldata, CD_MDISPS);
1028   BM_face_normal_flip_ex(bm, f, cd_loop_mdisp_offset, true);
1029 }
1030
1031 /**
1032  * BM POINT IN FACE
1033  *
1034  * Projects co onto face f, and returns true if it is inside
1035  * the face bounds.
1036  *
1037  * \note this uses a best-axis projection test,
1038  * instead of projecting co directly into f's orientation space,
1039  * so there might be accuracy issues.
1040  */
1041 bool BM_face_point_inside_test(const BMFace *f, const float co[3])
1042 {
1043   float axis_mat[3][3];
1044   float(*projverts)[2] = BLI_array_alloca(projverts, f->len);
1045
1046   float co_2d[2];
1047   BMLoop *l_iter;
1048   int i;
1049
1050   BLI_assert(BM_face_is_normal_valid(f));
1051
1052   axis_dominant_v3_to_m3(axis_mat, f->no);
1053
1054   mul_v2_m3v3(co_2d, axis_mat, co);
1055
1056   for (i = 0, l_iter = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l_iter = l_iter->next) {
1057     mul_v2_m3v3(projverts[i], axis_mat, l_iter->v->co);
1058   }
1059
1060   return isect_point_poly_v2(co_2d, projverts, f->len, false);
1061 }
1062
1063 /**
1064  * \brief BMESH TRIANGULATE FACE
1065  *
1066  * Breaks all quads and ngons down to triangles.
1067  * It uses polyfill for the ngons splitting, and
1068  * the beautify operator when use_beauty is true.
1069  *
1070  * \param r_faces_new: if non-null, must be an array of BMFace pointers,
1071  * with a length equal to (f->len - 3). It will be filled with the new
1072  * triangles (not including the original triangle).
1073  *
1074  * \param r_faces_double: When newly created faces are duplicates of existing faces,
1075  * they're added to this list. Caller must handle de-duplication.
1076  * This is done because its possible _all_ faces exist already,
1077  * and in that case we would have to remove all faces including the one passed,
1078  * which causes complications adding/removing faces while looking over them.
1079  *
1080  * \note The number of faces is _almost_ always (f->len - 3),
1081  *       However there may be faces that already occupying the
1082  *       triangles we would make, so the caller must check \a r_faces_new_tot.
1083  *
1084  * \note use_tag tags new flags and edges.
1085  */
1086 void BM_face_triangulate(BMesh *bm,
1087                          BMFace *f,
1088                          BMFace **r_faces_new,
1089                          int *r_faces_new_tot,
1090                          BMEdge **r_edges_new,
1091                          int *r_edges_new_tot,
1092                          LinkNode **r_faces_double,
1093                          const int quad_method,
1094                          const int ngon_method,
1095                          const bool use_tag,
1096                          /* use for ngons only! */
1097                          MemArena *pf_arena,
1098
1099                          /* use for MOD_TRIANGULATE_NGON_BEAUTY only! */
1100                          struct Heap *pf_heap)
1101 {
1102   const int cd_loop_mdisp_offset = CustomData_get_offset(&bm->ldata, CD_MDISPS);
1103   const bool use_beauty = (ngon_method == MOD_TRIANGULATE_NGON_BEAUTY);
1104   BMLoop *l_first, *l_new;
1105   BMFace *f_new;
1106   int nf_i = 0;
1107   int ne_i = 0;
1108
1109   BLI_assert(BM_face_is_normal_valid(f));
1110
1111   /* ensure both are valid or NULL */
1112   BLI_assert((r_faces_new == NULL) == (r_faces_new_tot == NULL));
1113
1114   BLI_assert(f->len > 3);
1115
1116   {
1117     BMLoop **loops = BLI_array_alloca(loops, f->len);
1118     uint(*tris)[3] = BLI_array_alloca(tris, f->len);
1119     const int totfilltri = f->len - 2;
1120     const int last_tri = f->len - 3;
1121     int i;
1122     /* for mdisps */
1123     float f_center[3];
1124
1125     if (f->len == 4) {
1126       /* even though we're not using BLI_polyfill, fill in 'tris' and 'loops'
1127        * so we can share code to handle face creation afterwards. */
1128       BMLoop *l_v1, *l_v2;
1129
1130       l_first = BM_FACE_FIRST_LOOP(f);
1131
1132       switch (quad_method) {
1133         case MOD_TRIANGULATE_QUAD_FIXED: {
1134           l_v1 = l_first;
1135           l_v2 = l_first->next->next;
1136           break;
1137         }
1138         case MOD_TRIANGULATE_QUAD_ALTERNATE: {
1139           l_v1 = l_first->next;
1140           l_v2 = l_first->prev;
1141           break;
1142         }
1143         case MOD_TRIANGULATE_QUAD_SHORTEDGE:
1144         case MOD_TRIANGULATE_QUAD_BEAUTY:
1145         default: {
1146           BMLoop *l_v3, *l_v4;
1147           bool split_24;
1148
1149           l_v1 = l_first->next;
1150           l_v2 = l_first->next->next;
1151           l_v3 = l_first->prev;
1152           l_v4 = l_first;
1153
1154           if (quad_method == MOD_TRIANGULATE_QUAD_SHORTEDGE) {
1155             float d1, d2;
1156             d1 = len_squared_v3v3(l_v4->v->co, l_v2->v->co);
1157             d2 = len_squared_v3v3(l_v1->v->co, l_v3->v->co);
1158             split_24 = ((d2 - d1) > 0.0f);
1159           }
1160           else {
1161             /* first check if the quad is concave on either diagonal */
1162             const int flip_flag = is_quad_flip_v3(
1163                 l_v1->v->co, l_v2->v->co, l_v3->v->co, l_v4->v->co);
1164             if (UNLIKELY(flip_flag & (1 << 0))) {
1165               split_24 = true;
1166             }
1167             else if (UNLIKELY(flip_flag & (1 << 1))) {
1168               split_24 = false;
1169             }
1170             else {
1171               split_24 = (BM_verts_calc_rotate_beauty(l_v1->v, l_v2->v, l_v3->v, l_v4->v, 0, 0) >
1172                           0.0f);
1173             }
1174           }
1175
1176           /* named confusingly, l_v1 is in fact the second vertex */
1177           if (split_24) {
1178             l_v1 = l_v4;
1179             // l_v2 = l_v2;
1180           }
1181           else {
1182             // l_v1 = l_v1;
1183             l_v2 = l_v3;
1184           }
1185           break;
1186         }
1187       }
1188
1189       loops[0] = l_v1;
1190       loops[1] = l_v1->next;
1191       loops[2] = l_v2;
1192       loops[3] = l_v2->next;
1193
1194       ARRAY_SET_ITEMS(tris[0], 0, 1, 2);
1195       ARRAY_SET_ITEMS(tris[1], 0, 2, 3);
1196     }
1197     else {
1198       BMLoop *l_iter;
1199       float axis_mat[3][3];
1200       float(*projverts)[2] = BLI_array_alloca(projverts, f->len);
1201
1202       axis_dominant_v3_to_m3_negate(axis_mat, f->no);
1203
1204       for (i = 0, l_iter = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l_iter = l_iter->next) {
1205         loops[i] = l_iter;
1206         mul_v2_m3v3(projverts[i], axis_mat, l_iter->v->co);
1207       }
1208
1209       BLI_polyfill_calc_arena(projverts, f->len, 1, tris, pf_arena);
1210
1211       if (use_beauty) {
1212         BLI_polyfill_beautify(projverts, f->len, tris, pf_arena, pf_heap);
1213       }
1214
1215       BLI_memarena_clear(pf_arena);
1216     }
1217
1218     if (cd_loop_mdisp_offset != -1) {
1219       BM_face_calc_center_median(f, f_center);
1220     }
1221
1222     /* loop over calculated triangles and create new geometry */
1223     for (i = 0; i < totfilltri; i++) {
1224       BMLoop *l_tri[3] = {loops[tris[i][0]], loops[tris[i][1]], loops[tris[i][2]]};
1225
1226       BMVert *v_tri[3] = {l_tri[0]->v, l_tri[1]->v, l_tri[2]->v};
1227
1228       f_new = BM_face_create_verts(bm, v_tri, 3, f, BM_CREATE_NOP, true);
1229       l_new = BM_FACE_FIRST_LOOP(f_new);
1230
1231       BLI_assert(v_tri[0] == l_new->v);
1232
1233       /* check for duplicate */
1234       if (l_new->radial_next != l_new) {
1235         BMLoop *l_iter = l_new->radial_next;
1236         do {
1237           if (UNLIKELY((l_iter->f->len == 3) && (l_new->prev->v == l_iter->prev->v))) {
1238             /* Check the last tri because we swap last f_new with f at the end... */
1239             BLI_linklist_prepend(r_faces_double, (i != last_tri) ? f_new : f);
1240             break;
1241           }
1242         } while ((l_iter = l_iter->radial_next) != l_new);
1243       }
1244
1245       /* copy CD data */
1246       BM_elem_attrs_copy(bm, bm, l_tri[0], l_new);
1247       BM_elem_attrs_copy(bm, bm, l_tri[1], l_new->next);
1248       BM_elem_attrs_copy(bm, bm, l_tri[2], l_new->prev);
1249
1250       /* add all but the last face which is swapped and removed (below) */
1251       if (i != last_tri) {
1252         if (use_tag) {
1253           BM_elem_flag_enable(f_new, BM_ELEM_TAG);
1254         }
1255         if (r_faces_new) {
1256           r_faces_new[nf_i++] = f_new;
1257         }
1258       }
1259
1260       if (use_tag || r_edges_new) {
1261         /* new faces loops */
1262         BMLoop *l_iter;
1263
1264         l_iter = l_first = l_new;
1265         do {
1266           BMEdge *e = l_iter->e;
1267           /* Confusing! if its not a boundary now, we know it will be later since this will be an
1268            * edge of one of the new faces which we're in the middle of creating. */
1269           bool is_new_edge = (l_iter == l_iter->radial_next);
1270
1271           if (is_new_edge) {
1272             if (use_tag) {
1273               BM_elem_flag_enable(e, BM_ELEM_TAG);
1274             }
1275             if (r_edges_new) {
1276               r_edges_new[ne_i++] = e;
1277             }
1278           }
1279           /* note, never disable tag's */
1280         } while ((l_iter = l_iter->next) != l_first);
1281       }
1282
1283       if (cd_loop_mdisp_offset != -1) {
1284         float f_new_center[3];
1285         BM_face_calc_center_median(f_new, f_new_center);
1286         BM_face_interp_multires_ex(bm, f_new, f, f_new_center, f_center, cd_loop_mdisp_offset);
1287       }
1288     }
1289
1290     {
1291       /* we can't delete the real face, because some of the callers expect it to remain valid.
1292        * so swap data and delete the last created tri */
1293       bmesh_face_swap_data(f, f_new);
1294       BM_face_kill(bm, f_new);
1295     }
1296   }
1297   bm->elem_index_dirty |= BM_FACE;
1298
1299   if (r_faces_new_tot) {
1300     *r_faces_new_tot = nf_i;
1301   }
1302
1303   if (r_edges_new_tot) {
1304     *r_edges_new_tot = ne_i;
1305   }
1306 }
1307
1308 /**
1309  * each pair of loops defines a new edge, a split.  this function goes
1310  * through and sets pairs that are geometrically invalid to null.  a
1311  * split is invalid, if it forms a concave angle or it intersects other
1312  * edges in the face, or it intersects another split.  in the case of
1313  * intersecting splits, only the first of the set of intersecting
1314  * splits survives
1315  */
1316 void BM_face_splits_check_legal(BMesh *bm, BMFace *f, BMLoop *(*loops)[2], int len)
1317 {
1318   float out[2] = {-FLT_MAX, -FLT_MAX};
1319   float center[2] = {0.0f, 0.0f};
1320   float axis_mat[3][3];
1321   float(*projverts)[2] = BLI_array_alloca(projverts, f->len);
1322   const float *(*edgeverts)[2] = BLI_array_alloca(edgeverts, len);
1323   BMLoop *l;
1324   int i, i_prev, j;
1325
1326   BLI_assert(BM_face_is_normal_valid(f));
1327
1328   axis_dominant_v3_to_m3(axis_mat, f->no);
1329
1330   for (i = 0, l = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l = l->next) {
1331     mul_v2_m3v3(projverts[i], axis_mat, l->v->co);
1332     add_v2_v2(center, projverts[i]);
1333   }
1334
1335   /* first test for completely convex face */
1336   if (is_poly_convex_v2(projverts, f->len)) {
1337     return;
1338   }
1339
1340   mul_v2_fl(center, 1.0f / f->len);
1341
1342   for (i = 0, l = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l = l->next) {
1343     BM_elem_index_set(l, i); /* set_dirty */
1344
1345     /* center the projection for maximum accuracy */
1346     sub_v2_v2(projverts[i], center);
1347
1348     out[0] = max_ff(out[0], projverts[i][0]);
1349     out[1] = max_ff(out[1], projverts[i][1]);
1350   }
1351   bm->elem_index_dirty |= BM_LOOP;
1352
1353   /* ensure we are well outside the face bounds (value is arbitrary) */
1354   add_v2_fl(out, 1.0f);
1355
1356   for (i = 0; i < len; i++) {
1357     edgeverts[i][0] = projverts[BM_elem_index_get(loops[i][0])];
1358     edgeverts[i][1] = projverts[BM_elem_index_get(loops[i][1])];
1359   }
1360
1361   /* do convexity test */
1362   for (i = 0; i < len; i++) {
1363     float mid[2];
1364     mid_v2_v2v2(mid, edgeverts[i][0], edgeverts[i][1]);
1365
1366     int isect = 0;
1367     int j_prev;
1368     for (j = 0, j_prev = f->len - 1; j < f->len; j_prev = j++) {
1369       const float *f_edge[2] = {projverts[j_prev], projverts[j]};
1370       if (isect_seg_seg_v2(UNPACK2(f_edge), mid, out) == ISECT_LINE_LINE_CROSS) {
1371         isect++;
1372       }
1373     }
1374
1375     if (isect % 2 == 0) {
1376       loops[i][0] = NULL;
1377     }
1378   }
1379
1380 #define EDGE_SHARE_VERT(e1, e2) \
1381   ((ELEM((e1)[0], (e2)[0], (e2)[1])) || (ELEM((e1)[1], (e2)[0], (e2)[1])))
1382
1383   /* do line crossing tests */
1384   for (i = 0, i_prev = f->len - 1; i < f->len; i_prev = i++) {
1385     const float *f_edge[2] = {projverts[i_prev], projverts[i]};
1386     for (j = 0; j < len; j++) {
1387       if ((loops[j][0] != NULL) && !EDGE_SHARE_VERT(f_edge, edgeverts[j])) {
1388         if (isect_seg_seg_v2(UNPACK2(f_edge), UNPACK2(edgeverts[j])) == ISECT_LINE_LINE_CROSS) {
1389           loops[j][0] = NULL;
1390         }
1391       }
1392     }
1393   }
1394
1395   /* self intersect tests */
1396   for (i = 0; i < len; i++) {
1397     if (loops[i][0]) {
1398       for (j = i + 1; j < len; j++) {
1399         if ((loops[j][0] != NULL) && !EDGE_SHARE_VERT(edgeverts[i], edgeverts[j])) {
1400           if (isect_seg_seg_v2(UNPACK2(edgeverts[i]), UNPACK2(edgeverts[j])) ==
1401               ISECT_LINE_LINE_CROSS) {
1402             loops[i][0] = NULL;
1403             break;
1404           }
1405         }
1406       }
1407     }
1408   }
1409
1410 #undef EDGE_SHARE_VERT
1411 }
1412
1413 /**
1414  * This simply checks that the verts don't connect faces which would have more optimal splits.
1415  * but _not_ check for correctness.
1416  */
1417 void BM_face_splits_check_optimal(BMFace *f, BMLoop *(*loops)[2], int len)
1418 {
1419   int i;
1420
1421   for (i = 0; i < len; i++) {
1422     BMLoop *l_a_dummy, *l_b_dummy;
1423     if (f != BM_vert_pair_share_face_by_angle(
1424                  loops[i][0]->v, loops[i][1]->v, &l_a_dummy, &l_b_dummy, false)) {
1425       loops[i][0] = NULL;
1426     }
1427   }
1428 }
1429
1430 /**
1431  * Small utility functions for fast access
1432  *
1433  * faster alternative to:
1434  * BM_iter_as_array(bm, BM_VERTS_OF_FACE, f, (void **)v, 3);
1435  */
1436 void BM_face_as_array_vert_tri(BMFace *f, BMVert *r_verts[3])
1437 {
1438   BMLoop *l = BM_FACE_FIRST_LOOP(f);
1439
1440   BLI_assert(f->len == 3);
1441
1442   r_verts[0] = l->v;
1443   l = l->next;
1444   r_verts[1] = l->v;
1445   l = l->next;
1446   r_verts[2] = l->v;
1447 }
1448
1449 /**
1450  * faster alternative to:
1451  * BM_iter_as_array(bm, BM_VERTS_OF_FACE, f, (void **)v, 4);
1452  */
1453 void BM_face_as_array_vert_quad(BMFace *f, BMVert *r_verts[4])
1454 {
1455   BMLoop *l = BM_FACE_FIRST_LOOP(f);
1456
1457   BLI_assert(f->len == 4);
1458
1459   r_verts[0] = l->v;
1460   l = l->next;
1461   r_verts[1] = l->v;
1462   l = l->next;
1463   r_verts[2] = l->v;
1464   l = l->next;
1465   r_verts[3] = l->v;
1466 }
1467
1468 /**
1469  * Small utility functions for fast access
1470  *
1471  * faster alternative to:
1472  * BM_iter_as_array(bm, BM_LOOPS_OF_FACE, f, (void **)l, 3);
1473  */
1474 void BM_face_as_array_loop_tri(BMFace *f, BMLoop *r_loops[3])
1475 {
1476   BMLoop *l = BM_FACE_FIRST_LOOP(f);
1477
1478   BLI_assert(f->len == 3);
1479
1480   r_loops[0] = l;
1481   l = l->next;
1482   r_loops[1] = l;
1483   l = l->next;
1484   r_loops[2] = l;
1485 }
1486
1487 /**
1488  * faster alternative to:
1489  * BM_iter_as_array(bm, BM_LOOPS_OF_FACE, f, (void **)l, 4);
1490  */
1491 void BM_face_as_array_loop_quad(BMFace *f, BMLoop *r_loops[4])
1492 {
1493   BMLoop *l = BM_FACE_FIRST_LOOP(f);
1494
1495   BLI_assert(f->len == 4);
1496
1497   r_loops[0] = l;
1498   l = l->next;
1499   r_loops[1] = l;
1500   l = l->next;
1501   r_loops[2] = l;
1502   l = l->next;
1503   r_loops[3] = l;
1504 }
1505
1506 /**
1507  * \brief BM_mesh_calc_tessellation get the looptris and its number from a certain bmesh
1508  * \param looptris:
1509  *
1510  * \note \a looptris Must be pre-allocated to at least the size of given by: poly_to_tri_count
1511  */
1512 void BM_mesh_calc_tessellation(BMesh *bm, BMLoop *(*looptris)[3], int *r_looptris_tot)
1513 {
1514   /* use this to avoid locking pthread for _every_ polygon
1515    * and calling the fill function */
1516 #define USE_TESSFACE_SPEEDUP
1517
1518   /* this assumes all faces can be scan-filled, which isn't always true,
1519    * worst case we over alloc a little which is acceptable */
1520 #ifndef NDEBUG
1521   const int looptris_tot = poly_to_tri_count(bm->totface, bm->totloop);
1522 #endif
1523
1524   BMIter iter;
1525   BMFace *efa;
1526   int i = 0;
1527
1528   MemArena *arena = NULL;
1529
1530   BM_ITER_MESH (efa, &iter, bm, BM_FACES_OF_MESH) {
1531     /* don't consider two-edged faces */
1532     if (UNLIKELY(efa->len < 3)) {
1533       /* do nothing */
1534     }
1535
1536 #ifdef USE_TESSFACE_SPEEDUP
1537
1538     /* no need to ensure the loop order, we know its ok */
1539
1540     else if (efa->len == 3) {
1541 #  if 0
1542       int j;
1543       BM_ITER_ELEM_INDEX(l, &liter, efa, BM_LOOPS_OF_FACE, j) {
1544         looptris[i][j] = l;
1545       }
1546       i += 1;
1547 #  else
1548       /* more cryptic but faster */
1549       BMLoop *l;
1550       BMLoop **l_ptr = looptris[i++];
1551       l_ptr[0] = l = BM_FACE_FIRST_LOOP(efa);
1552       l_ptr[1] = l = l->next;
1553       l_ptr[2] = l->next;
1554 #  endif
1555     }
1556     else if (efa->len == 4) {
1557 #  if 0
1558       BMLoop *ltmp[4];
1559       int j;
1560       BLI_array_grow_items(looptris, 2);
1561       BM_ITER_ELEM_INDEX(l, &liter, efa, BM_LOOPS_OF_FACE, j) {
1562         ltmp[j] = l;
1563       }
1564
1565       looptris[i][0] = ltmp[0];
1566       looptris[i][1] = ltmp[1];
1567       looptris[i][2] = ltmp[2];
1568       i += 1;
1569
1570       looptris[i][0] = ltmp[0];
1571       looptris[i][1] = ltmp[2];
1572       looptris[i][2] = ltmp[3];
1573       i += 1;
1574 #  else
1575       /* more cryptic but faster */
1576       BMLoop *l;
1577       BMLoop **l_ptr_a = looptris[i++];
1578       BMLoop **l_ptr_b = looptris[i++];
1579       (l_ptr_a[0] = l_ptr_b[0] = l = BM_FACE_FIRST_LOOP(efa));
1580       (l_ptr_a[1] = l = l->next);
1581       (l_ptr_a[2] = l_ptr_b[1] = l = l->next);
1582       (l_ptr_b[2] = l->next);
1583 #  endif
1584
1585       if (UNLIKELY(is_quad_flip_v3_first_third_fast(
1586               l_ptr_a[0]->v->co, l_ptr_a[1]->v->co, l_ptr_a[2]->v->co, l_ptr_b[2]->v->co))) {
1587         /* flip out of degenerate 0-2 state. */
1588         l_ptr_a[2] = l_ptr_b[2];
1589         l_ptr_b[0] = l_ptr_a[1];
1590       }
1591     }
1592
1593 #endif /* USE_TESSFACE_SPEEDUP */
1594
1595     else {
1596       int j;
1597
1598       BMLoop *l_iter;
1599       BMLoop *l_first;
1600       BMLoop **l_arr;
1601
1602       float axis_mat[3][3];
1603       float(*projverts)[2];
1604       uint(*tris)[3];
1605
1606       const int totfilltri = efa->len - 2;
1607
1608       if (UNLIKELY(arena == NULL)) {
1609         arena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, __func__);
1610       }
1611
1612       tris = BLI_memarena_alloc(arena, sizeof(*tris) * totfilltri);
1613       l_arr = BLI_memarena_alloc(arena, sizeof(*l_arr) * efa->len);
1614       projverts = BLI_memarena_alloc(arena, sizeof(*projverts) * efa->len);
1615
1616       axis_dominant_v3_to_m3_negate(axis_mat, efa->no);
1617
1618       j = 0;
1619       l_iter = l_first = BM_FACE_FIRST_LOOP(efa);
1620       do {
1621         l_arr[j] = l_iter;
1622         mul_v2_m3v3(projverts[j], axis_mat, l_iter->v->co);
1623         j++;
1624       } while ((l_iter = l_iter->next) != l_first);
1625
1626       BLI_polyfill_calc_arena(projverts, efa->len, 1, tris, arena);
1627
1628       for (j = 0; j < totfilltri; j++) {
1629         BMLoop **l_ptr = looptris[i++];
1630         uint *tri = tris[j];
1631
1632         l_ptr[0] = l_arr[tri[0]];
1633         l_ptr[1] = l_arr[tri[1]];
1634         l_ptr[2] = l_arr[tri[2]];
1635       }
1636
1637       BLI_memarena_clear(arena);
1638     }
1639   }
1640
1641   if (arena) {
1642     BLI_memarena_free(arena);
1643     arena = NULL;
1644   }
1645
1646   *r_looptris_tot = i;
1647
1648   BLI_assert(i <= looptris_tot);
1649
1650 #undef USE_TESSFACE_SPEEDUP
1651 }
1652
1653 /**
1654  * A version of #BM_mesh_calc_tessellation that avoids degenerate triangles.
1655  */
1656 void BM_mesh_calc_tessellation_beauty(BMesh *bm, BMLoop *(*looptris)[3], int *r_looptris_tot)
1657 {
1658   /* this assumes all faces can be scan-filled, which isn't always true,
1659    * worst case we over alloc a little which is acceptable */
1660 #ifndef NDEBUG
1661   const int looptris_tot = poly_to_tri_count(bm->totface, bm->totloop);
1662 #endif
1663
1664   BMIter iter;
1665   BMFace *efa;
1666   int i = 0;
1667
1668   MemArena *pf_arena = NULL;
1669
1670   /* use_beauty */
1671   Heap *pf_heap = NULL;
1672
1673   BM_ITER_MESH (efa, &iter, bm, BM_FACES_OF_MESH) {
1674     /* don't consider two-edged faces */
1675     if (UNLIKELY(efa->len < 3)) {
1676       /* do nothing */
1677     }
1678     else if (efa->len == 3) {
1679       BMLoop *l;
1680       BMLoop **l_ptr = looptris[i++];
1681       l_ptr[0] = l = BM_FACE_FIRST_LOOP(efa);
1682       l_ptr[1] = l = l->next;
1683       l_ptr[2] = l->next;
1684     }
1685     else if (efa->len == 4) {
1686       BMLoop *l_v1 = BM_FACE_FIRST_LOOP(efa);
1687       BMLoop *l_v2 = l_v1->next;
1688       BMLoop *l_v3 = l_v2->next;
1689       BMLoop *l_v4 = l_v1->prev;
1690
1691       /* #BM_verts_calc_rotate_beauty performs excessive checks we don't need!
1692        * It's meant for rotating edges, it also calculates a new normal.
1693        *
1694        * Use #BLI_polyfill_beautify_quad_rotate_calc since we have the normal.
1695        */
1696 #if 0
1697       const bool split_13 = (BM_verts_calc_rotate_beauty(
1698                                  l_v1->v, l_v2->v, l_v3->v, l_v4->v, 0, 0) < 0.0f);
1699 #else
1700       float axis_mat[3][3], v_quad[4][2];
1701       axis_dominant_v3_to_m3(axis_mat, efa->no);
1702       mul_v2_m3v3(v_quad[0], axis_mat, l_v1->v->co);
1703       mul_v2_m3v3(v_quad[1], axis_mat, l_v2->v->co);
1704       mul_v2_m3v3(v_quad[2], axis_mat, l_v3->v->co);
1705       mul_v2_m3v3(v_quad[3], axis_mat, l_v4->v->co);
1706
1707       const bool split_13 = BLI_polyfill_beautify_quad_rotate_calc(
1708                                 v_quad[0], v_quad[1], v_quad[2], v_quad[3]) < 0.0f;
1709 #endif
1710
1711       BMLoop **l_ptr_a = looptris[i++];
1712       BMLoop **l_ptr_b = looptris[i++];
1713       if (split_13) {
1714         l_ptr_a[0] = l_v1;
1715         l_ptr_a[1] = l_v2;
1716         l_ptr_a[2] = l_v3;
1717
1718         l_ptr_b[0] = l_v1;
1719         l_ptr_b[1] = l_v3;
1720         l_ptr_b[2] = l_v4;
1721       }
1722       else {
1723         l_ptr_a[0] = l_v1;
1724         l_ptr_a[1] = l_v2;
1725         l_ptr_a[2] = l_v4;
1726
1727         l_ptr_b[0] = l_v2;
1728         l_ptr_b[1] = l_v3;
1729         l_ptr_b[2] = l_v4;
1730       }
1731     }
1732     else {
1733       int j;
1734
1735       BMLoop *l_iter;
1736       BMLoop *l_first;
1737       BMLoop **l_arr;
1738
1739       float axis_mat[3][3];
1740       float(*projverts)[2];
1741       unsigned int(*tris)[3];
1742
1743       const int totfilltri = efa->len - 2;
1744
1745       if (UNLIKELY(pf_arena == NULL)) {
1746         pf_arena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, __func__);
1747         pf_heap = BLI_heap_new_ex(BLI_POLYFILL_ALLOC_NGON_RESERVE);
1748       }
1749
1750       tris = BLI_memarena_alloc(pf_arena, sizeof(*tris) * totfilltri);
1751       l_arr = BLI_memarena_alloc(pf_arena, sizeof(*l_arr) * efa->len);
1752       projverts = BLI_memarena_alloc(pf_arena, sizeof(*projverts) * efa->len);
1753
1754       axis_dominant_v3_to_m3_negate(axis_mat, efa->no);
1755
1756       j = 0;
1757       l_iter = l_first = BM_FACE_FIRST_LOOP(efa);
1758       do {
1759         l_arr[j] = l_iter;
1760         mul_v2_m3v3(projverts[j], axis_mat, l_iter->v->co);
1761         j++;
1762       } while ((l_iter = l_iter->next) != l_first);
1763
1764       BLI_polyfill_calc_arena(projverts, efa->len, 1, tris, pf_arena);
1765
1766       BLI_polyfill_beautify(projverts, efa->len, tris, pf_arena, pf_heap);
1767
1768       for (j = 0; j < totfilltri; j++) {
1769         BMLoop **l_ptr = looptris[i++];
1770         unsigned int *tri = tris[j];
1771
1772         l_ptr[0] = l_arr[tri[0]];
1773         l_ptr[1] = l_arr[tri[1]];
1774         l_ptr[2] = l_arr[tri[2]];
1775       }
1776
1777       BLI_memarena_clear(pf_arena);
1778     }
1779   }
1780
1781   if (pf_arena) {
1782     BLI_memarena_free(pf_arena);
1783
1784     BLI_heap_free(pf_heap, NULL);
1785   }
1786
1787   *r_looptris_tot = i;
1788
1789   BLI_assert(i <= looptris_tot);
1790 }