Curve: Remove duplicated and confusing argument
[blender.git] / source / blender / blenkernel / intern / curve.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  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
17  * All rights reserved.
18  */
19
20 /** \file
21  * \ingroup bke
22  */
23
24 #include <math.h>  // floor
25 #include <string.h>
26 #include <stdlib.h>
27
28 #include "MEM_guardedalloc.h"
29
30 #include "BLI_blenlib.h"
31 #include "BLI_math.h"
32 #include "BLI_utildefines.h"
33 #include "BLI_ghash.h"
34 #include "BLI_linklist.h"
35
36 #include "DNA_anim_types.h"
37 #include "DNA_curve_types.h"
38 #include "DNA_material_types.h"
39
40 /* for dereferencing pointers */
41 #include "DNA_key_types.h"
42 #include "DNA_scene_types.h"
43 #include "DNA_vfont_types.h"
44 #include "DNA_object_types.h"
45
46 #include "BKE_animsys.h"
47 #include "BKE_curve.h"
48 #include "BKE_displist.h"
49 #include "BKE_font.h"
50 #include "BKE_key.h"
51 #include "BKE_library.h"
52 #include "BKE_main.h"
53 #include "BKE_object.h"
54 #include "BKE_material.h"
55
56 #include "DEG_depsgraph.h"
57
58 #include "CLG_log.h"
59
60 /* globals */
61
62 /* local */
63 static CLG_LogRef LOG = {"bke.curve"};
64
65 static int cu_isectLL(const float v1[3],
66                       const float v2[3],
67                       const float v3[3],
68                       const float v4[3],
69                       short cox,
70                       short coy,
71                       float *lambda,
72                       float *mu,
73                       float vec[3]);
74
75 /* frees editcurve entirely */
76 void BKE_curve_editfont_free(Curve *cu)
77 {
78   if (cu->editfont) {
79     EditFont *ef = cu->editfont;
80
81     if (ef->textbuf) {
82       MEM_freeN(ef->textbuf);
83     }
84     if (ef->textbufinfo) {
85       MEM_freeN(ef->textbufinfo);
86     }
87     if (ef->selboxes) {
88       MEM_freeN(ef->selboxes);
89     }
90
91     MEM_freeN(ef);
92     cu->editfont = NULL;
93   }
94 }
95
96 static void curve_editNurb_keyIndex_cv_free_cb(void *val)
97 {
98   CVKeyIndex *index = val;
99   MEM_freeN(index->orig_cv);
100   MEM_freeN(val);
101 }
102
103 void BKE_curve_editNurb_keyIndex_delCV(GHash *keyindex, const void *cv)
104 {
105   BLI_assert(keyindex != NULL);
106   BLI_ghash_remove(keyindex, cv, NULL, curve_editNurb_keyIndex_cv_free_cb);
107 }
108
109 void BKE_curve_editNurb_keyIndex_free(GHash **keyindex)
110 {
111   if (!(*keyindex)) {
112     return;
113   }
114   BLI_ghash_free(*keyindex, NULL, curve_editNurb_keyIndex_cv_free_cb);
115   *keyindex = NULL;
116 }
117
118 void BKE_curve_editNurb_free(Curve *cu)
119 {
120   if (cu->editnurb) {
121     BKE_nurbList_free(&cu->editnurb->nurbs);
122     BKE_curve_editNurb_keyIndex_free(&cu->editnurb->keyindex);
123     MEM_freeN(cu->editnurb);
124     cu->editnurb = NULL;
125   }
126 }
127
128 /** Free (or release) any data used by this curve (does not free the curve itself). */
129 void BKE_curve_free(Curve *cu)
130 {
131   BKE_animdata_free((ID *)cu, false);
132
133   BKE_curve_batch_cache_free(cu);
134
135   BKE_nurbList_free(&cu->nurb);
136   BKE_curve_editfont_free(cu);
137
138   BKE_curve_editNurb_free(cu);
139
140   MEM_SAFE_FREE(cu->mat);
141   MEM_SAFE_FREE(cu->str);
142   MEM_SAFE_FREE(cu->strinfo);
143   MEM_SAFE_FREE(cu->bb);
144   MEM_SAFE_FREE(cu->tb);
145 }
146
147 void BKE_curve_init(Curve *cu)
148 {
149   /* BLI_assert(MEMCMP_STRUCT_AFTER_IS_ZERO(cu, id)); */ /* cu->type is already initialized... */
150
151   copy_v3_fl(cu->size, 1.0f);
152   cu->flag = CU_DEFORM_BOUNDS_OFF | CU_PATH_RADIUS;
153   cu->pathlen = 100;
154   cu->resolu = cu->resolv = (cu->type == OB_SURF) ? 4 : 12;
155   cu->width = 1.0;
156   cu->wordspace = 1.0;
157   cu->spacing = cu->linedist = 1.0;
158   cu->fsize = 1.0;
159   cu->ulheight = 0.05;
160   cu->texflag = CU_AUTOSPACE;
161   cu->smallcaps_scale = 0.75f;
162   /* XXX: this one seems to be the best one in most cases, at least for curve deform... */
163   cu->twist_mode = CU_TWIST_MINIMUM;
164   cu->bevfac1 = 0.0f;
165   cu->bevfac2 = 1.0f;
166   cu->bevfac1_mapping = CU_BEVFAC_MAP_RESOLU;
167   cu->bevfac2_mapping = CU_BEVFAC_MAP_RESOLU;
168   cu->bevresol = 4;
169
170   cu->bb = BKE_boundbox_alloc_unit();
171
172   if (cu->type == OB_FONT) {
173     cu->flag |= CU_FRONT | CU_BACK;
174     cu->vfont = cu->vfontb = cu->vfonti = cu->vfontbi = BKE_vfont_builtin_get();
175     cu->vfont->id.us += 4;
176     cu->str = MEM_malloc_arrayN(12, sizeof(unsigned char), "str");
177     BLI_strncpy(cu->str, "Text", 12);
178     cu->len = cu->len_wchar = cu->pos = 4;
179     cu->strinfo = MEM_calloc_arrayN(12, sizeof(CharInfo), "strinfo new");
180     cu->totbox = cu->actbox = 1;
181     cu->tb = MEM_calloc_arrayN(MAXTEXTBOX, sizeof(TextBox), "textbox");
182     cu->tb[0].w = cu->tb[0].h = 0.0;
183   }
184 }
185
186 Curve *BKE_curve_add(Main *bmain, const char *name, int type)
187 {
188   Curve *cu;
189
190   cu = BKE_libblock_alloc(bmain, ID_CU, name, 0);
191   cu->type = type;
192
193   BKE_curve_init(cu);
194
195   return cu;
196 }
197
198 /**
199  * Only copy internal data of Curve ID from source
200  * to already allocated/initialized destination.
201  * You probably never want to use that directly,
202  * use #BKE_id_copy or #BKE_id_copy_ex for typical needs.
203  *
204  * WARNING! This function will not handle ID user count!
205  *
206  * \param flag: Copying options (see BKE_library.h's LIB_ID_COPY_... flags for more).
207  */
208 void BKE_curve_copy_data(Main *bmain, Curve *cu_dst, const Curve *cu_src, const int flag)
209 {
210   BLI_listbase_clear(&cu_dst->nurb);
211   BKE_nurbList_duplicate(&(cu_dst->nurb), &(cu_src->nurb));
212
213   cu_dst->mat = MEM_dupallocN(cu_src->mat);
214
215   cu_dst->str = MEM_dupallocN(cu_src->str);
216   cu_dst->strinfo = MEM_dupallocN(cu_src->strinfo);
217   cu_dst->tb = MEM_dupallocN(cu_src->tb);
218   cu_dst->bb = MEM_dupallocN(cu_src->bb);
219   cu_dst->batch_cache = NULL;
220
221   if (cu_src->key && (flag & LIB_ID_COPY_SHAPEKEY)) {
222     BKE_id_copy_ex(bmain, &cu_src->key->id, (ID **)&cu_dst->key, flag);
223   }
224
225   cu_dst->editnurb = NULL;
226   cu_dst->editfont = NULL;
227 }
228
229 Curve *BKE_curve_copy(Main *bmain, const Curve *cu)
230 {
231   Curve *cu_copy;
232   BKE_id_copy(bmain, &cu->id, (ID **)&cu_copy);
233   return cu_copy;
234 }
235
236 void BKE_curve_make_local(Main *bmain, Curve *cu, const bool lib_local)
237 {
238   BKE_id_make_local_generic(bmain, &cu->id, true, lib_local);
239 }
240
241 /* Get list of nurbs from editnurbs structure */
242 ListBase *BKE_curve_editNurbs_get(Curve *cu)
243 {
244   if (cu->editnurb) {
245     return &cu->editnurb->nurbs;
246   }
247
248   return NULL;
249 }
250
251 short BKE_curve_type_get(Curve *cu)
252 {
253   Nurb *nu;
254   int type = cu->type;
255
256   if (cu->vfont) {
257     return OB_FONT;
258   }
259
260   if (!cu->type) {
261     type = OB_CURVE;
262
263     for (nu = cu->nurb.first; nu; nu = nu->next) {
264       if (nu->pntsv > 1) {
265         type = OB_SURF;
266       }
267     }
268   }
269
270   return type;
271 }
272
273 void BKE_curve_curve_dimension_update(Curve *cu)
274 {
275   ListBase *nurbs = BKE_curve_nurbs_get(cu);
276   Nurb *nu = nurbs->first;
277
278   if (cu->flag & CU_3D) {
279     for (; nu; nu = nu->next) {
280       nu->flag &= ~CU_2D;
281     }
282   }
283   else {
284     for (; nu; nu = nu->next) {
285       nu->flag |= CU_2D;
286       BKE_nurb_test_2d(nu);
287
288       /* since the handles are moved they need to be auto-located again */
289       if (nu->type == CU_BEZIER) {
290         BKE_nurb_handles_calc(nu);
291       }
292     }
293   }
294 }
295
296 void BKE_curve_type_test(Object *ob)
297 {
298   ob->type = BKE_curve_type_get(ob->data);
299
300   if (ob->type == OB_CURVE) {
301     BKE_curve_curve_dimension_update((Curve *)ob->data);
302   }
303 }
304
305 void BKE_curve_boundbox_calc(Curve *cu, float r_loc[3], float r_size[3])
306 {
307   BoundBox *bb;
308   float min[3], max[3];
309   float mloc[3], msize[3];
310
311   if (cu->bb == NULL) {
312     cu->bb = MEM_callocN(sizeof(BoundBox), "boundbox");
313   }
314   bb = cu->bb;
315
316   if (!r_loc) {
317     r_loc = mloc;
318   }
319   if (!r_size) {
320     r_size = msize;
321   }
322
323   INIT_MINMAX(min, max);
324   if (!BKE_curve_minmax(cu, true, min, max)) {
325     min[0] = min[1] = min[2] = -1.0f;
326     max[0] = max[1] = max[2] = 1.0f;
327   }
328
329   mid_v3_v3v3(r_loc, min, max);
330
331   r_size[0] = (max[0] - min[0]) / 2.0f;
332   r_size[1] = (max[1] - min[1]) / 2.0f;
333   r_size[2] = (max[2] - min[2]) / 2.0f;
334
335   BKE_boundbox_init_from_minmax(bb, min, max);
336
337   bb->flag &= ~BOUNDBOX_DIRTY;
338 }
339
340 BoundBox *BKE_curve_boundbox_get(Object *ob)
341 {
342   /* This is Object-level data access,
343    * DO NOT touch to Mesh's bb, would be totally thread-unsafe. */
344   if (ob->runtime.bb == NULL || ob->runtime.bb->flag & BOUNDBOX_DIRTY) {
345     Curve *cu = ob->data;
346     float min[3], max[3];
347
348     INIT_MINMAX(min, max);
349     BKE_curve_minmax(cu, true, min, max);
350
351     if (ob->runtime.bb == NULL) {
352       ob->runtime.bb = MEM_mallocN(sizeof(*ob->runtime.bb), __func__);
353     }
354     BKE_boundbox_init_from_minmax(ob->runtime.bb, min, max);
355     ob->runtime.bb->flag &= ~BOUNDBOX_DIRTY;
356   }
357
358   return ob->runtime.bb;
359 }
360
361 void BKE_curve_texspace_calc(Curve *cu)
362 {
363   float loc[3], size[3];
364   int a;
365
366   BKE_curve_boundbox_calc(cu, loc, size);
367
368   if (cu->texflag & CU_AUTOSPACE) {
369     for (a = 0; a < 3; a++) {
370       if (size[a] == 0.0f) {
371         size[a] = 1.0f;
372       }
373       else if (size[a] > 0.0f && size[a] < 0.00001f) {
374         size[a] = 0.00001f;
375       }
376       else if (size[a] < 0.0f && size[a] > -0.00001f) {
377         size[a] = -0.00001f;
378       }
379     }
380
381     copy_v3_v3(cu->loc, loc);
382     copy_v3_v3(cu->size, size);
383     zero_v3(cu->rot);
384   }
385 }
386
387 BoundBox *BKE_curve_texspace_get(Curve *cu, float r_loc[3], float r_rot[3], float r_size[3])
388 {
389   if (cu->bb == NULL || (cu->bb->flag & BOUNDBOX_DIRTY)) {
390     BKE_curve_texspace_calc(cu);
391   }
392
393   if (r_loc) {
394     copy_v3_v3(r_loc, cu->loc);
395   }
396   if (r_rot) {
397     copy_v3_v3(r_rot, cu->rot);
398   }
399   if (r_size) {
400     copy_v3_v3(r_size, cu->size);
401   }
402
403   return cu->bb;
404 }
405
406 bool BKE_nurbList_index_get_co(ListBase *nurb, const int index, float r_co[3])
407 {
408   Nurb *nu;
409   int tot = 0;
410
411   for (nu = nurb->first; nu; nu = nu->next) {
412     int tot_nu;
413     if (nu->type == CU_BEZIER) {
414       tot_nu = nu->pntsu;
415       if (index - tot < tot_nu) {
416         copy_v3_v3(r_co, nu->bezt[index - tot].vec[1]);
417         return true;
418       }
419     }
420     else {
421       tot_nu = nu->pntsu * nu->pntsv;
422       if (index - tot < tot_nu) {
423         copy_v3_v3(r_co, nu->bp[index - tot].vec);
424         return true;
425       }
426     }
427     tot += tot_nu;
428   }
429
430   return false;
431 }
432
433 int BKE_nurbList_verts_count(ListBase *nurb)
434 {
435   Nurb *nu;
436   int tot = 0;
437
438   nu = nurb->first;
439   while (nu) {
440     if (nu->bezt) {
441       tot += 3 * nu->pntsu;
442     }
443     else if (nu->bp) {
444       tot += nu->pntsu * nu->pntsv;
445     }
446
447     nu = nu->next;
448   }
449   return tot;
450 }
451
452 int BKE_nurbList_verts_count_without_handles(ListBase *nurb)
453 {
454   Nurb *nu;
455   int tot = 0;
456
457   nu = nurb->first;
458   while (nu) {
459     if (nu->bezt) {
460       tot += nu->pntsu;
461     }
462     else if (nu->bp) {
463       tot += nu->pntsu * nu->pntsv;
464     }
465
466     nu = nu->next;
467   }
468   return tot;
469 }
470
471 /* **************** NURBS ROUTINES ******************** */
472
473 void BKE_nurb_free(Nurb *nu)
474 {
475
476   if (nu == NULL) {
477     return;
478   }
479
480   if (nu->bezt) {
481     MEM_freeN(nu->bezt);
482   }
483   nu->bezt = NULL;
484   if (nu->bp) {
485     MEM_freeN(nu->bp);
486   }
487   nu->bp = NULL;
488   if (nu->knotsu) {
489     MEM_freeN(nu->knotsu);
490   }
491   nu->knotsu = NULL;
492   if (nu->knotsv) {
493     MEM_freeN(nu->knotsv);
494   }
495   nu->knotsv = NULL;
496   /* if (nu->trim.first) freeNurblist(&(nu->trim)); */
497
498   MEM_freeN(nu);
499 }
500
501 void BKE_nurbList_free(ListBase *lb)
502 {
503   Nurb *nu, *next;
504
505   if (lb == NULL) {
506     return;
507   }
508
509   nu = lb->first;
510   while (nu) {
511     next = nu->next;
512     BKE_nurb_free(nu);
513     nu = next;
514   }
515   BLI_listbase_clear(lb);
516 }
517
518 Nurb *BKE_nurb_duplicate(const Nurb *nu)
519 {
520   Nurb *newnu;
521   int len;
522
523   newnu = (Nurb *)MEM_mallocN(sizeof(Nurb), "duplicateNurb");
524   if (newnu == NULL) {
525     return NULL;
526   }
527   memcpy(newnu, nu, sizeof(Nurb));
528
529   if (nu->bezt) {
530     newnu->bezt = (BezTriple *)MEM_malloc_arrayN(nu->pntsu, sizeof(BezTriple), "duplicateNurb2");
531     memcpy(newnu->bezt, nu->bezt, nu->pntsu * sizeof(BezTriple));
532   }
533   else {
534     len = nu->pntsu * nu->pntsv;
535     newnu->bp = (BPoint *)MEM_malloc_arrayN(len, sizeof(BPoint), "duplicateNurb3");
536     memcpy(newnu->bp, nu->bp, len * sizeof(BPoint));
537
538     newnu->knotsu = newnu->knotsv = NULL;
539
540     if (nu->knotsu) {
541       len = KNOTSU(nu);
542       if (len) {
543         newnu->knotsu = MEM_malloc_arrayN(len, sizeof(float), "duplicateNurb4");
544         memcpy(newnu->knotsu, nu->knotsu, sizeof(float) * len);
545       }
546     }
547     if (nu->pntsv > 1 && nu->knotsv) {
548       len = KNOTSV(nu);
549       if (len) {
550         newnu->knotsv = MEM_malloc_arrayN(len, sizeof(float), "duplicateNurb5");
551         memcpy(newnu->knotsv, nu->knotsv, sizeof(float) * len);
552       }
553     }
554   }
555   return newnu;
556 }
557
558 /* copy the nurb but allow for different number of points (to be copied after this) */
559 Nurb *BKE_nurb_copy(Nurb *src, int pntsu, int pntsv)
560 {
561   Nurb *newnu = (Nurb *)MEM_mallocN(sizeof(Nurb), "copyNurb");
562   memcpy(newnu, src, sizeof(Nurb));
563
564   if (pntsu == 1) {
565     SWAP(int, pntsu, pntsv);
566   }
567   newnu->pntsu = pntsu;
568   newnu->pntsv = pntsv;
569
570   /* caller can manually handle these arrays */
571   newnu->knotsu = NULL;
572   newnu->knotsv = NULL;
573
574   if (src->bezt) {
575     newnu->bezt = (BezTriple *)MEM_malloc_arrayN(pntsu * pntsv, sizeof(BezTriple), "copyNurb2");
576   }
577   else {
578     newnu->bp = (BPoint *)MEM_malloc_arrayN(pntsu * pntsv, sizeof(BPoint), "copyNurb3");
579   }
580
581   return newnu;
582 }
583
584 void BKE_nurbList_duplicate(ListBase *lb1, const ListBase *lb2)
585 {
586   Nurb *nu, *nun;
587
588   BKE_nurbList_free(lb1);
589
590   nu = lb2->first;
591   while (nu) {
592     nun = BKE_nurb_duplicate(nu);
593     BLI_addtail(lb1, nun);
594
595     nu = nu->next;
596   }
597 }
598
599 void BKE_nurb_test_2d(Nurb *nu)
600 {
601   BezTriple *bezt;
602   BPoint *bp;
603   int a;
604
605   if ((nu->flag & CU_2D) == 0) {
606     return;
607   }
608
609   if (nu->type == CU_BEZIER) {
610     a = nu->pntsu;
611     bezt = nu->bezt;
612     while (a--) {
613       bezt->vec[0][2] = 0.0;
614       bezt->vec[1][2] = 0.0;
615       bezt->vec[2][2] = 0.0;
616       bezt++;
617     }
618   }
619   else {
620     a = nu->pntsu * nu->pntsv;
621     bp = nu->bp;
622     while (a--) {
623       bp->vec[2] = 0.0;
624       bp++;
625     }
626   }
627 }
628
629 /**
630  * if use_radius is truth, minmax will take points' radius into account,
631  * which will make boundbox closer to beveled curve.
632  */
633 void BKE_nurb_minmax(Nurb *nu, bool use_radius, float min[3], float max[3])
634 {
635   BezTriple *bezt;
636   BPoint *bp;
637   int a;
638   float point[3];
639
640   if (nu->type == CU_BEZIER) {
641     a = nu->pntsu;
642     bezt = nu->bezt;
643     while (a--) {
644       if (use_radius) {
645         float radius_vector[3];
646         radius_vector[0] = radius_vector[1] = radius_vector[2] = bezt->radius;
647
648         add_v3_v3v3(point, bezt->vec[1], radius_vector);
649         minmax_v3v3_v3(min, max, point);
650
651         sub_v3_v3v3(point, bezt->vec[1], radius_vector);
652         minmax_v3v3_v3(min, max, point);
653       }
654       else {
655         minmax_v3v3_v3(min, max, bezt->vec[1]);
656       }
657       minmax_v3v3_v3(min, max, bezt->vec[0]);
658       minmax_v3v3_v3(min, max, bezt->vec[2]);
659       bezt++;
660     }
661   }
662   else {
663     a = nu->pntsu * nu->pntsv;
664     bp = nu->bp;
665     while (a--) {
666       if (nu->pntsv == 1 && use_radius) {
667         float radius_vector[3];
668         radius_vector[0] = radius_vector[1] = radius_vector[2] = bp->radius;
669
670         add_v3_v3v3(point, bp->vec, radius_vector);
671         minmax_v3v3_v3(min, max, point);
672
673         sub_v3_v3v3(point, bp->vec, radius_vector);
674         minmax_v3v3_v3(min, max, point);
675       }
676       else {
677         /* Surfaces doesn't use bevel, so no need to take radius into account. */
678         minmax_v3v3_v3(min, max, bp->vec);
679       }
680       bp++;
681     }
682   }
683 }
684
685 float BKE_nurb_calc_length(const Nurb *nu, int resolution)
686 {
687   BezTriple *bezt, *prevbezt;
688   BPoint *bp, *prevbp;
689   int a, b;
690   float length = 0.0f;
691   int resolu = resolution ? resolution : nu->resolu;
692   int pntsu = nu->pntsu;
693   float *points, *pntsit, *prevpntsit;
694
695   if (nu->type == CU_POLY) {
696     a = nu->pntsu - 1;
697     bp = nu->bp;
698     if (nu->flagu & CU_NURB_CYCLIC) {
699       ++a;
700       prevbp = nu->bp + (nu->pntsu - 1);
701     }
702     else {
703       prevbp = bp;
704       bp++;
705     }
706
707     while (a--) {
708       length += len_v3v3(prevbp->vec, bp->vec);
709       prevbp = bp;
710       ++bp;
711     }
712   }
713   else if (nu->type == CU_BEZIER) {
714     points = MEM_mallocN(sizeof(float[3]) * (resolu + 1), "getLength_bezier");
715     a = nu->pntsu - 1;
716     bezt = nu->bezt;
717     if (nu->flagu & CU_NURB_CYCLIC) {
718       ++a;
719       prevbezt = nu->bezt + (nu->pntsu - 1);
720     }
721     else {
722       prevbezt = bezt;
723       ++bezt;
724     }
725
726     while (a--) {
727       if (prevbezt->h2 == HD_VECT && bezt->h1 == HD_VECT) {
728         length += len_v3v3(prevbezt->vec[1], bezt->vec[1]);
729       }
730       else {
731         for (int j = 0; j < 3; j++) {
732           BKE_curve_forward_diff_bezier(prevbezt->vec[1][j],
733                                         prevbezt->vec[2][j],
734                                         bezt->vec[0][j],
735                                         bezt->vec[1][j],
736                                         points + j,
737                                         resolu,
738                                         3 * sizeof(float));
739         }
740
741         prevpntsit = pntsit = points;
742         b = resolu;
743         while (b--) {
744           pntsit += 3;
745           length += len_v3v3(prevpntsit, pntsit);
746           prevpntsit = pntsit;
747         }
748       }
749       prevbezt = bezt;
750       ++bezt;
751     }
752
753     MEM_freeN(points);
754   }
755   else if (nu->type == CU_NURBS) {
756     if (nu->pntsv == 1) {
757       /* important to zero for BKE_nurb_makeCurve. */
758       points = MEM_callocN(sizeof(float[3]) * pntsu * resolu, "getLength_nurbs");
759
760       BKE_nurb_makeCurve(nu, points, NULL, NULL, NULL, resolu, sizeof(float[3]));
761
762       if (nu->flagu & CU_NURB_CYCLIC) {
763         b = pntsu * resolu + 1;
764         prevpntsit = points + 3 * (pntsu * resolu - 1);
765         pntsit = points;
766       }
767       else {
768         b = (pntsu - 1) * resolu;
769         prevpntsit = points;
770         pntsit = points + 3;
771       }
772
773       while (--b) {
774         length += len_v3v3(prevpntsit, pntsit);
775         prevpntsit = pntsit;
776         pntsit += 3;
777       }
778
779       MEM_freeN(points);
780     }
781   }
782
783   return length;
784 }
785
786 /* be sure to call makeknots after this */
787 void BKE_nurb_points_add(Nurb *nu, int number)
788 {
789   BPoint *bp;
790   int i;
791
792   nu->bp = MEM_recallocN(nu->bp, (nu->pntsu + number) * sizeof(BPoint));
793
794   for (i = 0, bp = &nu->bp[nu->pntsu]; i < number; i++, bp++) {
795     bp->radius = 1.0f;
796   }
797
798   nu->pntsu += number;
799 }
800
801 void BKE_nurb_bezierPoints_add(Nurb *nu, int number)
802 {
803   BezTriple *bezt;
804   int i;
805
806   nu->bezt = MEM_recallocN(nu->bezt, (nu->pntsu + number) * sizeof(BezTriple));
807
808   for (i = 0, bezt = &nu->bezt[nu->pntsu]; i < number; i++, bezt++) {
809     bezt->radius = 1.0f;
810   }
811
812   nu->pntsu += number;
813 }
814
815 int BKE_nurb_index_from_uv(Nurb *nu, int u, int v)
816 {
817   const int totu = nu->pntsu;
818   const int totv = nu->pntsv;
819
820   if (nu->flagu & CU_NURB_CYCLIC) {
821     u = mod_i(u, totu);
822   }
823   else if (u < 0 || u >= totu) {
824     return -1;
825   }
826
827   if (nu->flagv & CU_NURB_CYCLIC) {
828     v = mod_i(v, totv);
829   }
830   else if (v < 0 || v >= totv) {
831     return -1;
832   }
833
834   return (v * totu) + u;
835 }
836
837 void BKE_nurb_index_to_uv(Nurb *nu, int index, int *r_u, int *r_v)
838 {
839   const int totu = nu->pntsu;
840   const int totv = nu->pntsv;
841   BLI_assert(index >= 0 && index < (nu->pntsu * nu->pntsv));
842   *r_u = (index % totu);
843   *r_v = (index / totu) % totv;
844 }
845
846 BezTriple *BKE_nurb_bezt_get_next(Nurb *nu, BezTriple *bezt)
847 {
848   BezTriple *bezt_next;
849
850   BLI_assert(ARRAY_HAS_ITEM(bezt, nu->bezt, nu->pntsu));
851
852   if (bezt == &nu->bezt[nu->pntsu - 1]) {
853     if (nu->flagu & CU_NURB_CYCLIC) {
854       bezt_next = nu->bezt;
855     }
856     else {
857       bezt_next = NULL;
858     }
859   }
860   else {
861     bezt_next = bezt + 1;
862   }
863
864   return bezt_next;
865 }
866
867 BPoint *BKE_nurb_bpoint_get_next(Nurb *nu, BPoint *bp)
868 {
869   BPoint *bp_next;
870
871   BLI_assert(ARRAY_HAS_ITEM(bp, nu->bp, nu->pntsu));
872
873   if (bp == &nu->bp[nu->pntsu - 1]) {
874     if (nu->flagu & CU_NURB_CYCLIC) {
875       bp_next = nu->bp;
876     }
877     else {
878       bp_next = NULL;
879     }
880   }
881   else {
882     bp_next = bp + 1;
883   }
884
885   return bp_next;
886 }
887
888 BezTriple *BKE_nurb_bezt_get_prev(Nurb *nu, BezTriple *bezt)
889 {
890   BezTriple *bezt_prev;
891
892   BLI_assert(ARRAY_HAS_ITEM(bezt, nu->bezt, nu->pntsu));
893   BLI_assert(nu->pntsv <= 1);
894
895   if (bezt == nu->bezt) {
896     if (nu->flagu & CU_NURB_CYCLIC) {
897       bezt_prev = &nu->bezt[nu->pntsu - 1];
898     }
899     else {
900       bezt_prev = NULL;
901     }
902   }
903   else {
904     bezt_prev = bezt - 1;
905   }
906
907   return bezt_prev;
908 }
909
910 BPoint *BKE_nurb_bpoint_get_prev(Nurb *nu, BPoint *bp)
911 {
912   BPoint *bp_prev;
913
914   BLI_assert(ARRAY_HAS_ITEM(bp, nu->bp, nu->pntsu));
915   BLI_assert(nu->pntsv == 1);
916
917   if (bp == nu->bp) {
918     if (nu->flagu & CU_NURB_CYCLIC) {
919       bp_prev = &nu->bp[nu->pntsu - 1];
920     }
921     else {
922       bp_prev = NULL;
923     }
924   }
925   else {
926     bp_prev = bp - 1;
927   }
928
929   return bp_prev;
930 }
931
932 void BKE_nurb_bezt_calc_normal(struct Nurb *UNUSED(nu), BezTriple *bezt, float r_normal[3])
933 {
934   /* calculate the axis matrix from the spline */
935   float dir_prev[3], dir_next[3];
936
937   sub_v3_v3v3(dir_prev, bezt->vec[0], bezt->vec[1]);
938   sub_v3_v3v3(dir_next, bezt->vec[1], bezt->vec[2]);
939
940   normalize_v3(dir_prev);
941   normalize_v3(dir_next);
942
943   add_v3_v3v3(r_normal, dir_prev, dir_next);
944   normalize_v3(r_normal);
945 }
946
947 void BKE_nurb_bezt_calc_plane(struct Nurb *nu, BezTriple *bezt, float r_plane[3])
948 {
949   float dir_prev[3], dir_next[3];
950
951   sub_v3_v3v3(dir_prev, bezt->vec[0], bezt->vec[1]);
952   sub_v3_v3v3(dir_next, bezt->vec[1], bezt->vec[2]);
953
954   normalize_v3(dir_prev);
955   normalize_v3(dir_next);
956
957   cross_v3_v3v3(r_plane, dir_prev, dir_next);
958   if (normalize_v3(r_plane) < FLT_EPSILON) {
959     BezTriple *bezt_prev = BKE_nurb_bezt_get_prev(nu, bezt);
960     BezTriple *bezt_next = BKE_nurb_bezt_get_next(nu, bezt);
961
962     if (bezt_prev) {
963       sub_v3_v3v3(dir_prev, bezt_prev->vec[1], bezt->vec[1]);
964       normalize_v3(dir_prev);
965     }
966     if (bezt_next) {
967       sub_v3_v3v3(dir_next, bezt->vec[1], bezt_next->vec[1]);
968       normalize_v3(dir_next);
969     }
970     cross_v3_v3v3(r_plane, dir_prev, dir_next);
971   }
972
973   /* matches with bones more closely */
974   {
975     float dir_mid[3], tvec[3];
976     add_v3_v3v3(dir_mid, dir_prev, dir_next);
977     cross_v3_v3v3(tvec, r_plane, dir_mid);
978     copy_v3_v3(r_plane, tvec);
979   }
980
981   normalize_v3(r_plane);
982 }
983
984 void BKE_nurb_bpoint_calc_normal(struct Nurb *nu, BPoint *bp, float r_normal[3])
985 {
986   BPoint *bp_prev = BKE_nurb_bpoint_get_prev(nu, bp);
987   BPoint *bp_next = BKE_nurb_bpoint_get_next(nu, bp);
988
989   zero_v3(r_normal);
990
991   if (bp_prev) {
992     float dir_prev[3];
993     sub_v3_v3v3(dir_prev, bp_prev->vec, bp->vec);
994     normalize_v3(dir_prev);
995     add_v3_v3(r_normal, dir_prev);
996   }
997   if (bp_next) {
998     float dir_next[3];
999     sub_v3_v3v3(dir_next, bp->vec, bp_next->vec);
1000     normalize_v3(dir_next);
1001     add_v3_v3(r_normal, dir_next);
1002   }
1003
1004   normalize_v3(r_normal);
1005 }
1006
1007 void BKE_nurb_bpoint_calc_plane(struct Nurb *nu, BPoint *bp, float r_plane[3])
1008 {
1009   BPoint *bp_prev = BKE_nurb_bpoint_get_prev(nu, bp);
1010   BPoint *bp_next = BKE_nurb_bpoint_get_next(nu, bp);
1011
1012   float dir_prev[3] = {0.0f}, dir_next[3] = {0.0f};
1013
1014   if (bp_prev) {
1015     sub_v3_v3v3(dir_prev, bp_prev->vec, bp->vec);
1016     normalize_v3(dir_prev);
1017   }
1018   if (bp_next) {
1019     sub_v3_v3v3(dir_next, bp->vec, bp_next->vec);
1020     normalize_v3(dir_next);
1021   }
1022   cross_v3_v3v3(r_plane, dir_prev, dir_next);
1023
1024   /* matches with bones more closely */
1025   {
1026     float dir_mid[3], tvec[3];
1027     add_v3_v3v3(dir_mid, dir_prev, dir_next);
1028     cross_v3_v3v3(tvec, r_plane, dir_mid);
1029     copy_v3_v3(r_plane, tvec);
1030   }
1031
1032   normalize_v3(r_plane);
1033 }
1034
1035 /* ~~~~~~~~~~~~~~~~~~~~Non Uniform Rational B Spline calculations ~~~~~~~~~~~ */
1036
1037 static void calcknots(float *knots, const int pnts, const short order, const short flag)
1038 {
1039   /* knots: number of pnts NOT corrected for cyclic */
1040   const int pnts_order = pnts + order;
1041   float k;
1042   int a;
1043
1044   switch (flag & (CU_NURB_ENDPOINT | CU_NURB_BEZIER)) {
1045     case CU_NURB_ENDPOINT:
1046       k = 0.0;
1047       for (a = 1; a <= pnts_order; a++) {
1048         knots[a - 1] = k;
1049         if (a >= order && a <= pnts) {
1050           k += 1.0f;
1051         }
1052       }
1053       break;
1054     case CU_NURB_BEZIER:
1055       /* Warning, the order MUST be 2 or 4,
1056        * if this is not enforced, the displist will be corrupt */
1057       if (order == 4) {
1058         k = 0.34;
1059         for (a = 0; a < pnts_order; a++) {
1060           knots[a] = floorf(k);
1061           k += (1.0f / 3.0f);
1062         }
1063       }
1064       else if (order == 3) {
1065         k = 0.6f;
1066         for (a = 0; a < pnts_order; a++) {
1067           if (a >= order && a <= pnts) {
1068             k += 0.5f;
1069           }
1070           knots[a] = floorf(k);
1071         }
1072       }
1073       else {
1074         CLOG_ERROR(&LOG, "bez nurb curve order is not 3 or 4, should never happen");
1075       }
1076       break;
1077     default:
1078       for (a = 0; a < pnts_order; a++) {
1079         knots[a] = (float)a;
1080       }
1081       break;
1082   }
1083 }
1084
1085 static void makecyclicknots(float *knots, int pnts, short order)
1086 /* pnts, order: number of pnts NOT corrected for cyclic */
1087 {
1088   int a, b, order2, c;
1089
1090   if (knots == NULL) {
1091     return;
1092   }
1093
1094   order2 = order - 1;
1095
1096   /* do first long rows (order -1), remove identical knots at endpoints */
1097   if (order > 2) {
1098     b = pnts + order2;
1099     for (a = 1; a < order2; a++) {
1100       if (knots[b] != knots[b - a]) {
1101         break;
1102       }
1103     }
1104     if (a == order2) {
1105       knots[pnts + order - 2] += 1.0f;
1106     }
1107   }
1108
1109   b = order;
1110   c = pnts + order + order2;
1111   for (a = pnts + order2; a < c; a++) {
1112     knots[a] = knots[a - 1] + (knots[b] - knots[b - 1]);
1113     b--;
1114   }
1115 }
1116
1117 static void makeknots(Nurb *nu, short uv)
1118 {
1119   if (nu->type == CU_NURBS) {
1120     if (uv == 1) {
1121       if (nu->knotsu) {
1122         MEM_freeN(nu->knotsu);
1123       }
1124       if (BKE_nurb_check_valid_u(nu)) {
1125         nu->knotsu = MEM_calloc_arrayN(KNOTSU(nu) + 1, sizeof(float), "makeknots");
1126         if (nu->flagu & CU_NURB_CYCLIC) {
1127           calcknots(nu->knotsu, nu->pntsu, nu->orderu, 0); /* cyclic should be uniform */
1128           makecyclicknots(nu->knotsu, nu->pntsu, nu->orderu);
1129         }
1130         else {
1131           calcknots(nu->knotsu, nu->pntsu, nu->orderu, nu->flagu);
1132         }
1133       }
1134       else {
1135         nu->knotsu = NULL;
1136       }
1137     }
1138     else if (uv == 2) {
1139       if (nu->knotsv) {
1140         MEM_freeN(nu->knotsv);
1141       }
1142       if (BKE_nurb_check_valid_v(nu)) {
1143         nu->knotsv = MEM_calloc_arrayN(KNOTSV(nu) + 1, sizeof(float), "makeknots");
1144         if (nu->flagv & CU_NURB_CYCLIC) {
1145           calcknots(nu->knotsv, nu->pntsv, nu->orderv, 0); /* cyclic should be uniform */
1146           makecyclicknots(nu->knotsv, nu->pntsv, nu->orderv);
1147         }
1148         else {
1149           calcknots(nu->knotsv, nu->pntsv, nu->orderv, nu->flagv);
1150         }
1151       }
1152       else {
1153         nu->knotsv = NULL;
1154       }
1155     }
1156   }
1157 }
1158
1159 void BKE_nurb_knot_calc_u(Nurb *nu)
1160 {
1161   makeknots(nu, 1);
1162 }
1163
1164 void BKE_nurb_knot_calc_v(Nurb *nu)
1165 {
1166   makeknots(nu, 2);
1167 }
1168
1169 static void basisNurb(
1170     float t, short order, int pnts, float *knots, float *basis, int *start, int *end)
1171 {
1172   float d, e;
1173   int i, i1 = 0, i2 = 0, j, orderpluspnts, opp2, o2;
1174
1175   orderpluspnts = order + pnts;
1176   opp2 = orderpluspnts - 1;
1177
1178   /* this is for float inaccuracy */
1179   if (t < knots[0]) {
1180     t = knots[0];
1181   }
1182   else if (t > knots[opp2]) {
1183     t = knots[opp2];
1184   }
1185
1186   /* this part is order '1' */
1187   o2 = order + 1;
1188   for (i = 0; i < opp2; i++) {
1189     if (knots[i] != knots[i + 1] && t >= knots[i] && t <= knots[i + 1]) {
1190       basis[i] = 1.0;
1191       i1 = i - o2;
1192       if (i1 < 0) {
1193         i1 = 0;
1194       }
1195       i2 = i;
1196       i++;
1197       while (i < opp2) {
1198         basis[i] = 0.0;
1199         i++;
1200       }
1201       break;
1202     }
1203     else {
1204       basis[i] = 0.0;
1205     }
1206   }
1207   basis[i] = 0.0;
1208
1209   /* this is order 2, 3, ... */
1210   for (j = 2; j <= order; j++) {
1211
1212     if (i2 + j >= orderpluspnts) {
1213       i2 = opp2 - j;
1214     }
1215
1216     for (i = i1; i <= i2; i++) {
1217       if (basis[i] != 0.0f) {
1218         d = ((t - knots[i]) * basis[i]) / (knots[i + j - 1] - knots[i]);
1219       }
1220       else {
1221         d = 0.0f;
1222       }
1223
1224       if (basis[i + 1] != 0.0f) {
1225         e = ((knots[i + j] - t) * basis[i + 1]) / (knots[i + j] - knots[i + 1]);
1226       }
1227       else {
1228         e = 0.0;
1229       }
1230
1231       basis[i] = d + e;
1232     }
1233   }
1234
1235   *start = 1000;
1236   *end = 0;
1237
1238   for (i = i1; i <= i2; i++) {
1239     if (basis[i] > 0.0f) {
1240       *end = i;
1241       if (*start == 1000) {
1242         *start = i;
1243       }
1244     }
1245   }
1246 }
1247
1248 /**
1249  * \param coord_array: has to be (3 * 4 * resolu * resolv) in size, and zero-ed.
1250  */
1251 void BKE_nurb_makeFaces(const Nurb *nu, float *coord_array, int rowstride, int resolu, int resolv)
1252 {
1253   BPoint *bp;
1254   float *basisu, *basis, *basisv, *sum, *fp, *in;
1255   float u, v, ustart, uend, ustep, vstart, vend, vstep, sumdiv;
1256   int i, j, iofs, jofs, cycl, len, curu, curv;
1257   int istart, iend, jsta, jen, *jstart, *jend, ratcomp;
1258
1259   int totu = nu->pntsu * resolu, totv = nu->pntsv * resolv;
1260
1261   if (nu->knotsu == NULL || nu->knotsv == NULL) {
1262     return;
1263   }
1264   if (nu->orderu > nu->pntsu) {
1265     return;
1266   }
1267   if (nu->orderv > nu->pntsv) {
1268     return;
1269   }
1270   if (coord_array == NULL) {
1271     return;
1272   }
1273
1274   /* allocate and initialize */
1275   len = totu * totv;
1276   if (len == 0) {
1277     return;
1278   }
1279
1280   sum = (float *)MEM_calloc_arrayN(len, sizeof(float), "makeNurbfaces1");
1281
1282   bp = nu->bp;
1283   i = nu->pntsu * nu->pntsv;
1284   ratcomp = 0;
1285   while (i--) {
1286     if (bp->vec[3] != 1.0f) {
1287       ratcomp = 1;
1288       break;
1289     }
1290     bp++;
1291   }
1292
1293   fp = nu->knotsu;
1294   ustart = fp[nu->orderu - 1];
1295   if (nu->flagu & CU_NURB_CYCLIC) {
1296     uend = fp[nu->pntsu + nu->orderu - 1];
1297   }
1298   else {
1299     uend = fp[nu->pntsu];
1300   }
1301   ustep = (uend - ustart) / ((nu->flagu & CU_NURB_CYCLIC) ? totu : totu - 1);
1302
1303   basisu = (float *)MEM_malloc_arrayN(KNOTSU(nu), sizeof(float), "makeNurbfaces3");
1304
1305   fp = nu->knotsv;
1306   vstart = fp[nu->orderv - 1];
1307
1308   if (nu->flagv & CU_NURB_CYCLIC) {
1309     vend = fp[nu->pntsv + nu->orderv - 1];
1310   }
1311   else {
1312     vend = fp[nu->pntsv];
1313   }
1314   vstep = (vend - vstart) / ((nu->flagv & CU_NURB_CYCLIC) ? totv : totv - 1);
1315
1316   len = KNOTSV(nu);
1317   basisv = (float *)MEM_malloc_arrayN(len * totv, sizeof(float), "makeNurbfaces3");
1318   jstart = (int *)MEM_malloc_arrayN(totv, sizeof(float), "makeNurbfaces4");
1319   jend = (int *)MEM_malloc_arrayN(totv, sizeof(float), "makeNurbfaces5");
1320
1321   /* precalculation of basisv and jstart, jend */
1322   if (nu->flagv & CU_NURB_CYCLIC) {
1323     cycl = nu->orderv - 1;
1324   }
1325   else {
1326     cycl = 0;
1327   }
1328   v = vstart;
1329   basis = basisv;
1330   curv = totv;
1331   while (curv--) {
1332     basisNurb(v, nu->orderv, nu->pntsv + cycl, nu->knotsv, basis, jstart + curv, jend + curv);
1333     basis += KNOTSV(nu);
1334     v += vstep;
1335   }
1336
1337   if (nu->flagu & CU_NURB_CYCLIC) {
1338     cycl = nu->orderu - 1;
1339   }
1340   else {
1341     cycl = 0;
1342   }
1343   in = coord_array;
1344   u = ustart;
1345   curu = totu;
1346   while (curu--) {
1347     basisNurb(u, nu->orderu, nu->pntsu + cycl, nu->knotsu, basisu, &istart, &iend);
1348
1349     basis = basisv;
1350     curv = totv;
1351     while (curv--) {
1352       jsta = jstart[curv];
1353       jen = jend[curv];
1354
1355       /* calculate sum */
1356       sumdiv = 0.0;
1357       fp = sum;
1358
1359       for (j = jsta; j <= jen; j++) {
1360
1361         if (j >= nu->pntsv) {
1362           jofs = (j - nu->pntsv);
1363         }
1364         else {
1365           jofs = j;
1366         }
1367         bp = nu->bp + nu->pntsu * jofs + istart - 1;
1368
1369         for (i = istart; i <= iend; i++, fp++) {
1370           if (i >= nu->pntsu) {
1371             iofs = i - nu->pntsu;
1372             bp = nu->bp + nu->pntsu * jofs + iofs;
1373           }
1374           else {
1375             bp++;
1376           }
1377
1378           if (ratcomp) {
1379             *fp = basisu[i] * basis[j] * bp->vec[3];
1380             sumdiv += *fp;
1381           }
1382           else {
1383             *fp = basisu[i] * basis[j];
1384           }
1385         }
1386       }
1387
1388       if (ratcomp) {
1389         fp = sum;
1390         for (j = jsta; j <= jen; j++) {
1391           for (i = istart; i <= iend; i++, fp++) {
1392             *fp /= sumdiv;
1393           }
1394         }
1395       }
1396
1397       zero_v3(in);
1398
1399       /* one! (1.0) real point now */
1400       fp = sum;
1401       for (j = jsta; j <= jen; j++) {
1402
1403         if (j >= nu->pntsv) {
1404           jofs = (j - nu->pntsv);
1405         }
1406         else {
1407           jofs = j;
1408         }
1409         bp = nu->bp + nu->pntsu * jofs + istart - 1;
1410
1411         for (i = istart; i <= iend; i++, fp++) {
1412           if (i >= nu->pntsu) {
1413             iofs = i - nu->pntsu;
1414             bp = nu->bp + nu->pntsu * jofs + iofs;
1415           }
1416           else {
1417             bp++;
1418           }
1419
1420           if (*fp != 0.0f) {
1421             madd_v3_v3fl(in, bp->vec, *fp);
1422           }
1423         }
1424       }
1425
1426       in += 3;
1427       basis += KNOTSV(nu);
1428     }
1429     u += ustep;
1430     if (rowstride != 0) {
1431       in = (float *)(((unsigned char *)in) + (rowstride - 3 * totv * sizeof(*in)));
1432     }
1433   }
1434
1435   /* free */
1436   MEM_freeN(sum);
1437   MEM_freeN(basisu);
1438   MEM_freeN(basisv);
1439   MEM_freeN(jstart);
1440   MEM_freeN(jend);
1441 }
1442
1443 /**
1444  * \param coord_array: Has to be 3 * 4 * pntsu * resolu in size and zero-ed
1445  * \param tilt_array: set when non-NULL
1446  * \param radius_array: set when non-NULL
1447  */
1448 void BKE_nurb_makeCurve(const Nurb *nu,
1449                         float *coord_array,
1450                         float *tilt_array,
1451                         float *radius_array,
1452                         float *weight_array,
1453                         int resolu,
1454                         int stride)
1455 {
1456   const float eps = 1e-6f;
1457   BPoint *bp;
1458   float u, ustart, uend, ustep, sumdiv;
1459   float *basisu, *sum, *fp;
1460   float *coord_fp = coord_array, *tilt_fp = tilt_array, *radius_fp = radius_array,
1461         *weight_fp = weight_array;
1462   int i, len, istart, iend, cycl;
1463
1464   if (nu->knotsu == NULL) {
1465     return;
1466   }
1467   if (nu->orderu > nu->pntsu) {
1468     return;
1469   }
1470   if (coord_array == NULL) {
1471     return;
1472   }
1473
1474   /* allocate and initialize */
1475   len = nu->pntsu;
1476   if (len == 0) {
1477     return;
1478   }
1479   sum = (float *)MEM_calloc_arrayN(len, sizeof(float), "makeNurbcurve1");
1480
1481   resolu = (resolu * SEGMENTSU(nu));
1482
1483   if (resolu == 0) {
1484     MEM_freeN(sum);
1485     return;
1486   }
1487
1488   fp = nu->knotsu;
1489   ustart = fp[nu->orderu - 1];
1490   if (nu->flagu & CU_NURB_CYCLIC) {
1491     uend = fp[nu->pntsu + nu->orderu - 1];
1492   }
1493   else {
1494     uend = fp[nu->pntsu];
1495   }
1496   ustep = (uend - ustart) / (resolu - ((nu->flagu & CU_NURB_CYCLIC) ? 0 : 1));
1497
1498   basisu = (float *)MEM_malloc_arrayN(KNOTSU(nu), sizeof(float), "makeNurbcurve3");
1499
1500   if (nu->flagu & CU_NURB_CYCLIC) {
1501     cycl = nu->orderu - 1;
1502   }
1503   else {
1504     cycl = 0;
1505   }
1506
1507   u = ustart;
1508   while (resolu--) {
1509     basisNurb(u, nu->orderu, nu->pntsu + cycl, nu->knotsu, basisu, &istart, &iend);
1510
1511     /* calc sum */
1512     sumdiv = 0.0;
1513     fp = sum;
1514     bp = nu->bp + istart - 1;
1515     for (i = istart; i <= iend; i++, fp++) {
1516       if (i >= nu->pntsu) {
1517         bp = nu->bp + (i - nu->pntsu);
1518       }
1519       else {
1520         bp++;
1521       }
1522
1523       *fp = basisu[i] * bp->vec[3];
1524       sumdiv += *fp;
1525     }
1526     if ((sumdiv != 0.0f) && (sumdiv < 1.0f - eps || sumdiv > 1.0f + eps)) {
1527       /* is normalizing needed? */
1528       fp = sum;
1529       for (i = istart; i <= iend; i++, fp++) {
1530         *fp /= sumdiv;
1531       }
1532     }
1533
1534     zero_v3(coord_fp);
1535
1536     /* one! (1.0) real point */
1537     fp = sum;
1538     bp = nu->bp + istart - 1;
1539     for (i = istart; i <= iend; i++, fp++) {
1540       if (i >= nu->pntsu) {
1541         bp = nu->bp + (i - nu->pntsu);
1542       }
1543       else {
1544         bp++;
1545       }
1546
1547       if (*fp != 0.0f) {
1548         madd_v3_v3fl(coord_fp, bp->vec, *fp);
1549
1550         if (tilt_fp) {
1551           (*tilt_fp) += (*fp) * bp->tilt;
1552         }
1553
1554         if (radius_fp) {
1555           (*radius_fp) += (*fp) * bp->radius;
1556         }
1557
1558         if (weight_fp) {
1559           (*weight_fp) += (*fp) * bp->weight;
1560         }
1561       }
1562     }
1563
1564     coord_fp = POINTER_OFFSET(coord_fp, stride);
1565
1566     if (tilt_fp) {
1567       tilt_fp = POINTER_OFFSET(tilt_fp, stride);
1568     }
1569     if (radius_fp) {
1570       radius_fp = POINTER_OFFSET(radius_fp, stride);
1571     }
1572     if (weight_fp) {
1573       weight_fp = POINTER_OFFSET(weight_fp, stride);
1574     }
1575
1576     u += ustep;
1577   }
1578
1579   /* free */
1580   MEM_freeN(sum);
1581   MEM_freeN(basisu);
1582 }
1583
1584 /**
1585  * Calculate the length for arrays filled in by #BKE_curve_calc_coords_axis.
1586  */
1587 unsigned int BKE_curve_calc_coords_axis_len(const unsigned int bezt_array_len,
1588                                             const unsigned int resolu,
1589                                             const bool is_cyclic,
1590                                             const bool use_cyclic_duplicate_endpoint)
1591 {
1592   const unsigned int segments = bezt_array_len - (is_cyclic ? 0 : 1);
1593   const unsigned int points_len = (segments * resolu) +
1594                                   (is_cyclic ? (use_cyclic_duplicate_endpoint) : 1);
1595   return points_len;
1596 }
1597
1598 /**
1599  * Calculate an array for the entire curve (cyclic or non-cyclic).
1600  * \note Call for each axis.
1601  *
1602  * \param use_cyclic_duplicate_endpoint: Duplicate values at the beginning & end of the array.
1603  */
1604 void BKE_curve_calc_coords_axis(const BezTriple *bezt_array,
1605                                 const unsigned int bezt_array_len,
1606                                 const unsigned int resolu,
1607                                 const bool is_cyclic,
1608                                 const bool use_cyclic_duplicate_endpoint,
1609                                 /* array params */
1610                                 const unsigned int axis,
1611                                 const unsigned int stride,
1612                                 float *r_points)
1613 {
1614   const unsigned int points_len = BKE_curve_calc_coords_axis_len(
1615       bezt_array_len, resolu, is_cyclic, use_cyclic_duplicate_endpoint);
1616   float *r_points_offset = r_points;
1617
1618   const unsigned int resolu_stride = resolu * stride;
1619   const unsigned int bezt_array_last = bezt_array_len - 1;
1620
1621   for (unsigned int i = 0; i < bezt_array_last; i++) {
1622     const BezTriple *bezt_curr = &bezt_array[i];
1623     const BezTriple *bezt_next = &bezt_array[i + 1];
1624     BKE_curve_forward_diff_bezier(bezt_curr->vec[1][axis],
1625                                   bezt_curr->vec[2][axis],
1626                                   bezt_next->vec[0][axis],
1627                                   bezt_next->vec[1][axis],
1628                                   r_points_offset,
1629                                   (int)resolu,
1630                                   stride);
1631     r_points_offset = POINTER_OFFSET(r_points_offset, resolu_stride);
1632   }
1633
1634   if (is_cyclic) {
1635     const BezTriple *bezt_curr = &bezt_array[bezt_array_last];
1636     const BezTriple *bezt_next = &bezt_array[0];
1637     BKE_curve_forward_diff_bezier(bezt_curr->vec[1][axis],
1638                                   bezt_curr->vec[2][axis],
1639                                   bezt_next->vec[0][axis],
1640                                   bezt_next->vec[1][axis],
1641                                   r_points_offset,
1642                                   (int)resolu,
1643                                   stride);
1644     r_points_offset = POINTER_OFFSET(r_points_offset, resolu_stride);
1645     if (use_cyclic_duplicate_endpoint) {
1646       *r_points_offset = *r_points;
1647       r_points_offset = POINTER_OFFSET(r_points_offset, stride);
1648     }
1649   }
1650   else {
1651     float *r_points_last = POINTER_OFFSET(r_points, bezt_array_last * resolu_stride);
1652     *r_points_last = bezt_array[bezt_array_last].vec[1][axis];
1653     r_points_offset = POINTER_OFFSET(r_points_offset, stride);
1654   }
1655
1656   BLI_assert(POINTER_OFFSET(r_points, points_len * stride) == r_points_offset);
1657   UNUSED_VARS_NDEBUG(points_len);
1658 }
1659
1660 /* forward differencing method for bezier curve */
1661 void BKE_curve_forward_diff_bezier(
1662     float q0, float q1, float q2, float q3, float *p, int it, int stride)
1663 {
1664   float rt0, rt1, rt2, rt3, f;
1665   int a;
1666
1667   f = (float)it;
1668   rt0 = q0;
1669   rt1 = 3.0f * (q1 - q0) / f;
1670   f *= f;
1671   rt2 = 3.0f * (q0 - 2.0f * q1 + q2) / f;
1672   f *= it;
1673   rt3 = (q3 - q0 + 3.0f * (q1 - q2)) / f;
1674
1675   q0 = rt0;
1676   q1 = rt1 + rt2 + rt3;
1677   q2 = 2 * rt2 + 6 * rt3;
1678   q3 = 6 * rt3;
1679
1680   for (a = 0; a <= it; a++) {
1681     *p = q0;
1682     p = POINTER_OFFSET(p, stride);
1683     q0 += q1;
1684     q1 += q2;
1685     q2 += q3;
1686   }
1687 }
1688
1689 /* forward differencing method for first derivative of cubic bezier curve */
1690 void BKE_curve_forward_diff_tangent_bezier(
1691     float q0, float q1, float q2, float q3, float *p, int it, int stride)
1692 {
1693   float rt0, rt1, rt2, f;
1694   int a;
1695
1696   f = 1.0f / (float)it;
1697
1698   rt0 = 3.0f * (q1 - q0);
1699   rt1 = f * (3.0f * (q3 - q0) + 9.0f * (q1 - q2));
1700   rt2 = 6.0f * (q0 + q2) - 12.0f * q1;
1701
1702   q0 = rt0;
1703   q1 = f * (rt1 + rt2);
1704   q2 = 2.0f * f * rt1;
1705
1706   for (a = 0; a <= it; a++) {
1707     *p = q0;
1708     p = POINTER_OFFSET(p, stride);
1709     q0 += q1;
1710     q1 += q2;
1711   }
1712 }
1713
1714 static void forward_diff_bezier_cotangent(const float p0[3],
1715                                           const float p1[3],
1716                                           const float p2[3],
1717                                           const float p3[3],
1718                                           float p[3],
1719                                           int it,
1720                                           int stride)
1721 {
1722   /* note that these are not perpendicular to the curve
1723    * they need to be rotated for this,
1724    *
1725    * This could also be optimized like BKE_curve_forward_diff_bezier */
1726   int a;
1727   for (a = 0; a <= it; a++) {
1728     float t = (float)a / (float)it;
1729
1730     int i;
1731     for (i = 0; i < 3; i++) {
1732       p[i] = (-6.0f * t + 6.0f) * p0[i] + (18.0f * t - 12.0f) * p1[i] +
1733              (-18.0f * t + 6.0f) * p2[i] + (6.0f * t) * p3[i];
1734     }
1735     normalize_v3(p);
1736     p = POINTER_OFFSET(p, stride);
1737   }
1738 }
1739
1740 /* ***************** BEVEL ****************** */
1741
1742 void BKE_curve_bevel_make(Depsgraph *depsgraph,
1743                           Scene *scene,
1744                           Object *ob,
1745                           ListBase *disp,
1746                           const bool for_render,
1747                           LinkNode *ob_cyclic_list)
1748 {
1749   DispList *dl, *dlnew;
1750   Curve *bevcu, *cu;
1751   float *fp, facx, facy, angle, dangle;
1752   int nr, a;
1753
1754   cu = ob->data;
1755   BLI_listbase_clear(disp);
1756
1757   /* if a font object is being edited, then do nothing */
1758   // XXX  if ( ob == obedit && ob->type == OB_FONT ) return;
1759
1760   if (cu->bevobj) {
1761     if (cu->bevobj->type != OB_CURVE) {
1762       return;
1763     }
1764
1765     bevcu = cu->bevobj->data;
1766     if (bevcu->ext1 == 0.0f && bevcu->ext2 == 0.0f) {
1767       ListBase bevdisp = {NULL, NULL};
1768       facx = cu->bevobj->scale[0];
1769       facy = cu->bevobj->scale[1];
1770
1771       if (for_render) {
1772         if (BLI_linklist_index(ob_cyclic_list, cu->bevobj) == -1) {
1773           BKE_displist_make_curveTypes_forRender(depsgraph,
1774                                                  scene,
1775                                                  cu->bevobj,
1776                                                  &bevdisp,
1777                                                  NULL,
1778                                                  false,
1779                                                  &(LinkNode){
1780                                                      .link = ob,
1781                                                      .next = ob_cyclic_list,
1782                                                  });
1783           dl = bevdisp.first;
1784         }
1785         else {
1786           dl = NULL;
1787         }
1788       }
1789       else if (cu->bevobj->runtime.curve_cache) {
1790         dl = cu->bevobj->runtime.curve_cache->disp.first;
1791       }
1792       else {
1793         BLI_assert(cu->bevobj->runtime.curve_cache != NULL);
1794         dl = NULL;
1795       }
1796
1797       while (dl) {
1798         if (ELEM(dl->type, DL_POLY, DL_SEGM)) {
1799           dlnew = MEM_mallocN(sizeof(DispList), "makebevelcurve1");
1800           *dlnew = *dl;
1801           dlnew->verts = MEM_malloc_arrayN(
1802               dl->parts * dl->nr, 3 * sizeof(float), "makebevelcurve1");
1803           memcpy(dlnew->verts, dl->verts, 3 * sizeof(float) * dl->parts * dl->nr);
1804
1805           if (dlnew->type == DL_SEGM) {
1806             dlnew->flag |= (DL_FRONT_CURVE | DL_BACK_CURVE);
1807           }
1808
1809           BLI_addtail(disp, dlnew);
1810           fp = dlnew->verts;
1811           nr = dlnew->parts * dlnew->nr;
1812           while (nr--) {
1813             fp[2] = fp[1] * facy;
1814             fp[1] = -fp[0] * facx;
1815             fp[0] = 0.0;
1816             fp += 3;
1817           }
1818         }
1819         dl = dl->next;
1820       }
1821
1822       BKE_displist_free(&bevdisp);
1823     }
1824   }
1825   else if (cu->ext1 == 0.0f && cu->ext2 == 0.0f) {
1826     /* pass */
1827   }
1828   else if (cu->ext2 == 0.0f) {
1829     dl = MEM_callocN(sizeof(DispList), "makebevelcurve2");
1830     dl->verts = MEM_malloc_arrayN(2, sizeof(float[3]), "makebevelcurve2");
1831     BLI_addtail(disp, dl);
1832     dl->type = DL_SEGM;
1833     dl->parts = 1;
1834     dl->flag = DL_FRONT_CURVE | DL_BACK_CURVE;
1835     dl->nr = 2;
1836
1837     fp = dl->verts;
1838     fp[0] = fp[1] = 0.0;
1839     fp[2] = -cu->ext1;
1840     fp[3] = fp[4] = 0.0;
1841     fp[5] = cu->ext1;
1842   }
1843   else if ((cu->flag & (CU_FRONT | CU_BACK)) == 0 &&
1844            cu->ext1 == 0.0f) { /* we make a full round bevel in that case */
1845     nr = 4 + 2 * cu->bevresol;
1846
1847     dl = MEM_callocN(sizeof(DispList), "makebevelcurve p1");
1848     dl->verts = MEM_malloc_arrayN(nr, sizeof(float[3]), "makebevelcurve p1");
1849     BLI_addtail(disp, dl);
1850     dl->type = DL_POLY;
1851     dl->parts = 1;
1852     dl->flag = DL_BACK_CURVE;
1853     dl->nr = nr;
1854
1855     /* a circle */
1856     fp = dl->verts;
1857     dangle = (2.0f * (float)M_PI / (nr));
1858     angle = -(nr - 1) * dangle;
1859
1860     for (a = 0; a < nr; a++) {
1861       fp[0] = 0.0;
1862       fp[1] = (cosf(angle) * (cu->ext2));
1863       fp[2] = (sinf(angle) * (cu->ext2)) - cu->ext1;
1864       angle += dangle;
1865       fp += 3;
1866     }
1867   }
1868   else {
1869     short dnr;
1870
1871     /* bevel now in three parts, for proper vertex normals */
1872     /* part 1, back */
1873
1874     if ((cu->flag & CU_BACK) || !(cu->flag & CU_FRONT)) {
1875       dnr = nr = 2 + cu->bevresol;
1876       if ((cu->flag & (CU_FRONT | CU_BACK)) == 0) {
1877         nr = 3 + 2 * cu->bevresol;
1878       }
1879       dl = MEM_callocN(sizeof(DispList), "makebevelcurve p1");
1880       dl->verts = MEM_malloc_arrayN(nr, sizeof(float[3]), "makebevelcurve p1");
1881       BLI_addtail(disp, dl);
1882       dl->type = DL_SEGM;
1883       dl->parts = 1;
1884       dl->flag = DL_BACK_CURVE;
1885       dl->nr = nr;
1886
1887       /* half a circle */
1888       fp = dl->verts;
1889       dangle = ((float)M_PI_2 / (dnr - 1));
1890       angle = -(nr - 1) * dangle;
1891
1892       for (a = 0; a < nr; a++) {
1893         fp[0] = 0.0;
1894         fp[1] = (float)(cosf(angle) * (cu->ext2));
1895         fp[2] = (float)(sinf(angle) * (cu->ext2)) - cu->ext1;
1896         angle += dangle;
1897         fp += 3;
1898       }
1899     }
1900
1901     /* part 2, sidefaces */
1902     if (cu->ext1 != 0.0f) {
1903       nr = 2;
1904
1905       dl = MEM_callocN(sizeof(DispList), "makebevelcurve p2");
1906       dl->verts = MEM_malloc_arrayN(nr, sizeof(float[3]), "makebevelcurve p2");
1907       BLI_addtail(disp, dl);
1908       dl->type = DL_SEGM;
1909       dl->parts = 1;
1910       dl->nr = nr;
1911
1912       fp = dl->verts;
1913       fp[1] = cu->ext2;
1914       fp[2] = -cu->ext1;
1915       fp[4] = cu->ext2;
1916       fp[5] = cu->ext1;
1917
1918       if ((cu->flag & (CU_FRONT | CU_BACK)) == 0) {
1919         dl = MEM_dupallocN(dl);
1920         dl->verts = MEM_dupallocN(dl->verts);
1921         BLI_addtail(disp, dl);
1922
1923         fp = dl->verts;
1924         fp[1] = -fp[1];
1925         fp[2] = -fp[2];
1926         fp[4] = -fp[4];
1927         fp[5] = -fp[5];
1928       }
1929     }
1930
1931     /* part 3, front */
1932     if ((cu->flag & CU_FRONT) || !(cu->flag & CU_BACK)) {
1933       dnr = nr = 2 + cu->bevresol;
1934       if ((cu->flag & (CU_FRONT | CU_BACK)) == 0) {
1935         nr = 3 + 2 * cu->bevresol;
1936       }
1937       dl = MEM_callocN(sizeof(DispList), "makebevelcurve p3");
1938       dl->verts = MEM_malloc_arrayN(nr, sizeof(float[3]), "makebevelcurve p3");
1939       BLI_addtail(disp, dl);
1940       dl->type = DL_SEGM;
1941       dl->flag = DL_FRONT_CURVE;
1942       dl->parts = 1;
1943       dl->nr = nr;
1944
1945       /* half a circle */
1946       fp = dl->verts;
1947       angle = 0.0;
1948       dangle = ((float)M_PI_2 / (dnr - 1));
1949
1950       for (a = 0; a < nr; a++) {
1951         fp[0] = 0.0;
1952         fp[1] = (float)(cosf(angle) * (cu->ext2));
1953         fp[2] = (float)(sinf(angle) * (cu->ext2)) + cu->ext1;
1954         angle += dangle;
1955         fp += 3;
1956       }
1957     }
1958   }
1959 }
1960
1961 static int cu_isectLL(const float v1[3],
1962                       const float v2[3],
1963                       const float v3[3],
1964                       const float v4[3],
1965                       short cox,
1966                       short coy,
1967                       float *lambda,
1968                       float *mu,
1969                       float vec[3])
1970 {
1971   /* return:
1972    * -1: collinear
1973    *  0: no intersection of segments
1974    *  1: exact intersection of segments
1975    *  2: cross-intersection of segments
1976    */
1977   float deler;
1978
1979   deler = (v1[cox] - v2[cox]) * (v3[coy] - v4[coy]) - (v3[cox] - v4[cox]) * (v1[coy] - v2[coy]);
1980   if (deler == 0.0f) {
1981     return -1;
1982   }
1983
1984   *lambda = (v1[coy] - v3[coy]) * (v3[cox] - v4[cox]) - (v1[cox] - v3[cox]) * (v3[coy] - v4[coy]);
1985   *lambda = -(*lambda / deler);
1986
1987   deler = v3[coy] - v4[coy];
1988   if (deler == 0) {
1989     deler = v3[cox] - v4[cox];
1990     *mu = -(*lambda * (v2[cox] - v1[cox]) + v1[cox] - v3[cox]) / deler;
1991   }
1992   else {
1993     *mu = -(*lambda * (v2[coy] - v1[coy]) + v1[coy] - v3[coy]) / deler;
1994   }
1995   vec[cox] = *lambda * (v2[cox] - v1[cox]) + v1[cox];
1996   vec[coy] = *lambda * (v2[coy] - v1[coy]) + v1[coy];
1997
1998   if (*lambda >= 0.0f && *lambda <= 1.0f && *mu >= 0.0f && *mu <= 1.0f) {
1999     if (*lambda == 0.0f || *lambda == 1.0f || *mu == 0.0f || *mu == 1.0f) {
2000       return 1;
2001     }
2002     return 2;
2003   }
2004   return 0;
2005 }
2006
2007 static bool bevelinside(BevList *bl1, BevList *bl2)
2008 {
2009   /* is bl2 INSIDE bl1 ? with left-right method and "lambda's" */
2010   /* returns '1' if correct hole  */
2011   BevPoint *bevp, *prevbevp;
2012   float min, max, vec[3], hvec1[3], hvec2[3], lab, mu;
2013   int nr, links = 0, rechts = 0, mode;
2014
2015   /* take first vertex of possible hole */
2016
2017   bevp = bl2->bevpoints;
2018   hvec1[0] = bevp->vec[0];
2019   hvec1[1] = bevp->vec[1];
2020   hvec1[2] = 0.0;
2021   copy_v3_v3(hvec2, hvec1);
2022   hvec2[0] += 1000;
2023
2024   /* test it with all edges of potential surrounding poly */
2025   /* count number of transitions left-right  */
2026
2027   bevp = bl1->bevpoints;
2028   nr = bl1->nr;
2029   prevbevp = bevp + (nr - 1);
2030
2031   while (nr--) {
2032     min = prevbevp->vec[1];
2033     max = bevp->vec[1];
2034     if (max < min) {
2035       min = max;
2036       max = prevbevp->vec[1];
2037     }
2038     if (min != max) {
2039       if (min <= hvec1[1] && max >= hvec1[1]) {
2040         /* there's a transition, calc intersection point */
2041         mode = cu_isectLL(prevbevp->vec, bevp->vec, hvec1, hvec2, 0, 1, &lab, &mu, vec);
2042         /* if lab==0.0 or lab==1.0 then the edge intersects exactly a transition
2043          * only allow for one situation: we choose lab= 1.0
2044          */
2045         if (mode >= 0 && lab != 0.0f) {
2046           if (vec[0] < hvec1[0]) {
2047             links++;
2048           }
2049           else {
2050             rechts++;
2051           }
2052         }
2053       }
2054     }
2055     prevbevp = bevp;
2056     bevp++;
2057   }
2058
2059   return (links & 1) && (rechts & 1);
2060 }
2061
2062 struct BevelSort {
2063   BevList *bl;
2064   float left;
2065   int dir;
2066 };
2067
2068 static int vergxcobev(const void *a1, const void *a2)
2069 {
2070   const struct BevelSort *x1 = a1, *x2 = a2;
2071
2072   if (x1->left > x2->left) {
2073     return 1;
2074   }
2075   else if (x1->left < x2->left) {
2076     return -1;
2077   }
2078   return 0;
2079 }
2080
2081 /* this function cannot be replaced with atan2, but why? */
2082
2083 static void calc_bevel_sin_cos(
2084     float x1, float y1, float x2, float y2, float *r_sina, float *r_cosa)
2085 {
2086   float t01, t02, x3, y3;
2087
2088   t01 = sqrtf(x1 * x1 + y1 * y1);
2089   t02 = sqrtf(x2 * x2 + y2 * y2);
2090   if (t01 == 0.0f) {
2091     t01 = 1.0f;
2092   }
2093   if (t02 == 0.0f) {
2094     t02 = 1.0f;
2095   }
2096
2097   x1 /= t01;
2098   y1 /= t01;
2099   x2 /= t02;
2100   y2 /= t02;
2101
2102   t02 = x1 * x2 + y1 * y2;
2103   if (fabsf(t02) >= 1.0f) {
2104     t02 = M_PI_2;
2105   }
2106   else {
2107     t02 = (saacos(t02)) / 2.0f;
2108   }
2109
2110   t02 = sinf(t02);
2111   if (t02 == 0.0f) {
2112     t02 = 1.0f;
2113   }
2114
2115   x3 = x1 - x2;
2116   y3 = y1 - y2;
2117   if (x3 == 0 && y3 == 0) {
2118     x3 = y1;
2119     y3 = -x1;
2120   }
2121   else {
2122     t01 = sqrtf(x3 * x3 + y3 * y3);
2123     x3 /= t01;
2124     y3 /= t01;
2125   }
2126
2127   *r_sina = -y3 / t02;
2128   *r_cosa = x3 / t02;
2129 }
2130
2131 static void tilt_bezpart(BezTriple *prevbezt,
2132                          BezTriple *bezt,
2133                          Nurb *nu,
2134                          float *tilt_array,
2135                          float *radius_array,
2136                          float *weight_array,
2137                          int resolu,
2138                          int stride)
2139 {
2140   BezTriple *pprev, *next, *last;
2141   float fac, dfac, t[4];
2142   int a;
2143
2144   if (tilt_array == NULL && radius_array == NULL) {
2145     return;
2146   }
2147
2148   last = nu->bezt + (nu->pntsu - 1);
2149
2150   /* returns a point */
2151   if (prevbezt == nu->bezt) {
2152     if (nu->flagu & CU_NURB_CYCLIC) {
2153       pprev = last;
2154     }
2155     else {
2156       pprev = prevbezt;
2157     }
2158   }
2159   else {
2160     pprev = prevbezt - 1;
2161   }
2162
2163   /* next point */
2164   if (bezt == last) {
2165     if (nu->flagu & CU_NURB_CYCLIC) {
2166       next = nu->bezt;
2167     }
2168     else {
2169       next = bezt;
2170     }
2171   }
2172   else {
2173     next = bezt + 1;
2174   }
2175
2176   fac = 0.0;
2177   dfac = 1.0f / (float)resolu;
2178
2179   for (a = 0; a < resolu; a++, fac += dfac) {
2180     if (tilt_array) {
2181       if (nu->tilt_interp ==
2182           KEY_CU_EASE) { /* May as well support for tilt also 2.47 ease interp */
2183         *tilt_array = prevbezt->tilt +
2184                       (bezt->tilt - prevbezt->tilt) * (3.0f * fac * fac - 2.0f * fac * fac * fac);
2185       }
2186       else {
2187         key_curve_position_weights(fac, t, nu->tilt_interp);
2188         *tilt_array = t[0] * pprev->tilt + t[1] * prevbezt->tilt + t[2] * bezt->tilt +
2189                       t[3] * next->tilt;
2190       }
2191
2192       tilt_array = POINTER_OFFSET(tilt_array, stride);
2193     }
2194
2195     if (radius_array) {
2196       if (nu->radius_interp == KEY_CU_EASE) {
2197         /* Support 2.47 ease interp
2198          * Note! - this only takes the 2 points into account,
2199          * giving much more localized results to changes in radius, sometimes you want that */
2200         *radius_array = prevbezt->radius + (bezt->radius - prevbezt->radius) *
2201                                                (3.0f * fac * fac - 2.0f * fac * fac * fac);
2202       }
2203       else {
2204
2205         /* reuse interpolation from tilt if we can */
2206         if (tilt_array == NULL || nu->tilt_interp != nu->radius_interp) {
2207           key_curve_position_weights(fac, t, nu->radius_interp);
2208         }
2209         *radius_array = t[0] * pprev->radius + t[1] * prevbezt->radius + t[2] * bezt->radius +
2210                         t[3] * next->radius;
2211       }
2212
2213       radius_array = POINTER_OFFSET(radius_array, stride);
2214     }
2215
2216     if (weight_array) {
2217       /* basic interpolation for now, could copy tilt interp too  */
2218       *weight_array = prevbezt->weight + (bezt->weight - prevbezt->weight) *
2219                                              (3.0f * fac * fac - 2.0f * fac * fac * fac);
2220
2221       weight_array = POINTER_OFFSET(weight_array, stride);
2222     }
2223   }
2224 }
2225
2226 /* make_bevel_list_3D_* funcs, at a minimum these must
2227  * fill in the bezp->quat and bezp->dir values */
2228
2229 /* utility for make_bevel_list_3D_* funcs */
2230 static void bevel_list_calc_bisect(BevList *bl)
2231 {
2232   BevPoint *bevp2, *bevp1, *bevp0;
2233   int nr;
2234   bool is_cyclic = bl->poly != -1;
2235
2236   if (is_cyclic) {
2237     bevp2 = bl->bevpoints;
2238     bevp1 = bevp2 + (bl->nr - 1);
2239     bevp0 = bevp1 - 1;
2240     nr = bl->nr;
2241   }
2242   else {
2243     /* If spline is not cyclic, direction of first and
2244      * last bevel points matches direction of CV handle.
2245      *
2246      * This is getting calculated earlier when we know
2247      * CV's handles and here we might simply skip evaluation
2248      * of direction for this guys.
2249      */
2250
2251     bevp0 = bl->bevpoints;
2252     bevp1 = bevp0 + 1;
2253     bevp2 = bevp1 + 1;
2254
2255     nr = bl->nr - 2;
2256   }
2257
2258   while (nr--) {
2259     /* totally simple */
2260     bisect_v3_v3v3v3(bevp1->dir, bevp0->vec, bevp1->vec, bevp2->vec);
2261
2262     bevp0 = bevp1;
2263     bevp1 = bevp2;
2264     bevp2++;
2265   }
2266 }
2267 static void bevel_list_flip_tangents(BevList *bl)
2268 {
2269   BevPoint *bevp2, *bevp1, *bevp0;
2270   int nr;
2271
2272   bevp2 = bl->bevpoints;
2273   bevp1 = bevp2 + (bl->nr - 1);
2274   bevp0 = bevp1 - 1;
2275
2276   nr = bl->nr;
2277   while (nr--) {
2278     if (angle_normalized_v3v3(bevp0->tan, bevp1->tan) > DEG2RADF(90.0f)) {
2279       negate_v3(bevp1->tan);
2280     }
2281
2282     bevp0 = bevp1;
2283     bevp1 = bevp2;
2284     bevp2++;
2285   }
2286 }
2287 /* apply user tilt */
2288 static void bevel_list_apply_tilt(BevList *bl)
2289 {
2290   BevPoint *bevp2, *bevp1;
2291   int nr;
2292   float q[4];
2293
2294   bevp2 = bl->bevpoints;
2295   bevp1 = bevp2 + (bl->nr - 1);
2296
2297   nr = bl->nr;
2298   while (nr--) {
2299     axis_angle_to_quat(q, bevp1->dir, bevp1->tilt);
2300     mul_qt_qtqt(bevp1->quat, q, bevp1->quat);
2301     normalize_qt(bevp1->quat);
2302
2303     bevp1 = bevp2;
2304     bevp2++;
2305   }
2306 }
2307 /* smooth quats, this function should be optimized, it can get slow with many iterations. */
2308 static void bevel_list_smooth(BevList *bl, int smooth_iter)
2309 {
2310   BevPoint *bevp2, *bevp1, *bevp0;
2311   int nr;
2312
2313   float q[4];
2314   float bevp0_quat[4];
2315   int a;
2316
2317   for (a = 0; a < smooth_iter; a++) {
2318     bevp2 = bl->bevpoints;
2319     bevp1 = bevp2 + (bl->nr - 1);
2320     bevp0 = bevp1 - 1;
2321
2322     nr = bl->nr;
2323
2324     if (bl->poly == -1) { /* check its not cyclic */
2325       /* skip the first point */
2326       /* bevp0 = bevp1; */
2327       bevp1 = bevp2;
2328       bevp2++;
2329       nr--;
2330
2331       bevp0 = bevp1;
2332       bevp1 = bevp2;
2333       bevp2++;
2334       nr--;
2335     }
2336
2337     copy_qt_qt(bevp0_quat, bevp0->quat);
2338
2339     while (nr--) {
2340       /* interpolate quats */
2341       float zaxis[3] = {0, 0, 1}, cross[3], q2[4];
2342       interp_qt_qtqt(q, bevp0_quat, bevp2->quat, 0.5);
2343       normalize_qt(q);
2344
2345       mul_qt_v3(q, zaxis);
2346       cross_v3_v3v3(cross, zaxis, bevp1->dir);
2347       axis_angle_to_quat(q2, cross, angle_normalized_v3v3(zaxis, bevp1->dir));
2348       normalize_qt(q2);
2349
2350       copy_qt_qt(bevp0_quat, bevp1->quat);
2351       mul_qt_qtqt(q, q2, q);
2352       interp_qt_qtqt(bevp1->quat, bevp1->quat, q, 0.5);
2353       normalize_qt(bevp1->quat);
2354
2355       /* bevp0 = bevp1; */ /* UNUSED */
2356       bevp1 = bevp2;
2357       bevp2++;
2358     }
2359   }
2360 }
2361
2362 static void make_bevel_list_3D_zup(BevList *bl)
2363 {
2364   BevPoint *bevp = bl->bevpoints;
2365   int nr = bl->nr;
2366
2367   bevel_list_calc_bisect(bl);
2368
2369   while (nr--) {
2370     vec_to_quat(bevp->quat, bevp->dir, 5, 1);
2371     bevp++;
2372   }
2373 }
2374
2375 static void minimum_twist_between_two_points(BevPoint *current_point, BevPoint *previous_point)
2376 {
2377   float angle = angle_normalized_v3v3(previous_point->dir, current_point->dir);
2378   float q[4];
2379
2380   if (angle > 0.0f) { /* otherwise we can keep as is */
2381     float cross_tmp[3];
2382     cross_v3_v3v3(cross_tmp, previous_point->dir, current_point->dir);
2383     axis_angle_to_quat(q, cross_tmp, angle);
2384     mul_qt_qtqt(current_point->quat, q, previous_point->quat);
2385   }
2386   else {
2387     copy_qt_qt(current_point->quat, previous_point->quat);
2388   }
2389 }
2390
2391 static void make_bevel_list_3D_minimum_twist(BevList *bl)
2392 {
2393   BevPoint *bevp2, *bevp1, *bevp0; /* standard for all make_bevel_list_3D_* funcs */
2394   int nr;
2395   float q[4];
2396
2397   bevel_list_calc_bisect(bl);
2398
2399   bevp2 = bl->bevpoints;
2400   bevp1 = bevp2 + (bl->nr - 1);
2401   bevp0 = bevp1 - 1;
2402
2403   nr = bl->nr;
2404   while (nr--) {
2405
2406     if (nr + 4 > bl->nr) { /* first time and second time, otherwise first point adjusts last */
2407       vec_to_quat(bevp1->quat, bevp1->dir, 5, 1);
2408     }
2409     else {
2410       minimum_twist_between_two_points(bevp1, bevp0);
2411     }
2412
2413     bevp0 = bevp1;
2414     bevp1 = bevp2;
2415     bevp2++;
2416   }
2417
2418   if (bl->poly != -1) { /* check for cyclic */
2419
2420     /* Need to correct for the start/end points not matching
2421      * do this by calculating the tilt angle difference, then apply
2422      * the rotation gradually over the entire curve
2423      *
2424      * note that the split is between last and second last, rather than first/last as youd expect.
2425      *
2426      * real order is like this
2427      * 0,1,2,3,4 --> 1,2,3,4,0
2428      *
2429      * this is why we compare last with second last
2430      * */
2431     float vec_1[3] = {0, 1, 0}, vec_2[3] = {0, 1, 0}, angle, ang_fac, cross_tmp[3];
2432
2433     BevPoint *bevp_first;
2434     BevPoint *bevp_last;
2435
2436     bevp_first = bl->bevpoints;
2437     bevp_first += bl->nr - 1;
2438     bevp_last = bevp_first;
2439     bevp_last--;
2440
2441     /* quats and vec's are normalized, should not need to re-normalize */
2442     mul_qt_v3(bevp_first->quat, vec_1);
2443     mul_qt_v3(bevp_last->quat, vec_2);
2444     normalize_v3(vec_1);
2445     normalize_v3(vec_2);
2446
2447     /* align the vector, can avoid this and it looks 98% OK but
2448      * better to align the angle quat roll's before comparing */
2449     {
2450       cross_v3_v3v3(cross_tmp, bevp_last->dir, bevp_first->dir);
2451       angle = angle_normalized_v3v3(bevp_first->dir, bevp_last->dir);
2452       axis_angle_to_quat(q, cross_tmp, angle);
2453       mul_qt_v3(q, vec_2);
2454     }
2455
2456     angle = angle_normalized_v3v3(vec_1, vec_2);
2457
2458     /* flip rotation if needs be */
2459     cross_v3_v3v3(cross_tmp, vec_1, vec_2);
2460     normalize_v3(cross_tmp);
2461     if (angle_normalized_v3v3(bevp_first->dir, cross_tmp) < DEG2RADF(90.0f)) {
2462       angle = -angle;
2463     }
2464
2465     bevp2 = bl->bevpoints;
2466     bevp1 = bevp2 + (bl->nr - 1);
2467     bevp0 = bevp1 - 1;
2468
2469     nr = bl->nr;
2470     while (nr--) {
2471       ang_fac = angle * (1.0f - ((float)nr / bl->nr)); /* also works */
2472
2473       axis_angle_to_quat(q, bevp1->dir, ang_fac);
2474       mul_qt_qtqt(bevp1->quat, q, bevp1->quat);
2475
2476       bevp0 = bevp1;
2477       bevp1 = bevp2;
2478       bevp2++;
2479     }
2480   }
2481   else {
2482     /* Need to correct quat for the first/last point,
2483      * this is so because previously it was only calculated
2484      * using it's own direction, which might not correspond
2485      * the twist of neighbor point.
2486      */
2487     bevp1 = bl->bevpoints;
2488     bevp0 = bevp1 + 1;
2489     minimum_twist_between_two_points(bevp1, bevp0);
2490
2491     bevp2 = bl->bevpoints;
2492     bevp1 = bevp2 + (bl->nr - 1);
2493     bevp0 = bevp1 - 1;
2494     minimum_twist_between_two_points(bevp1, bevp0);
2495   }
2496 }
2497
2498 static void make_bevel_list_3D_tangent(BevList *bl)
2499 {
2500   BevPoint *bevp2, *bevp1, *bevp0; /* standard for all make_bevel_list_3D_* funcs */
2501   int nr;
2502
2503   float bevp0_tan[3];
2504
2505   bevel_list_calc_bisect(bl);
2506   bevel_list_flip_tangents(bl);
2507
2508   /* correct the tangents */
2509   bevp2 = bl->bevpoints;
2510   bevp1 = bevp2 + (bl->nr - 1);
2511   bevp0 = bevp1 - 1;
2512
2513   nr = bl->nr;
2514   while (nr--) {
2515     float cross_tmp[3];
2516     cross_v3_v3v3(cross_tmp, bevp1->tan, bevp1->dir);
2517     cross_v3_v3v3(bevp1->tan, cross_tmp, bevp1->dir);
2518     normalize_v3(bevp1->tan);
2519
2520     bevp0 = bevp1;
2521     bevp1 = bevp2;
2522     bevp2++;
2523   }
2524
2525   /* now for the real twist calc */
2526   bevp2 = bl->bevpoints;
2527   bevp1 = bevp2 + (bl->nr - 1);
2528   bevp0 = bevp1 - 1;
2529
2530   copy_v3_v3(bevp0_tan, bevp0->tan);
2531
2532   nr = bl->nr;
2533   while (nr--) {
2534     /* make perpendicular, modify tan in place, is ok */
2535     float cross_tmp[3];
2536     float zero[3] = {0, 0, 0};
2537
2538     cross_v3_v3v3(cross_tmp, bevp1->tan, bevp1->dir);
2539     normalize_v3(cross_tmp);
2540     tri_to_quat(bevp1->quat, zero, cross_tmp, bevp1->tan); /* XXX - could be faster */
2541
2542     /* bevp0 = bevp1; */ /* UNUSED */
2543     bevp1 = bevp2;
2544     bevp2++;
2545   }
2546 }
2547
2548 static void make_bevel_list_3D(BevList *bl, int smooth_iter, int twist_mode)
2549 {
2550   switch (twist_mode) {
2551     case CU_TWIST_TANGENT:
2552       make_bevel_list_3D_tangent(bl);
2553       break;
2554     case CU_TWIST_MINIMUM:
2555       make_bevel_list_3D_minimum_twist(bl);
2556       break;
2557     default: /* CU_TWIST_Z_UP default, pre 2.49c */
2558       make_bevel_list_3D_zup(bl);
2559       break;
2560   }
2561
2562   if (smooth_iter) {
2563     bevel_list_smooth(bl, smooth_iter);
2564   }
2565
2566   bevel_list_apply_tilt(bl);
2567 }
2568
2569 /* only for 2 points */
2570 static void make_bevel_list_segment_3D(BevList *bl)
2571 {
2572   float q[4];
2573
2574   BevPoint *bevp2 = bl->bevpoints;
2575   BevPoint *bevp1 = bevp2 + 1;
2576
2577   /* simple quat/dir */
2578   sub_v3_v3v3(bevp1->dir, bevp1->vec, bevp2->vec);
2579   normalize_v3(bevp1->dir);
2580
2581   vec_to_quat(bevp1->quat, bevp1->dir, 5, 1);
2582
2583   axis_angle_to_quat(q, bevp1->dir, bevp1->tilt);
2584   mul_qt_qtqt(bevp1->quat, q, bevp1->quat);
2585   normalize_qt(bevp1->quat);
2586   copy_v3_v3(bevp2->dir, bevp1->dir);
2587   copy_qt_qt(bevp2->quat, bevp1->quat);
2588 }
2589
2590 /* only for 2 points */
2591 static void make_bevel_list_segment_2D(BevList *bl)
2592 {
2593   BevPoint *bevp2 = bl->bevpoints;
2594   BevPoint *bevp1 = bevp2 + 1;
2595
2596   const float x1 = bevp1->vec[0] - bevp2->vec[0];
2597   const float y1 = bevp1->vec[1] - bevp2->vec[1];
2598
2599   calc_bevel_sin_cos(x1, y1, -x1, -y1, &(bevp1->sina), &(bevp1->cosa));
2600   bevp2->sina = bevp1->sina;
2601   bevp2->cosa = bevp1->cosa;
2602
2603   /* fill in dir & quat */
2604   make_bevel_list_segment_3D(bl);
2605 }
2606
2607 static void make_bevel_list_2D(BevList *bl)
2608 {
2609   /* note: bevp->dir and bevp->quat are not needed for beveling but are
2610    * used when making a path from a 2D curve, therefore they need to be set - Campbell */
2611
2612   BevPoint *bevp0, *bevp1, *bevp2;
2613   int nr;
2614
2615   if (bl->poly != -1) {
2616     bevp2 = bl->bevpoints;
2617     bevp1 = bevp2 + (bl->nr - 1);
2618     bevp0 = bevp1 - 1;
2619     nr = bl->nr;
2620   }
2621   else {
2622     bevp0 = bl->bevpoints;
2623     bevp1 = bevp0 + 1;
2624     bevp2 = bevp1 + 1;
2625
2626     nr = bl->nr - 2;
2627   }
2628
2629   while (nr--) {
2630     const float x1 = bevp1->vec[0] - bevp0->vec[0];
2631     const float x2 = bevp1->vec[0] - bevp2->vec[0];
2632     const float y1 = bevp1->vec[1] - bevp0->vec[1];
2633     const float y2 = bevp1->vec[1] - bevp2->vec[1];
2634
2635     calc_bevel_sin_cos(x1, y1, x2, y2, &(bevp1->sina), &(bevp1->cosa));
2636
2637     /* from: make_bevel_list_3D_zup, could call but avoid a second loop.
2638      * no need for tricky tilt calculation as with 3D curves */
2639     bisect_v3_v3v3v3(bevp1->dir, bevp0->vec, bevp1->vec, bevp2->vec);
2640     vec_to_quat(bevp1->quat, bevp1->dir, 5, 1);
2641     /* done with inline make_bevel_list_3D_zup */
2642
2643     bevp0 = bevp1;
2644     bevp1 = bevp2;
2645     bevp2++;
2646   }
2647
2648   /* correct non-cyclic cases */
2649   if (bl->poly == -1) {
2650     BevPoint *bevp;
2651     float angle;
2652
2653     /* first */
2654     bevp = bl->bevpoints;
2655     angle = atan2f(bevp->dir[0], bevp->dir[1]) - (float)M_PI_2;
2656     bevp->sina = sinf(angle);
2657     bevp->cosa = cosf(angle);
2658     vec_to_quat(bevp->quat, bevp->dir, 5, 1);
2659
2660     /* last */
2661     bevp = bl->bevpoints;
2662     bevp += (bl->nr - 1);
2663     angle = atan2f(bevp->dir[0], bevp->dir[1]) - (float)M_PI_2;
2664     bevp->sina = sinf(angle);
2665     bevp->cosa = cosf(angle);
2666     vec_to_quat(bevp->quat, bevp->dir, 5, 1);
2667   }
2668 }
2669
2670 static void bevlist_firstlast_direction_calc_from_bpoint(Nurb *nu, BevList *bl)
2671 {
2672   if (nu->pntsu > 1) {
2673     BPoint *first_bp = nu->bp, *last_bp = nu->bp + (nu->pntsu - 1);
2674     BevPoint *first_bevp, *last_bevp;
2675
2676     first_bevp = bl->bevpoints;
2677     last_bevp = first_bevp + (bl->nr - 1);
2678
2679     sub_v3_v3v3(first_bevp->dir, (first_bp + 1)->vec, first_bp->vec);
2680     normalize_v3(first_bevp->dir);
2681
2682     sub_v3_v3v3(last_bevp->dir, last_bp->vec, (last_bp - 1)->vec);
2683     normalize_v3(last_bevp->dir);
2684   }
2685 }
2686
2687 void BKE_curve_bevelList_free(ListBase *bev)
2688 {
2689   BevList *bl, *blnext;
2690   for (bl = bev->first; bl != NULL; bl = blnext) {
2691     blnext = bl->next;
2692     if (bl->seglen != NULL) {
2693       MEM_freeN(bl->seglen);
2694     }
2695     if (bl->segbevcount != NULL) {
2696       MEM_freeN(bl->segbevcount);
2697     }
2698     if (bl->bevpoints != NULL) {
2699       MEM_freeN(bl->bevpoints);
2700     }
2701     MEM_freeN(bl);
2702   }
2703
2704   BLI_listbase_clear(bev);
2705 }
2706
2707 void BKE_curve_bevelList_make(Object *ob, ListBase *nurbs, bool for_render)
2708 {
2709   /*
2710    * - convert all curves to polys, with indication of resol and flags for double-vertices
2711    * - possibly; do a smart vertice removal (in case Nurb)
2712    * - separate in individual blocks with BoundBox
2713    * - AutoHole detection
2714    */
2715
2716   /* this function needs an object, because of tflag and upflag */
2717   Curve *cu = ob->data;
2718   Nurb *nu;
2719   BezTriple *bezt, *prevbezt;
2720   BPoint *bp;
2721   BevList *bl, *blnew, *blnext;
2722   BevPoint *bevp2, *bevp1 = NULL, *bevp0;
2723   const float treshold = 0.00001f;
2724   float min, inp;
2725   float *seglen = NULL;
2726   struct BevelSort *sortdata, *sd, *sd1;
2727   int a, b, nr, poly, resolu = 0, len = 0, segcount;
2728   int *segbevcount;
2729   bool do_tilt, do_radius, do_weight;
2730   bool is_editmode = false;
2731   ListBase *bev;
2732
2733   /* segbevcount alsp requires seglen. */
2734   const bool need_seglen = ELEM(
2735                                cu->bevfac1_mapping, CU_BEVFAC_MAP_SEGMENT, CU_BEVFAC_MAP_SPLINE) ||
2736                            ELEM(cu->bevfac2_mapping, CU_BEVFAC_MAP_SEGMENT, CU_BEVFAC_MAP_SPLINE);
2737
2738   bev = &ob->runtime.curve_cache->bev;
2739
2740 #if 0
2741   /* do we need to calculate the radius for each point? */
2742   do_radius = (cu->bevobj || cu->taperobj || (cu->flag & CU_FRONT) || (cu->flag & CU_BACK)) ? 0 :
2743                                                                                               1;
2744 #endif
2745
2746   /* STEP 1: MAKE POLYS  */
2747
2748   BKE_curve_bevelList_free(&ob->runtime.curve_cache->bev);
2749   nu = nurbs->first;
2750   if (cu->editnurb && ob->type != OB_FONT) {
2751     is_editmode = 1;
2752   }
2753
2754   for (; nu; nu = nu->next) {
2755     if (nu->hide && is_editmode) {
2756       continue;
2757     }
2758
2759     /* check if we will calculate tilt data */
2760     do_tilt = CU_DO_TILT(cu, nu);
2761     do_radius = CU_DO_RADIUS(
2762         cu, nu); /* normal display uses the radius, better just to calculate them */
2763     do_weight = true;
2764
2765     /* check we are a single point? also check we are not a surface and that the orderu is sane,
2766      * enforced in the UI but can go wrong possibly */
2767     if (!BKE_nurb_check_valid_u(nu)) {
2768       bl = MEM_callocN(sizeof(BevList), "makeBevelList1");
2769       bl->bevpoints = MEM_calloc_arrayN(1, sizeof(BevPoint), "makeBevelPoints1");
2770       BLI_addtail(bev, bl);
2771       bl->nr = 0;
2772       bl->charidx = nu->charidx;
2773     }
2774     else {
2775       BevPoint *bevp;
2776
2777       if (for_render && cu->resolu_ren != 0) {
2778         resolu = cu->resolu_ren;
2779       }
2780       else {
2781         resolu = nu->resolu;
2782       }
2783
2784       segcount = SEGMENTSU(nu);
2785
2786       if (nu->type == CU_POLY) {
2787         len = nu->pntsu;
2788         bl = MEM_callocN(sizeof(BevList), "makeBevelList2");
2789         bl->bevpoints = MEM_calloc_arrayN(len, sizeof(BevPoint), "makeBevelPoints2");
2790         if (need_seglen && (nu->flagu & CU_NURB_CYCLIC) == 0) {
2791           bl->seglen = MEM_malloc_arrayN(segcount, sizeof(float), "makeBevelList2_seglen");
2792           bl->segbevcount = MEM_malloc_arrayN(segcount, sizeof(int), "makeBevelList2_segbevcount");
2793         }
2794         BLI_addtail(bev, bl);
2795
2796         bl->poly = (nu->flagu & CU_NURB_CYCLIC) ? 0 : -1;
2797         bl->nr = len;
2798         bl->dupe_nr = 0;
2799         bl->charidx = nu->charidx;
2800         bevp = bl->bevpoints;
2801         bevp->offset = 0;
2802         bp = nu->bp;
2803         seglen = bl->seglen;
2804         segbevcount = bl->segbevcount;
2805
2806         while (len--) {
2807           copy_v3_v3(bevp->vec, bp->vec);
2808           bevp->tilt = bp->tilt;
2809           bevp->radius = bp->radius;
2810           bevp->weight = bp->weight;
2811           bevp->split_tag = true;
2812           bp++;
2813           if (seglen != NULL && len != 0) {
2814             *seglen = len_v3v3(bevp->vec, bp->vec);
2815             bevp++;
2816             bevp->offset = *seglen;
2817             if (*seglen > treshold) {
2818               *segbevcount = 1;
2819             }
2820             else {
2821               *segbevcount = 0;
2822             }
2823             seglen++;
2824             segbevcount++;
2825           }
2826           else {
2827             bevp++;
2828           }
2829         }
2830
2831         if ((nu->flagu & CU_NURB_CYCLIC) == 0) {
2832           bevlist_firstlast_direction_calc_from_bpoint(nu, bl);
2833         }
2834       }
2835       else if (nu->type == CU_BEZIER) {
2836         /* in case last point is not cyclic */
2837         len = segcount * resolu + 1;
2838
2839         bl = MEM_callocN(sizeof(BevList), "makeBevelBPoints");
2840         bl->bevpoints = MEM_calloc_arrayN(len, sizeof(BevPoint), "makeBevelBPointsPoints");
2841         if (need_seglen && (nu->flagu & CU_NURB_CYCLIC) == 0) {
2842           bl->seglen = MEM_malloc_arrayN(segcount, sizeof(float), "makeBevelBPoints_seglen");
2843           bl->segbevcount = MEM_malloc_arrayN(
2844               segcount, sizeof(int), "makeBevelBPoints_segbevcount");
2845         }
2846         BLI_addtail(bev, bl);
2847
2848         bl->poly = (nu->flagu & CU_NURB_CYCLIC) ? 0 : -1;
2849         bl->charidx = nu->charidx;
2850
2851         bevp = bl->bevpoints;
2852         seglen = bl->seglen;
2853         segbevcount = bl->segbevcount;
2854
2855         bevp->offset = 0;
2856         if (seglen != NULL) {
2857           *seglen = 0;
2858           *segbevcount = 0;
2859         }
2860
2861         a = nu->pntsu - 1;
2862         bezt = nu->bezt;
2863         if (nu->flagu & CU_NURB_CYCLIC) {
2864           a++;
2865           prevbezt = nu->bezt + (nu->pntsu - 1);
2866         }
2867         else {
2868           prevbezt = bezt;
2869           bezt++;
2870         }
2871
2872         sub_v3_v3v3(bevp->dir, prevbezt->vec[2], prevbezt->vec[1]);
2873         normalize_v3(bevp->dir);
2874
2875         BLI_assert(segcount >= a);
2876
2877         while (a--) {
2878           if (prevbezt->h2 == HD_VECT && bezt->h1 == HD_VECT) {
2879
2880             copy_v3_v3(bevp->vec, prevbezt->vec[1]);
2881             bevp->tilt = prevbezt->tilt;
2882             bevp->radius = prevbezt->radius;
2883             bevp->weight = prevbezt->weight;
2884             bevp->split_tag = true;
2885             bevp->dupe_tag = false;
2886             bevp++;
2887             bl->nr++;
2888             bl->dupe_nr = 1;
2889             if (seglen != NULL) {
2890               *seglen = len_v3v3(prevbezt->vec[1], bezt->vec[1]);
2891               bevp->offset = *seglen;
2892               seglen++;
2893               /* match segbevcount to the cleaned up bevel lists (see STEP 2) */
2894               if (bevp->offset > treshold) {
2895                 *segbevcount = 1;
2896               }
2897               segbevcount++;
2898             }
2899           }
2900           else {
2901             /* always do all three, to prevent data hanging around */
2902             int j;
2903
2904             /* BevPoint must stay aligned to 4 so sizeof(BevPoint)/sizeof(float) works */
2905             for (j = 0; j < 3; j++) {
2906               BKE_curve_forward_diff_bezier(prevbezt->vec[1][j],
2907                                             prevbezt->vec[2][j],
2908                                             bezt->vec[0][j],
2909                                             bezt->vec[1][j],
2910                                             &(bevp->vec[j]),
2911                                             resolu,
2912                                             sizeof(BevPoint));
2913             }
2914
2915             /* if both arrays are NULL do nothiong */
2916             tilt_bezpart(prevbezt,
2917                          bezt,
2918                          nu,
2919                          do_tilt ? &bevp->tilt : NULL,
2920                          do_radius ? &bevp->radius : NULL,
2921                          do_weight ? &bevp->weight : NULL,
2922                          resolu,
2923                          sizeof(BevPoint));
2924
2925             if (cu->twist_mode == CU_TWIST_TANGENT) {
2926               forward_diff_bezier_cotangent(prevbezt->vec[1],
2927                                             prevbezt->vec[2],
2928                                             bezt->vec[0],
2929                                             bezt->vec[1],
2930                                             bevp->tan,
2931                                             resolu,
2932                                             sizeof(BevPoint));
2933             }
2934
2935             /* indicate with handlecodes double points */
2936             if (prevbezt->h1 == prevbezt->h2) {
2937               if (prevbezt->h1 == 0 || prevbezt->h1 == HD_VECT) {
2938                 bevp->split_tag = true;
2939               }
2940             }
2941             else {
2942               if (prevbezt->h1 == 0 || prevbezt->h1 == HD_VECT) {
2943                 bevp->split_tag = true;
2944               }
2945               else if (prevbezt->h2 == 0 || prevbezt->h2 == HD_VECT) {
2946                 bevp->split_tag = true;
2947               }
2948             }
2949
2950             /* seglen */
2951             if (seglen != NULL) {
2952               *seglen = 0;
2953               *segbevcount = 0;
2954               for (j = 0; j < resolu; j++) {
2955                 bevp0 = bevp;
2956                 bevp++;
2957                 bevp->offset = len_v3v3(bevp0->vec, bevp->vec);
2958                 /* match seglen and segbevcount to the cleaned up bevel lists (see STEP 2) */
2959                 if (bevp->offset > treshold) {
2960                   *seglen += bevp->offset;
2961                   *segbevcount += 1;
2962                 }
2963               }
2964               seglen++;
2965               segbevcount++;
2966             }
2967             else {
2968               bevp += resolu;
2969             }
2970             bl->nr += resolu;
2971           }
2972           prevbezt = bezt;
2973           bezt++;
2974         }
2975
2976         if ((nu->flagu & CU_NURB_CYCLIC) == 0) { /* not cyclic: endpoint */
2977           copy_v3_v3(bevp->vec, prevbezt->vec[1]);
2978           bevp->tilt = prevbezt->tilt;
2979           bevp->radius = prevbezt->radius;
2980           bevp->weight = prevbezt->weight;
2981
2982           sub_v3_v3v3(bevp->dir, prevbezt->vec[1], prevbezt->vec[0]);
2983           normalize_v3(bevp->dir);
2984
2985           bl->nr++;
2986         }
2987       }
2988       else if (nu->type == CU_NURBS) {
2989         if (nu->pntsv == 1) {
2990           len = (resolu * segcount);
2991
2992           bl = MEM_callocN(sizeof(BevList), "makeBevelList3");
2993           bl->bevpoints = MEM_calloc_arrayN(len, sizeof(BevPoint), "makeBevelPoints3");
2994           if (need_seglen && (nu->flagu & CU_NURB_CYCLIC) == 0) {
2995             bl->seglen = MEM_malloc_arrayN(segcount, sizeof(float), "makeBevelList3_seglen");
2996             bl->segbevcount = MEM_malloc_arrayN(
2997                 segcount, sizeof(int), "makeBevelList3_segbevcount");
2998           }
2999           BLI_addtail(bev, bl);
3000           bl->nr = len;
3001           bl->dupe_nr = 0;
3002           bl->poly = (nu->flagu & CU_NURB_CYCLIC) ? 0 : -1;
3003           bl->charidx = nu->charidx;
3004
3005           bevp = bl->bevpoints;
3006           seglen = bl->seglen;
3007           segbevcount = bl->segbevcount;
3008
3009           BKE_nurb_makeCurve(nu,
3010                              &bevp->vec[0],
3011                              do_tilt ? &bevp->tilt : NULL,
3012                              do_radius ? &bevp->radius : NULL,
3013                              do_weight ? &bevp->weight : NULL,
3014                              resolu,
3015                              sizeof(BevPoint));
3016
3017           /* match seglen and segbevcount to the cleaned up bevel lists (see STEP 2) */
3018           if (seglen != NULL) {
3019             nr = segcount;
3020             bevp0 = bevp;
3021             bevp++;
3022             while (nr) {
3023               int j;
3024               *seglen = 0;
3025               *segbevcount = 0;
3026               /* We keep last bevel segment zero-length. */
3027               for (j = 0; j < ((nr == 1) ? (resolu - 1) : resolu); j++) {
3028                 bevp->offset = len_v3v3(bevp0->vec, bevp->vec);
3029                 if (bevp->offset > treshold) {
3030                   *seglen += bevp->offset;
3031                   *segbevcount += 1;
3032                 }
3033                 bevp0 = bevp;
3034                 bevp++;
3035               }
3036               seglen++;
3037               segbevcount++;
3038               nr--;
3039             }
3040           }
3041
3042           if ((nu->flagu & CU_NURB_CYCLIC) == 0) {
3043             bevlist_firstlast_direction_calc_from_bpoint(nu, bl);
3044           }
3045         }
3046       }
3047     }
3048   }
3049
3050   /* STEP 2: DOUBLE POINTS AND AUTOMATIC RESOLUTION, REDUCE DATABLOCKS */
3051   bl = bev->first;
3052   while (bl) {
3053     if (bl->nr) { /* null bevel items come from single points */
3054       bool is_cyclic = bl->poly != -1;
3055       nr = bl->nr;
3056       if (is_cyclic) {
3057         bevp1 = bl->bevpoints;
3058         bevp0 = bevp1 + (nr - 1);
3059       }
3060       else {
3061         bevp0 = bl->bevpoints;
3062         bevp0->offset = 0;
3063         bevp1 = bevp0 + 1;
3064       }
3065       nr--;
3066       while (nr--) {
3067         if (seglen != NULL) {
3068           if (fabsf(bevp1->offset) < treshold) {
3069             bevp0->dupe_tag = true;
3070             bl->dupe_nr++;
3071           }
3072         }
3073         else {
3074           if (fabsf(bevp0->vec[0] - bevp1->vec[0]) < 0.00001f) {
3075             if (fabsf(bevp0->vec[1] - bevp1->vec[1]) < 0.00001f) {
3076               if (fabsf(bevp0->vec[2] - bevp1->vec[2]) < 0.00001f) {
3077                 bevp0->dupe_tag = true;
3078                 bl->dupe_nr++;
3079               }
3080             }
3081           }
3082         }
3083         bevp0 = bevp1;
3084         bevp1++;
3085       }
3086     }
3087     bl = bl->next;
3088   }
3089   bl = bev->first;
3090   while (bl) {
3091     blnext = bl->next;
3092     if (bl->nr && bl->dupe_nr) {
3093       nr = bl->nr - bl->dupe_nr + 1; /* +1 because vectorbezier sets flag too */
3094       blnew = MEM_callocN(sizeof(BevList), "makeBevelList4");
3095       memcpy(blnew, bl, sizeof(BevList));
3096       blnew->bevpoints = MEM_calloc_arrayN(nr, sizeof(BevPoint), "makeBevelPoints4");
3097       if (!blnew->bevpoints) {
3098         MEM_freeN(blnew);
3099         break;
3100       }
3101       blnew->segbevcount = bl->segbevcount;
3102       blnew->seglen = bl->seglen;
3103       blnew->nr = 0;
3104       BLI_remlink(bev, bl);
3105       BLI_insertlinkbefore(bev, blnext, blnew); /* to make sure bevlijst is tuned with nurblist */
3106       bevp0 = bl->bevpoints;
3107       bevp1 = blnew->bevpoints;
3108       nr = bl->nr;
3109       while (nr--) {
3110         if (bevp0->dupe_tag == 0) {
3111           memcpy(bevp1, bevp0, sizeof(BevPoint));
3112           bevp1++;
3113           blnew->nr++;
3114         }
3115         bevp0++;
3116       }
3117       if (bl->bevpoints != NULL) {
3118         MEM_freeN(bl->bevpoints);
3119       }
3120       MEM_freeN(bl);
3121       blnew->dupe_nr = 0;
3122     }
3123     bl = blnext;
3124   }
3125
3126   /* STEP 3: POLYS COUNT AND AUTOHOLE */
3127   bl = bev->first;
3128   poly = 0;
3129   while (bl) {
3130     if (bl->nr && bl->poly >= 0) {
3131       poly++;
3132       bl->poly = poly;
3133       bl->hole = 0;
3134     }
3135     bl = bl->next;
3136   }
3137
3138   /* find extreme left points, also test (turning) direction */
3139   if (poly > 0) {
3140     sd = sortdata = MEM_malloc_arrayN(poly, sizeof(struct BevelSort), "makeBevelList5");
3141     bl = bev->first;
3142     while (bl) {
3143       if (bl->poly > 0) {
3144         BevPoint *bevp;
3145
3146         min = 300000.0;
3147         bevp = bl->bevpoints;
3148         nr = bl->nr;
3149         while (nr--) {
3150           if (min > bevp->vec[0]) {
3151             min = bevp->vec[0];
3152             bevp1 = bevp;
3153           }
3154           bevp++;
3155         }
3156         sd->bl = bl;
3157         sd->left = min;
3158
3159         bevp = bl->bevpoints;
3160         if (bevp1 == bevp) {
3161           bevp0 = bevp + (bl->nr - 1);
3162         }
3163         else {
3164           bevp0 = bevp1 - 1;
3165         }
3166         bevp = bevp + (bl->nr - 1);
3167         if (bevp1 == bevp) {
3168           bevp2 = bl->bevpoints;
3169         }
3170         else {
3171           bevp2 = bevp1 + 1;
3172         }
3173
3174         inp = ((bevp1->vec[0] - bevp0->vec[0]) * (bevp0->vec[1] - bevp2->vec[1]) +
3175                (bevp0->vec[1] - bevp1->vec[1]) * (bevp0->vec[0] - bevp2->vec[0]));
3176
3177         if (inp > 0.0f) {
3178           sd->dir = 1;
3179         }
3180         else {
3181           sd->dir = 0;
3182         }
3183
3184         sd++;
3185       }
3186
3187       bl = bl->next;
3188     }
3189     qsort(sortdata, poly, sizeof(struct BevelSort), vergxcobev);
3190
3191     sd = sortdata + 1;
3192     for (a = 1; a < poly; a++, sd++) {
3193       bl = sd->bl; /* is bl a hole? */
3194       sd1 = sortdata + (a - 1);
3195       for (b = a - 1; b >= 0; b--, sd1--) {    /* all polys to the left */
3196         if (sd1->bl->charidx == bl->charidx) { /* for text, only check matching char */
3197           if (bevelinside(sd1->bl, bl)) {
3198             bl->hole = 1 - sd1->bl->hole;
3199             break;
3200           }
3201         }
3202       }
3203     }
3204
3205     /* turning direction */
3206     if ((cu->flag & CU_3D) == 0) {
3207       sd = sortdata;
3208       for (a = 0; a < poly; a++, sd++) {
3209         if (sd->bl->hole == sd->dir) {
3210           bl = sd->bl;
3211           bevp1 = bl->bevpoints;
3212           bevp2 = bevp1 + (bl->nr - 1);
3213           nr = bl->nr / 2;
3214           while (nr--) {
3215             SWAP(BevPoint, *bevp1, *bevp2);
3216             bevp1++;
3217             bevp2--;
3218           }
3219         }
3220       }
3221     }
3222     MEM_freeN(sortdata);
3223   }
3224
3225   /* STEP 4: 2D-COSINES or 3D ORIENTATION */
3226   if ((cu->flag & CU_3D) == 0) {
3227     /* 2D Curves */
3228     for (bl = bev->first; bl; bl = bl->next) {
3229       if (bl->nr < 2) {
3230         BevPoint *bevp = bl->bevpoints;
3231         unit_qt(bevp->quat);
3232       }
3233       else if (bl->nr == 2) { /* 2 pnt, treat separate */
3234         make_bevel_list_segment_2D(bl);
3235       }
3236       else {
3237         make_bevel_list_2D(bl);
3238       }
3239     }
3240   }
3241   else {
3242     /* 3D Curves */
3243     for (bl = bev->first; bl; bl = bl->next) {
3244       if (bl->nr < 2) {
3245         BevPoint *bevp = bl->bevpoints;
3246         unit_qt(bevp->quat);
3247       }
3248       else if (bl->nr == 2) { /* 2 pnt, treat separate */
3249         make_bevel_list_segment_3D(bl);
3250       }
3251       else {
3252         make_bevel_list_3D(bl, (int)(resolu * cu->twist_smooth), cu->twist_mode);
3253       }
3254     }
3255   }
3256 }
3257
3258 /* ****************** HANDLES ************** */
3259
3260 static void calchandleNurb_intern(BezTriple *bezt,
3261                                   const BezTriple *prev,
3262                                   const BezTriple *next,
3263                                   bool is_fcurve,
3264                                   bool skip_align,
3265                                   char fcurve_smoothing)
3266 {
3267   /* defines to avoid confusion */
3268 #define p2_h1 ((p2)-3)
3269 #define p2_h2 ((p2) + 3)
3270
3271   const float *p1, *p3;
3272   float *p2;
3273   float pt[3];
3274   float dvec_a[3], dvec_b[3];
3275   float len, len_a, len_b;
3276   float len_ratio;
3277   const float eps = 1e-5;
3278
3279   /* assume normal handle until we check */
3280   bezt->f5 = HD_AUTOTYPE_NORMAL;
3281
3282   if (bezt->h1 == 0 && bezt->h2 == 0) {
3283     return;
3284   }
3285
3286   p2 = bezt->vec[1];
3287
3288   if (prev == NULL) {
3289     p3 = next->vec[1];
3290     pt[0] = 2.0f * p2[0] - p3[0];
3291     pt[1] = 2.0f * p2[1] - p3[1];
3292     pt[2] = 2.0f * p2[2] - p3[2];
3293     p1 = pt;
3294   }
3295   else {
3296     p1 = prev->vec[1];
3297   }
3298
3299   if (next == NULL) {
3300     pt[0] = 2.0f * p2[0] - p1[0];
3301     pt[1] = 2.0f * p2[1] - p1[1];
3302     pt[2] = 2.0f * p2[2] - p1[2];
3303     p3 = pt;
3304   }
3305   else {
3306     p3 = next->vec[1];
3307   }
3308
3309   sub_v3_v3v3(dvec_a, p2, p1);
3310   sub_v3_v3v3(dvec_b, p3, p2);
3311
3312   if (is_fcurve) {
3313     len_a = dvec_a[0];
3314     len_b = dvec_b[0];
3315   }
3316   else {
3317     len_a = len_v3(dvec_a);
3318     len_b = len_v3(dvec_b);
3319   }
3320
3321   if (len_a == 0.0f) {
3322     len_a = 1.0f;
3323   }
3324   if (len_b == 0.0f) {
3325     len_b = 1.0f;
3326   }
3327
3328   len_ratio = len_a / len_b;
3329
3330   if (ELEM(bezt->h1, HD_AUTO, HD_AUTO_ANIM) || ELEM(bezt->h2, HD_AUTO, HD_AUTO_ANIM)) { /* auto */
3331     float tvec[3];
3332     tvec[0] = dvec_b[0] / len_b + dvec_a[0] / len_a;
3333     tvec[1] = dvec_b[1] / len_b + dvec_a[1] / len_a;
3334     tvec[2] = dvec_b[2] / len_b + dvec_a[2] / len_a;
3335
3336     if (is_fcurve) {
3337       if (fcurve_smoothing != FCURVE_SMOOTH_NONE) {
3338         /* force the horizontal handle size to be 1/3 of the key interval so that
3339          * the X component of the parametric bezier curve is a linear spline */
3340         len = 6.0f / 2.5614f;
3341       }
3342       else {
3343         len = tvec[0];
3344       }
3345     }
3346     else {
3347       len = len_v3(tvec);
3348     }
3349     len *= 2.5614f;
3350
3351     if (len != 0.0f) {
3352       /* only for fcurves */
3353       bool leftviolate = false, rightviolate = false;
3354
3355       if (!is_fcurve || fcurve_smoothing == FCURVE_SMOOTH_NONE) {
3356         if (len_a > 5.0f * len_b) {
3357           len_a = 5.0f * len_b;
3358         }
3359         if (len_b > 5.0f * len_a) {
3360           len_b = 5.0f * len_a;
3361         }
3362       }
3363
3364       if (ELEM(bezt->h1, HD_AUTO, HD_AUTO_ANIM)) {
3365         len_a /= len;
3366         madd_v3_v3v3fl(p2_h1, p2, tvec, -len_a);
3367
3368         if ((bezt->h1 == HD_AUTO_ANIM) && next && prev) { /* keep horizontal if extrema */
3369           float ydiff1 = prev->vec[1][1] - bezt->vec[1][1];
3370           float ydiff2 = next->vec[1][1] - bezt->vec[1][1];
3371           if ((ydiff1 <= 0.0f && ydiff2 <= 0.0f) || (ydiff1 >= 0.0f && ydiff2 >= 0.0f)) {
3372             bezt->vec[0][1] = bezt->vec[1][1];
3373             bezt->f5 = HD_AUTOTYPE_SPECIAL;
3374           }
3375           else { /* handles should not be beyond y coord of two others */
3376             if (ydiff1 <= 0.0f) {
3377               if (prev->vec[1][1] > bezt->vec[0][1]) {
3378                 bezt->vec[0][1] = prev->vec[1][1];
3379                 leftviolate = 1;
3380               }
3381             }
3382             else {
3383               if (prev->vec[1][1] < bezt->vec[0][1]) {
3384                 bezt->vec[0][1] = prev->vec[1][1];
3385                 leftviolate = 1;
3386               }
3387             }
3388           }
3389         }
3390       }
3391       if (ELEM(bezt->h2, HD_AUTO, HD_AUTO_ANIM)) {
3392         len_b /= len;
3393         madd_v3_v3v3fl(p2_h2, p2, tvec, len_b);
3394
3395         if ((bezt->h2 == HD_AUTO_ANIM) && next && prev) { /* keep horizontal if extrema */
3396           float ydiff1 = prev->vec[1][1] - bezt->vec[1][1];
3397           float ydiff2 = next->vec[1][1] - bezt->vec[1][1];
3398           if ((ydiff1 <= 0.0f && ydiff2 <= 0.0f) || (ydiff1 >= 0.0f && ydiff2 >= 0.0f)) {
3399             bezt->vec[2][1] = bezt->vec[1][1];
3400             bezt->f5 = HD_AUTOTYPE_SPECIAL;
3401           }
3402           else { /* handles should not be beyond y coord of two others */
3403             if (ydiff1 <= 0.0f) {
3404               if (next->vec[1][1] < bezt->vec[2][1]) {
3405                 bezt->vec[2][1] = next->vec[1][1];
3406                 rightviolate = 1;
3407               }
3408             }
3409             else {
3410               if (next->vec[1][1] > bezt->vec[2][1]) {
3411                 bezt->vec[2][1] = next->vec[1][1];
3412                 rightviolate = 1;
3413               }
3414             }
3415           }
3416         }
3417       }
3418       if (leftviolate || rightviolate) { /* align left handle */
3419         BLI_assert(is_fcurve);
3420         /* simple 2d calculation */
3421         float h1_x = p2_h1[0] - p2[0];
3422         float h2_x = p2[0] - p2_h2[0];
3423
3424         if (leftviolate) {
3425           p2_h2[1] = p2[1] + ((p2[1] - p2_h1[1]) / h1_x) * h2_x;
3426         }
3427         else {
3428           p2_h1[1] = p2[1] + ((p2[1] - p2_h2[1]) / h2_x) * h1_x;
3429         }
3430       }
3431     }
3432   }
3433
3434   if (bezt->h1 == HD_VECT) { /* vector */
3435     madd_v3_v3v3fl(p2_h1, p2, dvec_a, -1.0f / 3.0f);
3436   }
3437   if (bezt->h2 == HD_VECT) {
3438     madd_v3_v3v3fl(p2_h2, p2, dvec_b, 1.0f / 3.0f);
3439   }
3440
3441   if (skip_align ||
3442       /* when one handle is free, alignming makes no sense, see: T35952 */
3443       (ELEM(HD_FREE, bezt->h1, bezt->h2)) ||
3444       /* also when no handles are aligned, skip this step */
3445       (!ELEM(HD_ALIGN, bezt->h1, bezt->h2) && !ELEM(HD_ALIGN_DOUBLESIDE, bezt->h1, bezt->h2))) {
3446     /* handles need to be updated during animation and applying stuff like hooks,
3447      * but in such situations it's quite difficult to distinguish in which order
3448      * align handles should be aligned so skip them for now */
3449     return;
3450   }
3451
3452   len_a = len_v3v3(p2, p2_h1);
3453   len_b = len_v3v3(p2, p2_h2);
3454
3455   if (len_a == 0.0f) {
3456     len_a = 1.0f;
3457   }
3458   if (len_b == 0.0f) {
3459     len_b = 1.0f;
3460   }
3461
3462   len_ratio = len_a / len_b;
3463
3464   if (bezt->f1 & SELECT) {                               /* order of calculation */
3465     if (ELEM(bezt->h2, HD_ALIGN, HD_ALIGN_DOUBLESIDE)) { /* aligned */
3466       if (len_a > eps) {
3467         len = 1.0f / len_ratio;
3468         p2_h2[0] = p2[0] + len * (p2[0] - p2_h1[0]);
3469         p2_h2[1] = p2[1] + len * (p2[1] - p2_h1[1]);
3470         p2_h2[2] = p2[2] + len * (p2[2] - p2_h1[2]);
3471       }
3472     }
3473     if (ELEM(bezt->h1, HD_ALIGN, HD_ALIGN_DOUBLESIDE)) {
3474       if (len_b > eps) {
3475         len = len_ratio;
3476         p2_h1[0] = p2[0] + len * (p2[0] - p2_h2[0]);
3477         p2_h1[1] = p2[1] + len * (p2[1] - p2_h2[1]);
3478         p2_h1[2] = p2[2] + len * (p2[2] - p2_h2[2]);
3479       }
3480     }
3481   }
3482   else {
3483     if (ELEM(bezt->h1, HD_ALIGN, HD_ALIGN_DOUBLESIDE)) {
3484       if (len_b > eps) {
3485         len = len_ratio;
3486         p2_h1[0] = p2[0] + len * (p2[0] - p2_h2[0]);
3487         p2_h1[1] = p2[1] + len * (p2[1] - p2_h2[1]);
3488         p2_h1[2] = p2[2] + len * (p2[2] - p2_h2[2]);
3489       }
3490     }
3491     if (ELEM(bezt->h2, HD_ALIGN, HD_ALIGN_DOUBLESIDE)) { /* aligned */
3492       if (len_a > eps) {
3493         len = 1.0f / len_ratio;
3494         p2_h2[0] = p2[0] + len * (p2[0] - p2_h1[0]);
3495         p2_h2[1] = p2[1] + len * (p2[1] - p2_h1[1]);
3496         p2_h2[2] = p2[2] + len * (p2[2] - p2_h1[2]);
3497       }
3498     }
3499   }
3500
3501 #undef p2_h1
3502 #undef p2_h2
3503 }
3504
3505 static void calchandlesNurb_intern(Nurb *nu, bool skip_align)
3506 {
3507   BezTriple *bezt, *prev, *next;
3508   int a;
3509
3510   if (nu->type != CU_BEZIER) {
3511     return;
3512   }
3513   if (nu->pntsu < 2) {
3514     return;
3515   }
3516
3517   a = nu->pntsu;
3518   bezt = nu->bezt;
3519   if (nu->flagu & CU_NURB_CYCLIC) {
3520     prev = bezt + (a - 1);
3521   }
3522   else {
3523     prev = NULL;
3524   }
3525   next = bezt + 1;
3526
3527   while (a--) {
3528     calchandleNurb_intern(bezt, prev, next, 0, skip_align, 0);
3529     prev = bezt;
3530     if (a == 1) {
3531       if (nu->flagu & CU_NURB_CYCLIC) {
3532         next = nu->bezt;
3533       }
3534       else {
3535         next = NULL;
3536       }
3537     }
3538     else {
3539       next++;
3540     }
3541
3542     bezt++;
3543   }
3544 }
3545
3546 /* A utility function for allocating a number of arrays of the same length
3547  * with easy error checking and deallocation, and an easy way to add or remove
3548  * arrays that are processed in this way when changing code.
3549  *
3550  * floats, chars: NULL-terminated arrays of pointers to array pointers that need to be allocated.
3551  *
3552  * Returns: pointer to the buffer that contains all of the arrays.
3553  */
3554 static void *allocate_arrays(int count, float ***floats, char ***chars, const char *name)
3555 {
3556   size_t num_floats = 0, num_chars = 0;
3557
3558   while (floats && floats[num_floats]) {
3559     num_floats++;
3560   }
3561
3562   while (chars && chars[num_chars]) {
3563     num_chars++;
3564   }
3565
3566   void *buffer = (float *)MEM_malloc_arrayN(count, (sizeof(float) * num_floats + num_chars), name);
3567
3568   if (!buffer) {
3569     return NULL;
3570   }
3571
3572   float *fptr = buffer;
3573
3574   for (int i = 0; i < num_floats; i++, fptr += count) {
3575     *floats[i] = fptr;
3576   }
3577
3578   char *cptr = (char *)fptr;
3579
3580   for (int i = 0; i < num_chars; i++, cptr += count) {
3581     *chars[i] = cptr;
3582   }
3583
3584   return buffer;
3585 }
3586
3587 static void free_arrays(void *buffer)
3588 {
3589   MEM_freeN(buffer);
3590 }
3591
3592 /* computes in which direction to change h[i] to satisfy conditions better */
3593 static float bezier_relax_direction(
3594     float *a, float *b, float *c, float *d, float *h, int i, int count)
3595 {
3596   /* current deviation between sides of the equation */
3597   float state = a[i] * h[(i + count - 1) % count] + b[i] * h[i] + c[i] * h[(i + 1) % count] - d[i];
3598
3599   /* only the sign is meaningful */
3600   return -state * b[i];
3601 }
3602
3603 static void bezier_lock_unknown(float *a, float *b, float *c, float *d, int i, float value)
3604 {
3605   a[i] = c[i] = 0.0f;
3606   b[i] = 1.0f;
3607   d[i] = value;
3608 }
3609
3610 static void bezier_restore_equation(
3611     float *a, float *b, float *c, float *d, float *a0, float *b0, float *c0, float *d0, int i)
3612 {
3613   a[i] = a0[i];
3614   b[i] = b0[i];
3615   c[i] = c0[i];
3616   d[i] = d0[i];
3617 }
3618
3619 static bool tridiagonal_solve_with_limits(
3620     float *a, float *b, float *c, float *d, float *h, float *hmin, float *hmax, int solve_count)
3621 {
3622   float *a0, *b0, *c0, *d0;
3623   float **arrays[] = {&a0, &b0, &c0, &d0, NULL};
3624   char *is_locked, *num_unlocks;
3625   char **flagarrays[] = {&is_locked, &num_unlocks, NULL};
3626
3627   void *tmps = allocate_arrays(solve_count, arrays, flagarrays, "tridiagonal_solve_with_limits");
3628   if (!tmps) {
3629     return false;
3630   }
3631
3632   memcpy(a0, a, sizeof(float) * solve_count);
3633   memcpy(b0, b, sizeof(float) * solve_count);
3634   memcpy(c0, c, sizeof(float) * solve_count);
3635   memcpy(d0, d, sizeof(float) * solve_count);
3636
3637   memset(is_locked, 0, solve_count);
3638   memset(num_unlocks, 0, solve_count);
3639
3640   bool overshoot, unlocked;
3641
3642   do {
3643     if (!BLI_tridiagonal_solve_cyclic(a, b, c, d, h, solve_count)) {
3644       free_arrays(tmps);
3645       return false;
3646     }
3647
3648     /* first check if any handles overshoot the limits, and lock them */
3649     bool all = false, locked = false;
3650
3651     overshoot = unlocked = false;
3652
3653     do {
3654       for (int i = 0; i < solve_count; i++) {
3655         if (h[i] >= hmin[i] && h[i] <= hmax[i]) {
3656           continue;
3657         }
3658
3659         overshoot = true;
3660
3661         float target = h[i] > hmax[i] ? hmax[i] : hmin[i];
3662
3663         /* heuristically only lock handles that go in the right direction if there are such ones */
3664         if (target != 0.0f || all) {
3665           /* mark item locked */
3666           is_locked[i] = 1;
3667
3668           bezier_lock_unknown(a, b, c, d, i, target);
3669           locked = true;
3670         }
3671       }
3672
3673       all = true;
3674     } while (overshoot && !locked);
3675
3676     /* If no handles overshot and were locked,
3677      * see if it may be a good idea to unlock some handles. */
3678     if (!locked) {
3679       for (int i = 0; i < solve_count; i++) {
3680         // to definitely avoid infinite loops limit this to 2 times
3681         if (!is_locked[i] || num_unlocks[i] >= 2) {
3682           continue;
3683         }
3684
3685         /* if the handle wants to move in allowable direction, release it */
3686         float relax = bezier_relax_direction(a0, b0, c0, d0, h, i, solve_count);
3687
3688         if ((relax > 0 && h[i] < hmax[i]) || (relax < 0 && h[i] > hmin[i])) {
3689           bezier_restore_equation(a, b, c, d, a0, b0, c0, d0, i);
3690
3691           is_locked[i] = 0;
3692           num_unlocks[i]++;
3693           unlocked = true;
3694         }
3695       }
3696     }
3697   } while (overshoot || unlocked);
3698
3699   free_arrays(tmps);
3700   return true;
3701 }
3702
3703 /* Keep ascii art. */
3704 /* clang-format off */
3705 /*
3706  * This function computes the handles of a series of auto bezier points
3707  * on the basis of 'no acceleration discontinuities' at the points.
3708  * The first and last bezier points are considered 'fixed' (their handles are not touched)
3709  * The result is the smoothest possible trajectory going through intermediate points.
3710  * The difficulty is that the handles depends on their neighbors.
3711  *
3712  * The exact solution is found by solving a tridiagonal matrix equation formed
3713  * by the continuity and boundary conditions. Although theoretically handle position
3714  * is affected by all other points of the curve segment, in practice the influence
3715  * decreases exponentially with distance.
3716  *
3717  * Note: this algorithm assumes that the handle horizontal size if always 1/3 of the
3718  * of the interval to the next point. This rule ensures linear interpolation of time.
3719  *
3720  * ^ height (co 1)
3721  * |                                            yN
3722  * |                                   yN-1     |
3723  * |                      y2           |        |
3724  * |           y1         |            |        |
3725  * |    y0     |          |            |        |
3726  * |    |      |          |            |        |
3727  * |    |      |          |            |        |
3728  * |    |      |          |            |        |
3729  * |-------t1---------t2--------- ~ --------tN-------------------> time (co 0)
3730  * Mathematical basis:
3731  *
3732  * 1. Handle lengths on either side of each point are connected by a factor
3733  *    ensuring continuity of the first derivative:
3734  *
3735  *    l[i] = t[i+1]/t[i]
3736  *
3737  * 2. The tridiagonal system is formed by the following equation, which is derived
3738  *    by differentiating the bezier curve and specifies second derivative continuity
3739  *    at every point:
3740  *
3741  *    l[i]^2 * h[i-1] + (2*l[i]+2) * h[i] + 1/l[i+1] * h[i+1] = (y[i]-y[i-1])*l[i]^2 + y[i+1]-y[i]
3742  *
3743  * 3. If this point is adjacent to a manually set handle with X size not equal to 1/3
3744  *    of the horizontal interval, this equation becomes slightly more complex: