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