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