Fix T64842: crash rendering files with bevel curves
[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(Object *ob, ListBase *disp)
1743 {
1744   DispList *dl, *dlnew;
1745   Curve *bevcu, *cu;
1746   float *fp, facx, facy, angle, dangle;
1747   int nr, a;
1748
1749   cu = ob->data;
1750   BLI_listbase_clear(disp);
1751
1752   /* if a font object is being edited, then do nothing */
1753   // XXX  if ( ob == obedit && ob->type == OB_FONT ) return;
1754
1755   if (cu->bevobj) {
1756     if (cu->bevobj->type != OB_CURVE) {
1757       return;
1758     }
1759
1760     bevcu = cu->bevobj->data;
1761     if (bevcu->ext1 == 0.0f && bevcu->ext2 == 0.0f) {
1762       ListBase bevdisp = {NULL, NULL};
1763       facx = cu->bevobj->scale[0];
1764       facy = cu->bevobj->scale[1];
1765
1766       if (cu->bevobj->runtime.curve_cache) {
1767         dl = cu->bevobj->runtime.curve_cache->disp.first;
1768       }
1769       else {
1770         BLI_assert(cu->bevobj->runtime.curve_cache != NULL);
1771         dl = NULL;
1772       }
1773
1774       while (dl) {
1775         if (ELEM(dl->type, DL_POLY, DL_SEGM)) {
1776           dlnew = MEM_mallocN(sizeof(DispList), "makebevelcurve1");
1777           *dlnew = *dl;
1778           dlnew->verts = MEM_malloc_arrayN(
1779               dl->parts * dl->nr, 3 * sizeof(float), "makebevelcurve1");
1780           memcpy(dlnew->verts, dl->verts, 3 * sizeof(float) * dl->parts * dl->nr);
1781
1782           if (dlnew->type == DL_SEGM) {
1783             dlnew->flag |= (DL_FRONT_CURVE | DL_BACK_CURVE);
1784           }
1785
1786           BLI_addtail(disp, dlnew);
1787           fp = dlnew->verts;
1788           nr = dlnew->parts * dlnew->nr;
1789           while (nr--) {
1790             fp[2] = fp[1] * facy;
1791             fp[1] = -fp[0] * facx;
1792             fp[0] = 0.0;
1793             fp += 3;
1794           }
1795         }
1796         dl = dl->next;
1797       }
1798
1799       BKE_displist_free(&bevdisp);
1800     }
1801   }
1802   else if (cu->ext1 == 0.0f && cu->ext2 == 0.0f) {
1803     /* pass */
1804   }
1805   else if (cu->ext2 == 0.0f) {
1806     dl = MEM_callocN(sizeof(DispList), "makebevelcurve2");
1807     dl->verts = MEM_malloc_arrayN(2, sizeof(float[3]), "makebevelcurve2");
1808     BLI_addtail(disp, dl);
1809     dl->type = DL_SEGM;
1810     dl->parts = 1;
1811     dl->flag = DL_FRONT_CURVE | DL_BACK_CURVE;
1812     dl->nr = 2;
1813
1814     fp = dl->verts;
1815     fp[0] = fp[1] = 0.0;
1816     fp[2] = -cu->ext1;
1817     fp[3] = fp[4] = 0.0;
1818     fp[5] = cu->ext1;
1819   }
1820   else if ((cu->flag & (CU_FRONT | CU_BACK)) == 0 &&
1821            cu->ext1 == 0.0f) { /* we make a full round bevel in that case */
1822     nr = 4 + 2 * cu->bevresol;
1823
1824     dl = MEM_callocN(sizeof(DispList), "makebevelcurve p1");
1825     dl->verts = MEM_malloc_arrayN(nr, sizeof(float[3]), "makebevelcurve p1");
1826     BLI_addtail(disp, dl);
1827     dl->type = DL_POLY;
1828     dl->parts = 1;
1829     dl->flag = DL_BACK_CURVE;
1830     dl->nr = nr;
1831
1832     /* a circle */
1833     fp = dl->verts;
1834     dangle = (2.0f * (float)M_PI / (nr));
1835     angle = -(nr - 1) * dangle;
1836
1837     for (a = 0; a < nr; a++) {
1838       fp[0] = 0.0;
1839       fp[1] = (cosf(angle) * (cu->ext2));
1840       fp[2] = (sinf(angle) * (cu->ext2)) - cu->ext1;
1841       angle += dangle;
1842       fp += 3;
1843     }
1844   }
1845   else {
1846     short dnr;
1847
1848     /* bevel now in three parts, for proper vertex normals */
1849     /* part 1, back */
1850
1851     if ((cu->flag & CU_BACK) || !(cu->flag & CU_FRONT)) {
1852       dnr = nr = 2 + cu->bevresol;
1853       if ((cu->flag & (CU_FRONT | CU_BACK)) == 0) {
1854         nr = 3 + 2 * cu->bevresol;
1855       }
1856       dl = MEM_callocN(sizeof(DispList), "makebevelcurve p1");
1857       dl->verts = MEM_malloc_arrayN(nr, sizeof(float[3]), "makebevelcurve p1");
1858       BLI_addtail(disp, dl);
1859       dl->type = DL_SEGM;
1860       dl->parts = 1;
1861       dl->flag = DL_BACK_CURVE;
1862       dl->nr = nr;
1863
1864       /* half a circle */
1865       fp = dl->verts;
1866       dangle = ((float)M_PI_2 / (dnr - 1));
1867       angle = -(nr - 1) * dangle;
1868
1869       for (a = 0; a < nr; a++) {
1870         fp[0] = 0.0;
1871         fp[1] = (float)(cosf(angle) * (cu->ext2));
1872         fp[2] = (float)(sinf(angle) * (cu->ext2)) - cu->ext1;
1873         angle += dangle;
1874         fp += 3;
1875       }
1876     }
1877
1878     /* part 2, sidefaces */
1879     if (cu->ext1 != 0.0f) {
1880       nr = 2;
1881
1882       dl = MEM_callocN(sizeof(DispList), "makebevelcurve p2");
1883       dl->verts = MEM_malloc_arrayN(nr, sizeof(float[3]), "makebevelcurve p2");
1884       BLI_addtail(disp, dl);
1885       dl->type = DL_SEGM;
1886       dl->parts = 1;
1887       dl->nr = nr;
1888
1889       fp = dl->verts;
1890       fp[1] = cu->ext2;
1891       fp[2] = -cu->ext1;
1892       fp[4] = cu->ext2;
1893       fp[5] = cu->ext1;
1894
1895       if ((cu->flag & (CU_FRONT | CU_BACK)) == 0) {
1896         dl = MEM_dupallocN(dl);
1897         dl->verts = MEM_dupallocN(dl->verts);
1898         BLI_addtail(disp, dl);
1899
1900         fp = dl->verts;
1901         fp[1] = -fp[1];
1902         fp[2] = -fp[2];
1903         fp[4] = -fp[4];
1904         fp[5] = -fp[5];
1905       }
1906     }
1907
1908     /* part 3, front */
1909     if ((cu->flag & CU_FRONT) || !(cu->flag & CU_BACK)) {
1910       dnr = nr = 2 + cu->bevresol;
1911       if ((cu->flag & (CU_FRONT | CU_BACK)) == 0) {
1912         nr = 3 + 2 * cu->bevresol;
1913       }
1914       dl = MEM_callocN(sizeof(DispList), "makebevelcurve p3");
1915       dl->verts = MEM_malloc_arrayN(nr, sizeof(float[3]), "makebevelcurve p3");
1916       BLI_addtail(disp, dl);
1917       dl->type = DL_SEGM;
1918       dl->flag = DL_FRONT_CURVE;
1919       dl->parts = 1;
1920       dl->nr = nr;
1921
1922       /* half a circle */
1923       fp = dl->verts;
1924       angle = 0.0;
1925       dangle = ((float)M_PI_2 / (dnr - 1));
1926
1927       for (a = 0; a < nr; a++) {
1928         fp[0] = 0.0;
1929         fp[1] = (float)(cosf(angle) * (cu->ext2));
1930         fp[2] = (float)(sinf(angle) * (cu->ext2)) + cu->ext1;
1931         angle += dangle;
1932         fp += 3;
1933       }
1934     }
1935   }
1936 }
1937
1938 static int cu_isectLL(const float v1[3],
1939                       const float v2[3],
1940                       const float v3[3],
1941                       const float v4[3],
1942                       short cox,
1943                       short coy,
1944                       float *lambda,
1945                       float *mu,
1946                       float vec[3])
1947 {
1948   /* return:
1949    * -1: collinear
1950    *  0: no intersection of segments
1951    *  1: exact intersection of segments
1952    *  2: cross-intersection of segments
1953    */
1954   float deler;
1955
1956   deler = (v1[cox] - v2[cox]) * (v3[coy] - v4[coy]) - (v3[cox] - v4[cox]) * (v1[coy] - v2[coy]);
1957   if (deler == 0.0f) {
1958     return -1;
1959   }
1960
1961   *lambda = (v1[coy] - v3[coy]) * (v3[cox] - v4[cox]) - (v1[cox] - v3[cox]) * (v3[coy] - v4[coy]);
1962   *lambda = -(*lambda / deler);
1963
1964   deler = v3[coy] - v4[coy];
1965   if (deler == 0) {
1966     deler = v3[cox] - v4[cox];
1967     *mu = -(*lambda * (v2[cox] - v1[cox]) + v1[cox] - v3[cox]) / deler;
1968   }
1969   else {
1970     *mu = -(*lambda * (v2[coy] - v1[coy]) + v1[coy] - v3[coy]) / deler;
1971   }
1972   vec[cox] = *lambda * (v2[cox] - v1[cox]) + v1[cox];
1973   vec[coy] = *lambda * (v2[coy] - v1[coy]) + v1[coy];
1974
1975   if (*lambda >= 0.0f && *lambda <= 1.0f && *mu >= 0.0f && *mu <= 1.0f) {
1976     if (*lambda == 0.0f || *lambda == 1.0f || *mu == 0.0f || *mu == 1.0f) {
1977       return 1;
1978     }
1979     return 2;
1980   }
1981   return 0;
1982 }
1983
1984 static bool bevelinside(BevList *bl1, BevList *bl2)
1985 {
1986   /* is bl2 INSIDE bl1 ? with left-right method and "lambda's" */
1987   /* returns '1' if correct hole  */
1988   BevPoint *bevp, *prevbevp;
1989   float min, max, vec[3], hvec1[3], hvec2[3], lab, mu;
1990   int nr, links = 0, rechts = 0, mode;
1991
1992   /* take first vertex of possible hole */
1993
1994   bevp = bl2->bevpoints;
1995   hvec1[0] = bevp->vec[0];
1996   hvec1[1] = bevp->vec[1];
1997   hvec1[2] = 0.0;
1998   copy_v3_v3(hvec2, hvec1);
1999   hvec2[0] += 1000;
2000
2001   /* test it with all edges of potential surrounding poly */
2002   /* count number of transitions left-right  */
2003
2004   bevp = bl1->bevpoints;
2005   nr = bl1->nr;
2006   prevbevp = bevp + (nr - 1);
2007
2008   while (nr--) {
2009     min = prevbevp->vec[1];
2010     max = bevp->vec[1];
2011     if (max < min) {
2012       min = max;
2013       max = prevbevp->vec[1];
2014     }
2015     if (min != max) {
2016       if (min <= hvec1[1] && max >= hvec1[1]) {
2017         /* there's a transition, calc intersection point */
2018         mode = cu_isectLL(prevbevp->vec, bevp->vec, hvec1, hvec2, 0, 1, &lab, &mu, vec);
2019         /* if lab==0.0 or lab==1.0 then the edge intersects exactly a transition
2020          * only allow for one situation: we choose lab= 1.0
2021          */
2022         if (mode >= 0 && lab != 0.0f) {
2023           if (vec[0] < hvec1[0]) {
2024             links++;
2025           }
2026           else {
2027             rechts++;
2028           }
2029         }
2030       }
2031     }
2032     prevbevp = bevp;
2033     bevp++;
2034   }
2035
2036   return (links & 1) && (rechts & 1);
2037 }
2038
2039 struct BevelSort {
2040   BevList *bl;
2041   float left;
2042   int dir;
2043 };
2044
2045 static int vergxcobev(const void *a1, const void *a2)
2046 {
2047   const struct BevelSort *x1 = a1, *x2 = a2;
2048
2049   if (x1->left > x2->left) {
2050     return 1;
2051   }
2052   else if (x1->left < x2->left) {
2053     return -1;
2054   }
2055   return 0;
2056 }
2057
2058 /* this function cannot be replaced with atan2, but why? */
2059
2060 static void calc_bevel_sin_cos(
2061     float x1, float y1, float x2, float y2, float *r_sina, float *r_cosa)
2062 {
2063   float t01, t02, x3, y3;
2064
2065   t01 = sqrtf(x1 * x1 + y1 * y1);
2066   t02 = sqrtf(x2 * x2 + y2 * y2);
2067   if (t01 == 0.0f) {
2068     t01 = 1.0f;
2069   }
2070   if (t02 == 0.0f) {
2071     t02 = 1.0f;
2072   }
2073
2074   x1 /= t01;
2075   y1 /= t01;
2076   x2 /= t02;
2077   y2 /= t02;
2078
2079   t02 = x1 * x2 + y1 * y2;
2080   if (fabsf(t02) >= 1.0f) {
2081     t02 = M_PI_2;
2082   }
2083   else {
2084     t02 = (saacos(t02)) / 2.0f;
2085   }
2086
2087   t02 = sinf(t02);
2088   if (t02 == 0.0f) {
2089     t02 = 1.0f;
2090   }
2091
2092   x3 = x1 - x2;
2093   y3 = y1 - y2;
2094   if (x3 == 0 && y3 == 0) {
2095     x3 = y1;
2096     y3 = -x1;
2097   }
2098   else {
2099     t01 = sqrtf(x3 * x3 + y3 * y3);
2100     x3 /= t01;
2101     y3 /= t01;
2102   }
2103
2104   *r_sina = -y3 / t02;
2105   *r_cosa = x3 / t02;
2106 }
2107
2108 static void tilt_bezpart(BezTriple *prevbezt,
2109                          BezTriple *bezt,
2110                          Nurb *nu,
2111                          float *tilt_array,
2112                          float *radius_array,
2113                          float *weight_array,
2114                          int resolu,
2115                          int stride)
2116 {
2117   BezTriple *pprev, *next, *last;
2118   float fac, dfac, t[4];
2119   int a;
2120
2121   if (tilt_array == NULL && radius_array == NULL) {
2122     return;
2123   }
2124
2125   last = nu->bezt + (nu->pntsu - 1);
2126
2127   /* returns a point */
2128   if (prevbezt == nu->bezt) {
2129     if (nu->flagu & CU_NURB_CYCLIC) {
2130       pprev = last;
2131     }
2132     else {
2133       pprev = prevbezt;
2134     }
2135   }
2136   else {
2137     pprev = prevbezt - 1;
2138   }
2139
2140   /* next point */
2141   if (bezt == last) {
2142     if (nu->flagu & CU_NURB_CYCLIC) {
2143       next = nu->bezt;
2144     }
2145     else {
2146       next = bezt;
2147     }
2148   }
2149   else {
2150     next = bezt + 1;
2151   }
2152
2153   fac = 0.0;
2154   dfac = 1.0f / (float)resolu;
2155
2156   for (a = 0; a < resolu; a++, fac += dfac) {
2157     if (tilt_array) {
2158       if (nu->tilt_interp ==
2159           KEY_CU_EASE) { /* May as well support for tilt also 2.47 ease interp */
2160         *tilt_array = prevbezt->tilt +
2161                       (bezt->tilt - prevbezt->tilt) * (3.0f * fac * fac - 2.0f * fac * fac * fac);
2162       }
2163       else {
2164         key_curve_position_weights(fac, t, nu->tilt_interp);
2165         *tilt_array = t[0] * pprev->tilt + t[1] * prevbezt->tilt + t[2] * bezt->tilt +
2166                       t[3] * next->tilt;
2167       }
2168
2169       tilt_array = POINTER_OFFSET(tilt_array, stride);
2170     }
2171
2172     if (radius_array) {
2173       if (nu->radius_interp == KEY_CU_EASE) {
2174         /* Support 2.47 ease interp
2175          * Note! - this only takes the 2 points into account,
2176          * giving much more localized results to changes in radius, sometimes you want that */
2177         *radius_array = prevbezt->radius + (bezt->radius - prevbezt->radius) *
2178                                                (3.0f * fac * fac - 2.0f * fac * fac * fac);
2179       }
2180       else {
2181
2182         /* reuse interpolation from tilt if we can */
2183         if (tilt_array == NULL || nu->tilt_interp != nu->radius_interp) {
2184           key_curve_position_weights(fac, t, nu->radius_interp);
2185         }
2186         *radius_array = t[0] * pprev->radius + t[1] * prevbezt->radius + t[2] * bezt->radius +
2187                         t[3] * next->radius;
2188       }
2189
2190       radius_array = POINTER_OFFSET(radius_array, stride);
2191     }
2192
2193     if (weight_array) {
2194       /* basic interpolation for now, could copy tilt interp too  */
2195       *weight_array = prevbezt->weight + (bezt->weight - prevbezt->weight) *
2196                                              (3.0f * fac * fac - 2.0f * fac * fac * fac);
2197
2198       weight_array = POINTER_OFFSET(weight_array, stride);
2199     }
2200   }
2201 }
2202
2203 /* make_bevel_list_3D_* funcs, at a minimum these must
2204  * fill in the bezp->quat and bezp->dir values */
2205
2206 /* utility for make_bevel_list_3D_* funcs */
2207 static void bevel_list_calc_bisect(BevList *bl)
2208 {
2209   BevPoint *bevp2, *bevp1, *bevp0;
2210   int nr;
2211   bool is_cyclic = bl->poly != -1;
2212
2213   if (is_cyclic) {
2214     bevp2 = bl->bevpoints;
2215     bevp1 = bevp2 + (bl->nr - 1);
2216     bevp0 = bevp1 - 1;
2217     nr = bl->nr;
2218   }
2219   else {
2220     /* If spline is not cyclic, direction of first and
2221      * last bevel points matches direction of CV handle.
2222      *
2223      * This is getting calculated earlier when we know
2224      * CV's handles and here we might simply skip evaluation
2225      * of direction for this guys.
2226      */
2227
2228     bevp0 = bl->bevpoints;
2229     bevp1 = bevp0 + 1;
2230     bevp2 = bevp1 + 1;
2231
2232     nr = bl->nr - 2;
2233   }
2234
2235   while (nr--) {
2236     /* totally simple */
2237     bisect_v3_v3v3v3(bevp1->dir, bevp0->vec, bevp1->vec, bevp2->vec);
2238
2239     bevp0 = bevp1;
2240     bevp1 = bevp2;
2241     bevp2++;
2242   }
2243 }
2244 static void bevel_list_flip_tangents(BevList *bl)
2245 {
2246   BevPoint *bevp2, *bevp1, *bevp0;
2247   int nr;
2248
2249   bevp2 = bl->bevpoints;
2250   bevp1 = bevp2 + (bl->nr - 1);
2251   bevp0 = bevp1 - 1;
2252
2253   nr = bl->nr;
2254   while (nr--) {
2255     if (angle_normalized_v3v3(bevp0->tan, bevp1->tan) > DEG2RADF(90.0f)) {
2256       negate_v3(bevp1->tan);
2257     }
2258
2259     bevp0 = bevp1;
2260     bevp1 = bevp2;
2261     bevp2++;
2262   }
2263 }
2264 /* apply user tilt */
2265 static void bevel_list_apply_tilt(BevList *bl)
2266 {
2267   BevPoint *bevp2, *bevp1;
2268   int nr;
2269   float q[4];
2270
2271   bevp2 = bl->bevpoints;
2272   bevp1 = bevp2 + (bl->nr - 1);
2273
2274   nr = bl->nr;
2275   while (nr--) {
2276     axis_angle_to_quat(q, bevp1->dir, bevp1->tilt);
2277     mul_qt_qtqt(bevp1->quat, q, bevp1->quat);
2278     normalize_qt(bevp1->quat);
2279
2280     bevp1 = bevp2;
2281     bevp2++;
2282   }
2283 }
2284 /* smooth quats, this function should be optimized, it can get slow with many iterations. */
2285 static void bevel_list_smooth(BevList *bl, int smooth_iter)
2286 {
2287   BevPoint *bevp2, *bevp1, *bevp0;
2288   int nr;
2289
2290   float q[4];
2291   float bevp0_quat[4];
2292   int a;
2293
2294   for (a = 0; a < smooth_iter; a++) {
2295     bevp2 = bl->bevpoints;
2296     bevp1 = bevp2 + (bl->nr - 1);
2297     bevp0 = bevp1 - 1;
2298
2299     nr = bl->nr;
2300
2301     if (bl->poly == -1) { /* check its not cyclic */
2302       /* skip the first point */
2303       /* bevp0 = bevp1; */
2304       bevp1 = bevp2;
2305       bevp2++;
2306       nr--;
2307
2308       bevp0 = bevp1;
2309       bevp1 = bevp2;
2310       bevp2++;
2311       nr--;
2312     }
2313
2314     copy_qt_qt(bevp0_quat, bevp0->quat);
2315
2316     while (nr--) {
2317       /* interpolate quats */
2318       float zaxis[3] = {0, 0, 1}, cross[3], q2[4];
2319       interp_qt_qtqt(q, bevp0_quat, bevp2->quat, 0.5);
2320       normalize_qt(q);
2321
2322       mul_qt_v3(q, zaxis);
2323       cross_v3_v3v3(cross, zaxis, bevp1->dir);
2324       axis_angle_to_quat(q2, cross, angle_normalized_v3v3(zaxis, bevp1->dir));
2325       normalize_qt(q2);
2326
2327       copy_qt_qt(bevp0_quat, bevp1->quat);
2328       mul_qt_qtqt(q, q2, q);
2329       interp_qt_qtqt(bevp1->quat, bevp1->quat, q, 0.5);
2330       normalize_qt(bevp1->quat);
2331
2332       /* bevp0 = bevp1; */ /* UNUSED */
2333       bevp1 = bevp2;
2334       bevp2++;
2335     }
2336   }
2337 }
2338
2339 static void make_bevel_list_3D_zup(BevList *bl)
2340 {
2341   BevPoint *bevp = bl->bevpoints;
2342   int nr = bl->nr;
2343
2344   bevel_list_calc_bisect(bl);
2345
2346   while (nr--) {
2347     vec_to_quat(bevp->quat, bevp->dir, 5, 1);
2348     bevp++;
2349   }
2350 }
2351
2352 static void minimum_twist_between_two_points(BevPoint *current_point, BevPoint *previous_point)
2353 {
2354   float angle = angle_normalized_v3v3(previous_point->dir, current_point->dir);
2355   float q[4];
2356
2357   if (angle > 0.0f) { /* otherwise we can keep as is */
2358     float cross_tmp[3];
2359     cross_v3_v3v3(cross_tmp, previous_point->dir, current_point->dir);
2360     axis_angle_to_quat(q, cross_tmp, angle);
2361     mul_qt_qtqt(current_point->quat, q, previous_point->quat);
2362   }
2363   else {
2364     copy_qt_qt(current_point->quat, previous_point->quat);
2365   }
2366 }
2367
2368 static void make_bevel_list_3D_minimum_twist(BevList *bl)
2369 {
2370   BevPoint *bevp2, *bevp1, *bevp0; /* standard for all make_bevel_list_3D_* funcs */
2371   int nr;
2372   float q[4];
2373
2374   bevel_list_calc_bisect(bl);
2375
2376   bevp2 = bl->bevpoints;
2377   bevp1 = bevp2 + (bl->nr - 1);
2378   bevp0 = bevp1 - 1;
2379
2380   nr = bl->nr;
2381   while (nr--) {
2382
2383     if (nr + 4 > bl->nr) { /* first time and second time, otherwise first point adjusts last */
2384       vec_to_quat(bevp1->quat, bevp1->dir, 5, 1);
2385     }
2386     else {
2387       minimum_twist_between_two_points(bevp1, bevp0);
2388     }
2389
2390     bevp0 = bevp1;
2391     bevp1 = bevp2;
2392     bevp2++;
2393   }
2394
2395   if (bl->poly != -1) { /* check for cyclic */
2396
2397     /* Need to correct for the start/end points not matching
2398      * do this by calculating the tilt angle difference, then apply
2399      * the rotation gradually over the entire curve
2400      *
2401      * note that the split is between last and second last, rather than first/last as youd expect.
2402      *
2403      * real order is like this
2404      * 0,1,2,3,4 --> 1,2,3,4,0
2405      *
2406      * this is why we compare last with second last
2407      * */
2408     float vec_1[3] = {0, 1, 0}, vec_2[3] = {0, 1, 0}, angle, ang_fac, cross_tmp[3];
2409
2410     BevPoint *bevp_first;
2411     BevPoint *bevp_last;
2412
2413     bevp_first = bl->bevpoints;
2414     bevp_first += bl->nr - 1;
2415     bevp_last = bevp_first;
2416     bevp_last--;
2417
2418     /* quats and vec's are normalized, should not need to re-normalize */
2419     mul_qt_v3(bevp_first->quat, vec_1);
2420     mul_qt_v3(bevp_last->quat, vec_2);
2421     normalize_v3(vec_1);
2422     normalize_v3(vec_2);
2423
2424     /* align the vector, can avoid this and it looks 98% OK but
2425      * better to align the angle quat roll's before comparing */
2426     {
2427       cross_v3_v3v3(cross_tmp, bevp_last->dir, bevp_first->dir);
2428       angle = angle_normalized_v3v3(bevp_first->dir, bevp_last->dir);
2429       axis_angle_to_quat(q, cross_tmp, angle);
2430       mul_qt_v3(q, vec_2);
2431     }
2432
2433     angle = angle_normalized_v3v3(vec_1, vec_2);
2434
2435     /* flip rotation if needs be */
2436     cross_v3_v3v3(cross_tmp, vec_1, vec_2);
2437     normalize_v3(cross_tmp);
2438     if (angle_normalized_v3v3(bevp_first->dir, cross_tmp) < DEG2RADF(90.0f)) {
2439       angle = -angle;
2440     }
2441
2442     bevp2 = bl->bevpoints;
2443     bevp1 = bevp2 + (bl->nr - 1);
2444     bevp0 = bevp1 - 1;
2445
2446     nr = bl->nr;
2447     while (nr--) {
2448       ang_fac = angle * (1.0f - ((float)nr / bl->nr)); /* also works */
2449
2450       axis_angle_to_quat(q, bevp1->dir, ang_fac);
2451       mul_qt_qtqt(bevp1->quat, q, bevp1->quat);
2452
2453       bevp0 = bevp1;
2454       bevp1 = bevp2;
2455       bevp2++;
2456     }
2457   }
2458   else {
2459     /* Need to correct quat for the first/last point,
2460      * this is so because previously it was only calculated
2461      * using it's own direction, which might not correspond
2462      * the twist of neighbor point.
2463      */
2464     bevp1 = bl->bevpoints;
2465     bevp0 = bevp1 + 1;
2466     minimum_twist_between_two_points(bevp1, bevp0);
2467
2468     bevp2 = bl->bevpoints;
2469     bevp1 = bevp2 + (bl->nr - 1);
2470     bevp0 = bevp1 - 1;
2471     minimum_twist_between_two_points(bevp1, bevp0);
2472   }
2473 }
2474
2475 static void make_bevel_list_3D_tangent(BevList *bl)
2476 {
2477   BevPoint *bevp2, *bevp1, *bevp0; /* standard for all make_bevel_list_3D_* funcs */
2478   int nr;
2479
2480   float bevp0_tan[3];
2481
2482   bevel_list_calc_bisect(bl);
2483   bevel_list_flip_tangents(bl);
2484
2485   /* correct the tangents */
2486   bevp2 = bl->bevpoints;
2487   bevp1 = bevp2 + (bl->nr - 1);
2488   bevp0 = bevp1 - 1;
2489
2490   nr = bl->nr;
2491   while (nr--) {
2492     float cross_tmp[3];
2493     cross_v3_v3v3(cross_tmp, bevp1->tan, bevp1->dir);
2494     cross_v3_v3v3(bevp1->tan, cross_tmp, bevp1->dir);
2495     normalize_v3(bevp1->tan);
2496
2497     bevp0 = bevp1;
2498     bevp1 = bevp2;
2499     bevp2++;
2500   }
2501
2502   /* now for the real twist calc */
2503   bevp2 = bl->bevpoints;
2504   bevp1 = bevp2 + (bl->nr - 1);
2505   bevp0 = bevp1 - 1;
2506
2507   copy_v3_v3(bevp0_tan, bevp0->tan);
2508
2509   nr = bl->nr;
2510   while (nr--) {
2511     /* make perpendicular, modify tan in place, is ok */
2512     float cross_tmp[3];
2513     float zero[3] = {0, 0, 0};
2514
2515     cross_v3_v3v3(cross_tmp, bevp1->tan, bevp1->dir);
2516     normalize_v3(cross_tmp);
2517     tri_to_quat(bevp1->quat, zero, cross_tmp, bevp1->tan); /* XXX - could be faster */
2518
2519     /* bevp0 = bevp1; */ /* UNUSED */
2520     bevp1 = bevp2;
2521     bevp2++;
2522   }
2523 }
2524
2525 static void make_bevel_list_3D(BevList *bl, int smooth_iter, int twist_mode)
2526 {
2527   switch (twist_mode) {
2528     case CU_TWIST_TANGENT:
2529       make_bevel_list_3D_tangent(bl);
2530       break;
2531     case CU_TWIST_MINIMUM:
2532       make_bevel_list_3D_minimum_twist(bl);
2533       break;
2534     default: /* CU_TWIST_Z_UP default, pre 2.49c */
2535       make_bevel_list_3D_zup(bl);
2536       break;
2537   }
2538
2539   if (smooth_iter) {
2540     bevel_list_smooth(bl, smooth_iter);
2541   }
2542
2543   bevel_list_apply_tilt(bl);
2544 }
2545
2546 /* only for 2 points */
2547 static void make_bevel_list_segment_3D(BevList *bl)
2548 {
2549   float q[4];
2550
2551   BevPoint *bevp2 = bl->bevpoints;
2552   BevPoint *bevp1 = bevp2 + 1;
2553
2554   /* simple quat/dir */
2555   sub_v3_v3v3(bevp1->dir, bevp1->vec, bevp2->vec);
2556   normalize_v3(bevp1->dir);
2557
2558   vec_to_quat(bevp1->quat, bevp1->dir, 5, 1);
2559
2560   axis_angle_to_quat(q, bevp1->dir, bevp1->tilt);
2561   mul_qt_qtqt(bevp1->quat, q, bevp1->quat);
2562   normalize_qt(bevp1->quat);
2563   copy_v3_v3(bevp2->dir, bevp1->dir);
2564   copy_qt_qt(bevp2->quat, bevp1->quat);
2565 }
2566
2567 /* only for 2 points */
2568 static void make_bevel_list_segment_2D(BevList *bl)
2569 {
2570   BevPoint *bevp2 = bl->bevpoints;
2571   BevPoint *bevp1 = bevp2 + 1;
2572
2573   const float x1 = bevp1->vec[0] - bevp2->vec[0];
2574   const float y1 = bevp1->vec[1] - bevp2->vec[1];
2575
2576   calc_bevel_sin_cos(x1, y1, -x1, -y1, &(bevp1->sina), &(bevp1->cosa));
2577   bevp2->sina = bevp1->sina;
2578   bevp2->cosa = bevp1->cosa;
2579
2580   /* fill in dir & quat */
2581   make_bevel_list_segment_3D(bl);
2582 }
2583
2584 static void make_bevel_list_2D(BevList *bl)
2585 {
2586   /* note: bevp->dir and bevp->quat are not needed for beveling but are
2587    * used when making a path from a 2D curve, therefore they need to be set - Campbell */
2588
2589   BevPoint *bevp0, *bevp1, *bevp2;
2590   int nr;
2591
2592   if (bl->poly != -1) {
2593     bevp2 = bl->bevpoints;
2594     bevp1 = bevp2 + (bl->nr - 1);
2595     bevp0 = bevp1 - 1;
2596     nr = bl->nr;
2597   }
2598   else {
2599     bevp0 = bl->bevpoints;
2600     bevp1 = bevp0 + 1;
2601     bevp2 = bevp1 + 1;
2602
2603     nr = bl->nr - 2;
2604   }
2605
2606   while (nr--) {
2607     const float x1 = bevp1->vec[0] - bevp0->vec[0];
2608     const float x2 = bevp1->vec[0] - bevp2->vec[0];
2609     const float y1 = bevp1->vec[1] - bevp0->vec[1];
2610     const float y2 = bevp1->vec[1] - bevp2->vec[1];
2611
2612     calc_bevel_sin_cos(x1, y1, x2, y2, &(bevp1->sina), &(bevp1->cosa));
2613
2614     /* from: make_bevel_list_3D_zup, could call but avoid a second loop.
2615      * no need for tricky tilt calculation as with 3D curves */
2616     bisect_v3_v3v3v3(bevp1->dir, bevp0->vec, bevp1->vec, bevp2->vec);
2617     vec_to_quat(bevp1->quat, bevp1->dir, 5, 1);
2618     /* done with inline make_bevel_list_3D_zup */
2619
2620     bevp0 = bevp1;
2621     bevp1 = bevp2;
2622     bevp2++;
2623   }
2624
2625   /* correct non-cyclic cases */
2626   if (bl->poly == -1) {
2627     BevPoint *bevp;
2628     float angle;
2629
2630     /* first */
2631     bevp = bl->bevpoints;
2632     angle = atan2f(bevp->dir[0], bevp->dir[1]) - (float)M_PI_2;
2633     bevp->sina = sinf(angle);
2634     bevp->cosa = cosf(angle);
2635     vec_to_quat(bevp->quat, bevp->dir, 5, 1);
2636
2637     /* last */
2638     bevp = bl->bevpoints;
2639     bevp += (bl->nr - 1);
2640     angle = atan2f(bevp->dir[0], bevp->dir[1]) - (float)M_PI_2;
2641     bevp->sina = sinf(angle);
2642     bevp->cosa = cosf(angle);
2643     vec_to_quat(bevp->quat, bevp->dir, 5, 1);
2644   }
2645 }
2646
2647 static void bevlist_firstlast_direction_calc_from_bpoint(Nurb *nu, BevList *bl)
2648 {
2649   if (nu->pntsu > 1) {
2650     BPoint *first_bp = nu->bp, *last_bp = nu->bp + (nu->pntsu - 1);
2651     BevPoint *first_bevp, *last_bevp;
2652
2653     first_bevp = bl->bevpoints;
2654     last_bevp = first_bevp + (bl->nr - 1);
2655
2656     sub_v3_v3v3(first_bevp->dir, (first_bp + 1)->vec, first_bp->vec);
2657     normalize_v3(first_bevp->dir);
2658
2659     sub_v3_v3v3(last_bevp->dir, last_bp->vec, (last_bp - 1)->vec);
2660     normalize_v3(last_bevp->dir);
2661   }
2662 }
2663
2664 void BKE_curve_bevelList_free(ListBase *bev)
2665 {
2666   BevList *bl, *blnext;
2667   for (bl = bev->first; bl != NULL; bl = blnext) {
2668     blnext = bl->next;
2669     if (bl->seglen != NULL) {
2670       MEM_freeN(bl->seglen);
2671     }
2672     if (bl->segbevcount != NULL) {
2673       MEM_freeN(bl->segbevcount);
2674     }
2675     if (bl->bevpoints != NULL) {
2676       MEM_freeN(bl->bevpoints);
2677     }
2678     MEM_freeN(bl);
2679   }
2680
2681   BLI_listbase_clear(bev);
2682 }
2683
2684 void BKE_curve_bevelList_make(Object *ob, ListBase *nurbs, bool for_render)
2685 {
2686   /*
2687    * - convert all curves to polys, with indication of resol and flags for double-vertices
2688    * - possibly; do a smart vertice removal (in case Nurb)
2689    * - separate in individual blocks with BoundBox
2690    * - AutoHole detection
2691    */
2692
2693   /* this function needs an object, because of tflag and upflag */
2694   Curve *cu = ob->data;
2695   Nurb *nu;
2696   BezTriple *bezt, *prevbezt;
2697   BPoint *bp;
2698   BevList *bl, *blnew, *blnext;
2699   BevPoint *bevp2, *bevp1 = NULL, *bevp0;
2700   const float treshold = 0.00001f;
2701   float min, inp;
2702   float *seglen = NULL;
2703   struct BevelSort *sortdata, *sd, *sd1;
2704   int a, b, nr, poly, resolu = 0, len = 0, segcount;
2705   int *segbevcount;
2706   bool do_tilt, do_radius, do_weight;
2707   bool is_editmode = false;
2708   ListBase *bev;
2709
2710   /* segbevcount alsp requires seglen. */
2711   const bool need_seglen = ELEM(
2712                                cu->bevfac1_mapping, CU_BEVFAC_MAP_SEGMENT, CU_BEVFAC_MAP_SPLINE) ||
2713                            ELEM(cu->bevfac2_mapping, CU_BEVFAC_MAP_SEGMENT, CU_BEVFAC_MAP_SPLINE);
2714
2715   bev = &ob->runtime.curve_cache->bev;
2716
2717 #if 0
2718   /* do we need to calculate the radius for each point? */
2719   do_radius = (cu->bevobj || cu->taperobj || (cu->flag & CU_FRONT) || (cu->flag & CU_BACK)) ? 0 :
2720                                                                                               1;
2721 #endif
2722
2723   /* STEP 1: MAKE POLYS  */
2724
2725   BKE_curve_bevelList_free(&ob->runtime.curve_cache->bev);
2726   nu = nurbs->first;
2727   if (cu->editnurb && ob->type != OB_FONT) {
2728     is_editmode = 1;
2729   }
2730
2731   for (; nu; nu = nu->next) {
2732     if (nu->hide && is_editmode) {
2733       continue;
2734     }
2735
2736     /* check if we will calculate tilt data */
2737     do_tilt = CU_DO_TILT(cu, nu);
2738     do_radius = CU_DO_RADIUS(
2739         cu, nu); /* normal display uses the radius, better just to calculate them */
2740     do_weight = true;
2741
2742     /* check we are a single point? also check we are not a surface and that the orderu is sane,
2743      * enforced in the UI but can go wrong possibly */
2744     if (!BKE_nurb_check_valid_u(nu)) {
2745       bl = MEM_callocN(sizeof(BevList), "makeBevelList1");
2746       bl->bevpoints = MEM_calloc_arrayN(1, sizeof(BevPoint), "makeBevelPoints1");
2747       BLI_addtail(bev, bl);
2748       bl->nr = 0;
2749       bl->charidx = nu->charidx;
2750     }
2751     else {
2752       BevPoint *bevp;
2753
2754       if (for_render && cu->resolu_ren != 0) {
2755         resolu = cu->resolu_ren;
2756       }
2757       else {
2758         resolu = nu->resolu;
2759       }
2760
2761       segcount = SEGMENTSU(nu);
2762
2763       if (nu->type == CU_POLY) {
2764         len = nu->pntsu;
2765         bl = MEM_callocN(sizeof(BevList), "makeBevelList2");
2766         bl->bevpoints = MEM_calloc_arrayN(len, sizeof(BevPoint), "makeBevelPoints2");
2767         if (need_seglen && (nu->flagu & CU_NURB_CYCLIC) == 0) {
2768           bl->seglen = MEM_malloc_arrayN(segcount, sizeof(float), "makeBevelList2_seglen");
2769           bl->segbevcount = MEM_malloc_arrayN(segcount, sizeof(int), "makeBevelList2_segbevcount");
2770         }
2771         BLI_addtail(bev, bl);
2772
2773         bl->poly = (nu->flagu & CU_NURB_CYCLIC) ? 0 : -1;
2774         bl->nr = len;
2775         bl->dupe_nr = 0;
2776         bl->charidx = nu->charidx;
2777         bevp = bl->bevpoints;
2778         bevp->offset = 0;
2779         bp = nu->bp;
2780         seglen = bl->seglen;
2781         segbevcount = bl->segbevcount;
2782
2783         while (len--) {
2784           copy_v3_v3(bevp->vec, bp->vec);
2785           bevp->tilt = bp->tilt;
2786           bevp->radius = bp->radius;
2787           bevp->weight = bp->weight;
2788           bevp->split_tag = true;
2789           bp++;
2790           if (seglen != NULL && len != 0) {
2791             *seglen = len_v3v3(bevp->vec, bp->vec);
2792             bevp++;
2793             bevp->offset = *seglen;
2794             if (*seglen > treshold) {
2795               *segbevcount = 1;
2796             }
2797             else {
2798               *segbevcount = 0;
2799             }
2800             seglen++;
2801             segbevcount++;
2802           }
2803           else {
2804             bevp++;
2805           }
2806         }
2807
2808         if ((nu->flagu & CU_NURB_CYCLIC) == 0) {
2809           bevlist_firstlast_direction_calc_from_bpoint(nu, bl);
2810         }
2811       }
2812       else if (nu->type == CU_BEZIER) {
2813         /* in case last point is not cyclic */
2814         len = segcount * resolu + 1;
2815
2816         bl = MEM_callocN(sizeof(BevList), "makeBevelBPoints");
2817         bl->bevpoints = MEM_calloc_arrayN(len, sizeof(BevPoint), "makeBevelBPointsPoints");
2818         if (need_seglen && (nu->flagu & CU_NURB_CYCLIC) == 0) {
2819           bl->seglen = MEM_malloc_arrayN(segcount, sizeof(float), "makeBevelBPoints_seglen");
2820           bl->segbevcount = MEM_malloc_arrayN(
2821               segcount, sizeof(int), "makeBevelBPoints_segbevcount");
2822         }
2823         BLI_addtail(bev, bl);
2824
2825         bl->poly = (nu->flagu & CU_NURB_CYCLIC) ? 0 : -1;
2826         bl->charidx = nu->charidx;
2827
2828         bevp = bl->bevpoints;
2829         seglen = bl->seglen;
2830         segbevcount = bl->segbevcount;
2831
2832         bevp->offset = 0;
2833         if (seglen != NULL) {
2834           *seglen = 0;
2835           *segbevcount = 0;
2836         }
2837
2838         a = nu->pntsu - 1;
2839         bezt = nu->bezt;
2840         if (nu->flagu & CU_NURB_CYCLIC) {
2841           a++;
2842           prevbezt = nu->bezt + (nu->pntsu - 1);
2843         }
2844         else {
2845           prevbezt = bezt;
2846           bezt++;
2847         }
2848
2849         sub_v3_v3v3(bevp->dir, prevbezt->vec[2], prevbezt->vec[1]);
2850         normalize_v3(bevp->dir);
2851
2852         BLI_assert(segcount >= a);
2853
2854         while (a--) {
2855           if (prevbezt->h2 == HD_VECT && bezt->h1 == HD_VECT) {
2856
2857             copy_v3_v3(bevp->vec, prevbezt->vec[1]);
2858             bevp->tilt = prevbezt->tilt;
2859             bevp->radius = prevbezt->radius;
2860             bevp->weight = prevbezt->weight;
2861             bevp->split_tag = true;
2862             bevp->dupe_tag = false;
2863             bevp++;
2864             bl->nr++;
2865             bl->dupe_nr = 1;
2866             if (seglen != NULL) {
2867               *seglen = len_v3v3(prevbezt->vec[1], bezt->vec[1]);
2868               bevp->offset = *seglen;
2869               seglen++;
2870               /* match segbevcount to the cleaned up bevel lists (see STEP 2) */
2871               if (bevp->offset > treshold) {
2872                 *segbevcount = 1;
2873               }
2874               segbevcount++;
2875             }
2876           }
2877           else {
2878             /* always do all three, to prevent data hanging around */
2879             int j;
2880
2881             /* BevPoint must stay aligned to 4 so sizeof(BevPoint)/sizeof(float) works */
2882             for (j = 0; j < 3; j++) {
2883               BKE_curve_forward_diff_bezier(prevbezt->vec[1][j],
2884                                             prevbezt->vec[2][j],
2885                                             bezt->vec[0][j],
2886                                             bezt->vec[1][j],
2887                                             &(bevp->vec[j]),
2888                                             resolu,
2889                                             sizeof(BevPoint));
2890             }
2891
2892             /* if both arrays are NULL do nothiong */
2893             tilt_bezpart(prevbezt,
2894                          bezt,
2895                          nu,
2896                          do_tilt ? &bevp->tilt : NULL,
2897                          do_radius ? &bevp->radius : NULL,
2898                          do_weight ? &bevp->weight : NULL,
2899                          resolu,
2900                          sizeof(BevPoint));
2901
2902             if (cu->twist_mode == CU_TWIST_TANGENT) {
2903               forward_diff_bezier_cotangent(prevbezt->vec[1],
2904                                             prevbezt->vec[2],
2905                                             bezt->vec[0],
2906                                             bezt->vec[1],
2907                                             bevp->tan,
2908                                             resolu,
2909                                             sizeof(BevPoint));
2910             }
2911
2912             /* indicate with handlecodes double points */
2913             if (prevbezt->h1 == prevbezt->h2) {
2914               if (prevbezt->h1 == 0 || prevbezt->h1 == HD_VECT) {
2915                 bevp->split_tag = true;
2916               }
2917             }
2918             else {
2919               if (prevbezt->h1 == 0 || prevbezt->h1 == HD_VECT) {
2920                 bevp->split_tag = true;
2921               }
2922               else if (prevbezt->h2 == 0 || prevbezt->h2 == HD_VECT) {
2923                 bevp->split_tag = true;
2924               }
2925             }
2926
2927             /* seglen */
2928             if (seglen != NULL) {
2929               *seglen = 0;
2930               *segbevcount = 0;
2931               for (j = 0; j < resolu; j++) {
2932                 bevp0 = bevp;
2933                 bevp++;
2934                 bevp->offset = len_v3v3(bevp0->vec, bevp->vec);
2935                 /* match seglen and segbevcount to the cleaned up bevel lists (see STEP 2) */
2936                 if (bevp->offset > treshold) {
2937                   *seglen += bevp->offset;
2938                   *segbevcount += 1;
2939                 }
2940               }
2941               seglen++;
2942               segbevcount++;
2943             }
2944             else {
2945               bevp += resolu;
2946             }
2947             bl->nr += resolu;
2948           }
2949           prevbezt = bezt;
2950           bezt++;
2951         }
2952
2953         if ((nu->flagu & CU_NURB_CYCLIC) == 0) { /* not cyclic: endpoint */
2954           copy_v3_v3(bevp->vec, prevbezt->vec[1]);
2955           bevp->tilt = prevbezt->tilt;
2956           bevp->radius = prevbezt->radius;
2957           bevp->weight = prevbezt->weight;
2958
2959           sub_v3_v3v3(bevp->dir, prevbezt->vec[1], prevbezt->vec[0]);
2960           normalize_v3(bevp->dir);
2961
2962           bl->nr++;
2963         }
2964       }
2965       else if (nu->type == CU_NURBS) {
2966         if (nu->pntsv == 1) {
2967           len = (resolu * segcount);
2968
2969           bl = MEM_callocN(sizeof(BevList), "makeBevelList3");
2970           bl->bevpoints = MEM_calloc_arrayN(len, sizeof(BevPoint), "makeBevelPoints3");
2971           if (need_seglen && (nu->flagu & CU_NURB_CYCLIC) == 0) {
2972             bl->seglen = MEM_malloc_arrayN(segcount, sizeof(float), "makeBevelList3_seglen");
2973             bl->segbevcount = MEM_malloc_arrayN(
2974                 segcount, sizeof(int), "makeBevelList3_segbevcount");
2975           }
2976           BLI_addtail(bev, bl);
2977           bl->nr = len;
2978           bl->dupe_nr = 0;
2979           bl->poly = (nu->flagu & CU_NURB_CYCLIC) ? 0 : -1;
2980           bl->charidx = nu->charidx;
2981
2982           bevp = bl->bevpoints;
2983           seglen = bl->seglen;
2984           segbevcount = bl->segbevcount;
2985
2986           BKE_nurb_makeCurve(nu,
2987                              &bevp->vec[0],
2988                              do_tilt ? &bevp->tilt : NULL,
2989                              do_radius ? &bevp->radius : NULL,
2990                              do_weight ? &bevp->weight : NULL,
2991                              resolu,
2992                              sizeof(BevPoint));
2993
2994           /* match seglen and segbevcount to the cleaned up bevel lists (see STEP 2) */
2995           if (seglen != NULL) {
2996             nr = segcount;
2997             bevp0 = bevp;
2998             bevp++;
2999             while (nr) {
3000               int j;
3001               *seglen = 0;
3002               *segbevcount = 0;
3003               /* We keep last bevel segment zero-length. */
3004               for (j = 0; j < ((nr == 1) ? (resolu - 1) : resolu); j++) {
3005                 bevp->offset = len_v3v3(bevp0->vec, bevp->vec);
3006                 if (bevp->offset > treshold) {
3007                   *seglen += bevp->offset;
3008                   *segbevcount += 1;
3009                 }
3010                 bevp0 = bevp;
3011                 bevp++;
3012               }
3013               seglen++;
3014               segbevcount++;
3015               nr--;
3016             }
3017           }
3018
3019           if ((nu->flagu & CU_NURB_CYCLIC) == 0) {
3020             bevlist_firstlast_direction_calc_from_bpoint(nu, bl);
3021           }
3022         }
3023       }
3024     }
3025   }
3026
3027   /* STEP 2: DOUBLE POINTS AND AUTOMATIC RESOLUTION, REDUCE DATABLOCKS */
3028   bl = bev->first;
3029   while (bl) {
3030     if (bl->nr) { /* null bevel items come from single points */
3031       bool is_cyclic = bl->poly != -1;
3032       nr = bl->nr;
3033       if (is_cyclic) {
3034         bevp1 = bl->bevpoints;
3035         bevp0 = bevp1 + (nr - 1);
3036       }
3037       else {
3038         bevp0 = bl->bevpoints;
3039         bevp0->offset = 0;
3040         bevp1 = bevp0 + 1;
3041       }
3042       nr--;
3043       while (nr--) {
3044         if (seglen != NULL) {
3045           if (fabsf(bevp1->offset) < treshold) {
3046             bevp0->dupe_tag = true;
3047             bl->dupe_nr++;
3048           }
3049         }
3050         else {
3051           if (fabsf(bevp0->vec[0] - bevp1->vec[0]) < 0.00001f) {
3052             if (fabsf(bevp0->vec[1] - bevp1->vec[1]) < 0.00001f) {
3053               if (fabsf(bevp0->vec[2] - bevp1->vec[2]) < 0.00001f) {
3054                 bevp0->dupe_tag = true;
3055                 bl->dupe_nr++;
3056               }
3057             }
3058           }
3059         }
3060         bevp0 = bevp1;
3061         bevp1++;
3062       }
3063     }
3064     bl = bl->next;
3065   }
3066   bl = bev->first;
3067   while (bl) {
3068     blnext = bl->next;
3069     if (bl->nr && bl->dupe_nr) {
3070       nr = bl->nr - bl->dupe_nr + 1; /* +1 because vectorbezier sets flag too */
3071       blnew = MEM_callocN(sizeof(BevList), "makeBevelList4");
3072       memcpy(blnew, bl, sizeof(BevList));
3073       blnew->bevpoints = MEM_calloc_arrayN(nr, sizeof(BevPoint), "makeBevelPoints4");
3074       if (!blnew->bevpoints) {
3075         MEM_freeN(blnew);
3076         break;
3077       }
3078       blnew->segbevcount = bl->segbevcount;
3079       blnew->seglen = bl->seglen;
3080       blnew->nr = 0;
3081       BLI_remlink(bev, bl);
3082       BLI_insertlinkbefore(bev, blnext, blnew); /* to make sure bevlijst is tuned with nurblist */
3083       bevp0 = bl->bevpoints;
3084       bevp1 = blnew->bevpoints;
3085       nr = bl->nr;
3086       while (nr--) {
3087         if (bevp0->dupe_tag == 0) {
3088           memcpy(bevp1, bevp0, sizeof(BevPoint));
3089           bevp1++;
3090           blnew->nr++;
3091         }
3092         bevp0++;
3093       }
3094       if (bl->bevpoints != NULL) {
3095         MEM_freeN(bl->bevpoints);
3096       }
3097       MEM_freeN(bl);
3098       blnew->dupe_nr = 0;
3099     }
3100     bl = blnext;
3101   }
3102
3103   /* STEP 3: POLYS COUNT AND AUTOHOLE */
3104   bl = bev->first;
3105   poly = 0;
3106   while (bl) {
3107     if (bl->nr && bl->poly >= 0) {
3108       poly++;
3109       bl->poly = poly;
3110       bl->hole = 0;
3111     }
3112     bl = bl->next;
3113   }
3114
3115   /* find extreme left points, also test (turning) direction */
3116   if (poly > 0) {
3117     sd = sortdata = MEM_malloc_arrayN(poly, sizeof(struct BevelSort), "makeBevelList5");
3118     bl = bev->first;
3119     while (bl) {
3120       if (bl->poly > 0) {
3121         BevPoint *bevp;
3122
3123         min = 300000.0;
3124         bevp = bl->bevpoints;
3125         nr = bl->nr;
3126         while (nr--) {
3127           if (min > bevp->vec[0]) {
3128             min = bevp->vec[0];
3129             bevp1 = bevp;
3130           }
3131           bevp++;
3132         }
3133         sd->bl = bl;
3134         sd->left = min;
3135
3136         bevp = bl->bevpoints;
3137         if (bevp1 == bevp) {
3138           bevp0 = bevp + (bl->nr - 1);
3139         }
3140         else {
3141           bevp0 = bevp1 - 1;
3142         }
3143         bevp = bevp + (bl->nr - 1);
3144         if (bevp1 == bevp) {
3145           bevp2 = bl->bevpoints;
3146         }
3147         else {
3148           bevp2 = bevp1 + 1;
3149         }
3150
3151         inp = ((bevp1->vec[0] - bevp0->vec[0]) * (bevp0->vec[1] - bevp2->vec[1]) +
3152                (bevp0->vec[1] - bevp1->vec[1]) * (bevp0->vec[0] - bevp2->vec[0]));
3153
3154         if (inp > 0.0f) {
3155           sd->dir = 1;
3156         }
3157         else {
3158           sd->dir = 0;
3159         }
3160
3161         sd++;
3162       }
3163
3164       bl = bl->next;
3165     }
3166     qsort(sortdata, poly, sizeof(struct BevelSort), vergxcobev);
3167
3168     sd = sortdata + 1;
3169     for (a = 1; a < poly; a++, sd++) {
3170       bl = sd->bl; /* is bl a hole? */
3171       sd1 = sortdata + (a - 1);
3172       for (b = a - 1; b >= 0; b--, sd1--) {    /* all polys to the left */
3173         if (sd1->bl->charidx == bl->charidx) { /* for text, only check matching char */
3174           if (bevelinside(sd1->bl, bl)) {
3175             bl->hole = 1 - sd1->bl->hole;
3176             break;
3177           }
3178         }
3179       }
3180     }
3181
3182     /* turning direction */
3183     if ((cu->flag & CU_3D) == 0) {
3184       sd = sortdata;
3185       for (a = 0; a < poly; a++, sd++) {
3186         if (sd->bl->hole == sd->dir) {
3187           bl = sd->bl;
3188           bevp1 = bl->bevpoints;
3189           bevp2 = bevp1 + (bl->nr - 1);
3190           nr = bl->nr / 2;
3191           while (nr--) {
3192             SWAP(BevPoint, *bevp1, *bevp2);
3193             bevp1++;
3194             bevp2--;
3195           }
3196         }
3197       }
3198     }
3199     MEM_freeN(sortdata);
3200   }
3201
3202   /* STEP 4: 2D-COSINES or 3D ORIENTATION */
3203   if ((cu->flag & CU_3D) == 0) {
3204     /* 2D Curves */
3205     for (bl = bev->first; bl; bl = bl->next) {
3206       if (bl->nr < 2) {
3207         BevPoint *bevp = bl->bevpoints;
3208         unit_qt(bevp->quat);
3209       }
3210       else if (bl->nr == 2) { /* 2 pnt, treat separate */
3211         make_bevel_list_segment_2D(bl);
3212       }
3213       else {
3214         make_bevel_list_2D(bl);
3215       }
3216     }
3217   }
3218   else {
3219     /* 3D Curves */
3220     for (bl = bev->first; bl; bl = bl->next) {
3221       if (bl->nr < 2) {
3222         BevPoint *bevp = bl->bevpoints;
3223         unit_qt(bevp->quat);
3224       }
3225       else if (bl->nr == 2) { /* 2 pnt, treat separate */
3226         make_bevel_list_segment_3D(bl);
3227       }
3228       else {
3229         make_bevel_list_3D(bl, (int)(resolu * cu->twist_smooth), cu->twist_mode);
3230       }
3231     }
3232   }
3233 }
3234
3235 /* ****************** HANDLES ************** */
3236
3237 static void calchandleNurb_intern(BezTriple *bezt,
3238                                   const BezTriple *prev,
3239                                   const BezTriple *next,
3240                                   bool is_fcurve,
3241                                   bool skip_align,
3242                                   char fcurve_smoothing)
3243 {
3244   /* defines to avoid confusion */
3245 #define p2_h1 ((p2)-3)
3246 #define p2_h2 ((p2) + 3)
3247
3248   const float *p1, *p3;
3249   float *p2;
3250   float pt[3];
3251   float dvec_a[3], dvec_b[3];
3252   float len, len_a, len_b;
3253   float len_ratio;
3254   const float eps = 1e-5;
3255
3256   /* assume normal handle until we check */
3257   bezt->f5 = HD_AUTOTYPE_NORMAL;
3258
3259   if (bezt->h1 == 0 && bezt->h2 == 0) {
3260     return;
3261   }
3262
3263   p2 = bezt->vec[1];
3264
3265   if (prev == NULL) {
3266     p3 = next->vec[1];
3267     pt[0] = 2.0f * p2[0] - p3[0];
3268     pt[1] = 2.0f * p2[1] - p3[1];
3269     pt[2] = 2.0f * p2[2] - p3[2];
3270     p1 = pt;
3271   }
3272   else {
3273     p1 = prev->vec[1];
3274   }
3275
3276   if (next == NULL) {
3277     pt[0] = 2.0f * p2[0] - p1[0];
3278     pt[1] = 2.0f * p2[1] - p1[1];
3279     pt[2] = 2.0f * p2[2] - p1[2];
3280     p3 = pt;
3281   }
3282   else {
3283     p3 = next->vec[1];
3284   }
3285
3286   sub_v3_v3v3(dvec_a, p2, p1);
3287   sub_v3_v3v3(dvec_b, p3, p2);
3288
3289   if (is_fcurve) {
3290     len_a = dvec_a[0];
3291     len_b = dvec_b[0];
3292   }
3293   else {
3294     len_a = len_v3(dvec_a);
3295     len_b = len_v3(dvec_b);
3296   }
3297
3298   if (len_a == 0.0f) {
3299     len_a = 1.0f;
3300   }
3301   if (len_b == 0.0f) {
3302     len_b = 1.0f;
3303   }
3304
3305   len_ratio = len_a / len_b;
3306
3307   if (ELEM(bezt->h1, HD_AUTO, HD_AUTO_ANIM) || ELEM(bezt->h2, HD_AUTO, HD_AUTO_ANIM)) { /* auto */
3308     float tvec[3];
3309     tvec[0] = dvec_b[0] / len_b + dvec_a[0] / len_a;
3310     tvec[1] = dvec_b[1] / len_b + dvec_a[1] / len_a;
3311     tvec[2] = dvec_b[2] / len_b + dvec_a[2] / len_a;
3312
3313     if (is_fcurve) {
3314       if (fcurve_smoothing != FCURVE_SMOOTH_NONE) {
3315         /* force the horizontal handle size to be 1/3 of the key interval so that
3316          * the X component of the parametric bezier curve is a linear spline */
3317         len = 6.0f / 2.5614f;
3318       }
3319       else {
3320         len = tvec[0];
3321       }
3322     }
3323     else {
3324       len = len_v3(tvec);
3325     }
3326     len *= 2.5614f;
3327
3328     if (len != 0.0f) {
3329       /* only for fcurves */
3330       bool leftviolate = false, rightviolate = false;
3331
3332       if (!is_fcurve || fcurve_smoothing == FCURVE_SMOOTH_NONE) {
3333         if (len_a > 5.0f * len_b) {
3334           len_a = 5.0f * len_b;
3335         }
3336         if (len_b > 5.0f * len_a) {
3337           len_b = 5.0f * len_a;
3338         }
3339       }
3340
3341       if (ELEM(bezt->h1, HD_AUTO, HD_AUTO_ANIM)) {
3342         len_a /= len;
3343         madd_v3_v3v3fl(p2_h1, p2, tvec, -len_a);
3344
3345         if ((bezt->h1 == HD_AUTO_ANIM) && next && prev) { /* keep horizontal if extrema */
3346           float ydiff1 = prev->vec[1][1] - bezt->vec[1][1];
3347           float ydiff2 = next->vec[1][1] - bezt->vec[1][1];
3348           if ((ydiff1 <= 0.0f && ydiff2 <= 0.0f) || (ydiff1 >= 0.0f && ydiff2 >= 0.0f)) {
3349             bezt->vec[0][1] = bezt->vec[1][1];
3350             bezt->f5 = HD_AUTOTYPE_SPECIAL;
3351           }
3352           else { /* handles should not be beyond y coord of two others */
3353             if (ydiff1 <= 0.0f) {
3354               if (prev->vec[1][1] > bezt->vec[0][1]) {
3355                 bezt->vec[0][1] = prev->vec[1][1];
3356                 leftviolate = 1;
3357               }
3358             }
3359             else {
3360               if (prev->vec[1][1] < bezt->vec[0][1]) {
3361                 bezt->vec[0][1] = prev->vec[1][1];
3362                 leftviolate = 1;
3363               }
3364             }
3365           }
3366         }
3367       }
3368       if (ELEM(bezt->h2, HD_AUTO, HD_AUTO_ANIM)) {
3369         len_b /= len;
3370         madd_v3_v3v3fl(p2_h2, p2, tvec, len_b);
3371
3372         if ((bezt->h2 == HD_AUTO_ANIM) && next && prev) { /* keep horizontal if extrema */
3373           float ydiff1 = prev->vec[1][1] - bezt->vec[1][1];
3374           float ydiff2 = next->vec[1][1] - bezt->vec[1][1];
3375           if ((ydiff1 <= 0.0f && ydiff2 <= 0.0f) || (ydiff1 >= 0.0f && ydiff2 >= 0.0f)) {
3376             bezt->vec[2][1] = bezt->vec[1][1];
3377             bezt->f5 = HD_AUTOTYPE_SPECIAL;
3378           }
3379           else { /* handles should not be beyond y coord of two others */
3380             if (ydiff1 <= 0.0f) {
3381               if (next->vec[1][1] < bezt->vec[2][1]) {
3382                 bezt->vec[2][1] = next->vec[1][1];
3383                 rightviolate = 1;
3384               }
3385             }
3386             else {
3387               if (next->vec[1][1] > bezt->vec[2][1]) {
3388                 bezt->vec[2][1] = next->vec[1][1];
3389                 rightviolate = 1;
3390               }
3391             }
3392           }
3393         }
3394       }
3395       if (leftviolate || rightviolate) { /* align left handle */
3396         BLI_assert(is_fcurve);
3397         /* simple 2d calculation */
3398         float h1_x = p2_h1[0] - p2[0];
3399         float h2_x = p2[0] - p2_h2[0];
3400
3401         if (leftviolate) {
3402           p2_h2[1] = p2[1] + ((p2[1] - p2_h1[1]) / h1_x) * h2_x;
3403         }
3404         else {
3405           p2_h1[1] = p2[1] + ((p2[1] - p2_h2[1]) / h2_x) * h1_x;
3406         }
3407       }
3408     }
3409   }
3410
3411   if (bezt->h1 == HD_VECT) { /* vector */
3412     madd_v3_v3v3fl(p2_h1, p2, dvec_a, -1.0f / 3.0f);
3413   }
3414   if (bezt->h2 == HD_VECT) {
3415     madd_v3_v3v3fl(p2_h2, p2, dvec_b, 1.0f / 3.0f);
3416   }
3417
3418   if (skip_align ||
3419       /* when one handle is free, alignming makes no sense, see: T35952 */
3420       (ELEM(HD_FREE, bezt->h1, bezt->h2)) ||
3421       /* also when no handles are aligned, skip this step */
3422       (!ELEM(HD_ALIGN, bezt->h1, bezt->h2) && !ELEM(HD_ALIGN_DOUBLESIDE, bezt->h1, bezt->h2))) {
3423     /* handles need to be updated during animation and applying stuff like hooks,
3424      * but in such situations it's quite difficult to distinguish in which order
3425      * align handles should be aligned so skip them for now */
3426     return;
3427   }
3428
3429   len_a = len_v3v3(p2, p2_h1);
3430   len_b = len_v3v3(p2, p2_h2);
3431
3432   if (len_a == 0.0f) {
3433     len_a = 1.0f;
3434   }
3435   if (len_b == 0.0f) {
3436     len_b = 1.0f;
3437   }
3438
3439   len_ratio = len_a / len_b;
3440
3441   if (bezt->f1 & SELECT) {                               /* order of calculation */
3442     if (ELEM(bezt->h2, HD_ALIGN, HD_ALIGN_DOUBLESIDE)) { /* aligned */
3443       if (len_a > eps) {
3444         len = 1.0f / len_ratio;
3445         p2_h2[0] = p2[0] + len * (p2[0] - p2_h1[0]);
3446         p2_h2[1] = p2[1] + len * (p2[1] - p2_h1[1]);
3447         p2_h2[2] = p2[2] + len * (p2[2] - p2_h1[2]);
3448       }
3449     }
3450     if (ELEM(bezt->h1, HD_ALIGN, HD_ALIGN_DOUBLESIDE)) {
3451       if (len_b > eps) {
3452         len = len_ratio;
3453         p2_h1[0] = p2[0] + len * (p2[0] - p2_h2[0]);
3454         p2_h1[1] = p2[1] + len * (p2[1] - p2_h2[1]);
3455         p2_h1[2] = p2[2] + len * (p2[2] - p2_h2[2]);
3456       }
3457     }
3458   }
3459   else {
3460     if (ELEM(bezt->h1, HD_ALIGN, HD_ALIGN_DOUBLESIDE)) {
3461       if (len_b > eps) {
3462         len = len_ratio;
3463         p2_h1[0] = p2[0] + len * (p2[0] - p2_h2[0]);
3464         p2_h1[1] = p2[1] + len * (p2[1] - p2_h2[1]);
3465         p2_h1[2] = p2[2] + len * (p2[2] - p2_h2[2]);
3466       }
3467     }
3468     if (ELEM(bezt->h2, HD_ALIGN, HD_ALIGN_DOUBLESIDE)) { /* aligned */
3469       if (len_a > eps) {
3470         len = 1.0f / len_ratio;
3471         p2_h2[0] = p2[0] + len * (p2[0] - p2_h1[0]);
3472         p2_h2[1] = p2[1] + len * (p2[1] - p2_h1[1]);
3473         p2_h2[2] = p2[2] + len * (p2[2] - p2_h1[2]);
3474       }
3475     }
3476   }
3477
3478 #undef p2_h1
3479 #undef p2_h2
3480 }
3481
3482 static void calchandlesNurb_intern(Nurb *nu, bool skip_align)
3483 {
3484   BezTriple *bezt, *prev, *next;
3485   int a;
3486
3487   if (nu->type != CU_BEZIER) {
3488     return;
3489   }
3490   if (nu->pntsu < 2) {
3491     return;
3492   }
3493
3494   a = nu->pntsu;
3495   bezt = nu->bezt;
3496   if (nu->flagu & CU_NURB_CYCLIC) {
3497     prev = bezt + (a - 1);
3498   }
3499   else {
3500     prev = NULL;
3501   }
3502   next = bezt + 1;
3503
3504   while (a--) {
3505     calchandleNurb_intern(bezt, prev, next, 0, skip_align, 0);
3506     prev = bezt;
3507     if (a == 1) {
3508       if (nu->flagu & CU_NURB_CYCLIC) {
3509         next = nu->bezt;
3510       }
3511       else {
3512         next = NULL;
3513       }
3514     }
3515     else {
3516       next++;
3517     }
3518
3519     bezt++;
3520   }
3521 }
3522
3523 /* A utility function for allocating a number of arrays of the same length
3524  * with easy error checking and deallocation, and an easy way to add or remove
3525  * arrays that are processed in this way when changing code.
3526  *
3527  * floats, chars: NULL-terminated arrays of pointers to array pointers that need to be allocated.
3528  *
3529  * Returns: pointer to the buffer that contains all of the arrays.
3530  */
3531 static void *allocate_arrays(int count, float ***floats, char ***chars, const char *name)
3532 {
3533   size_t num_floats = 0, num_chars = 0;
3534
3535   while (floats && floats[num_floats]) {
3536     num_floats++;
3537   }
3538
3539   while (chars && chars[num_chars]) {
3540     num_chars++;
3541   }
3542
3543   void *buffer = (float *)MEM_malloc_arrayN(count, (sizeof(float) * num_floats + num_chars), name);
3544
3545   if (!buffer) {
3546     return NULL;
3547   }
3548
3549   float *fptr = buffer;
3550
3551   for (int i = 0; i < num_floats; i++, fptr += count) {
3552     *floats[i] = fptr;
3553   }
3554
3555   char *cptr = (char *)fptr;
3556
3557   for (int i = 0; i < num_chars; i++, cptr += count) {
3558     *chars[i] = cptr;
3559   }
3560
3561   return buffer;
3562 }
3563
3564 static void free_arrays(void *buffer)
3565 {
3566   MEM_freeN(buffer);
3567 }
3568
3569 /* computes in which direction to change h[i] to satisfy conditions better */
3570 static float bezier_relax_direction(
3571     float *a, float *b, float *c, float *d, float *h, int i, int count)
3572 {
3573   /* current deviation between sides of the equation */
3574   float state = a[i] * h[(i + count - 1) % count] + b[i] * h[i] + c[i] * h[(i + 1) % count] - d[i];
3575
3576   /* only the sign is meaningful */
3577   return -state * b[i];
3578 }
3579
3580 static void bezier_lock_unknown(float *a, float *b, float *c, float *d, int i, float value)
3581 {
3582   a[i] = c[i] = 0.0f;
3583   b[i] = 1.0f;
3584   d[i] = value;
3585 }
3586
3587 static void bezier_restore_equation(
3588     float *a, float *b, float *c, float *d, float *a0, float *b0, float *c0, float *d0, int i)
3589 {
3590   a[i] = a0[i];
3591   b[i] = b0[i];
3592   c[i] = c0[i];
3593   d[i] = d0[i];
3594 }
3595
3596 static bool tridiagonal_solve_with_limits(
3597     float *a, float *b, float *c, float *d, float *h, float *hmin, float *hmax, int solve_count)
3598 {
3599   float *a0, *b0, *c0, *d0;
3600   float **arrays[] = {&a0, &b0, &c0, &d0, NULL};
3601   char *is_locked, *num_unlocks;
3602   char **flagarrays[] = {&is_locked, &num_unlocks, NULL};
3603
3604   void *tmps = allocate_arrays(solve_count, arrays, flagarrays, "tridiagonal_solve_with_limits");
3605   if (!tmps) {
3606     return false;
3607   }
3608
3609   memcpy(a0, a, sizeof(float) * solve_count);
3610   memcpy(b0, b, sizeof(float) * solve_count);
3611   memcpy(c0, c, sizeof(float) * solve_count);
3612   memcpy(d0, d, sizeof(float) * solve_count);
3613
3614   memset(is_locked, 0, solve_count);
3615   memset(num_unlocks, 0, solve_count);
3616
3617   bool overshoot, unlocked;
3618
3619   do {
3620     if (!BLI_tridiagonal_solve_cyclic(a, b, c, d, h, solve_count)) {
3621       free_arrays(tmps);
3622       return false;
3623     }
3624
3625     /* first check if any handles overshoot the limits, and lock them */
3626     bool all = false, locked = false;
3627
3628     overshoot = unlocked = false;
3629
3630     do {
3631       for (int i = 0; i < solve_count; i++) {
3632         if (h[i] >= hmin[i] && h[i] <= hmax[i]) {
3633           continue;
3634         }
3635
3636         overshoot = true;
3637
3638         float target = h[i] > hmax[i] ? hmax[i] : hmin[i];
3639
3640         /* heuristically only lock handles that go in the right direction if there are such ones */
3641         if (target != 0.0f || all) {
3642           /* mark item locked */
3643           is_locked[i] = 1;
3644
3645           bezier_lock_unknown(a, b, c, d, i, target);
3646           locked = true;
3647         }
3648       }
3649
3650       all = true;
3651     } while (overshoot && !locked);
3652
3653     /* If no handles overshot and were locked,
3654      * see if it may be a good idea to unlock some handles. */
3655     if (!locked) {
3656       for (int i = 0; i < solve_count; i++) {
3657         // to definitely avoid infinite loops limit this to 2 times
3658         if (!is_locked[i] || num_unlocks[i] >= 2) {
3659           continue;
3660         }
3661
3662         /* if the handle wants to move in allowable direction, release it */
3663         float relax = bezier_relax_direction(a0, b0, c0, d0, h, i, solve_count);
3664
3665         if ((relax > 0 && h[i] < hmax[i]) || (relax < 0 && h[i] > hmin[i])) {
3666           bezier_restore_equation(a, b, c, d, a0, b0, c0, d0, i);
3667
3668           is_locked[i] = 0;
3669           num_unlocks[i]++;
3670           unlocked = true;
3671         }
3672       }
3673     }
3674   } while (overshoot || unlocked);
3675
3676   free_arrays(tmps);
3677   return true;
3678 }
3679
3680 /* Keep ascii art. */
3681 /* clang-format off */
3682 /*
3683  * This function computes the handles of a series of auto bezier points
3684  * on the basis of 'no acceleration discontinuities' at the points.
3685  * The first and last bezier points are considered 'fixed' (their handles are not touched)
3686  * The result is the smoothest possible trajectory going through intermediate points.
3687  * The difficulty is that the handles depends on their neighbors.
3688  *
3689  * The exact solution is found by solving a tridiagonal matrix equation formed
3690  * by the continuity and boundary conditions. Although theoretically handle position
3691  * is affected by all other points of the curve segment, in practice the influence
3692  * decreases exponentially with distance.
3693  *
3694  * Note: this algorithm assumes that the handle horizontal size if always 1/3 of the
3695  * of the interval to the next point. This rule ensures linear interpolation of time.
3696  *
3697  * ^ height (co 1)
3698  * |                                            yN
3699  * |                                   yN-1     |
3700  * |                      y2           |        |
3701  * |           y1         |            |        |
3702  * |    y0     |          |            |        |
3703  * |    |      |          |            |        |
3704  * |    |      |          |            |        |
3705  * |    |      |          |            |        |
3706  * |-------t1---------t2--------- ~ --------tN-------------------> time (co 0)
3707  * Mathematical basis:
3708  *
3709  * 1. Handle lengths on either side of each point are connected by a factor
3710  *    ensuring continuity of the first derivative:
3711  *
3712  *    l[i] = t[i+1]/t[i]
3713  *
3714  * 2. The tridiagonal system is formed by the following equation, which is derived
3715  *    by differentiating the bezier curve and specifies second derivative continuity
3716  *    at every point:
3717  *
3718  *    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]
3719  *
3720  * 3. If this point is adjacent to a manually set handle with X size not equal to 1/3
3721  *    of the horizontal interval, this equation becomes slightly more complex:
3722  *
3723  *    l[i]^2 * h[i-1] + (3*(1-R[i-1])*l[i] + 3*(1-L[i+1])) * h[i] + 1/l[i+1] * h[i+1] = (y[i]-y[i-1])*l[i]^2 + y[i+1]-y[i]
3724  *
3725  *    The difference between equations amounts to this, and it's obvious that when R[i-1]
3726  *    and L[i+1] are both 1/3, it becomes zero:
3727  *
3728  *    ( (1-3*R[i-1])*l[i] + (1-3*L[i+1]) ) * h[i]
3729  *
3730  * 4. The equations for zero acceleration border conditions are basically the above
3731  *    equation with parts omitted, so the handle size correction also applies.
3732  */
3733 /* clang-format on */
3734
3735 static void bezier_eq_continuous(
3736     float *a, float *b, float *c, float *d, float *dy, float *l, int i)
3737 {
3738   a[i] = l[i] * l[i];
3739   b[i] = 2.0f * (l[i] + 1);
3740   c[i] = 1.0f / l[i + 1];
3741   d[i] = dy[i] * l[i] * l[i] + dy[i + 1];
3742 }
3743
3744 static void bezier_eq_noaccel_right(
3745     float *a, float *b, float *c, float *d, float *dy, float *l, int i)
3746 {
3747   a[i] = 0.0f;
3748   b[i] = 2.0f;
3749   c[i] = 1.0f / l[i + 1];
3750   d[i] = dy[i + 1];
3751 }
3752
3753 static void bezier_eq_noaccel_left(
3754     float *a, float *b, float *c, float *d, float *dy, float *l, int i)
3755 {