More text object fancyness, and fixes:
[blender.git] / source / blender / blenkernel / intern / curve.c
1
2 /*  curve.c 
3  * 
4  *  
5  * $Id$
6  *
7  * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
8  *
9  * This program is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU General Public License
11  * as published by the Free Software Foundation; either version 2
12  * of the License, or (at your option) any later version. The Blender
13  * Foundation also sells licenses for use in proprietary software under
14  * the Blender License.  See http://www.blender.org/BL/ for information
15  * about this.
16  *
17  * This program is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with this program; if not, write to the Free Software Foundation,
24  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
25  *
26  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
27  * All rights reserved.
28  *
29  * The Original Code is: all of this file.
30  *
31  * Contributor(s): none yet.
32  *
33  * ***** END GPL/BL DUAL LICENSE BLOCK *****
34  */
35
36 #include <math.h>  // floor
37 #include <string.h>
38 #include <stdlib.h>  
39
40 #ifdef HAVE_CONFIG_H
41 #include <config.h>
42 #endif
43
44 #include "MEM_guardedalloc.h"
45 #include "BLI_blenlib.h"  
46 #include "BLI_arithb.h"  
47
48 #include "DNA_object_types.h"  
49 #include "DNA_curve_types.h"  
50 #include "DNA_material_types.h"  
51
52 /* for dereferencing pointers */
53 #include "DNA_ID.h"  
54 #include "DNA_vfont_types.h"  
55 #include "DNA_key_types.h"  
56 #include "DNA_ipo_types.h"  
57
58 #include "BKE_global.h" 
59 #include "BKE_main.h"  
60 #include "BKE_utildefines.h"  // VECCOPY
61 #include "BKE_object.h"  
62 #include "BKE_mesh.h" 
63 #include "BKE_curve.h"  
64 #include "BKE_displist.h"  
65 #include "BKE_ipo.h"  
66 #include "BKE_anim.h"  
67 #include "BKE_library.h"  
68 #include "BKE_key.h"  
69
70
71 /* globals */
72
73 extern ListBase editNurb;  /* editcurve.c */
74
75 /* local */
76 int cu_isectLL(float *v1, float *v2, float *v3, float *v4, 
77                            short cox, short coy, 
78                            float *labda, float *mu, float *vec);
79
80 void unlink_curve(Curve *cu)
81 {
82         int a;
83         
84         for(a=0; a<cu->totcol; a++) {
85                 if(cu->mat[a]) cu->mat[a]->id.us--;
86                 cu->mat[a]= 0;
87         }
88         if(cu->vfont) cu->vfont->id.us--; 
89         cu->vfont= 0;
90         if(cu->key) cu->key->id.us--;
91         cu->key= 0;
92         if(cu->ipo) cu->ipo->id.us--;
93         cu->ipo= 0;
94 }
95
96
97 /* niet curve zelf vrijgeven */
98 void free_curve(Curve *cu)
99 {
100
101         freeNurblist(&cu->nurb);
102         BLI_freelistN(&cu->bev);
103         freedisplist(&cu->disp);
104         
105         unlink_curve(cu);
106         
107         if(cu->mat) MEM_freeN(cu->mat);
108         if(cu->str) MEM_freeN(cu->str);
109         if(cu->strinfo) MEM_freeN(cu->strinfo);
110         if(cu->bb) MEM_freeN(cu->bb);
111         if(cu->path) free_path(cu->path);
112         if(cu->tb) MEM_freeN(cu->tb);
113 }
114
115 Curve *add_curve(int type)
116 {
117         Curve *cu;
118         char *str;
119         
120         if(type==OB_CURVE) str= "Curve";
121         else if(type==OB_SURF) str= "Surf";
122         else str= "Text";
123
124         cu= alloc_libblock(&G.main->curve, ID_CU, str);
125         
126         cu->size[0]= cu->size[1]= cu->size[2]= 1.0;
127         cu->flag= CU_FRONT+CU_BACK;
128         cu->pathlen= 100;
129         cu->resolu= cu->resolv= 6;
130         cu->width= 1.0;
131         cu->wordspace = 1.0;
132         cu->spacing= cu->linedist= 1.0;
133         cu->fsize= 1.0;
134         cu->ulheight = 0.05;    
135         cu->texflag= CU_AUTOSPACE;
136         
137         cu->bb= unit_boundbox();
138         
139         return cu;
140 }
141
142 Curve *copy_curve(Curve *cu)
143 {
144         Curve *cun;
145         int a;
146         
147         cun= copy_libblock(cu);
148         cun->nurb.first= cun->nurb.last= 0;
149         duplicateNurblist( &(cun->nurb), &(cu->nurb));
150
151         cun->mat= MEM_dupallocN(cu->mat);
152         for(a=0; a<cun->totcol; a++) {
153                 id_us_plus((ID *)cun->mat[a]);
154         }
155         
156         cun->str= MEM_dupallocN(cu->str);
157         cun->strinfo= MEM_dupallocN(cu->strinfo);       
158         cun->tb= MEM_dupallocN(cu->tb);
159         cun->bb= MEM_dupallocN(cu->bb);
160         
161         cun->key= copy_key(cu->key);
162         if(cun->key) cun->key->from= (ID *)cun;
163         
164         cun->disp.first= cun->disp.last= 0;
165         cun->bev.first= cun->bev.last= 0;
166         cun->path= 0;
167
168         /* single user ipo too */
169         if(cun->ipo) cun->ipo= copy_ipo(cun->ipo);
170
171         id_us_plus((ID *)cun->vfont);
172         id_us_plus((ID *)cun->vfontb);  
173         id_us_plus((ID *)cun->vfonti);
174         id_us_plus((ID *)cun->vfontbi);
175         
176         return cun;
177 }
178
179 void make_local_curve(Curve *cu)
180 {
181         Object *ob = 0;
182         Curve *cun;
183         int local=0, lib=0;
184         
185         /* - when there are only lib users: don't do
186          * - when there are only local users: set flag
187          * - mixed: do a copy
188          */
189         
190         if(cu->id.lib==0) return;
191         
192         if(cu->vfont) cu->vfont->id.lib= 0;
193         
194         if(cu->id.us==1) {
195                 cu->id.lib= 0;
196                 cu->id.flag= LIB_LOCAL;
197                 new_id(0, (ID *)cu, 0);
198                 return;
199         }
200         
201         ob= G.main->object.first;
202         while(ob) {
203                 if(ob->data==cu) {
204                         if(ob->id.lib) lib= 1;
205                         else local= 1;
206                 }
207                 ob= ob->id.next;
208         }
209         
210         if(local && lib==0) {
211                 cu->id.lib= 0;
212                 cu->id.flag= LIB_LOCAL;
213                 new_id(0, (ID *)cu, 0);
214         }
215         else if(local && lib) {
216                 cun= copy_curve(cu);
217                 cun->id.us= 0;
218                 
219                 ob= G.main->object.first;
220                 while(ob) {
221                         if(ob->data==cu) {
222                                 
223                                 if(ob->id.lib==0) {
224                                         ob->data= cun;
225                                         cun->id.us++;
226                                         cu->id.us--;
227                                 }
228                         }
229                         ob= ob->id.next;
230                 }
231         }
232 }
233
234
235 void test_curve_type(Object *ob)
236 {
237         Nurb *nu;
238         Curve *cu;
239         
240         cu= ob->data;
241         if(cu->vfont) {
242                 ob->type= OB_FONT;
243                 return;
244         }
245         else {
246                 nu= cu->nurb.first;
247                 while(nu) {
248                         if(nu->pntsv>1) {
249                                 ob->type= OB_SURF;
250                                 return;
251                         }
252                         nu= nu->next;
253                 }
254         }
255         ob->type= OB_CURVE;
256 }
257
258 void tex_space_curve(Curve *cu)
259 {
260         DispList *dl;
261         BoundBox *bb;
262         float *data, min[3], max[3], loc[3], size[3];
263         int tot, doit= 0;
264         
265         if(cu->bb==0) cu->bb= MEM_callocN(sizeof(BoundBox), "boundbox");
266         bb= cu->bb;
267         
268         INIT_MINMAX(min, max);
269
270         dl= cu->disp.first;
271         while(dl) {
272                 
273                 if(dl->type==DL_INDEX3 || dl->type==DL_INDEX3) tot= dl->nr;
274                 else tot= dl->nr*dl->parts;
275                 
276                 if(tot) doit= 1;
277                 data= dl->verts;
278                 while(tot--) {
279                         DO_MINMAX(data, min, max);
280                         data+= 3;
281                 }
282                 dl= dl->next;
283         }
284
285         if(!doit) {
286                 min[0] = min[1] = min[2] = -1.0f;
287                 max[0] = max[1] = max[2] = 1.0f;
288         }
289         
290         loc[0]= (min[0]+max[0])/2.0f;
291         loc[1]= (min[1]+max[1])/2.0f;
292         loc[2]= (min[2]+max[2])/2.0f;
293         
294         size[0]= (max[0]-min[0])/2.0f;
295         size[1]= (max[1]-min[1])/2.0f;
296         size[2]= (max[2]-min[2])/2.0f;
297
298         boundbox_set_from_min_max(bb, min, max);
299
300         if(cu->texflag & CU_AUTOSPACE) {
301                 VECCOPY(cu->loc, loc);
302                 VECCOPY(cu->size, size);
303                 cu->rot[0]= cu->rot[1]= cu->rot[2]= 0.0;
304
305                 if(cu->size[0]==0.0) cu->size[0]= 1.0;
306                 else if(cu->size[0]>0.0 && cu->size[0]<0.00001) cu->size[0]= 0.00001;
307                 else if(cu->size[0]<0.0 && cu->size[0]> -0.00001) cu->size[0]= -0.00001;
308         
309                 if(cu->size[1]==0.0) cu->size[1]= 1.0;
310                 else if(cu->size[1]>0.0 && cu->size[1]<0.00001) cu->size[1]= 0.00001;
311                 else if(cu->size[1]<0.0 && cu->size[1]> -0.00001) cu->size[1]= -0.00001;
312         
313                 if(cu->size[2]==0.0) cu->size[2]= 1.0;
314                 else if(cu->size[2]>0.0 && cu->size[2]<0.00001) cu->size[2]= 0.00001;
315                 else if(cu->size[2]<0.0 && cu->size[2]> -0.00001) cu->size[2]= -0.00001;
316
317         }
318 }
319
320
321 int count_curveverts(ListBase *nurb)
322 {
323         Nurb *nu;
324         int tot=0;
325         
326         nu= nurb->first;
327         while(nu) {
328                 if(nu->bezt) tot+= 3*nu->pntsu;
329                 else if(nu->bp) tot+= nu->pntsu*nu->pntsv;
330                 
331                 nu= nu->next;
332         }
333         return tot;
334 }
335
336
337
338 /* **************** NURBS ROUTINES ******************** */
339
340 void freeNurb(Nurb *nu)
341 {
342
343         if(nu==0) return;
344
345         if(nu->bezt) MEM_freeN(nu->bezt);
346         nu->bezt= 0;
347         if(nu->bp) MEM_freeN(nu->bp);
348         nu->bp= 0;
349         if(nu->knotsu) MEM_freeN(nu->knotsu);
350         nu->knotsu= 0;
351         if(nu->knotsv) MEM_freeN(nu->knotsv);
352         nu->knotsv= 0;
353         /* if(nu->trim.first) freeNurblist(&(nu->trim)); */
354
355         MEM_freeN(nu);
356
357 }
358
359
360 void freeNurblist(ListBase *lb)
361 {
362         Nurb *nu, *next;
363
364         if(lb==0) return;
365
366         nu= lb->first;
367         while(nu) {
368                 next= nu->next;
369                 freeNurb(nu);
370                 nu= next;
371         }
372         lb->first= lb->last= 0;
373 }
374
375 Nurb *duplicateNurb(Nurb *nu)
376 {
377         Nurb *newnu;
378         int len;
379
380         newnu= (Nurb*)MEM_mallocN(sizeof(Nurb),"duplicateNurb");
381         if(newnu==0) return 0;
382         memcpy(newnu, nu, sizeof(Nurb));
383
384         if(nu->bezt) {
385                 newnu->bezt=
386                         (BezTriple*)MEM_mallocN((nu->pntsu)* sizeof(BezTriple),"duplicateNurb2");
387                 memcpy(newnu->bezt, nu->bezt, nu->pntsu*sizeof(BezTriple));
388         }
389         else {
390                 len= nu->pntsu*nu->pntsv;
391                 newnu->bp=
392                         (BPoint*)MEM_mallocN((len)* sizeof(BPoint),"duplicateNurb3");
393                 memcpy(newnu->bp, nu->bp, len*sizeof(BPoint));
394                 
395                 newnu->knotsu=newnu->knotsv= 0;
396                 
397                 if(nu->knotsu) {
398                         len= KNOTSU(nu);
399                         if(len) {
400                                 newnu->knotsu= MEM_mallocN(len*sizeof(float), "duplicateNurb4");
401                                 memcpy(newnu->knotsu, nu->knotsu, sizeof(float)*len);
402                         }
403                 }
404                 if(nu->pntsv>1 && nu->knotsv) {
405                         len= KNOTSV(nu);
406                         if(len) {
407                                 newnu->knotsv= MEM_mallocN(len*sizeof(float), "duplicateNurb5");
408                                 memcpy(newnu->knotsv, nu->knotsv, sizeof(float)*len);
409                         }
410                 }
411         }
412         return newnu;
413 }
414
415 void duplicateNurblist(ListBase *lb1, ListBase *lb2)
416 {
417         Nurb *nu, *nun;
418         
419         freeNurblist(lb1);
420         
421         nu= lb2->first;
422         while(nu) {
423                 nun= duplicateNurb(nu);
424                 BLI_addtail(lb1, nun);
425                 
426                 nu= nu->next;
427         }
428 }
429
430 void test2DNurb(Nurb *nu)
431 {
432         BezTriple *bezt;
433         BPoint *bp;
434         int a;
435
436         if( nu->type== CU_BEZIER+CU_2D ) {
437                 a= nu->pntsu;
438                 bezt= nu->bezt;
439                 while(a--) {
440                         bezt->vec[0][2]= 0.0; 
441                         bezt->vec[1][2]= 0.0; 
442                         bezt->vec[2][2]= 0.0;
443                         bezt++;
444                 }
445         }
446         else if(nu->type & CU_2D) {
447                 a= nu->pntsu*nu->pntsv;
448                 bp= nu->bp;
449                 while(a--) {
450                         bp->vec[2]= 0.0;
451                         bp++;
452                 }
453         }
454 }
455
456 void minmaxNurb(Nurb *nu, float *min, float *max)
457 {
458         BezTriple *bezt;
459         BPoint *bp;
460         int a;
461
462         if( (nu->type & 7)==CU_BEZIER ) {
463                 a= nu->pntsu;
464                 bezt= nu->bezt;
465                 while(a--) {
466                         DO_MINMAX(bezt->vec[0], min, max);
467                         DO_MINMAX(bezt->vec[1], min, max);
468                         DO_MINMAX(bezt->vec[2], min, max);
469                         bezt++;
470                 }
471         }
472         else {
473                 a= nu->pntsu*nu->pntsv;
474                 bp= nu->bp;
475                 while(a--) {
476                         DO_MINMAX(bp->vec, min, max);
477                         bp++;
478                 }
479         }
480
481 }
482
483 /* ~~~~~~~~~~~~~~~~~~~~Non Uniform Rational B Spline calculations ~~~~~~~~~~~ */
484
485
486 /* actually, doubles should be used here as much as possible */
487
488 void extend_spline(float * pnts, int in, int out)
489 {
490         float *_pnts;
491         double * add;
492         int i, j, k, in2;
493
494         _pnts = pnts;
495         add = (double*)MEM_mallocN((in)* sizeof(double), "extend_spline");
496
497         in2 = in -1;
498
499         for (k = 3; k > 0; k--){
500                 pnts = _pnts;
501
502                 /* copy points to 'add' */
503                 for (i = 0; i < in; i++){
504                         add[i] = *pnts;
505                         pnts += 3;
506                 }
507
508                 /* inverse forward differencing */
509                 for (i = 0; i < in2; i++){
510                         for (j = in2; j > i; j--){
511                                 add[j] -= add[j - 1];
512                         }
513                 }
514
515                 pnts = _pnts;
516                 for (i = out; i > 0; i--){
517                         *pnts = (float)(add[0]);
518                         pnts += 3;
519                         for (j = 0; j < in2; j++){
520                                 add[j] += add[j+1];
521                         }
522                 }
523
524                 _pnts++;
525         }
526
527         MEM_freeN(add);
528 }
529
530
531 void calcknots(float *knots, short aantal, short order, short type)
532 /* knots: number of pnts NOT corrected for cyclic */
533 /* type;         0: uniform, 1: endpoints, 2: bezier */
534 {
535         float k;
536         int a, t;
537
538         t = aantal+order;
539         if(type==0) {
540
541                 for(a=0;a<t;a++) {
542                         knots[a]= (float)a;
543                 }
544         }
545         else if(type==1) {
546                 k= 0.0;
547                 for(a=1;a<=t;a++) {
548                         knots[a-1]= k;
549                         if(a>=order && a<=aantal) k+= 1.0;
550                 }
551         }
552         else if(type==2) {
553                 if(order==4) {
554                         k= 0.34;
555                         for(a=0;a<t;a++) {
556                                 knots[a]= (float)floor(k);
557                                 k+= (1.0/3.0);
558                         }
559                 }
560                 else if(order==3) {
561                         k= 0.6;
562                         for(a=0;a<t;a++) {
563                                 if(a>=order && a<=aantal) k+= (0.5);
564                                 knots[a]= (float)floor(k);
565                         }
566                 }
567         }
568 }
569
570 void makecyclicknots(float *knots, short pnts, short order)
571 /* pnts, order: number of pnts NOT corrected for cyclic */
572 {
573         int a, b, order2, c;
574
575         if(knots==0) return;
576         order2=order-1;
577
578         /* do first long rows (order -1), remove identical knots at endpoints */
579         if(order>2) {
580                 b= pnts+order2;
581                 for(a=1; a<order2; a++) {
582                         if(knots[b]!= knots[b-a]) break;
583                 }
584                 if(a==order2) knots[pnts+order-2]+= 1.0;
585         }
586
587         b= order;
588         c=pnts + order + order2;
589         for(a=pnts+order2; a<c; a++) {
590                 knots[a]= knots[a-1]+ (knots[b]-knots[b-1]);
591                 b--;
592         }
593 }
594
595
596 void makeknots(Nurb *nu, short uv, short type)  /* 0: uniform, 1: endpoints, 2: bezier */
597 {
598         if( (nu->type & 7)==CU_NURBS ) {
599                 if(uv & 1) {
600                         if(nu->knotsu) MEM_freeN(nu->knotsu);
601                         if(nu->pntsu>1) {
602                                 nu->knotsu= MEM_callocN(4+sizeof(float)*KNOTSU(nu), "makeknots");
603                                 calcknots(nu->knotsu, nu->pntsu, nu->orderu, type);
604                                 if(nu->flagu & 1) makecyclicknots(nu->knotsu, nu->pntsu, nu->orderu);
605                         }
606                         else nu->knotsu= 0;
607                 }
608                 if(uv & 2) {
609                         if(nu->knotsv) MEM_freeN(nu->knotsv);
610                         if(nu->pntsv>1) {
611                                 nu->knotsv= MEM_callocN(4+sizeof(float)*KNOTSV(nu), "makeknots");
612                                 calcknots(nu->knotsv, nu->pntsv, nu->orderv, type);
613                                 if(nu->flagv & 1) makecyclicknots(nu->knotsv, nu->pntsv, nu->orderv);
614                         }
615                         else nu->knotsv= 0;
616                 }
617         }
618 }
619
620 void basisNurb(float t, short order, short pnts, float *knots, float *basis, int *start, int *end)
621 {
622         float d, e;
623         int i, i1 = 0, i2 = 0 ,j, orderpluspnts, opp2, o2;
624
625         orderpluspnts= order+pnts;
626         opp2 = orderpluspnts-1;
627
628         /* this is for float inaccuracy */
629         if(t < knots[0]) t= knots[0];
630         else if(t > knots[opp2]) t= knots[opp2];
631
632         /* this part is order '1' */
633         o2 = order + 1;
634         for(i=0;i<opp2;i++) {
635                 if(knots[i]!=knots[i+1] && t>= knots[i] && t<=knots[i+1]) {
636                         basis[i]= 1.0;
637                         i1= i-o2;
638                         if(i1<0) i1= 0;
639                         i2= i;
640                         i++;
641                         while(i<opp2) {
642                                 basis[i]= 0.0;
643                                 i++;
644                         }
645                         break;
646                 }
647                 else basis[i]= 0.0;
648         }
649         basis[i]= 0.0;
650         
651         /* this is order 2,3,... */
652         for(j=2; j<=order; j++) {
653
654                 if(i2+j>= orderpluspnts) i2= opp2-j;
655
656                 for(i= i1; i<=i2; i++) {
657                         if(basis[i]!=0.0)
658                                 d= ((t-knots[i])*basis[i]) / (knots[i+j-1]-knots[i]);
659                         else
660                                 d= 0.0;
661
662                         if(basis[i+1]!=0.0)
663                                 e= ((knots[i+j]-t)*basis[i+1]) / (knots[i+j]-knots[i+1]);
664                         else
665                                 e= 0.0;
666
667                         basis[i]= d+e;
668                 }
669         }
670
671         *start= 1000;
672         *end= 0;
673
674         for(i=i1; i<=i2; i++) {
675                 if(basis[i]>0.0) {
676                         *end= i;
677                         if(*start==1000) *start= i;
678                 }
679         }
680 }
681
682
683 void makeNurbfaces(Nurb *nu, float *data, int rowstride) 
684 /* data  has to be 3*4*resolu*resolv in size, and zero-ed */
685 {
686         BPoint *bp;
687         float *basisu, *basis, *basisv, *sum, *fp, *in;
688         float u, v, ustart, uend, ustep, vstart, vend, vstep, sumdiv;
689         int i, j, iofs, jofs, cycl, len, resolu, resolv;
690         int istart, iend, jsta, jen, *jstart, *jend, ratcomp;
691
692         if(nu->knotsu==0 || nu->knotsv==0) return;
693         if(nu->orderu>nu->pntsu) return;
694         if(nu->orderv>nu->pntsv) return;
695         if(data==0) return;
696
697         /* allocate and initialize */
698         len= nu->pntsu*nu->pntsv;
699         if(len==0) return;
700         
701
702         
703         sum= (float *)MEM_callocN(sizeof(float)*len, "makeNurbfaces1");
704
705         resolu= nu->resolu;
706         resolv= nu->resolv;
707         len= resolu*resolv;
708         if(len==0) {
709                 MEM_freeN(sum);
710                 return;
711         }
712
713         bp= nu->bp;
714         i= nu->pntsu*nu->pntsv;
715         ratcomp=0;
716         while(i--) {
717                 if(bp->vec[3]!=1.0) {
718                         ratcomp= 1;
719                         break;
720                 }
721                 bp++;
722         }
723
724         fp= nu->knotsu;
725         ustart= fp[nu->orderu-1];
726         if(nu->flagu & 1) uend= fp[nu->pntsu+nu->orderu-1];
727         else uend= fp[nu->pntsu];
728         ustep= (uend-ustart)/(resolu-1+(nu->flagu & 1));
729         basisu= (float *)MEM_mallocN(sizeof(float)*KNOTSU(nu), "makeNurbfaces3");
730
731         fp= nu->knotsv;
732         vstart= fp[nu->orderv-1];
733         
734         if(nu->flagv & 1) vend= fp[nu->pntsv+nu->orderv-1];
735         else vend= fp[nu->pntsv];
736         vstep= (vend-vstart)/(resolv-1+(nu->flagv & 1));
737         len= KNOTSV(nu);
738         basisv= (float *)MEM_mallocN(sizeof(float)*len*resolv, "makeNurbfaces3");
739         jstart= (int *)MEM_mallocN(sizeof(float)*resolv, "makeNurbfaces4");
740         jend= (int *)MEM_mallocN(sizeof(float)*resolv, "makeNurbfaces5");
741
742         /* precalculation of basisv and jstart,jend */
743         if(nu->flagv & 1) cycl= nu->orderv-1; 
744         else cycl= 0;
745         v= vstart;
746         basis= basisv;
747         while(resolv--) {
748                 basisNurb(v, nu->orderv, (short)(nu->pntsv+cycl), nu->knotsv, basis, jstart+resolv, jend+resolv);
749                 basis+= KNOTSV(nu);
750                 v+= vstep;
751         }
752
753         if(nu->flagu & 1) cycl= nu->orderu-1; 
754         else cycl= 0;
755         in= data;
756         u= ustart;
757         while(resolu--) {
758
759                 basisNurb(u, nu->orderu, (short)(nu->pntsu+cycl), nu->knotsu, basisu, &istart, &iend);
760
761                 basis= basisv;
762                 resolv= nu->resolv;
763                 while(resolv--) {
764
765                         jsta= jstart[resolv];
766                         jen= jend[resolv];
767
768                         /* calculate sum */
769                         sumdiv= 0.0;
770                         fp= sum;
771
772                         for(j= jsta; j<=jen; j++) {
773
774                                 if(j>=nu->pntsv) jofs= (j - nu->pntsv);
775                                 else jofs= j;
776                                 bp= nu->bp+ nu->pntsu*jofs+istart-1;
777
778                                 for(i= istart; i<=iend; i++, fp++) {
779
780                                         if(i>= nu->pntsu) {
781                                                 iofs= i- nu->pntsu;
782                                                 bp= nu->bp+ nu->pntsu*jofs+iofs;
783                                         }
784                                         else bp++;
785
786                                         if(ratcomp) {
787                                                 *fp= basisu[i]*basis[j]*bp->vec[3];
788                                                 sumdiv+= *fp;
789                                         }
790                                         else *fp= basisu[i]*basis[j];
791                                 }
792                         }
793                 
794                         if(ratcomp) {
795                                 fp= sum;
796                                 for(j= jsta; j<=jen; j++) {
797                                         for(i= istart; i<=iend; i++, fp++) {
798                                                 *fp/= sumdiv;
799                                         }
800                                 }
801                         }
802
803                         /* one! (1.0) real point now */
804                         fp= sum;
805                         for(j= jsta; j<=jen; j++) {
806
807                                 if(j>=nu->pntsv) jofs= (j - nu->pntsv);
808                                 else jofs= j;
809                                 bp= nu->bp+ nu->pntsu*jofs+istart-1;
810
811                                 for(i= istart; i<=iend; i++, fp++) {
812
813                                         if(i>= nu->pntsu) {
814                                                 iofs= i- nu->pntsu;
815                                                 bp= nu->bp+ nu->pntsu*jofs+iofs;
816                                         }
817                                         else bp++;
818
819                                         if(*fp!=0.0) {
820                                                 in[0]+= (*fp) * bp->vec[0];
821                                                 in[1]+= (*fp) * bp->vec[1];
822                                                 in[2]+= (*fp) * bp->vec[2];
823                                         }
824                                 }
825                         }
826
827                         in+=3;
828                         basis+= KNOTSV(nu);
829                 }
830                 u+= ustep;
831                 if (rowstride!=0) in = (float*) (((unsigned char*) in) + (rowstride - 3*nu->resolv*sizeof(*in)));
832         }
833
834         for (i=0; i<144*3; i++) {
835 //              fprintf(stderr, "%f %f %f\n", nu->bp[i].vec[0], nu->bp[i].vec[1], nu->bp[i].vec[2]);
836                 fprintf(stderr, "%f ", data[i]);
837         }
838
839         /* free */
840         MEM_freeN(sum);
841         MEM_freeN(basisu);
842         MEM_freeN(basisv);
843         MEM_freeN(jstart);
844         MEM_freeN(jend);
845 }
846
847
848 void makeNurbcurve_forw(Nurb *nu, float *data)
849 /* *data: has to be 3*4*pntsu*resolu in size and zero-ed */
850 {
851         BPoint *bp;
852         float *basisu, *sum, *fp,  *in;
853         float u, ustart, uend, ustep, sumdiv;
854         int i, j, k, len, resolu, istart, iend;
855         int wanted, org;
856
857         if(nu->knotsu==0) return;
858         if(data==0) return;
859
860         /* allocate and init */
861         len= nu->pntsu;
862         if(len==0) return;
863         sum= (float *)MEM_callocN(sizeof(float)*len, "makeNurbcurve1");
864
865         resolu= nu->resolu*nu->pntsu;
866         if(resolu==0) {
867                 MEM_freeN(sum);
868                 return;
869         }
870
871         fp= nu->knotsu;
872         ustart= fp[nu->orderu-1];
873         uend= fp[nu->pntsu];
874         ustep= (uend-ustart)/(resolu-1);
875         basisu= (float *)MEM_mallocN(sizeof(float)*(nu->orderu+nu->pntsu), "makeNurbcurve3");
876
877         in= data;
878         u= ustart;
879         for (k = nu->orderu - 1; k < nu->pntsu; k++){
880
881                 wanted = (int)((nu->knotsu[k+1] - nu->knotsu[k]) / ustep);
882                 org = 4;        /* equal to order */
883                 if (org > wanted) org = wanted;
884
885                 for (j = org; j > 0; j--){
886
887                         basisNurb(u, nu->orderu, nu->pntsu, nu->knotsu, basisu, &istart, &iend);
888                         /* calc sum */
889                         sumdiv= 0.0;
890                         fp= sum;
891                         for(i= istart; i<=iend; i++, fp++) {
892                                 /* do the rational component */
893                                 *fp= basisu[i];
894                                 sumdiv+= *fp;
895                         }
896                         if(sumdiv!=0.0) if(sumdiv<0.999 || sumdiv>1.001) {
897                                 /* is this normalizing needed? */
898                                 fp= sum;
899                                 for(i= istart; i<=iend; i++, fp++) {
900                                         *fp/= sumdiv;
901                                 }
902                         }
903
904                         /* one! (1.0) real point */
905                         fp= sum;
906                         bp= nu->bp+ istart;
907                         for(i= istart; i<=iend; i++, bp++, fp++) {
908
909                                 if(*fp!=0.0) {
910                                         in[0]+= (*fp) * bp->vec[0];
911                                         in[1]+= (*fp) * bp->vec[1];
912                                         in[2]+= (*fp) * bp->vec[2];
913                                 }
914                         }
915
916                         in+=3;
917
918                         u+= ustep;
919                 }
920
921                 if (wanted > org){
922                         extend_spline(in - 3 * org, org, wanted);
923                         in += 3 * (wanted - org);
924                         u += ustep * (wanted - org);
925                 }
926         }
927
928         /* free */
929         MEM_freeN(sum);
930         MEM_freeN(basisu);
931 }
932
933
934 void makeNurbcurve(Nurb *nu, float *data, int dim)
935 /* data has to be dim*4*pntsu*resolu in size and zero-ed */
936 {
937         BPoint *bp;
938         float u, ustart, uend, ustep, sumdiv;
939         float *basisu, *sum, *fp,  *in;
940         int i, len, resolu, istart, iend, cycl;
941
942         if(nu->knotsu==0) return;
943         if(nu->orderu>nu->pntsu) return;
944         if(data==0) return;
945
946         /* allocate and initialize */
947         len= nu->pntsu;
948         if(len==0) return;
949         sum= (float *)MEM_callocN(sizeof(float)*len, "makeNurbcurve1");
950
951         resolu= nu->resolu*nu->pntsu;
952         if(resolu==0) {
953                 MEM_freeN(sum);
954                 return;
955         }
956
957         fp= nu->knotsu;
958         ustart= fp[nu->orderu-1];
959         if(nu->flagu & 1) uend= fp[nu->pntsu+nu->orderu-1];
960         else uend= fp[nu->pntsu];
961         ustep= (uend-ustart)/(resolu-1+(nu->flagu & 1));
962         basisu= (float *)MEM_mallocN(sizeof(float)*KNOTSU(nu), "makeNurbcurve3");
963
964         if(nu->flagu & 1) cycl= nu->orderu-1; 
965         else cycl= 0;
966
967         in= data;
968         u= ustart;
969         while(resolu--) {
970
971                 basisNurb(u, nu->orderu, (short)(nu->pntsu+cycl), nu->knotsu, basisu, &istart, &iend);
972                 /* calc sum */
973                 sumdiv= 0.0;
974                 fp= sum;
975                 bp= nu->bp+ istart-1;
976                 for(i= istart; i<=iend; i++, fp++) {
977
978                         if(i>=nu->pntsu) bp= nu->bp+(i - nu->pntsu);
979                         else bp++;
980
981                         *fp= basisu[i]*bp->vec[3];
982                         sumdiv+= *fp;
983                 }
984                 if(sumdiv!=0.0) if(sumdiv<0.999 || sumdiv>1.001) {
985                         /* is normalizing needed? */
986                         fp= sum;
987                         for(i= istart; i<=iend; i++, fp++) {
988                                 *fp/= sumdiv;
989                         }
990                 }
991
992                 /* one! (1.0) real point */
993                 fp= sum;
994                 bp= nu->bp+ istart-1;
995                 for(i= istart; i<=iend; i++, fp++) {
996
997                         if(i>=nu->pntsu) bp= nu->bp+(i - nu->pntsu);
998                         else bp++;
999
1000                         if(*fp!=0.0) {
1001                                 
1002                                 in[0]+= (*fp) * bp->vec[0];
1003                                 in[1]+= (*fp) * bp->vec[1];
1004                                 if(dim>=3) {
1005                                         in[2]+= (*fp) * bp->vec[2];
1006                                         if(dim==4) in[3]+= (*fp) * bp->alfa;
1007                                 }
1008                         }
1009                 }
1010
1011                 in+= dim;
1012
1013                 u+= ustep;
1014         }
1015
1016         /* free */
1017         MEM_freeN(sum);
1018         MEM_freeN(basisu);
1019 }
1020
1021 /* forward differencing method for bezier curve */
1022 void forward_diff_bezier(float q0, float q1, float q2, float q3, float *p, int it, int stride)
1023 {
1024         float rt0,rt1,rt2,rt3,f;
1025         int a;
1026
1027         f= (float)it;
1028         rt0= q0;
1029         rt1= 3.0f*(q1-q0)/f;
1030         f*= f;
1031         rt2= 3.0f*(q0-2.0f*q1+q2)/f;
1032         f*= it;
1033         rt3= (q3-q0+3.0f*(q1-q2))/f;
1034         
1035         q0= rt0;
1036         q1= rt1+rt2+rt3;
1037         q2= 2*rt2+6*rt3;
1038         q3= 6*rt3;
1039   
1040         for(a=0; a<=it; a++) {
1041                 *p= q0;
1042                 p+= stride;
1043                 q0+= q1;
1044                 q1+= q2;
1045                 q2+= q3;
1046         }
1047 }       
1048
1049 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
1050
1051 float *make_orco_surf(Object *ob)
1052 {
1053         Curve *cu = ob->data;
1054         Nurb *nu;
1055         int a, b, tot=0;
1056         int sizeu, sizev;// ###
1057         float *data;
1058         float *orco;
1059
1060                 /* first calculate the size of the datablock */
1061         for (nu=cu->nurb.first; nu; nu=nu->next) {
1062                 sizeu = nu->resolu; sizev = nu->resolv;
1063                 if(nu->pntsv>1) tot+= sizeu * sizev;
1064         }
1065                                 /* makeNurbfaces wants zeros */
1066         data= orco= MEM_callocN(3*sizeof(float)*tot, "make_orco");
1067
1068         for (nu=cu->nurb.first; nu; nu=nu->next) {
1069                 if(nu->pntsv>1) {
1070                         sizeu = nu->resolu;
1071                         sizev = nu->resolv;
1072                         
1073                         if(cu->flag & CU_UV_ORCO) {
1074                                 for(b=0; b< sizeu; b++) {
1075                                         for(a=0; a< sizev; a++) {
1076                                         
1077                                                 if(sizev <2) data[0]= 0.0f;
1078                                                 else data[0]= -1.0f + 2.0f*((float)a)/(sizev - 1);
1079                                                 
1080                                                 if(sizeu <2) data[1]= 0.0f;
1081                                                 else data[1]= -1.0f + 2.0f*((float)b)/(sizeu - 1);
1082                                                 
1083                                                 data[2]= 0.0;
1084                 
1085                                                 data+= 3;
1086                                         }
1087                                 }
1088                         }
1089                         else {
1090                                 makeNurbfaces(nu, data, sizeof(*data)*sizev*3);
1091
1092                                 for(b=0; b<sizeu; b++) {
1093                                         for(a=0; a<sizev; a++) {
1094                                                 data = orco + 3 * (b * sizev + a);
1095                                                 data[0]= (data[0]-cu->loc[0])/cu->size[0];
1096                                                 data[1]= (data[1]-cu->loc[1])/cu->size[1];
1097                                                 data[2]= (data[2]-cu->loc[2])/cu->size[2];
1098                                         }
1099                                 }
1100                         }
1101                 }
1102         }
1103
1104         return orco;
1105 }
1106
1107
1108         /* NOTE: This routine is tied to the order of vertex
1109          * built by displist and as passed to the renderer.
1110          */
1111 float *make_orco_curve(Object *ob)
1112 {
1113         Curve *cu = ob->data;
1114         DispList *dl;
1115         int u, v, numVerts;
1116         float *fp, *orco;
1117         int remakeDisp = 0;
1118
1119         if (!(cu->flag&CU_UV_ORCO) && cu->key && cu->key->refkey) {
1120                 cp_cu_key(cu, cu->key->refkey, 0, count_curveverts(&cu->nurb));
1121                 makeDispListCurveTypes(ob, 1);
1122                 remakeDisp = 1;
1123         }
1124
1125                 /* Assumes displist has been built */
1126
1127         numVerts = 0;
1128         for (dl=cu->disp.first; dl; dl=dl->next) {
1129                 if (dl->type==DL_INDEX3) {
1130                         numVerts += dl->nr;
1131                 } else if (dl->type==DL_SURF) {
1132                         numVerts += dl->parts*dl->nr;
1133                 }
1134         }
1135
1136         fp= orco= MEM_mallocN(3*sizeof(float)*numVerts, "cu_orco");
1137
1138         for (dl=cu->disp.first; dl; dl=dl->next) {
1139                 if (dl->type==DL_INDEX3) {
1140                         for (u=0; u<dl->nr; u++,fp+=3) {
1141                                 if (cu->flag&CU_UV_ORCO) {
1142                                         fp[0]= 2.0f*u/(dl->nr-1) - 1.0f;
1143                                         fp[1]= 0.0;
1144                                         fp[2]= 0.0;
1145                                 } else {
1146                                         VECCOPY(fp, &dl->verts[u*3]);
1147
1148                                         fp[0]= (fp[0]-cu->loc[0])/cu->size[0];
1149                                         fp[1]= (fp[1]-cu->loc[1])/cu->size[1];
1150                                         fp[2]= (fp[2]-cu->loc[2])/cu->size[2];
1151                                 }
1152                         }
1153                 } else if (dl->type==DL_SURF) {
1154                         for (u=0; u<dl->parts; u++) {
1155                                 for (v=0; v<dl->nr; v++,fp+=3) {
1156                                         if (cu->flag&CU_UV_ORCO) {
1157                                                 fp[0]= 2.0f*u/(dl->parts-1) - 1.0f;
1158                                                 fp[1]= 2.0f*v/(dl->nr-1) - 1.0f;
1159                                                 fp[2]= 0.0;
1160                                         } else {
1161                                                 VECCOPY(fp, &dl->verts[(dl->nr*u + v)*3]);
1162
1163                                                 fp[0]= (fp[0]-cu->loc[0])/cu->size[0];
1164                                                 fp[1]= (fp[1]-cu->loc[1])/cu->size[1];
1165                                                 fp[2]= (fp[2]-cu->loc[2])/cu->size[2];
1166                                         }
1167                                 }
1168                         }
1169                 }
1170         }
1171
1172         if (remakeDisp) {
1173                 makeDispListCurveTypes(ob, 0);
1174         }
1175
1176         return orco;
1177 }
1178
1179
1180 /* ***************** BEVEL ****************** */
1181
1182 void makebevelcurve(Object *ob, ListBase *disp)
1183 {
1184         DispList *dl, *dlnew;
1185         Curve *bevcu, *cu;
1186         float *fp, facx, facy, hoek, dhoek;
1187         int nr, a;
1188
1189         cu= ob->data;
1190
1191         disp->first = disp->last = NULL;
1192         if(cu->bevobj && cu->bevobj!=ob) {
1193                 if(cu->bevobj->type==OB_CURVE) {
1194                         bevcu= cu->bevobj->data;
1195                         if(bevcu->ext1==0.0 && bevcu->ext2==0.0) {
1196                                 facx= cu->bevobj->size[0];
1197                                 facy= cu->bevobj->size[1];
1198
1199                                 dl= bevcu->disp.first;
1200                                 if(dl==0) {
1201                                         makeDispListCurveTypes(cu->bevobj, 0);
1202                                         dl= bevcu->disp.first;
1203                                 }
1204                                 while(dl) {
1205                                         if ELEM(dl->type, DL_POLY, DL_SEGM) {
1206                                                 dlnew= MEM_mallocN(sizeof(DispList), "makebevelcurve1");                                        
1207                                                 *dlnew= *dl;
1208                                                 dlnew->verts= MEM_mallocN(3*sizeof(float)*dl->parts*dl->nr, "makebevelcurve1");
1209                                                 memcpy(dlnew->verts, dl->verts, 3*sizeof(float)*dl->parts*dl->nr);
1210                                                 
1211                                                 if(dlnew->type==DL_SEGM) dlnew->flag |= (DL_FRONT_CURVE|DL_BACK_CURVE);
1212                                                 
1213                                                 BLI_addtail(disp, dlnew);
1214                                                 fp= dlnew->verts;
1215                                                 nr= dlnew->parts*dlnew->nr;
1216                                                 while(nr--) {
1217                                                         fp[2]= fp[1]*facy;
1218                                                         fp[1]= -fp[0]*facx;
1219                                                         fp[0]= 0.0;
1220                                                         fp+= 3;
1221                                                 }
1222                                         }
1223                                         dl= dl->next;
1224                                 }
1225                         }
1226                 }
1227         }
1228         else if(cu->ext1==0.0 && cu->ext2==0.0) {
1229                 ;
1230         }
1231         else if(cu->ext2==0.0) {
1232                 dl= MEM_callocN(sizeof(DispList), "makebevelcurve2");
1233                 dl->verts= MEM_mallocN(2*3*sizeof(float), "makebevelcurve2");
1234                 BLI_addtail(disp, dl);
1235                 dl->type= DL_SEGM;
1236                 dl->parts= 1;
1237                 dl->flag= DL_FRONT_CURVE|DL_BACK_CURVE;
1238                 dl->nr= 2;
1239                 
1240                 fp= dl->verts;
1241                 fp[0]= fp[1]= 0.0;
1242                 fp[2]= -cu->ext1;
1243                 fp[3]= fp[4]= 0.0;
1244                 fp[5]= cu->ext1;
1245         }
1246         else {
1247                 short dnr;
1248                 
1249                 /* bevel now in three parts, for proper vertex normals */
1250                 /* part 1 */
1251                 dnr= nr= 2+ cu->bevresol;
1252                 if( (cu->flag & (CU_FRONT|CU_BACK))==0) // we make a full round bevel in that case
1253                         nr= 3+ 2*cu->bevresol;
1254                    
1255                 dl= MEM_callocN(sizeof(DispList), "makebevelcurve p1");
1256                 dl->verts= MEM_mallocN(nr*3*sizeof(float), "makebevelcurve p1");
1257                 BLI_addtail(disp, dl);
1258                 dl->type= DL_SEGM;
1259                 dl->parts= 1;
1260                 dl->flag= DL_BACK_CURVE;
1261                 dl->nr= nr;
1262
1263                 /* half a circle */
1264                 fp= dl->verts;
1265                 dhoek= (0.5*M_PI/(dnr-1));
1266                 hoek= -(nr-1)*dhoek;
1267                 
1268                 for(a=0; a<nr; a++) {
1269                         fp[0]= 0.0;
1270                         fp[1]= (float)(cos(hoek)*(cu->ext2));
1271                         fp[2]= (float)(sin(hoek)*(cu->ext2)) - cu->ext1;
1272                         hoek+= dhoek;
1273                         fp+= 3;
1274                 }
1275                 
1276                 /* part 2, sidefaces */
1277                 if(cu->ext1!=0.0) {
1278                         nr= 2;
1279                         
1280                         dl= MEM_callocN(sizeof(DispList), "makebevelcurve p2");
1281                         dl->verts= MEM_callocN(nr*3*sizeof(float), "makebevelcurve p2");
1282                         BLI_addtail(disp, dl);
1283                         dl->type= DL_SEGM;
1284                         dl->parts= 1;
1285                         dl->nr= nr;
1286                         
1287                         fp= dl->verts;
1288                         fp[1]= cu->ext2;
1289                         fp[2]= -cu->ext1;
1290                         fp[4]= cu->ext2;
1291                         fp[5]= cu->ext1;
1292                         
1293                         if( (cu->flag & (CU_FRONT|CU_BACK))==0) {
1294                                 dl= MEM_dupallocN(dl);
1295                                 dl->verts= MEM_dupallocN(dl->verts);
1296                                 BLI_addtail(disp, dl);
1297                                 
1298                                 fp= dl->verts;
1299                                 fp[1]= -fp[1];
1300                                 fp[2]= -fp[2];
1301                                 fp[4]= -fp[4];
1302                                 fp[5]= -fp[5];
1303                         }
1304                 }
1305                 
1306                 /* part 3 */
1307                 dnr= nr= 2+ cu->bevresol;
1308                 if( (cu->flag & (CU_FRONT|CU_BACK))==0)
1309                         nr= 3+ 2*cu->bevresol;
1310                 
1311                 dl= MEM_callocN(sizeof(DispList), "makebevelcurve p3");
1312                 dl->verts= MEM_mallocN(nr*3*sizeof(float), "makebevelcurve p3");
1313                 BLI_addtail(disp, dl);
1314                 dl->type= DL_SEGM;
1315                 dl->flag= DL_FRONT_CURVE;
1316                 dl->parts= 1;
1317                 dl->nr= nr;
1318                 
1319                 /* half a circle */
1320                 fp= dl->verts;
1321                 hoek= 0.0;
1322                 dhoek= (0.5*M_PI/(dnr-1));
1323                 
1324                 for(a=0; a<nr; a++) {
1325                         fp[0]= 0.0;
1326                         fp[1]= (float)(cos(hoek)*(cu->ext2));
1327                         fp[2]= (float)(sin(hoek)*(cu->ext2)) + cu->ext1;
1328                         hoek+= dhoek;
1329                         fp+= 3;
1330                 }
1331         }
1332 }
1333
1334 int cu_isectLL(float *v1, float *v2, float *v3, float *v4, short cox, short coy, float *labda, float *mu, float *vec)
1335 {
1336         /* return:
1337                 -1: colliniar
1338                  0: no intersection of segments
1339                  1: exact intersection of segments
1340                  2: cross-intersection of segments
1341         */
1342         float deler;
1343
1344         deler= (v1[cox]-v2[cox])*(v3[coy]-v4[coy])-(v3[cox]-v4[cox])*(v1[coy]-v2[coy]);
1345         if(deler==0.0) return -1;
1346
1347         *labda= (v1[coy]-v3[coy])*(v3[cox]-v4[cox])-(v1[cox]-v3[cox])*(v3[coy]-v4[coy]);
1348         *labda= -(*labda/deler);
1349
1350         deler= v3[coy]-v4[coy];
1351         if(deler==0) {
1352                 deler=v3[cox]-v4[cox];
1353                 *mu= -(*labda*(v2[cox]-v1[cox])+v1[cox]-v3[cox])/deler;
1354         } else {
1355                 *mu= -(*labda*(v2[coy]-v1[coy])+v1[coy]-v3[coy])/deler;
1356         }
1357         vec[cox]= *labda*(v2[cox]-v1[cox])+v1[cox];
1358         vec[coy]= *labda*(v2[coy]-v1[coy])+v1[coy];
1359
1360         if(*labda>=0.0 && *labda<=1.0 && *mu>=0.0 && *mu<=1.0) {
1361                 if(*labda==0.0 || *labda==1.0 || *mu==0.0 || *mu==1.0) return 1;
1362                 return 2;
1363         }
1364         return 0;
1365 }
1366
1367
1368 short bevelinside(BevList *bl1,BevList *bl2)
1369 {
1370         /* is bl2 INSIDE bl1 ? with left-right method and "labda's" */
1371         /* returns '1' if correct hole  */
1372         BevPoint *bevp, *prevbevp;
1373         float min,max,vec[3],hvec1[3],hvec2[3],lab,mu;
1374         int nr, links=0,rechts=0,mode;
1375
1376         /* take first vertex of possible hole */
1377
1378         bevp= (BevPoint *)(bl2+1);
1379         hvec1[0]= bevp->x; 
1380         hvec1[1]= bevp->y; 
1381         hvec1[2]= 0.0;
1382         VECCOPY(hvec2,hvec1);
1383         hvec2[0]+=1000;
1384
1385         /* test it with all edges of potential surounding poly */
1386         /* count number of transitions left-right  */
1387
1388         bevp= (BevPoint *)(bl1+1);
1389         nr= bl1->nr;
1390         prevbevp= bevp+(nr-1);
1391
1392         while(nr--) {
1393                 min= prevbevp->y;
1394                 max= bevp->y;
1395                 if(max<min) {
1396                         min= max;
1397                         max= prevbevp->y;
1398                 }
1399                 if(min!=max) {
1400                         if(min<=hvec1[1] && max>=hvec1[1]) {
1401                                 /* there's a transition, calc intersection point */
1402                                 mode= cu_isectLL(&(prevbevp->x),&(bevp->x),hvec1,hvec2,0,1,&lab,&mu,vec);
1403                                 /* if lab==0.0 or lab==1.0 then the edge intersects exactly a transition
1404                                    only allow for one situation: we choose lab= 1.0
1405                                  */
1406                                 if(mode>=0 && lab!=0.0) {
1407                                         if(vec[0]<hvec1[0]) links++;
1408                                         else rechts++;
1409                                 }
1410                         }
1411                 }
1412                 prevbevp= bevp;
1413                 bevp++;
1414         }
1415         
1416         if( (links & 1) && (rechts & 1) ) return 1;
1417         return 0;
1418 }
1419
1420
1421 struct bevelsort {
1422         float left;
1423         BevList *bl;
1424         int dir;
1425 };
1426
1427 int vergxcobev(const void *a1, const void *a2)
1428 {
1429         const struct bevelsort *x1=a1,*x2=a2;
1430
1431         if( x1->left > x2->left ) return 1;
1432         else if( x1->left < x2->left) return -1;
1433         return 0;
1434 }
1435
1436 /* this function cannot be replaced with atan2, but why? */
1437
1438 void calc_bevel_sin_cos(float x1, float y1, float x2, float y2, float *sina, float *cosa)
1439 {
1440         float t01, t02, x3, y3;
1441
1442         t01= (float)sqrt(x1*x1+y1*y1);
1443         t02= (float)sqrt(x2*x2+y2*y2);
1444         if(t01==0.0) t01= 1.0;
1445         if(t02==0.0) t02= 1.0;
1446
1447         x1/=t01; 
1448         y1/=t01;
1449         x2/=t02; 
1450         y2/=t02;
1451
1452         t02= x1*x2+y1*y2;
1453         if(fabs(t02)>=1.0) t02= .5*M_PI;
1454         else t02= (saacos(t02))/2.0f;
1455
1456         t02= (float)sin(t02);
1457         if(t02==0.0) t02= 1.0;
1458
1459         x3= x1-x2;
1460         y3= y1-y2;
1461         if(x3==0 && y3==0) {
1462                 x3= y1;
1463                 y3= -x1;
1464         } else {
1465                 t01= (float)sqrt(x3*x3+y3*y3);
1466                 x3/=t01; 
1467                 y3/=t01;
1468         }
1469
1470         *sina= -y3/t02;
1471         *cosa= x3/t02;
1472
1473 }
1474
1475 void alfa_bezpart(BezTriple *prevbezt, BezTriple *bezt, Nurb *nu, float *data_a)
1476 {
1477         BezTriple *pprev, *next, *last;
1478         float fac, dfac, t[4];
1479         int a;
1480         
1481         last= nu->bezt+(nu->pntsu-1);
1482         
1483         /* returns a point */
1484         if(prevbezt==nu->bezt) {
1485                 if(nu->flagu & 1) pprev= last;
1486                 else pprev= prevbezt;
1487         }
1488         else pprev= prevbezt-1;
1489         
1490         /* next point */
1491         if(bezt==last) {
1492                 if(nu->flagu & 1) next= nu->bezt;
1493                 else next= bezt;
1494         }
1495         else next= bezt+1;
1496         
1497         fac= 0.0;
1498         dfac= 1.0f/(float)nu->resolu;
1499         
1500         for(a=0; a<nu->resolu; a++, fac+= dfac) {
1501                 
1502                 set_four_ipo(fac, t, KEY_BSPLINE);
1503                 
1504                 data_a[a]= t[0]*pprev->alfa + t[1]*prevbezt->alfa + t[2]*bezt->alfa + t[3]*next->alfa;
1505         }
1506 }
1507
1508 void makeBevelList(Object *ob)
1509 {
1510         /*
1511          - convert all curves to polys, with indication of resol and flags for double-vertices
1512          - possibly; do a smart vertice removal (in case Nurb)
1513          - separate in individual blicks with BoundBox
1514          - AutoHole detection
1515         */
1516         Curve *cu;
1517         Nurb *nu;
1518         BezTriple *bezt, *prevbezt;
1519         BPoint *bp;
1520         BevList *bl, *blnew, *blnext;
1521         BevPoint *bevp, *bevp2, *bevp1 = NULL, *bevp0;
1522         float  *data, *data_a, *v1, *v2, min, inp, x1, x2, y1, y2, vec[3];
1523         struct bevelsort *sortdata, *sd, *sd1;
1524         int a, b, len, nr, poly;
1525
1526         /* this function needs an object, because of tflag and upflag */
1527         cu= ob->data;
1528
1529         /* STEP 1: MAKE POLYS  */
1530
1531         BLI_freelistN(&(cu->bev));
1532         if(ob==G.obedit) nu= editNurb.first;
1533         else nu= cu->nurb.first;
1534         
1535         while(nu) {
1536                 if(nu->pntsu>1) {
1537                 
1538                         if((nu->type & 7)==CU_POLY) {
1539         
1540                                 len= nu->pntsu;
1541                                 bl= MEM_callocN(sizeof(BevList)+len*sizeof(BevPoint), "makeBevelList");
1542                                 BLI_addtail(&(cu->bev), bl);
1543         
1544                                 if(nu->flagu & 1) bl->poly= 0;
1545                                 else bl->poly= -1;
1546                                 bl->nr= len;
1547                                 bl->flag= 0;
1548                                 bevp= (BevPoint *)(bl+1);
1549                                 bp= nu->bp;
1550         
1551                                 while(len--) {
1552                                         bevp->x= bp->vec[0];
1553                                         bevp->y= bp->vec[1];
1554                                         bevp->z= bp->vec[2];
1555                                         bevp->alfa= bp->alfa;
1556                                         bevp->f1= 1;
1557                                         bevp++;
1558                                         bp++;
1559                                 }
1560                         }
1561                         else if((nu->type & 7)==CU_BEZIER) {
1562         
1563                                 len= nu->resolu*(nu->pntsu+ (nu->flagu & 1) -1)+1;      /* in case last point is not cyclic */
1564                                 bl= MEM_callocN(sizeof(BevList)+len*sizeof(BevPoint), "makeBevelList");
1565                                 BLI_addtail(&(cu->bev), bl);
1566         
1567                                 if(nu->flagu & 1) bl->poly= 0;
1568                                 else bl->poly= -1;
1569                                 bevp= (BevPoint *)(bl+1);
1570         
1571                                 a= nu->pntsu-1;
1572                                 bezt= nu->bezt;
1573                                 if(nu->flagu & 1) {
1574                                         a++;
1575                                         prevbezt= nu->bezt+(nu->pntsu-1);
1576                                 }
1577                                 else {
1578                                         prevbezt= bezt;
1579                                         bezt++;
1580                                 }
1581                                 
1582                                 data= MEM_mallocN(3*sizeof(float)*(nu->resolu+1), "makeBevelList2");
1583                                 data_a= MEM_callocN(sizeof(float)*(nu->resolu+1), "data_a");
1584                                 
1585                                 while(a--) {
1586                                         if(prevbezt->h2==HD_VECT && bezt->h1==HD_VECT) {
1587         
1588                                                 bevp->x= prevbezt->vec[1][0];
1589                                                 bevp->y= prevbezt->vec[1][1];
1590                                                 bevp->z= prevbezt->vec[1][2];
1591                                                 bevp->alfa= prevbezt->alfa;
1592                                                 bevp->f1= 1;
1593                                                 bevp->f2= 0;
1594                                                 bevp++;
1595                                                 bl->nr++;
1596                                                 bl->flag= 1;
1597                                         }
1598                                         else {
1599                                                 v1= prevbezt->vec[1];
1600                                                 v2= bezt->vec[0];
1601                                                 
1602                                                 /* always do all three, to prevent data hanging around */
1603                                                 forward_diff_bezier(v1[0], v1[3], v2[0], v2[3], data, nu->resolu, 3);
1604                                                 forward_diff_bezier(v1[1], v1[4], v2[1], v2[4], data+1, nu->resolu, 3);
1605                                                 forward_diff_bezier(v1[2], v1[5], v2[2], v2[5], data+2, nu->resolu, 3);
1606                                                 
1607                                                 if((nu->type & CU_2D)==0) {
1608                                                         if(cu->flag & CU_3D) {
1609                                                                 alfa_bezpart(prevbezt, bezt, nu, data_a);
1610                                                         }
1611                                                 }
1612                                                 
1613                                                 
1614                                                 /* indicate with handlecodes double points */
1615                                                 if(prevbezt->h1==prevbezt->h2) {
1616                                                         if(prevbezt->h1==0 || prevbezt->h1==HD_VECT) bevp->f1= 1;
1617                                                 }
1618                                                 else {
1619                                                         if(prevbezt->h1==0 || prevbezt->h1==HD_VECT) bevp->f1= 1;
1620                                                         else if(prevbezt->h2==0 || prevbezt->h2==HD_VECT) bevp->f1= 1;
1621                                                 }
1622                                                 
1623                                                 v1= data;
1624                                                 v2= data_a;
1625                                                 nr= nu->resolu;
1626                                                 
1627                                                 while(nr--) {
1628                                                         bevp->x= v1[0]; 
1629                                                         bevp->y= v1[1];
1630                                                         bevp->z= v1[2];
1631                                                         bevp->alfa= v2[0];
1632                                                         bevp++;
1633                                                         v1+=3;
1634                                                         v2++;
1635                                                 }
1636                                                 bl->nr+= nu->resolu;
1637         
1638                                         }
1639                                         prevbezt= bezt;
1640                                         bezt++;
1641                                 }
1642                                 
1643                                 MEM_freeN(data);
1644                                 MEM_freeN(data_a);
1645                                 
1646                                 if((nu->flagu & 1)==0) {            /* not cyclic: endpoint */
1647                                         bevp->x= prevbezt->vec[1][0];
1648                                         bevp->y= prevbezt->vec[1][1];
1649                                         bevp->z= prevbezt->vec[1][2];
1650                                         bl->nr++;
1651                                 }
1652         
1653                         }
1654                         else if((nu->type & 7)==CU_NURBS) {
1655                                 if(nu->pntsv==1) {
1656                                         len= nu->resolu*nu->pntsu;
1657                                         bl= MEM_mallocN(sizeof(BevList)+len*sizeof(BevPoint), "makeBevelList3");
1658                                         BLI_addtail(&(cu->bev), bl);
1659                                         bl->nr= len;
1660                                         bl->flag= 0;
1661                                         if(nu->flagu & 1) bl->poly= 0;
1662                                         else bl->poly= -1;
1663                                         bevp= (BevPoint *)(bl+1);
1664         
1665                                         data= MEM_callocN(4*sizeof(float)*len, "makeBevelList4");    /* has to be zero-ed */
1666                                         makeNurbcurve(nu, data, 4);
1667                                         
1668                                         v1= data;
1669                                         while(len--) {
1670                                                 bevp->x= v1[0]; 
1671                                                 bevp->y= v1[1];
1672                                                 bevp->z= v1[2];
1673                                                 bevp->alfa= v1[3];
1674                                                 
1675                                                 bevp->f1= bevp->f2= 0;
1676                                                 bevp++;
1677                                                 v1+=4;
1678                                         }
1679                                         MEM_freeN(data);
1680                                 }
1681                         }
1682                 }
1683                 nu= nu->next;
1684         }
1685
1686         /* STEP 2: DOUBLE POINTS AND AUTOMATIC RESOLUTION, REDUCE DATABLOCKS */
1687         bl= cu->bev.first;
1688         while(bl) {
1689                 nr= bl->nr;
1690                 bevp1= (BevPoint *)(bl+1);
1691                 bevp0= bevp1+(nr-1);
1692                 nr--;
1693                 while(nr--) {
1694                         if( fabs(bevp0->x-bevp1->x)<0.00001 ) {
1695                                 if( fabs(bevp0->y-bevp1->y)<0.00001 ) {
1696                                         if( fabs(bevp0->z-bevp1->z)<0.00001 ) {
1697                                                 bevp0->f2= 1;
1698                                                 bl->flag++;
1699                                         }
1700                                 }
1701                         }
1702                         bevp0= bevp1;
1703                         bevp1++;
1704                 }
1705                 bl= bl->next;
1706         }
1707         bl= cu->bev.first;
1708         while(bl) {
1709                 blnext= bl->next;
1710                 if(bl->flag) {
1711                         nr= bl->nr- bl->flag+1; /* +1 because vectorbezier sets flag too */
1712                         blnew= MEM_mallocN(sizeof(BevList)+nr*sizeof(BevPoint), "makeBevelList");
1713                         memcpy(blnew, bl, sizeof(BevList));
1714                         blnew->nr= 0;
1715                         BLI_remlink(&(cu->bev), bl);
1716                         BLI_insertlinkbefore(&(cu->bev),blnext,blnew);  /* to make sure bevlijst is tuned with nurblist */
1717                         bevp0= (BevPoint *)(bl+1);
1718                         bevp1= (BevPoint *)(blnew+1);
1719                         nr= bl->nr;
1720                         while(nr--) {
1721                                 if(bevp0->f2==0) {
1722                                         memcpy(bevp1, bevp0, sizeof(BevPoint));
1723                                         bevp1++;
1724                                         blnew->nr++;
1725                                 }
1726                                 bevp0++;
1727                         }
1728                         MEM_freeN(bl);
1729                         blnew->flag= 0;
1730                 }
1731                 bl= blnext;
1732         }
1733
1734         /* STEP 3: COUNT POLYS TELLEN AND AUTOHOLE */
1735         bl= cu->bev.first;
1736         poly= 0;
1737         while(bl) {
1738                 if(bl->poly>=0) {
1739                         poly++;
1740                         bl->poly= poly;
1741                         bl->gat= 0;     /* 'gat' is dutch for hole */
1742                 }
1743                 bl= bl->next;
1744         }
1745         
1746
1747         /* find extreme left points, also test (turning) direction */
1748         if(poly>0) {
1749                 sd= sortdata= MEM_mallocN(sizeof(struct bevelsort)*poly, "makeBevelList5");
1750                 bl= cu->bev.first;
1751                 while(bl) {
1752                         if(bl->poly>0) {
1753
1754                                 min= 300000.0;
1755                                 bevp= (BevPoint *)(bl+1);
1756                                 nr= bl->nr;
1757                                 while(nr--) {
1758                                         if(min>bevp->x) {
1759                                                 min= bevp->x;
1760                                                 bevp1= bevp;
1761                                         }
1762                                         bevp++;
1763                                 }
1764                                 sd->bl= bl;
1765                                 sd->left= min;
1766
1767                                 bevp= (BevPoint *)(bl+1);
1768                                 if(bevp1== bevp) bevp0= bevp+ (bl->nr-1);
1769                                 else bevp0= bevp1-1;
1770                                 bevp= bevp+ (bl->nr-1);
1771                                 if(bevp1== bevp) bevp2= (BevPoint *)(bl+1);
1772                                 else bevp2= bevp1+1;
1773
1774                                 inp= (bevp1->x- bevp0->x)*(bevp0->y- bevp2->y)
1775                                     +(bevp0->y- bevp1->y)*(bevp0->x- bevp2->x);
1776
1777                                 if(inp>0.0) sd->dir= 1;
1778                                 else sd->dir= 0;
1779
1780                                 sd++;
1781                         }
1782
1783                         bl= bl->next;
1784                 }
1785                 qsort(sortdata,poly,sizeof(struct bevelsort), vergxcobev);
1786
1787                 sd= sortdata+1;
1788                 for(a=1; a<poly; a++, sd++) {
1789                         bl= sd->bl;         /* is bl a hole? */
1790                         sd1= sortdata+ (a-1);
1791                         for(b=a-1; b>=0; b--, sd1--) {  /* all polys to the left */
1792                                 if(bevelinside(sd1->bl, bl)) {
1793                                         bl->gat= 1- sd1->bl->gat;
1794                                         break;
1795                                 }
1796                         }
1797                 }
1798
1799                 /* turning direction */
1800                 if((cu->flag & CU_3D)==0) {
1801                         sd= sortdata;
1802                         for(a=0; a<poly; a++, sd++) {
1803                                 if(sd->bl->gat==sd->dir) {
1804                                         bl= sd->bl;
1805                                         bevp1= (BevPoint *)(bl+1);
1806                                         bevp2= bevp1+ (bl->nr-1);
1807                                         nr= bl->nr/2;
1808                                         while(nr--) {
1809                                                 SWAP(BevPoint, *bevp1, *bevp2);
1810                                                 bevp1++;
1811                                                 bevp2--;
1812                                         }
1813                                 }
1814                         }
1815                 }
1816                 MEM_freeN(sortdata);
1817         }
1818
1819         /* STEP 4: COSINES */
1820         bl= cu->bev.first;
1821         while(bl) {
1822         
1823                 if(bl->nr==2) { /* 2 pnt, treat separate */
1824                         bevp2= (BevPoint *)(bl+1);
1825                         bevp1= bevp2+1;
1826
1827                         x1= bevp1->x- bevp2->x;
1828                         y1= bevp1->y- bevp2->y;
1829
1830                         calc_bevel_sin_cos(x1, y1, -x1, -y1, &(bevp1->sina), &(bevp1->cosa));
1831                         bevp2->sina= bevp1->sina;
1832                         bevp2->cosa= bevp1->cosa;
1833
1834                         if(cu->flag & CU_3D) {  /* 3D */
1835                                 float *quat, q[4];
1836                         
1837                                 vec[0]= bevp1->x - bevp2->x;
1838                                 vec[1]= bevp1->y - bevp2->y;
1839                                 vec[2]= bevp1->z - bevp2->z;
1840                                 
1841                                 quat= vectoquat(vec, 5, 1);
1842                                 
1843                                 Normalise(vec);
1844                                 q[0]= (float)cos(0.5*bevp1->alfa);
1845                                 x1= (float)sin(0.5*bevp1->alfa);
1846                                 q[1]= x1*vec[0];
1847                                 q[2]= x1*vec[1];
1848                                 q[3]= x1*vec[2];
1849                                 QuatMul(quat, q, quat);
1850                                 
1851                                 QuatToMat3(quat, bevp1->mat);
1852                                 Mat3CpyMat3(bevp2->mat, bevp1->mat);
1853                         }
1854
1855                 }
1856                 else if(bl->nr>2) {
1857                         bevp2= (BevPoint *)(bl+1);
1858                         bevp1= bevp2+(bl->nr-1);
1859                         bevp0= bevp1-1;
1860
1861                 
1862                         nr= bl->nr;
1863         
1864                         while(nr--) {
1865         
1866                                 if(cu->flag & CU_3D) {  /* 3D */
1867                                         float *quat, q[4];
1868                                 
1869                                         vec[0]= bevp2->x - bevp0->x;
1870                                         vec[1]= bevp2->y - bevp0->y;
1871                                         vec[2]= bevp2->z - bevp0->z;
1872                                         
1873                                         Normalise(vec);
1874
1875                                         quat= vectoquat(vec, 5, 1);
1876                                         
1877                                         q[0]= (float)cos(0.5*bevp1->alfa);
1878                                         x1= (float)sin(0.5*bevp1->alfa);
1879                                         q[1]= x1*vec[0];
1880                                         q[2]= x1*vec[1];
1881                                         q[3]= x1*vec[2];
1882                                         QuatMul(quat, q, quat);
1883                                         
1884                                         QuatToMat3(quat, bevp1->mat);
1885                                 }
1886                                 
1887                                 x1= bevp1->x- bevp0->x;
1888                                 x2= bevp1->x- bevp2->x;
1889                                 y1= bevp1->y- bevp0->y;
1890                                 y2= bevp1->y- bevp2->y;
1891                         
1892                                 calc_bevel_sin_cos(x1, y1, x2, y2, &(bevp1->sina), &(bevp1->cosa));
1893                                 
1894                                 
1895                                 bevp0= bevp1;
1896                                 bevp1= bevp2;
1897                                 bevp2++;
1898                         }
1899                         /* correct non-cyclic cases */
1900                         if(bl->poly== -1) {
1901                                 if(bl->nr>2) {
1902                                         bevp= (BevPoint *)(bl+1);
1903                                         bevp1= bevp+1;
1904                                         bevp->sina= bevp1->sina;
1905                                         bevp->cosa= bevp1->cosa;
1906                                         Mat3CpyMat3(bevp->mat, bevp1->mat);
1907                                         bevp= (BevPoint *)(bl+1);
1908                                         bevp+= (bl->nr-1);
1909                                         bevp1= bevp-1;
1910                                         bevp->sina= bevp1->sina;
1911                                         bevp->cosa= bevp1->cosa;
1912                                         Mat3CpyMat3(bevp->mat, bevp1->mat);
1913                                 }
1914                         }
1915                 }
1916                 bl= bl->next;
1917         }
1918 }
1919
1920 /* ****************** HANDLES ************** */
1921
1922 /*
1923  *   handlecodes:
1924  *              1: nothing,  1:auto,  2:vector,  3:aligned
1925  */
1926
1927 /* mode: is not zero when IpoCurve, is 2 when forced horizontal for autohandles */
1928 void calchandleNurb(BezTriple *bezt, BezTriple *prev, BezTriple *next, int mode)
1929 {
1930         float *p1,*p2,*p3,pt[3];
1931         float dx1,dy1,dz1,dx,dy,dz,vx,vy,vz,len,len1,len2;
1932
1933         if(bezt->h1==0 && bezt->h2==0) return;
1934
1935         p2= bezt->vec[1];
1936
1937         if(prev==0) {
1938                 p3= next->vec[1];
1939                 pt[0]= 2*p2[0]- p3[0];
1940                 pt[1]= 2*p2[1]- p3[1];
1941                 pt[2]= 2*p2[2]- p3[2];
1942                 p1= pt;
1943         }
1944         else p1= prev->vec[1];
1945
1946         if(next==0) {
1947                 pt[0]= 2*p2[0]- p1[0];
1948                 pt[1]= 2*p2[1]- p1[1];
1949                 pt[2]= 2*p2[2]- p1[2];
1950                 p3= pt;
1951         }
1952         else p3= next->vec[1];
1953
1954         if(mode && bezt->h1==HD_AUTO && prev) {
1955                 dx= p2[0] - (p1[0]+p1[3])/2.0f;
1956                 dy= p2[1] - (p1[1]+p1[4])/2.0f;
1957                 dz= p2[2] - (p1[2]+p1[5])/2.0f;
1958         }
1959         else {
1960                 dx= p2[0]- p1[0];
1961                 dy= p2[1]- p1[1];
1962                 dz= p2[2]- p1[2];
1963         }
1964         len1= (float)sqrt(dx*dx+dy*dy+dz*dz);
1965         
1966         if(mode && bezt->h2==HD_AUTO && next) {
1967                 dx1= (p3[0]+p3[-3])/2.0f - p2[0];
1968                 dy1= (p3[1]+p3[-2])/2.0f - p2[1];
1969                 dz1= (p3[2]+p3[-1])/2.0f - p2[2];
1970         }
1971         else {
1972                 dx1= p3[0]- p2[0];
1973                 dy1= p3[1]- p2[1];
1974                 dz1= p3[2]- p2[2];
1975         }
1976         len2= (float)sqrt(dx1*dx1+dy1*dy1+dz1*dz1);
1977
1978         if(len1==0.0f) len1=1.0f;
1979         if(len2==0.0f) len2=1.0f;
1980
1981
1982         if(bezt->h1==HD_AUTO || bezt->h2==HD_AUTO) {    /* auto */
1983                 vx= dx1/len2 + dx/len1;
1984                 vy= dy1/len2 + dy/len1;
1985                 vz= dz1/len2 + dz/len1;
1986                 len= 2.5614f*(float)sqrt(vx*vx + vy*vy + vz*vz);
1987                 if(len!=0.0f) {
1988                 
1989                         if(len1>5.0f*len2) len1= 5.0f*len2;     
1990                         if(len2>5.0f*len1) len2= 5.0f*len1;
1991                         
1992                         if(bezt->h1==HD_AUTO) {
1993                                 len1/=len;
1994                                 *(p2-3)= *p2-vx*len1;
1995                                 *(p2-2)= *(p2+1)-vy*len1;
1996                                 *(p2-1)= *(p2+2)-vz*len1;
1997                                 
1998                                 if(mode==2 && next && prev) {   // keep horizontal if extrema
1999                                         float ydiff1= prev->vec[1][1] - bezt->vec[1][1];
2000                                         float ydiff2= next->vec[1][1] - bezt->vec[1][1];
2001                                         if( (ydiff1<=0.0 && ydiff2<=0.0) || (ydiff1>=0.0 && ydiff2>=0.0) ) {
2002                                                 bezt->vec[0][1]= bezt->vec[1][1];
2003                                         }
2004                                 }
2005                         }
2006                         if(bezt->h2==HD_AUTO) {
2007                                 len2/=len;
2008                                 *(p2+3)= *p2+vx*len2;
2009                                 *(p2+4)= *(p2+1)+vy*len2;
2010                                 *(p2+5)= *(p2+2)+vz*len2;
2011                                 
2012                                 if(mode==2 && next && prev) {   // keep horizontal if extrema
2013                                         float ydiff1= prev->vec[1][1] - bezt->vec[1][1];
2014                                         float ydiff2= next->vec[1][1] - bezt->vec[1][1];
2015                                         if( (ydiff1<=0.0 && ydiff2<=0.0) || (ydiff1>=0.0 && ydiff2>=0.0) ) {
2016                                                 bezt->vec[2][1]= bezt->vec[1][1];
2017                                         }
2018                                 }
2019                         }
2020                         
2021                 }
2022         }
2023
2024         if(bezt->h1==HD_VECT) { /* vector */
2025                 dx/=3.0; 
2026                 dy/=3.0; 
2027                 dz/=3.0;
2028                 *(p2-3)= *p2-dx;
2029                 *(p2-2)= *(p2+1)-dy;
2030                 *(p2-1)= *(p2+2)-dz;
2031         }
2032         if(bezt->h2==HD_VECT) {
2033                 dx1/=3.0; 
2034                 dy1/=3.0; 
2035                 dz1/=3.0;
2036                 *(p2+3)= *p2+dx1;
2037                 *(p2+4)= *(p2+1)+dy1;
2038                 *(p2+5)= *(p2+2)+dz1;
2039         }
2040
2041         len2= VecLenf(p2, p2+3);
2042         len1= VecLenf(p2, p2-3);
2043         if(len1==0.0) len1=1.0;
2044         if(len2==0.0) len2=1.0;
2045
2046         if(bezt->f1 & 1) { /* order of calculation */
2047                 if(bezt->h2==HD_ALIGN) {        /* aligned */
2048                         len= len2/len1;
2049                         p2[3]= p2[0]+len*(p2[0]-p2[-3]);
2050                         p2[4]= p2[1]+len*(p2[1]-p2[-2]);
2051                         p2[5]= p2[2]+len*(p2[2]-p2[-1]);
2052                 }
2053                 if(bezt->h1==HD_ALIGN) {
2054                         len= len1/len2;
2055                         p2[-3]= p2[0]+len*(p2[0]-p2[3]);
2056                         p2[-2]= p2[1]+len*(p2[1]-p2[4]);
2057                         p2[-1]= p2[2]+len*(p2[2]-p2[5]);
2058                 }
2059         }
2060         else {
2061                 if(bezt->h1==HD_ALIGN) {
2062                         len= len1/len2;
2063                         p2[-3]= p2[0]+len*(p2[0]-p2[3]);
2064                         p2[-2]= p2[1]+len*(p2[1]-p2[4]);
2065                         p2[-1]= p2[2]+len*(p2[2]-p2[5]);
2066                 }
2067                 if(bezt->h2==HD_ALIGN) {        /* aligned */
2068                         len= len2/len1;
2069                         p2[3]= p2[0]+len*(p2[0]-p2[-3]);
2070                         p2[4]= p2[1]+len*(p2[1]-p2[-2]);
2071                         p2[5]= p2[2]+len*(p2[2]-p2[-1]);
2072                 }
2073         }
2074 }
2075
2076 void calchandlesNurb(Nurb *nu) /* first, if needed, set handle flags */
2077 {
2078         BezTriple *bezt, *prev, *next;
2079         short a;
2080
2081         if((nu->type & 7)!=1) return;
2082         if(nu->pntsu<2) return;
2083         
2084         a= nu->pntsu;
2085         bezt= nu->bezt;
2086         if(nu->flagu & 1) prev= bezt+(a-1);
2087         else prev= 0;
2088         next= bezt+1;
2089
2090         while(a--) {
2091                 calchandleNurb(bezt, prev, next, 0);
2092                 prev= bezt;
2093                 if(a==1) {
2094                         if(nu->flagu & 1) next= nu->bezt;
2095                         else next= 0;
2096                 }
2097                 else next++;
2098
2099                 bezt++;
2100         }
2101 }
2102
2103
2104 void testhandlesNurb(Nurb *nu)
2105 {
2106     /* use when something has changed with handles.
2107     it treats all BezTriples with the following rules:
2108     PHASE 1: do types have to be altered?
2109        Auto handles: become aligned when selection status is NOT(000 || 111)
2110        Vector handles: become 'nothing' when (one half selected AND other not)
2111     PHASE 2: recalculate handles
2112     */
2113         BezTriple *bezt;
2114         short flag, a;
2115
2116         if((nu->type & 7)!=CU_BEZIER) return;
2117
2118         bezt= nu->bezt;
2119         a= nu->pntsu;
2120         while(a--) {
2121                 flag= 0;
2122                 if(bezt->f1 & 1) flag++;
2123                 if(bezt->f2 & 1) flag += 2;
2124                 if(bezt->f3 & 1) flag += 4;
2125
2126                 if( !(flag==0 || flag==7) ) {
2127                         if(bezt->h1==HD_AUTO) {   /* auto */
2128                                 bezt->h1= HD_ALIGN;
2129                         }
2130                         if(bezt->h2==HD_AUTO) {   /* auto */
2131                                 bezt->h2= HD_ALIGN;
2132                         }
2133
2134                         if(bezt->h1==HD_VECT) {   /* vector */
2135                                 if(flag < 4) bezt->h1= 0;
2136                         }
2137                         if(bezt->h2==HD_VECT) {   /* vector */
2138                                 if( flag > 3) bezt->h2= 0;
2139                         }
2140                 }
2141                 bezt++;
2142         }
2143
2144         calchandlesNurb(nu);
2145 }
2146
2147 void autocalchandlesNurb(Nurb *nu, int flag)
2148 {
2149         /* checks handle coordinates and calculates type */
2150         
2151         BezTriple *bezt2, *bezt1, *bezt0;
2152         int i, align, leftsmall, rightsmall;
2153
2154         if(nu==0 || nu->bezt==0) return;
2155         
2156         bezt2 = nu->bezt;
2157         bezt1 = bezt2 + (nu->pntsu-1);
2158         bezt0 = bezt1 - 1;
2159         i = nu->pntsu;
2160
2161         while(i--) {
2162                 
2163                 align= leftsmall= rightsmall= 0;
2164                 
2165                 /* left handle: */
2166                 if(flag==0 || (bezt1->f1 & flag) ) {
2167                         bezt1->h1= 0;
2168                         /* distance too short: vectorhandle */
2169                         if( VecLenf( bezt1->vec[1], bezt0->vec[1] ) < 0.0001) {
2170                                 bezt1->h1= HD_VECT;
2171                                 leftsmall= 1;
2172                         }
2173                         else {
2174                                 /* aligned handle? */
2175                                 if(DistVL2Dfl(bezt1->vec[1], bezt1->vec[0], bezt1->vec[2]) < 0.0001) {
2176                                         align= 1;
2177                                         bezt1->h1= HD_ALIGN;
2178                                 }
2179                                 /* or vector handle? */
2180                                 if(DistVL2Dfl(bezt1->vec[0], bezt1->vec[1], bezt0->vec[1]) < 0.0001)
2181                                         bezt1->h1= HD_VECT;
2182                                 
2183                         }
2184                 }
2185                 /* right handle: */
2186                 if(flag==0 || (bezt1->f3 & flag) ) {
2187                         bezt1->h2= 0;
2188                         /* distance too short: vectorhandle */
2189                         if( VecLenf( bezt1->vec[1], bezt2->vec[1] ) < 0.0001) {
2190                                 bezt1->h2= HD_VECT;
2191                                 rightsmall= 1;
2192                         }
2193                         else {
2194                                 /* aligned handle? */
2195                                 if(align) bezt1->h2= HD_ALIGN;
2196
2197                                 /* or vector handle? */
2198                                 if(DistVL2Dfl(bezt1->vec[2], bezt1->vec[1], bezt2->vec[1]) < 0.0001)
2199                                         bezt1->h2= HD_VECT;
2200                                 
2201                         }
2202                 }
2203                 if(leftsmall && bezt1->h2==HD_ALIGN) bezt1->h2= 0;
2204                 if(rightsmall && bezt1->h1==HD_ALIGN) bezt1->h1= 0;
2205                 
2206                 /* undesired combination: */
2207                 if(bezt1->h1==HD_ALIGN && bezt1->h2==HD_VECT) bezt1->h1= 0;
2208                 if(bezt1->h2==HD_ALIGN && bezt1->h1==HD_VECT) bezt1->h2= 0;
2209                 
2210                 bezt0= bezt1;
2211                 bezt1= bezt2;
2212                 bezt2++;
2213         }
2214
2215         calchandlesNurb(nu);
2216 }
2217
2218 void autocalchandlesNurb_all(int flag)
2219 {
2220         Nurb *nu;
2221         
2222         nu= editNurb.first;
2223         while(nu) {
2224                 autocalchandlesNurb(nu, flag);
2225                 nu= nu->next;
2226         }
2227 }
2228
2229 void sethandlesNurb(short code)
2230 {
2231         /* code==1: set autohandle */
2232         /* code==2: set vectorhandle */
2233         /* code==3 (HD_ALIGN) it toggle, vectorhandles become HD_FREE */
2234         /* code==4: sets icu flag to become IPO_AUTO_HORIZ, horizontal extremes on auto-handles */
2235         Nurb *nu;
2236         BezTriple *bezt;
2237         short a, ok=0;
2238
2239         if(code==1 || code==2) {
2240                 nu= editNurb.first;
2241                 while(nu) {
2242                         if( (nu->type & 7)==1) {
2243                                 bezt= nu->bezt;
2244                                 a= nu->pntsu;
2245                                 while(a--) {
2246                                         if(bezt->f1 || bezt->f3) {
2247                                                 if(bezt->f1) bezt->h1= code;
2248                                                 if(bezt->f3) bezt->h2= code;
2249                                                 if(bezt->h1!=bezt->h2) {
2250                                                         if ELEM(bezt->h1, HD_ALIGN, HD_AUTO) bezt->h1= HD_FREE;
2251                                                         if ELEM(bezt->h2, HD_ALIGN, HD_AUTO) bezt->h2= HD_FREE;
2252                                                 }
2253                                         }
2254                                         bezt++;
2255                                 }
2256                                 calchandlesNurb(nu);
2257                         }
2258                         nu= nu->next;
2259                 }
2260         }
2261         else {
2262                 /* there is 1 handle not FREE: FREE it all, else make ALIGNED  */
2263                 
2264                 nu= editNurb.first;
2265                 while(nu) {
2266                         if( (nu->type & 7)==1) {
2267                                 bezt= nu->bezt;
2268                                 a= nu->pntsu;
2269                                 while(a--) {
2270                                         if(bezt->f1 && bezt->h1) ok= 1;
2271                                         if(bezt->f3 && bezt->h2) ok= 1;
2272                                         if(ok) break;
2273                                         bezt++;
2274                                 }
2275                         }
2276                         nu= nu->next;
2277                 }
2278                 if(ok) ok= HD_FREE;
2279                 else ok= HD_ALIGN;
2280                 
2281                 nu= editNurb.first;
2282                 while(nu) {
2283                         if( (nu->type & 7)==1) {
2284                                 bezt= nu->bezt;
2285                                 a= nu->pntsu;
2286                                 while(a--) {
2287                                         if(bezt->f1) bezt->h1= ok;
2288                                         if(bezt->f3 ) bezt->h2= ok;
2289         
2290                                         bezt++;
2291                                 }
2292                                 calchandlesNurb(nu);
2293                         }
2294                         nu= nu->next;
2295                 }
2296         }
2297 }
2298
2299 void swapdata(void *adr1, void *adr2, int len)
2300 {
2301
2302         if(len<=0) return;
2303
2304         if(len<65) {
2305                 char adr[64];
2306
2307                 memcpy(adr, adr1, len);
2308                 memcpy(adr1, adr2, len);
2309                 memcpy(adr2, adr, len);
2310         }
2311         else {
2312                 char *adr;
2313
2314                 adr= (char *)MEM_mallocN(len, "curve swap");
2315                 memcpy(adr, adr1, len);
2316                 memcpy(adr1, adr2, len);
2317                 memcpy(adr2, adr, len);
2318                 MEM_freeN(adr);
2319         }
2320 }
2321
2322 void switchdirectionNurb(Nurb *nu)
2323 {
2324         BezTriple *bezt1, *bezt2;
2325         BPoint *bp1, *bp2;
2326         float *fp1, *fp2, *tempf;
2327         int a, b;
2328
2329         if(nu->pntsu==1 && nu->pntsv==1) return;
2330
2331         if((nu->type & 7)==CU_BEZIER) {
2332                 a= nu->pntsu;
2333                 bezt1= nu->bezt;
2334                 bezt2= bezt1+(a-1);
2335                 if(a & 1) a+= 1;        /* if odd, also swap middle content */
2336                 a/= 2;
2337                 while(a>0) {
2338                         if(bezt1!=bezt2) SWAP(BezTriple, *bezt1, *bezt2);
2339
2340                         swapdata(bezt1->vec[0], bezt1->vec[2], 12);
2341                         if(bezt1!=bezt2) swapdata(bezt2->vec[0], bezt2->vec[2], 12);
2342
2343                         SWAP(char, bezt1->h1, bezt1->h2);
2344                         SWAP(short, bezt1->f1, bezt1->f3);
2345                         
2346                         if(bezt1!=bezt2) {
2347                                 SWAP(char, bezt2->h1, bezt2->h2);
2348                                 SWAP(short, bezt2->f1, bezt2->f3);
2349                                 bezt1->alfa= -bezt1->alfa;
2350                                 bezt2->alfa= -bezt2->alfa;
2351                         }
2352                         a--;
2353                         bezt1++; 
2354                         bezt2--;
2355                 }
2356         }
2357         else if(nu->pntsv==1) {
2358                 a= nu->pntsu;
2359                 bp1= nu->bp;
2360                 bp2= bp1+(a-1);
2361                 a/= 2;
2362                 while(bp1!=bp2 && a>0) {
2363                         SWAP(BPoint, *bp1, *bp2);
2364                         a--;
2365                         bp1->alfa= -bp1->alfa;
2366                         bp2->alfa= -bp2->alfa;
2367                         bp1++; 
2368                         bp2--;
2369                 }
2370                 if((nu->type & 7)==CU_NURBS) {
2371                         /* inverse knots */
2372                         a= KNOTSU(nu);
2373                         fp1= nu->knotsu;
2374                         fp2= fp1+(a-1);
2375                         a/= 2;
2376                         while(fp1!=fp2 && a>0) {
2377                                 SWAP(float, *fp1, *fp2);
2378                                 a--;
2379                                 fp1++; 
2380                                 fp2--;
2381                         }
2382                         /* and make in increasing order again */
2383                         a= KNOTSU(nu);
2384                         fp1= nu->knotsu;
2385                         fp2=tempf= MEM_mallocN(sizeof(float)*a, "switchdirect");
2386                         while(a--) {
2387                                 fp2[0]= fabs(fp1[1]-fp1[0]);
2388                                 fp1++;
2389                                 fp2++;
2390                         }
2391         
2392                         a= KNOTSU(nu)-1;
2393                         fp1= nu->knotsu;
2394                         fp2= tempf;
2395                         fp1[0]= 0.0;
2396                         fp1++;
2397                         while(a--) {
2398                                 fp1[0]= fp1[-1]+fp2[0];
2399                                 fp1++;
2400                                 fp2++;
2401                         }
2402                         MEM_freeN(tempf);
2403                 }
2404         }
2405         else {
2406                 
2407                 for(b=0; b<nu->pntsv; b++) {
2408                 
2409                         bp1= nu->bp+b*nu->pntsu;
2410                         a= nu->pntsu;
2411                         bp2= bp1+(a-1);
2412                         a/= 2;
2413                         
2414                         while(bp1!=bp2 && a>0) {
2415                                 SWAP(BPoint, *bp1, *bp2);
2416                                 a--;
2417                                 bp1++; 
2418                                 bp2--;
2419                         }
2420                 }
2421         }
2422 }
2423
2424
2425 float (*curve_getVertexCos(Curve *cu, ListBase *lb, int *numVerts_r))[3]
2426 {
2427         int i, numVerts = *numVerts_r = count_curveverts(lb);
2428         float *co, (*cos)[3] = MEM_mallocN(sizeof(*cos)*numVerts, "cu_vcos");
2429         Nurb *nu;
2430
2431         co = cos[0];
2432         for (nu=lb->first; nu; nu=nu->next) {
2433                 if ((nu->type & 7)==CU_BEZIER) {
2434                         BezTriple *bezt = nu->bezt;
2435
2436                         for (i=0; i<nu->pntsu; i++,bezt++) {
2437                                 VECCOPY(co, bezt->vec[0]); co+=3;
2438                                 VECCOPY(co, bezt->vec[1]); co+=3;
2439                                 VECCOPY(co, bezt->vec[2]); co+=3;
2440                         }
2441                 } else {
2442                         BPoint *bp = nu->bp;
2443
2444                         for (i=0; i<nu->pntsu*nu->pntsv; i++,bp++) {
2445                                 VECCOPY(co, bp->vec); co+=3;
2446                         }
2447                 }
2448         }
2449
2450         return cos;
2451 }
2452
2453 void curve_applyVertexCos(Curve *cu, ListBase *lb, float (*vertexCos)[3])
2454 {
2455         float *co = vertexCos[0];
2456         Nurb *nu;
2457         int i;
2458
2459         for (nu=lb->first; nu; nu=nu->next) {
2460                 if ((nu->type & 7)==CU_BEZIER) {
2461                         BezTriple *bezt = nu->bezt;
2462
2463                         for (i=0; i<nu->pntsu; i++,bezt++) {
2464                                 VECCOPY(bezt->vec[0], co); co+=3;
2465                                 VECCOPY(bezt->vec[1], co); co+=3;
2466                                 VECCOPY(bezt->vec[2], co); co+=3;
2467                         }
2468                 } else {
2469                         BPoint *bp = nu->bp;
2470
2471                         for (i=0; i<nu->pntsu*nu->pntsv; i++,bp++) {
2472                                 VECCOPY(bp->vec, co); co+=3;
2473                         }
2474                 }
2475         }
2476 }