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