[#21331] Pose snap not working in some cases
[blender.git] / source / blender / blenkernel / intern / armature.c
1 /**
2  * $Id$
3  *
4  * ***** BEGIN GPL LICENSE BLOCK *****
5  *
6  * This program is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU General Public License
8  * as published by the Free Software Foundation; either version 2
9  * of the License, or (at your option) any later version. 
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software Foundation,
18  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19  *
20  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
21  * All rights reserved.
22  *
23  * Contributor(s): Full recode, Ton Roosendaal, Crete 2005
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 #include <ctype.h>
29 #include <stdlib.h>
30 #include <math.h>
31 #include <string.h>
32 #include <stdio.h>
33 #include <float.h>
34
35 #include "MEM_guardedalloc.h"
36
37 #include "BLI_math.h"
38 #include "BLI_blenlib.h"
39
40 #include "DNA_anim_types.h"
41 #include "DNA_armature_types.h"
42 #include "DNA_constraint_types.h"
43 #include "DNA_mesh_types.h"
44 #include "DNA_lattice_types.h"
45 #include "DNA_meshdata_types.h"
46 #include "DNA_nla_types.h"
47 #include "DNA_scene_types.h"
48 #include "DNA_object_types.h"
49
50 #include "BKE_animsys.h"
51 #include "BKE_armature.h"
52 #include "BKE_action.h"
53 #include "BKE_anim.h"
54 #include "BKE_constraint.h"
55 #include "BKE_curve.h"
56 #include "BKE_depsgraph.h"
57 #include "BKE_DerivedMesh.h"
58 #include "BKE_deform.h"
59 #include "BKE_displist.h"
60 #include "BKE_global.h"
61 #include "BKE_idprop.h"
62 #include "BKE_library.h"
63 #include "BKE_lattice.h"
64 #include "BKE_main.h"
65 #include "BKE_object.h"
66 #include "BKE_utildefines.h"
67 #include "BIK_api.h"
68 #include "BKE_sketch.h"
69
70 /*      **************** Generic Functions, data level *************** */
71
72 bArmature *add_armature(char *name)
73 {
74         bArmature *arm;
75         
76         arm= alloc_libblock (&G.main->armature, ID_AR, name);
77         arm->deformflag = ARM_DEF_VGROUP|ARM_DEF_ENVELOPE;
78         arm->flag = ARM_COL_CUSTOM; /* custom bone-group colors */
79         arm->layer= 1;
80         return arm;
81 }
82
83 bArmature *get_armature(Object *ob)
84 {
85         if(ob->type==OB_ARMATURE)
86                 return (bArmature *)ob->data;
87         return NULL;
88 }
89
90 void free_bonelist (ListBase *lb)
91 {
92         Bone *bone;
93
94         for(bone=lb->first; bone; bone=bone->next) {
95                 if(bone->prop) {
96                         IDP_FreeProperty(bone->prop);
97                         MEM_freeN(bone->prop);
98                 }
99                 free_bonelist(&bone->childbase);
100         }
101         
102         BLI_freelistN(lb);
103 }
104
105 void free_armature(bArmature *arm)
106 {
107         if (arm) {
108                 free_bonelist(&arm->bonebase);
109                 
110                 /* free editmode data */
111                 if (arm->edbo) {
112                         BLI_freelistN(arm->edbo);
113                         
114                         MEM_freeN(arm->edbo);
115                         arm->edbo= NULL;
116                 }
117
118                 /* free sketch */
119                 if (arm->sketch) {
120                         freeSketch(arm->sketch);
121                         arm->sketch = NULL;
122                 }
123                 
124                 /* free animation data */
125                 if (arm->adt) {
126                         BKE_free_animdata(&arm->id);
127                         arm->adt= NULL;
128                 }
129         }
130 }
131
132 void make_local_armature(bArmature *arm)
133 {
134         int local=0, lib=0;
135         Object *ob;
136         bArmature *newArm;
137         
138         if (arm->id.lib==0)
139                 return;
140         if (arm->id.us==1) {
141                 arm->id.lib= 0;
142                 arm->id.flag= LIB_LOCAL;
143                 new_id(0, (ID*)arm, 0);
144                 return;
145         }
146         
147         if(local && lib==0) {
148                 arm->id.lib= 0;
149                 arm->id.flag= LIB_LOCAL;
150                 new_id(0, (ID *)arm, 0);
151         }
152         else if(local && lib) {
153                 newArm= copy_armature(arm);
154                 newArm->id.us= 0;
155                 
156                 ob= G.main->object.first;
157                 while(ob) {
158                         if(ob->data==arm) {
159                                 
160                                 if(ob->id.lib==0) {
161                                         ob->data= newArm;
162                                         newArm->id.us++;
163                                         arm->id.us--;
164                                 }
165                         }
166                         ob= ob->id.next;
167                 }
168         }
169 }
170
171 static void     copy_bonechildren (Bone* newBone, Bone* oldBone, Bone* actBone, Bone **newActBone)
172 {
173         Bone    *curBone, *newChildBone;
174         
175         if(oldBone == actBone)
176                 *newActBone= newBone;
177
178         if(oldBone->prop)
179                 newBone->prop= IDP_CopyProperty(oldBone->prop);
180
181         /*      Copy this bone's list*/
182         BLI_duplicatelist(&newBone->childbase, &oldBone->childbase);
183         
184         /*      For each child in the list, update it's children*/
185         newChildBone=newBone->childbase.first;
186         for (curBone=oldBone->childbase.first;curBone;curBone=curBone->next){
187                 newChildBone->parent=newBone;
188                 copy_bonechildren(newChildBone, curBone, actBone, newActBone);
189                 newChildBone=newChildBone->next;
190         }
191 }
192
193 bArmature *copy_armature(bArmature *arm)
194 {
195         bArmature *newArm;
196         Bone            *oldBone, *newBone;
197         Bone            *newActBone= NULL;
198         
199         newArm= copy_libblock (arm);
200         BLI_duplicatelist(&newArm->bonebase, &arm->bonebase);
201         
202         /*      Duplicate the childrens' lists*/
203         newBone=newArm->bonebase.first;
204         for (oldBone=arm->bonebase.first;oldBone;oldBone=oldBone->next){
205                 newBone->parent=NULL;
206                 copy_bonechildren (newBone, oldBone, arm->act_bone, &newActBone);
207                 newBone=newBone->next;
208         };
209         
210         newArm->act_bone= newActBone;
211         return newArm;
212 }
213
214 static Bone *get_named_bone_bonechildren (Bone *bone, const char *name)
215 {
216         Bone *curBone, *rbone;
217         
218         if (!strcmp (bone->name, name))
219                 return bone;
220         
221         for (curBone=bone->childbase.first; curBone; curBone=curBone->next){
222                 rbone=get_named_bone_bonechildren (curBone, name);
223                 if (rbone)
224                         return rbone;
225         }
226         
227         return NULL;
228 }
229
230
231 Bone *get_named_bone (bArmature *arm, const char *name)
232 /*
233         Walk the list until the bone is found
234  */
235 {
236         Bone *bone=NULL, *curBone;
237         
238         if (!arm) return NULL;
239         
240         for (curBone=arm->bonebase.first; curBone; curBone=curBone->next){
241                 bone = get_named_bone_bonechildren (curBone, name);
242                 if (bone)
243                         return bone;
244         }
245         
246         return bone;
247 }
248
249 /* Finds the best possible extension to the name on a particular axis. (For renaming, check for unique names afterwards)
250  * This assumes that bone names are at most 32 chars long!
251  *      strip_number: removes number extensions  (TODO: not used)
252  *      axis: the axis to name on
253  *      head/tail: the head/tail co-ordinate of the bone on the specified axis
254  */
255 int bone_autoside_name (char *name, int UNUSED(strip_number), short axis, float head, float tail)
256 {
257         unsigned int len;
258         char    basename[32]={""};
259         char    extension[5]={""};
260
261         len= strlen(name);
262         if (len == 0) return 0;
263         strcpy(basename, name);
264         
265         /* Figure out extension to append: 
266          *      - The extension to append is based upon the axis that we are working on.
267          *      - If head happens to be on 0, then we must consider the tail position as well to decide
268          *        which side the bone is on
269          *              -> If tail is 0, then it's bone is considered to be on axis, so no extension should be added
270          *              -> Otherwise, extension is added from perspective of object based on which side tail goes to
271          *      - If head is non-zero, extension is added from perspective of object based on side head is on
272          */
273         if (axis == 2) {
274                 /* z-axis - vertical (top/bottom) */
275                 if (IS_EQ(head, 0)) {
276                         if (tail < 0)
277                                 strcpy(extension, "Bot");
278                         else if (tail > 0)
279                                 strcpy(extension, "Top");
280                 }
281                 else {
282                         if (head < 0)
283                                 strcpy(extension, "Bot");
284                         else
285                                 strcpy(extension, "Top");
286                 }
287         }
288         else if (axis == 1) {
289                 /* y-axis - depth (front/back) */
290                 if (IS_EQ(head, 0)) {
291                         if (tail < 0)
292                                 strcpy(extension, "Fr");
293                         else if (tail > 0)
294                                 strcpy(extension, "Bk");
295                 }
296                 else {
297                         if (head < 0)
298                                 strcpy(extension, "Fr");
299                         else
300                                 strcpy(extension, "Bk");
301                 }
302         }
303         else {
304                 /* x-axis - horizontal (left/right) */
305                 if (IS_EQ(head, 0)) {
306                         if (tail < 0)
307                                 strcpy(extension, "R");
308                         else if (tail > 0)
309                                 strcpy(extension, "L");
310                 }
311                 else {
312                         if (head < 0)
313                                 strcpy(extension, "R");
314                         else if (head > 0)
315                                 strcpy(extension, "L");
316                 }
317         }
318
319         /* Simple name truncation 
320          *      - truncate if there is an extension and it wouldn't be able to fit
321          *      - otherwise, just append to end
322          */
323         if (extension[0]) {
324                 int change = 1;
325                 
326                 while (change) { /* remove extensions */
327                         change = 0;
328                         if (len > 2 && basename[len-2]=='.') {
329                                 if (basename[len-1]=='L' || basename[len-1] == 'R' ) { /* L R */
330                                         basename[len-2] = '\0';
331                                         len-=2;
332                                         change= 1;
333                                 }
334                         } else if (len > 3 && basename[len-3]=='.') {
335                                 if (    (basename[len-2]=='F' && basename[len-1] == 'r') ||     /* Fr */
336                                                 (basename[len-2]=='B' && basename[len-1] == 'k')        /* Bk */
337                                 ) {
338                                         basename[len-3] = '\0';
339                                         len-=3;
340                                         change= 1;
341                                 }
342                         } else if (len > 4 && basename[len-4]=='.') {
343                                 if (    (basename[len-3]=='T' && basename[len-2]=='o' && basename[len-1] == 'p') ||     /* Top */
344                                                 (basename[len-3]=='B' && basename[len-2]=='o' && basename[len-1] == 't')        /* Bot */
345                                 ) {
346                                         basename[len-4] = '\0';
347                                         len-=4;
348                                         change= 1;
349                                 }
350                         }
351                 }
352                 
353                 if ((32 - len) < strlen(extension) + 1) { /* add 1 for the '.' */
354                         strncpy(name, basename, len-strlen(extension));
355                 }
356                 
357                 sprintf(name, "%s.%s", basename, extension);
358                 
359                 return 1;
360         }
361
362         else {
363                 return 0;
364         }
365 }
366
367 /* ************* B-Bone support ******************* */
368
369 #define MAX_BBONE_SUBDIV        32
370
371 /* data has MAX_BBONE_SUBDIV+1 interpolated points, will become desired amount with equal distances */
372 static void equalize_bezier(float *data, int desired)
373 {
374         float *fp, totdist, ddist, dist, fac1, fac2;
375         float pdist[MAX_BBONE_SUBDIV+1];
376         float temp[MAX_BBONE_SUBDIV+1][4];
377         int a, nr;
378         
379         pdist[0]= 0.0f;
380         for(a=0, fp= data; a<MAX_BBONE_SUBDIV; a++, fp+=4) {
381                 QUATCOPY(temp[a], fp);
382                 pdist[a+1]= pdist[a]+len_v3v3(fp, fp+4);
383         }
384         /* do last point */
385         QUATCOPY(temp[a], fp);
386         totdist= pdist[a];
387         
388         /* go over distances and calculate new points */
389         ddist= totdist/((float)desired);
390         nr= 1;
391         for(a=1, fp= data+4; a<desired; a++, fp+=4) {
392                 
393                 dist= ((float)a)*ddist;
394                 
395                 /* we're looking for location (distance) 'dist' in the array */
396                 while((dist>= pdist[nr]) && nr<MAX_BBONE_SUBDIV) {
397                         nr++;
398                 }
399                 
400                 fac1= pdist[nr]- pdist[nr-1];
401                 fac2= pdist[nr]-dist;
402                 fac1= fac2/fac1;
403                 fac2= 1.0f-fac1;
404                 
405                 fp[0]= fac1*temp[nr-1][0]+ fac2*temp[nr][0];
406                 fp[1]= fac1*temp[nr-1][1]+ fac2*temp[nr][1];
407                 fp[2]= fac1*temp[nr-1][2]+ fac2*temp[nr][2];
408                 fp[3]= fac1*temp[nr-1][3]+ fac2*temp[nr][3];
409         }
410         /* set last point, needed for orientation calculus */
411         QUATCOPY(fp, temp[MAX_BBONE_SUBDIV]);
412 }
413
414 /* returns pointer to static array, filled with desired amount of bone->segments elements */
415 /* this calculation is done  within unit bone space */
416 Mat4 *b_bone_spline_setup(bPoseChannel *pchan, int rest)
417 {
418         static Mat4 bbone_array[MAX_BBONE_SUBDIV];
419         static Mat4 bbone_rest_array[MAX_BBONE_SUBDIV];
420         Mat4 *result_array= (rest)? bbone_rest_array: bbone_array;
421         bPoseChannel *next, *prev;
422         Bone *bone= pchan->bone;
423         float h1[3], h2[3], scale[3], length, hlength1, hlength2, roll1=0.0f, roll2;
424         float mat3[3][3], imat[4][4], posemat[4][4], scalemat[4][4], iscalemat[4][4];
425         float data[MAX_BBONE_SUBDIV+1][4], *fp;
426         int a, doscale= 0;
427
428         length= bone->length;
429
430         if(!rest) {
431                 /* check if we need to take non-uniform bone scaling into account */
432                 scale[0]= len_v3(pchan->pose_mat[0]);
433                 scale[1]= len_v3(pchan->pose_mat[1]);
434                 scale[2]= len_v3(pchan->pose_mat[2]);
435
436                 if(fabs(scale[0] - scale[1]) > 1e-6f || fabs(scale[1] - scale[2]) > 1e-6f) {
437                         unit_m4(scalemat);
438                         scalemat[0][0]= scale[0];
439                         scalemat[1][1]= scale[1];
440                         scalemat[2][2]= scale[2];
441                         invert_m4_m4(iscalemat, scalemat);
442
443                         length *= scale[1];
444                         doscale = 1;
445                 }
446         }
447         
448         hlength1= bone->ease1*length*0.390464f;         // 0.5*sqrt(2)*kappa, the handle length for near-perfect circles
449         hlength2= bone->ease2*length*0.390464f;
450         
451         /* evaluate next and prev bones */
452         if(bone->flag & BONE_CONNECTED)
453                 prev= pchan->parent;
454         else
455                 prev= NULL;
456         
457         next= pchan->child;
458         
459         /* find the handle points, since this is inside bone space, the 
460                 first point = (0,0,0)
461                 last point =  (0, length, 0) */
462         
463         if(rest) {
464                 invert_m4_m4(imat, pchan->bone->arm_mat);
465         }
466         else if(doscale) {
467                 copy_m4_m4(posemat, pchan->pose_mat);
468                 normalize_m4(posemat);
469                 invert_m4_m4(imat, posemat);
470         }
471         else
472                 invert_m4_m4(imat, pchan->pose_mat);
473         
474         if(prev) {
475                 float difmat[4][4], result[3][3], imat3[3][3];
476
477                 /* transform previous point inside this bone space */
478                 if(rest)
479                         VECCOPY(h1, prev->bone->arm_head)
480                 else
481                         VECCOPY(h1, prev->pose_head)
482                 mul_m4_v3(imat, h1);
483
484                 if(prev->bone->segments>1) {
485                         /* if previous bone is B-bone too, use average handle direction */
486                         h1[1]-= length;
487                         roll1= 0.0f;
488                 }
489
490                 normalize_v3(h1);
491                 mul_v3_fl(h1, -hlength1);
492
493                 if(prev->bone->segments==1) {
494                         /* find the previous roll to interpolate */
495                         if(rest)
496                                 mul_m4_m4m4(difmat, prev->bone->arm_mat, imat);
497                         else
498                                 mul_m4_m4m4(difmat, prev->pose_mat, imat);
499                         copy_m3_m4(result, difmat);                             // the desired rotation at beginning of next bone
500                         
501                         vec_roll_to_mat3(h1, 0.0f, mat3);                       // the result of vec_roll without roll
502                         
503                         invert_m3_m3(imat3, mat3);
504                         mul_m3_m3m3(mat3, result, imat3);                       // the matrix transforming vec_roll to desired roll
505                         
506                         roll1= (float)atan2(mat3[2][0], mat3[2][2]);
507                 }
508         }
509         else {
510                 h1[0]= 0.0f; h1[1]= hlength1; h1[2]= 0.0f;
511                 roll1= 0.0f;
512         }
513         if(next) {
514                 float difmat[4][4], result[3][3], imat3[3][3];
515                 
516                 /* transform next point inside this bone space */
517                 if(rest)
518                         VECCOPY(h2, next->bone->arm_tail)
519                 else
520                         VECCOPY(h2, next->pose_tail)
521                 mul_m4_v3(imat, h2);
522                 /* if next bone is B-bone too, use average handle direction */
523                 if(next->bone->segments>1);
524                 else h2[1]-= length;
525                 normalize_v3(h2);
526                 
527                 /* find the next roll to interpolate as well */
528                 if(rest)
529                         mul_m4_m4m4(difmat, next->bone->arm_mat, imat);
530                 else
531                         mul_m4_m4m4(difmat, next->pose_mat, imat);
532                 copy_m3_m4(result, difmat);                             // the desired rotation at beginning of next bone
533                 
534                 vec_roll_to_mat3(h2, 0.0f, mat3);                       // the result of vec_roll without roll
535                 
536                 invert_m3_m3(imat3, mat3);
537                 mul_m3_m3m3(mat3, imat3, result);                       // the matrix transforming vec_roll to desired roll
538                 
539                 roll2= (float)atan2(mat3[2][0], mat3[2][2]);
540                 
541                 /* and only now negate handle */
542                 mul_v3_fl(h2, -hlength2);
543         }
544         else {
545                 h2[0]= 0.0f; h2[1]= -hlength2; h2[2]= 0.0f;
546                 roll2= 0.0;
547         }
548
549         /* make curve */
550         if(bone->segments > MAX_BBONE_SUBDIV)
551                 bone->segments= MAX_BBONE_SUBDIV;
552         
553         forward_diff_bezier(0.0, h1[0],         h2[0],                  0.0,            data[0],        MAX_BBONE_SUBDIV, 4*sizeof(float));
554         forward_diff_bezier(0.0, h1[1],         length + h2[1], length,         data[0]+1,      MAX_BBONE_SUBDIV, 4*sizeof(float));
555         forward_diff_bezier(0.0, h1[2],         h2[2],                  0.0,            data[0]+2,      MAX_BBONE_SUBDIV, 4*sizeof(float));
556         forward_diff_bezier(roll1, roll1 + 0.390464f*(roll2-roll1), roll2 - 0.390464f*(roll2-roll1),    roll2,  data[0]+3,      MAX_BBONE_SUBDIV, 4*sizeof(float));
557         
558         equalize_bezier(data[0], bone->segments);       // note: does stride 4!
559         
560         /* make transformation matrices for the segments for drawing */
561         for(a=0, fp= data[0]; a<bone->segments; a++, fp+=4) {
562                 sub_v3_v3v3(h1, fp+4, fp);
563                 vec_roll_to_mat3(h1, fp[3], mat3);              // fp[3] is roll
564
565                 copy_m4_m3(result_array[a].mat, mat3);
566                 VECCOPY(result_array[a].mat[3], fp);
567
568                 if(doscale) {
569                         /* correct for scaling when this matrix is used in scaled space */
570                         mul_serie_m4(result_array[a].mat, iscalemat, result_array[a].mat,
571                                 scalemat, NULL, NULL, NULL, NULL, NULL);
572                 }
573         }
574         
575         return result_array;
576 }
577
578 /* ************ Armature Deform ******************* */
579
580 static void pchan_b_bone_defmats(bPoseChannel *pchan, int use_quaternion)
581 {
582         Bone *bone= pchan->bone;
583         Mat4 *b_bone= b_bone_spline_setup(pchan, 0);
584         Mat4 *b_bone_rest= b_bone_spline_setup(pchan, 1);
585         Mat4 *b_bone_mats;
586         DualQuat *b_bone_dual_quats= NULL;
587         float tmat[4][4];
588         int a;
589         
590         /* allocate b_bone matrices and dual quats */
591         b_bone_mats= MEM_mallocN((1+bone->segments)*sizeof(Mat4), "BBone defmats");
592         pchan->b_bone_mats= b_bone_mats;
593
594         if(use_quaternion) {
595                 b_bone_dual_quats= MEM_mallocN((bone->segments)*sizeof(DualQuat), "BBone dqs");
596                 pchan->b_bone_dual_quats= b_bone_dual_quats;
597         }
598         
599         /* first matrix is the inverse arm_mat, to bring points in local bone space
600            for finding out which segment it belongs to */
601         invert_m4_m4(b_bone_mats[0].mat, bone->arm_mat);
602
603         /* then we make the b_bone_mats:
604                 - first transform to local bone space
605                 - translate over the curve to the bbone mat space
606                 - transform with b_bone matrix
607                 - transform back into global space */
608         unit_m4(tmat);
609
610         for(a=0; a<bone->segments; a++) {
611                 invert_m4_m4(tmat, b_bone_rest[a].mat);
612
613                 mul_serie_m4(b_bone_mats[a+1].mat, pchan->chan_mat, bone->arm_mat,
614                         b_bone[a].mat, tmat, b_bone_mats[0].mat, NULL, NULL, NULL);
615
616                 if(use_quaternion)
617                         mat4_to_dquat( &b_bone_dual_quats[a],bone->arm_mat, b_bone_mats[a+1].mat);
618         }
619 }
620
621 static void b_bone_deform(bPoseChannel *pchan, Bone *bone, float *co, DualQuat *dq, float defmat[][3])
622 {
623         Mat4 *b_bone= pchan->b_bone_mats;
624         float (*mat)[4]= b_bone[0].mat;
625         float segment, y;
626         int a;
627         
628         /* need to transform co back to bonespace, only need y */
629         y= mat[0][1]*co[0] + mat[1][1]*co[1] + mat[2][1]*co[2] + mat[3][1];
630         
631         /* now calculate which of the b_bones are deforming this */
632         segment= bone->length/((float)bone->segments);
633         a= (int)(y/segment);
634         
635         /* note; by clamping it extends deform at endpoints, goes best with
636            straight joints in restpos. */
637         CLAMP(a, 0, bone->segments-1);
638
639         if(dq) {
640                 copy_dq_dq(dq, &((DualQuat*)pchan->b_bone_dual_quats)[a]);
641         }
642         else {
643                 mul_m4_v3(b_bone[a+1].mat, co);
644
645                 if(defmat)
646                         copy_m3_m4(defmat, b_bone[a+1].mat);
647         }
648 }
649
650 /* using vec with dist to bone b1 - b2 */
651 float distfactor_to_bone (float vec[3], float b1[3], float b2[3], float rad1, float rad2, float rdist)
652 {
653         float dist=0.0f; 
654         float bdelta[3];
655         float pdelta[3];
656         float hsqr, a, l, rad;
657         
658         sub_v3_v3v3(bdelta, b2, b1);
659         l = normalize_v3(bdelta);
660         
661         sub_v3_v3v3(pdelta, vec, b1);
662         
663         a = bdelta[0]*pdelta[0] + bdelta[1]*pdelta[1] + bdelta[2]*pdelta[2];
664         hsqr = ((pdelta[0]*pdelta[0]) + (pdelta[1]*pdelta[1]) + (pdelta[2]*pdelta[2]));
665         
666         if (a < 0.0F){
667                 /* If we're past the end of the bone, do a spherical field attenuation thing */
668                 dist= ((b1[0]-vec[0])*(b1[0]-vec[0]) +(b1[1]-vec[1])*(b1[1]-vec[1]) +(b1[2]-vec[2])*(b1[2]-vec[2])) ;
669                 rad= rad1;
670         }
671         else if (a > l){
672                 /* If we're past the end of the bone, do a spherical field attenuation thing */
673                 dist= ((b2[0]-vec[0])*(b2[0]-vec[0]) +(b2[1]-vec[1])*(b2[1]-vec[1]) +(b2[2]-vec[2])*(b2[2]-vec[2])) ;
674                 rad= rad2;
675         }
676         else {
677                 dist= (hsqr - (a*a));
678                 
679                 if(l!=0.0f) {
680                         rad= a/l;
681                         rad= rad*rad2 + (1.0f-rad)*rad1;
682                 }
683                 else rad= rad1;
684         }
685         
686         a= rad*rad;
687         if(dist < a) 
688                 return 1.0f;
689         else {
690                 l= rad+rdist;
691                 l*= l;
692                 if(rdist==0.0f || dist >= l) 
693                         return 0.0f;
694                 else {
695                         a= (float)sqrt(dist)-rad;
696                         return 1.0f-( a*a )/( rdist*rdist );
697                 }
698         }
699 }
700
701 static void pchan_deform_mat_add(bPoseChannel *pchan, float weight, float bbonemat[][3], float mat[][3])
702 {
703         float wmat[3][3];
704
705         if(pchan->bone->segments>1)
706                 copy_m3_m3(wmat, bbonemat);
707         else
708                 copy_m3_m4(wmat, pchan->chan_mat);
709
710         mul_m3_fl(wmat, weight);
711         add_m3_m3m3(mat, mat, wmat);
712 }
713
714 static float dist_bone_deform(bPoseChannel *pchan, float *vec, DualQuat *dq, float mat[][3], float *co)
715 {
716         Bone *bone= pchan->bone;
717         float fac, contrib=0.0;
718         float cop[3], bbonemat[3][3];
719         DualQuat bbonedq;
720
721         if(bone==NULL) return 0.0f;
722         
723         VECCOPY (cop, co);
724
725         fac= distfactor_to_bone(cop, bone->arm_head, bone->arm_tail, bone->rad_head, bone->rad_tail, bone->dist);
726         
727         if (fac>0.0) {
728                 
729                 fac*=bone->weight;
730                 contrib= fac;
731                 if(contrib>0.0) {
732                         if(vec) {
733                                 if(bone->segments>1)
734                                         // applies on cop and bbonemat
735                                         b_bone_deform(pchan, bone, cop, NULL, (mat)?bbonemat:NULL);
736                                 else
737                                         mul_m4_v3(pchan->chan_mat, cop);
738
739                                 //      Make this a delta from the base position
740                                 sub_v3_v3(cop, co);
741                                 madd_v3_v3fl(vec, cop, fac);
742
743                                 if(mat)
744                                         pchan_deform_mat_add(pchan, fac, bbonemat, mat);
745                         }
746                         else {
747                                 if(bone->segments>1) {
748                                         b_bone_deform(pchan, bone, cop, &bbonedq, NULL);
749                                         add_weighted_dq_dq(dq, &bbonedq, fac);
750                                 }
751                                 else
752                                         add_weighted_dq_dq(dq, pchan->dual_quat, fac);
753                         }
754                 }
755         }
756         
757         return contrib;
758 }
759
760 static void pchan_bone_deform(bPoseChannel *pchan, float weight, float *vec, DualQuat *dq, float mat[][3], float *co, float *contrib)
761 {
762         float cop[3], bbonemat[3][3];
763         DualQuat bbonedq;
764
765         if (!weight)
766                 return;
767
768         VECCOPY(cop, co);
769
770         if(vec) {
771                 if(pchan->bone->segments>1)
772                         // applies on cop and bbonemat
773                         b_bone_deform(pchan, pchan->bone, cop, NULL, (mat)?bbonemat:NULL);
774                 else
775                         mul_m4_v3(pchan->chan_mat, cop);
776                 
777                 vec[0]+=(cop[0]-co[0])*weight;
778                 vec[1]+=(cop[1]-co[1])*weight;
779                 vec[2]+=(cop[2]-co[2])*weight;
780
781                 if(mat)
782                         pchan_deform_mat_add(pchan, weight, bbonemat, mat);
783         }
784         else {
785                 if(pchan->bone->segments>1) {
786                         b_bone_deform(pchan, pchan->bone, cop, &bbonedq, NULL);
787                         add_weighted_dq_dq(dq, &bbonedq, weight);
788                 }
789                 else
790                         add_weighted_dq_dq(dq, pchan->dual_quat, weight);
791         }
792
793         (*contrib)+=weight;
794 }
795
796 void armature_deform_verts(Object *armOb, Object *target, DerivedMesh *dm,
797                                                    float (*vertexCos)[3], float (*defMats)[3][3],
798                                                    int numVerts, int deformflag, 
799                                                    float (*prevCos)[3], const char *defgrp_name)
800 {
801         bArmature *arm= armOb->data;
802         bPoseChannel *pchan, **defnrToPC = NULL;
803         MDeformVert *dverts = NULL;
804         bDeformGroup *dg;
805         DualQuat *dualquats= NULL;
806         float obinv[4][4], premat[4][4], postmat[4][4];
807         int use_envelope = deformflag & ARM_DEF_ENVELOPE;
808         int use_quaternion = deformflag & ARM_DEF_QUATERNION;
809         int invert_vgroup= deformflag & ARM_DEF_INVERT_VGROUP;
810         int numGroups = 0;              /* safety for vertexgroup index overflow */
811         int i, target_totvert = 0;      /* safety for vertexgroup overflow */
812         int use_dverts = 0;
813         int armature_def_nr;
814         int totchan;
815
816         if(arm->edbo) return;
817         
818         invert_m4_m4(obinv, target->obmat);
819         copy_m4_m4(premat, target->obmat);
820         mul_m4_m4m4(postmat, armOb->obmat, obinv);
821         invert_m4_m4(premat, postmat);
822
823         /* bone defmats are already in the channels, chan_mat */
824         
825         /* initialize B_bone matrices and dual quaternions */
826         if(use_quaternion) {
827                 totchan= BLI_countlist(&armOb->pose->chanbase);
828                 dualquats= MEM_callocN(sizeof(DualQuat)*totchan, "dualquats");
829         }
830
831         totchan= 0;
832         for(pchan = armOb->pose->chanbase.first; pchan; pchan = pchan->next) {
833                 if(!(pchan->bone->flag & BONE_NO_DEFORM)) {
834                         if(pchan->bone->segments > 1)
835                                 pchan_b_bone_defmats(pchan, use_quaternion);
836
837                         if(use_quaternion) {
838                                 pchan->dual_quat= &dualquats[totchan++];
839                                 mat4_to_dquat( pchan->dual_quat,pchan->bone->arm_mat, pchan->chan_mat);
840                         }
841                 }
842         }
843
844         /* get the def_nr for the overall armature vertex group if present */
845         armature_def_nr= defgroup_name_index(target, defgrp_name);
846         
847         if(ELEM(target->type, OB_MESH, OB_LATTICE)) {
848                 numGroups = BLI_countlist(&target->defbase);
849                 
850                 if(target->type==OB_MESH) {
851                         Mesh *me= target->data;
852                         dverts = me->dvert;
853                         if(dverts)
854                                 target_totvert = me->totvert;
855                 }
856                 else {
857                         Lattice *lt= target->data;
858                         dverts = lt->dvert;
859                         if(dverts)
860                                 target_totvert = lt->pntsu*lt->pntsv*lt->pntsw;
861                 }
862         }
863         
864         /* get a vertex-deform-index to posechannel array */
865         if(deformflag & ARM_DEF_VGROUP) {
866                 if(ELEM(target->type, OB_MESH, OB_LATTICE)) {
867                         /* if we have a DerivedMesh, only use dverts if it has them */
868                         if(dm)
869                                 if(dm->getVertData(dm, 0, CD_MDEFORMVERT))
870                                         use_dverts = 1;
871                                 else use_dverts = 0;
872                         else if(dverts) use_dverts = 1;
873
874                         if(use_dverts) {
875                                 defnrToPC = MEM_callocN(sizeof(*defnrToPC) * numGroups,
876                                                                                 "defnrToBone");
877                                 for(i = 0, dg = target->defbase.first; dg;
878                                         i++, dg = dg->next) {
879                                         defnrToPC[i] = get_pose_channel(armOb->pose, dg->name);
880                                         /* exclude non-deforming bones */
881                                         if(defnrToPC[i]) {
882                                                 if(defnrToPC[i]->bone->flag & BONE_NO_DEFORM)
883                                                         defnrToPC[i]= NULL;
884                                         }
885                                 }
886                         }
887                 }
888         }
889
890         for(i = 0; i < numVerts; i++) {
891                 MDeformVert *dvert;
892                 DualQuat sumdq, *dq = NULL;
893                 float *co, dco[3];
894                 float sumvec[3], summat[3][3];
895                 float *vec = NULL, (*smat)[3] = NULL;
896                 float contrib = 0.0f;
897                 float armature_weight = 1.0f;   /* default to 1 if no overall def group */
898                 float prevco_weight = 1.0f;             /* weight for optional cached vertexcos */
899                 int       j;
900
901                 if(use_quaternion) {
902                         memset(&sumdq, 0, sizeof(DualQuat));
903                         dq= &sumdq;
904                 }
905                 else {
906                         sumvec[0] = sumvec[1] = sumvec[2] = 0.0f;
907                         vec= sumvec;
908
909                         if(defMats) {
910                                 zero_m3(summat);
911                                 smat = summat;
912                         }
913                 }
914
915                 if(use_dverts || armature_def_nr >= 0) {
916                         if(dm) dvert = dm->getVertData(dm, i, CD_MDEFORMVERT);
917                         else if(dverts && i < target_totvert) dvert = dverts + i;
918                         else dvert = NULL;
919                 } else
920                         dvert = NULL;
921
922                 if(armature_def_nr >= 0 && dvert) {
923                         armature_weight = 0.0f; /* a def group was given, so default to 0 */
924                         for(j = 0; j < dvert->totweight; j++) {
925                                 if(dvert->dw[j].def_nr == armature_def_nr) {
926                                         armature_weight = dvert->dw[j].weight;
927                                         break;
928                                 }
929                         }
930                         /* hackish: the blending factor can be used for blending with prevCos too */
931                         if(prevCos) {
932                                 if(invert_vgroup)
933                                         prevco_weight= 1.0f-armature_weight;
934                                 else
935                                         prevco_weight= armature_weight;
936                                 armature_weight= 1.0f;
937                         }
938                 }
939
940                 /* check if there's any  point in calculating for this vert */
941                 if(armature_weight == 0.0f) continue;
942                 
943                 /* get the coord we work on */
944                 co= prevCos?prevCos[i]:vertexCos[i];
945                 
946                 /* Apply the object's matrix */
947                 mul_m4_v3(premat, co);
948                 
949                 if(use_dverts && dvert && dvert->totweight) { // use weight groups ?
950                         int deformed = 0;
951                         
952                         for(j = 0; j < dvert->totweight; j++){
953                                 int index = dvert->dw[j].def_nr;
954                                 pchan = index < numGroups?defnrToPC[index]:NULL;
955                                 if(pchan) {
956                                         float weight = dvert->dw[j].weight;
957                                         Bone *bone = pchan->bone;
958
959                                         deformed = 1;
960                                         
961                                         if(bone && bone->flag & BONE_MULT_VG_ENV) {
962                                                 weight *= distfactor_to_bone(co, bone->arm_head,
963                                                                                                          bone->arm_tail,
964                                                                                                          bone->rad_head,
965                                                                                                          bone->rad_tail,
966                                                                                                          bone->dist);
967                                         }
968                                         pchan_bone_deform(pchan, weight, vec, dq, smat, co, &contrib);
969                                 }
970                         }
971                         /* if there are vertexgroups but not groups with bones
972                          * (like for softbody groups)
973                          */
974                         if(deformed == 0 && use_envelope) {
975                                 for(pchan = armOb->pose->chanbase.first; pchan;
976                                         pchan = pchan->next) {
977                                         if(!(pchan->bone->flag & BONE_NO_DEFORM))
978                                                 contrib += dist_bone_deform(pchan, vec, dq, smat, co);
979                                 }
980                         }
981                 }
982                 else if(use_envelope) {
983                         for(pchan = armOb->pose->chanbase.first; pchan;
984                                 pchan = pchan->next) {
985                                 if(!(pchan->bone->flag & BONE_NO_DEFORM))
986                                         contrib += dist_bone_deform(pchan, vec, dq, smat, co);
987                         }
988                 }
989
990                 /* actually should be EPSILON? weight values and contrib can be like 10e-39 small */
991                 if(contrib > 0.0001f) {
992                         if(use_quaternion) {
993                                 normalize_dq(dq, contrib);
994
995                                 if(armature_weight != 1.0f) {
996                                         VECCOPY(dco, co);
997                                         mul_v3m3_dq( dco, (defMats)? summat: NULL,dq);
998                                         sub_v3_v3(dco, co);
999                                         mul_v3_fl(dco, armature_weight);
1000                                         add_v3_v3(co, dco);
1001                                 }
1002                                 else
1003                                         mul_v3m3_dq( co, (defMats)? summat: NULL,dq);
1004
1005                                 smat = summat;
1006                         }
1007                         else {
1008                                 mul_v3_fl(vec, armature_weight/contrib);
1009                                 add_v3_v3v3(co, vec, co);
1010                         }
1011
1012                         if(defMats) {
1013                                 float pre[3][3], post[3][3], tmpmat[3][3];
1014
1015                                 copy_m3_m4(pre, premat);
1016                                 copy_m3_m4(post, postmat);
1017                                 copy_m3_m3(tmpmat, defMats[i]);
1018
1019                                 if(!use_quaternion) /* quaternion already is scale corrected */
1020                                         mul_m3_fl(smat, armature_weight/contrib);
1021
1022                                 mul_serie_m3(defMats[i], tmpmat, pre, smat, post,
1023                                         NULL, NULL, NULL, NULL);
1024                         }
1025                 }
1026                 
1027                 /* always, check above code */
1028                 mul_m4_v3(postmat, co);
1029                 
1030                 
1031                 /* interpolate with previous modifier position using weight group */
1032                 if(prevCos) {
1033                         float mw= 1.0f - prevco_weight;
1034                         vertexCos[i][0]= prevco_weight*vertexCos[i][0] + mw*co[0];
1035                         vertexCos[i][1]= prevco_weight*vertexCos[i][1] + mw*co[1];
1036                         vertexCos[i][2]= prevco_weight*vertexCos[i][2] + mw*co[2];
1037                 }
1038         }
1039
1040         if(dualquats) MEM_freeN(dualquats);
1041         if(defnrToPC) MEM_freeN(defnrToPC);
1042         
1043         /* free B_bone matrices */
1044         for(pchan = armOb->pose->chanbase.first; pchan; pchan = pchan->next) {
1045                 if(pchan->b_bone_mats) {
1046                         MEM_freeN(pchan->b_bone_mats);
1047                         pchan->b_bone_mats = NULL;
1048                 }
1049                 if(pchan->b_bone_dual_quats) {
1050                         MEM_freeN(pchan->b_bone_dual_quats);
1051                         pchan->b_bone_dual_quats = NULL;
1052                 }
1053
1054                 pchan->dual_quat = NULL;
1055         }
1056 }
1057
1058 /* ************ END Armature Deform ******************* */
1059
1060 void get_objectspace_bone_matrix (struct Bone* bone, float M_accumulatedMatrix[][4], int UNUSED(root), int UNUSED(posed))
1061 {
1062         copy_m4_m4(M_accumulatedMatrix, bone->arm_mat);
1063 }
1064
1065 /* **************** Space to Space API ****************** */
1066
1067 /* Convert World-Space Matrix to Pose-Space Matrix */
1068 void armature_mat_world_to_pose(Object *ob, float inmat[][4], float outmat[][4]) 
1069 {
1070         float obmat[4][4];
1071         
1072         /* prevent crashes */
1073         if (ob==NULL) return;
1074         
1075         /* get inverse of (armature) object's matrix  */
1076         invert_m4_m4(obmat, ob->obmat);
1077         
1078         /* multiply given matrix by object's-inverse to find pose-space matrix */
1079         mul_m4_m4m4(outmat, obmat, inmat);
1080 }
1081
1082 /* Convert Wolrd-Space Location to Pose-Space Location
1083  * NOTE: this cannot be used to convert to pose-space location of the supplied
1084  *              pose-channel into its local space (i.e. 'visual'-keyframing) 
1085  */
1086 void armature_loc_world_to_pose(Object *ob, float *inloc, float *outloc) 
1087 {
1088         float xLocMat[4][4];
1089         float nLocMat[4][4];
1090         
1091         /* build matrix for location */
1092         unit_m4(xLocMat);
1093         VECCOPY(xLocMat[3], inloc);
1094
1095         /* get bone-space cursor matrix and extract location */
1096         armature_mat_world_to_pose(ob, xLocMat, nLocMat);
1097         VECCOPY(outloc, nLocMat[3]);
1098 }
1099
1100 /* Convert Pose-Space Matrix to Bone-Space Matrix 
1101  * NOTE: this cannot be used to convert to pose-space transforms of the supplied
1102  *              pose-channel into its local space (i.e. 'visual'-keyframing)
1103  */
1104 void armature_mat_pose_to_bone(bPoseChannel *pchan, float inmat[][4], float outmat[][4])
1105 {
1106         float pc_trans[4][4], inv_trans[4][4];
1107         float pc_posemat[4][4], inv_posemat[4][4];
1108         float pose_mat[4][4];
1109
1110         /* paranoia: prevent crashes with no pose-channel supplied */
1111         if (pchan==NULL) return;
1112
1113         /* default flag */
1114         if((pchan->bone->flag & BONE_NO_LOCAL_LOCATION)==0) {
1115                 /* get the inverse matrix of the pchan's transforms */
1116                 switch(pchan->rotmode) {
1117                 case ROT_MODE_QUAT:
1118                         loc_quat_size_to_mat4(pc_trans, pchan->loc, pchan->quat, pchan->size);
1119                         break;
1120                 case ROT_MODE_AXISANGLE:
1121                         loc_axisangle_size_to_mat4(pc_trans, pchan->loc, pchan->rotAxis, pchan->rotAngle, pchan->size);
1122                         break;
1123                 default: /* euler */
1124                         loc_eul_size_to_mat4(pc_trans, pchan->loc, pchan->eul, pchan->size);
1125                 }
1126
1127                 copy_m4_m4(pose_mat, pchan->pose_mat);
1128         }
1129         else {
1130                 /* local location, this is not default, different calculation
1131                  * note: only tested for location with pose bone snapping.
1132                  * If this is not useful in other cases the BONE_NO_LOCAL_LOCATION
1133                  * case may have to be split into its own function. */
1134                 unit_m4(pc_trans);
1135                 copy_v3_v3(pc_trans[3], pchan->loc);
1136
1137                 /* use parents rotation/scale space + own absolute position */
1138                 if(pchan->parent)       copy_m4_m4(pose_mat, pchan->parent->pose_mat);
1139                 else                            unit_m4(pose_mat);
1140
1141                 copy_v3_v3(pose_mat[3], pchan->pose_mat[3]);
1142         }
1143
1144
1145         invert_m4_m4(inv_trans, pc_trans);
1146         
1147         /* Remove the pchan's transforms from it's pose_mat.
1148          * This should leave behind the effects of restpose + 
1149          * parenting + constraints
1150          */
1151         mul_m4_m4m4(pc_posemat, inv_trans, pose_mat);
1152         
1153         /* get the inverse of the leftovers so that we can remove 
1154          * that component from the supplied matrix
1155          */
1156         invert_m4_m4(inv_posemat, pc_posemat);
1157         
1158         /* get the new matrix */
1159         mul_m4_m4m4(outmat, inmat, inv_posemat);
1160 }
1161
1162 /* Convert Pose-Space Location to Bone-Space Location
1163  * NOTE: this cannot be used to convert to pose-space location of the supplied
1164  *              pose-channel into its local space (i.e. 'visual'-keyframing) 
1165  */
1166 void armature_loc_pose_to_bone(bPoseChannel *pchan, float *inloc, float *outloc) 
1167 {
1168         float xLocMat[4][4];
1169         float nLocMat[4][4];
1170         
1171         /* build matrix for location */
1172         unit_m4(xLocMat);
1173         VECCOPY(xLocMat[3], inloc);
1174
1175         /* get bone-space cursor matrix and extract location */
1176         armature_mat_pose_to_bone(pchan, xLocMat, nLocMat);
1177         VECCOPY(outloc, nLocMat[3]);
1178 }
1179
1180
1181 /* Apply a 4x4 matrix to the pose bone,
1182  * similar to object_apply_mat4()
1183  */
1184 void pchan_apply_mat4(bPoseChannel *pchan, float mat[][4])
1185 {
1186         /* location */
1187         copy_v3_v3(pchan->loc, mat[3]);
1188
1189         /* scale */
1190         mat4_to_size(pchan->size, mat);
1191
1192         /* rotation */
1193         if (pchan->rotmode == ROT_MODE_AXISANGLE) {
1194                 float tmp_quat[4];
1195
1196                 /* need to convert to quat first (in temp var)... */
1197                 mat4_to_quat(tmp_quat, mat);
1198                 quat_to_axis_angle(pchan->rotAxis, &pchan->rotAngle, tmp_quat);
1199         }
1200         else if (pchan->rotmode == ROT_MODE_QUAT) {
1201                 mat4_to_quat(pchan->quat, mat);
1202         }
1203         else {
1204                 mat4_to_eulO(pchan->eul, pchan->rotmode, mat);
1205         }
1206 }
1207
1208 /* Remove rest-position effects from pose-transform for obtaining
1209  * 'visual' transformation of pose-channel.
1210  * (used by the Visual-Keyframing stuff)
1211  */
1212 void armature_mat_pose_to_delta(float delta_mat[][4], float pose_mat[][4], float arm_mat[][4])
1213 {
1214          float imat[4][4];
1215  
1216          invert_m4_m4(imat, arm_mat);
1217          mul_m4_m4m4(delta_mat, pose_mat, imat);
1218 }
1219
1220 /* **************** Rotation Mode Conversions ****************************** */
1221 /* Used for Objects and Pose Channels, since both can have multiple rotation representations */
1222
1223 /* Called from RNA when rotation mode changes 
1224  *      - the result should be that the rotations given in the provided pointers have had conversions 
1225  *        applied (as appropriate), such that the rotation of the element hasn't 'visually' changed 
1226  */
1227 void BKE_rotMode_change_values (float quat[4], float eul[3], float axis[3], float *angle, short oldMode, short newMode)
1228 {
1229         /* check if any change - if so, need to convert data */
1230         if (newMode > 0) { /* to euler */
1231                 if (oldMode == ROT_MODE_AXISANGLE) {
1232                         /* axis-angle to euler */
1233                         axis_angle_to_eulO( eul, newMode,axis, *angle);
1234                 }
1235                 else if (oldMode == ROT_MODE_QUAT) {
1236                         /* quat to euler */
1237                         quat_to_eulO( eul, newMode,quat);
1238                 }
1239                 /* else { no conversion needed } */
1240         }
1241         else if (newMode == ROT_MODE_QUAT) { /* to quat */
1242                 if (oldMode == ROT_MODE_AXISANGLE) {
1243                         /* axis angle to quat */
1244                         axis_angle_to_quat(quat, axis, *angle);
1245                 }
1246                 else if (oldMode > 0) {
1247                         /* euler to quat */
1248                         eulO_to_quat( quat,eul, oldMode);
1249                 }
1250                 /* else { no conversion needed } */
1251         }
1252         else if (newMode == ROT_MODE_AXISANGLE) { /* to axis-angle */
1253                 if (oldMode > 0) {
1254                         /* euler to axis angle */
1255                         eulO_to_axis_angle( axis, angle,eul, oldMode);
1256                 }
1257                 else if (oldMode == ROT_MODE_QUAT) {
1258                         /* quat to axis angle */
1259                         quat_to_axis_angle( axis, angle,quat);
1260                 }
1261                 
1262                 /* when converting to axis-angle, we need a special exception for the case when there is no axis */
1263                 if (IS_EQ(axis[0], axis[1]) && IS_EQ(axis[1], axis[2])) {
1264                         /* for now, rotate around y-axis then (so that it simply becomes the roll) */
1265                         axis[1]= 1.0f;
1266                 }
1267         }
1268 }
1269
1270 /* **************** The new & simple (but OK!) armature evaluation ********* */ 
1271
1272 /*  ****************** And how it works! ****************************************
1273
1274   This is the bone transformation trick; they're hierarchical so each bone(b)
1275   is in the coord system of bone(b-1):
1276
1277   arm_mat(b)= arm_mat(b-1) * yoffs(b-1) * d_root(b) * bone_mat(b) 
1278   
1279   -> yoffs is just the y axis translation in parent's coord system
1280   -> d_root is the translation of the bone root, also in parent's coord system
1281
1282   pose_mat(b)= pose_mat(b-1) * yoffs(b-1) * d_root(b) * bone_mat(b) * chan_mat(b)
1283
1284   we then - in init deform - store the deform in chan_mat, such that:
1285
1286   pose_mat(b)= arm_mat(b) * chan_mat(b)
1287   
1288   *************************************************************************** */
1289 /*  Computes vector and roll based on a rotation. "mat" must
1290          contain only a rotation, and no scaling. */ 
1291 void mat3_to_vec_roll(float mat[][3], float *vec, float *roll) 
1292 {
1293         if (vec)
1294                 copy_v3_v3(vec, mat[1]);
1295
1296         if (roll) {
1297                 float vecmat[3][3], vecmatinv[3][3], rollmat[3][3];
1298
1299                 vec_roll_to_mat3(mat[1], 0.0f, vecmat);
1300                 invert_m3_m3(vecmatinv, vecmat);
1301                 mul_m3_m3m3(rollmat, vecmatinv, mat);
1302
1303                 *roll= (float)atan2(rollmat[2][0], rollmat[2][2]);
1304         }
1305 }
1306
1307 /*      Calculates the rest matrix of a bone based
1308         On its vector and a roll around that vector */
1309 void vec_roll_to_mat3(float *vec, float roll, float mat[][3])
1310 {
1311         float   nor[3], axis[3], target[3]={0,1,0};
1312         float   theta;
1313         float   rMatrix[3][3], bMatrix[3][3];
1314
1315         normalize_v3_v3(nor, vec);
1316         
1317         /*      Find Axis & Amount for bone matrix*/
1318         cross_v3_v3v3(axis,target,nor);
1319
1320         /* was 0.0000000000001, caused bug [#23954], smaller values give unstable
1321          * roll when toggling editmode */
1322         if (dot_v3v3(axis,axis) > 0.00001) {
1323                 /* if nor is *not* a multiple of target ... */
1324                 normalize_v3(axis);
1325                 
1326                 theta= angle_normalized_v3v3(target, nor);
1327                 
1328                 /*      Make Bone matrix*/
1329                 vec_rot_to_mat3( bMatrix,axis, theta);
1330         }
1331         else {
1332                 /* if nor is a multiple of target ... */
1333                 float updown;
1334                 
1335                 /* point same direction, or opposite? */
1336                 updown = ( dot_v3v3(target,nor) > 0 ) ? 1.0f : -1.0f;
1337                 
1338                 /* I think this should work ... */
1339                 bMatrix[0][0]=updown; bMatrix[0][1]=0.0;    bMatrix[0][2]=0.0;
1340                 bMatrix[1][0]=0.0;    bMatrix[1][1]=updown; bMatrix[1][2]=0.0;
1341                 bMatrix[2][0]=0.0;    bMatrix[2][1]=0.0;    bMatrix[2][2]=1.0;
1342         }
1343         
1344         /*      Make Roll matrix*/
1345         vec_rot_to_mat3( rMatrix,nor, roll);
1346         
1347         /*      Combine and output result*/
1348         mul_m3_m3m3(mat, rMatrix, bMatrix);
1349 }
1350
1351
1352 /* recursive part, calculates restposition of entire tree of children */
1353 /* used by exiting editmode too */
1354 void where_is_armature_bone(Bone *bone, Bone *prevbone)
1355 {
1356         float vec[3];
1357         
1358         /* Bone Space */
1359         sub_v3_v3v3(vec, bone->tail, bone->head);
1360         vec_roll_to_mat3(vec, bone->roll, bone->bone_mat);
1361
1362         bone->length= len_v3v3(bone->head, bone->tail);
1363         
1364         /* this is called on old file reading too... */
1365         if(bone->xwidth==0.0) {
1366                 bone->xwidth= 0.1f;
1367                 bone->zwidth= 0.1f;
1368                 bone->segments= 1;
1369         }
1370         
1371         if(prevbone) {
1372                 float offs_bone[4][4];  // yoffs(b-1) + root(b) + bonemat(b)
1373                 
1374                 /* bone transform itself */
1375                 copy_m4_m3(offs_bone, bone->bone_mat);
1376                                 
1377                 /* The bone's root offset (is in the parent's coordinate system) */
1378                 VECCOPY(offs_bone[3], bone->head);
1379
1380                 /* Get the length translation of parent (length along y axis) */
1381                 offs_bone[3][1]+= prevbone->length;
1382                 
1383                 /* Compose the matrix for this bone  */
1384                 mul_m4_m4m4(bone->arm_mat, offs_bone, prevbone->arm_mat);
1385         }
1386         else {
1387                 copy_m4_m3(bone->arm_mat, bone->bone_mat);
1388                 VECCOPY(bone->arm_mat[3], bone->head);
1389         }
1390         
1391         /* head */
1392         VECCOPY(bone->arm_head, bone->arm_mat[3]);
1393         /* tail is in current local coord system */
1394         VECCOPY(vec, bone->arm_mat[1]);
1395         mul_v3_fl(vec, bone->length);
1396         add_v3_v3v3(bone->arm_tail, bone->arm_head, vec);
1397         
1398         /* and the kiddies */
1399         prevbone= bone;
1400         for(bone= bone->childbase.first; bone; bone= bone->next) {
1401                 where_is_armature_bone(bone, prevbone);
1402         }
1403 }
1404
1405 /* updates vectors and matrices on rest-position level, only needed 
1406    after editing armature itself, now only on reading file */
1407 void where_is_armature (bArmature *arm)
1408 {
1409         Bone *bone;
1410         
1411         /* hierarchical from root to children */
1412         for(bone= arm->bonebase.first; bone; bone= bone->next) {
1413                 where_is_armature_bone(bone, NULL);
1414         }
1415 }
1416
1417 /* if bone layer is protected, copy the data from from->pose
1418  * when used with linked libraries this copies from the linked pose into the local pose */
1419 static void pose_proxy_synchronize(Object *ob, Object *from, int layer_protected)
1420 {
1421         bPose *pose= ob->pose, *frompose= from->pose;
1422         bPoseChannel *pchan, *pchanp, pchanw;
1423         bConstraint *con;
1424         int error = 0;
1425         
1426         if (frompose==NULL) return;
1427
1428         /* in some cases when rigs change, we cant synchronize
1429          * to avoid crashing check for possible errors here */
1430         for (pchan= pose->chanbase.first; pchan; pchan= pchan->next) {
1431                 if (pchan->bone->layer & layer_protected) {
1432                         if(get_pose_channel(frompose, pchan->name) == NULL) {
1433                                 printf("failed to sync proxy armature because '%s' is missing pose channel '%s'\n", from->id.name, pchan->name);
1434                                 error = 1;
1435                         }
1436                 }
1437         }
1438
1439         if(error)
1440                 return;
1441         
1442         /* exception, armature local layer should be proxied too */
1443         if (pose->proxy_layer)
1444                 ((bArmature *)ob->data)->layer= pose->proxy_layer;
1445         
1446         /* clear all transformation values from library */
1447         rest_pose(frompose);
1448         
1449         /* copy over all of the proxy's bone groups */
1450                 /* TODO for later - implement 'local' bone groups as for constraints
1451                  *      Note: this isn't trivial, as bones reference groups by index not by pointer, 
1452                  *               so syncing things correctly needs careful attention
1453                  */
1454         BLI_freelistN(&pose->agroups);
1455         BLI_duplicatelist(&pose->agroups, &frompose->agroups);
1456         pose->active_group= frompose->active_group;
1457         
1458         for (pchan= pose->chanbase.first; pchan; pchan= pchan->next) {
1459                 pchanp= get_pose_channel(frompose, pchan->name);
1460
1461                 if (pchan->bone->layer & layer_protected) {
1462                         ListBase proxylocal_constraints = {NULL, NULL};
1463                         
1464                         /* copy posechannel to temp, but restore important pointers */
1465                         pchanw= *pchanp;
1466                         pchanw.prev= pchan->prev;
1467                         pchanw.next= pchan->next;
1468                         pchanw.parent= pchan->parent;
1469                         pchanw.child= pchan->child;
1470                         pchanw.path= NULL;
1471                         
1472                         /* this is freed so copy a copy, else undo crashes */
1473                         if(pchanw.prop) {
1474                                 pchanw.prop= IDP_CopyProperty(pchanw.prop);
1475
1476                                 /* use the values from the the existing props */
1477                                 if(pchan->prop) {
1478                                         IDP_SyncGroupValues(pchanw.prop, pchan->prop);
1479                                 }
1480                         }
1481
1482                         /* constraints - proxy constraints are flushed... local ones are added after 
1483                          *      1. extract constraints not from proxy (CONSTRAINT_PROXY_LOCAL) from pchan's constraints
1484                          *      2. copy proxy-pchan's constraints on-to new
1485                          *      3. add extracted local constraints back on top 
1486                          *
1487                          *  note for copy_constraints: when copying constraints, disable 'do_extern' otherwise we get the libs direct linked in this blend.
1488                          */
1489                         extract_proxylocal_constraints(&proxylocal_constraints, &pchan->constraints);
1490                         copy_constraints(&pchanw.constraints, &pchanp->constraints, FALSE);
1491                         addlisttolist(&pchanw.constraints, &proxylocal_constraints);
1492                         
1493                         /* constraints - set target ob pointer to own object */
1494                         for (con= pchanw.constraints.first; con; con= con->next) {
1495                                 bConstraintTypeInfo *cti= constraint_get_typeinfo(con);
1496                                 ListBase targets = {NULL, NULL};
1497                                 bConstraintTarget *ct;
1498                                 
1499                                 if (cti && cti->get_constraint_targets) {
1500                                         cti->get_constraint_targets(con, &targets);
1501                                         
1502                                         for (ct= targets.first; ct; ct= ct->next) {
1503                                                 if (ct->tar == from)
1504                                                         ct->tar = ob;
1505                                         }
1506                                         
1507                                         if (cti->flush_constraint_targets)
1508                                                 cti->flush_constraint_targets(con, &targets, 0);
1509                                 }
1510                         }
1511                         
1512                         /* free stuff from current channel */
1513                         free_pose_channel(pchan);
1514                         
1515                         /* the final copy */
1516                         *pchan= pchanw;
1517                 }
1518                 else {
1519                         /* always copy custom shape */
1520                         pchan->custom= pchanp->custom;
1521                         pchan->custom_tx= pchanp->custom_tx;
1522
1523                         /* ID-Property Syncing */
1524                         {
1525                                 IDProperty *prop_orig= pchan->prop;
1526                                 if(pchanp->prop) {
1527                                         pchan->prop= IDP_CopyProperty(pchanp->prop);
1528                                         if(prop_orig) {
1529                                                 /* copy existing values across when types match */
1530                                                 IDP_SyncGroupValues(pchan->prop, prop_orig);
1531                                         }
1532                                 }
1533                                 else {
1534                                         pchan->prop= NULL;
1535                                 }
1536                                 if (prop_orig) {
1537                                         IDP_FreeProperty(prop_orig);
1538                                         MEM_freeN(prop_orig);
1539                                 }
1540                         }
1541                 }
1542         }
1543 }
1544
1545 static int rebuild_pose_bone(bPose *pose, Bone *bone, bPoseChannel *parchan, int counter)
1546 {
1547         bPoseChannel *pchan = verify_pose_channel (pose, bone->name);   // verify checks and/or adds
1548
1549         pchan->bone= bone;
1550         pchan->parent= parchan;
1551         
1552         counter++;
1553         
1554         for(bone= bone->childbase.first; bone; bone= bone->next) {
1555                 counter= rebuild_pose_bone(pose, bone, pchan, counter);
1556                 /* for quick detecting of next bone in chain, only b-bone uses it now */
1557                 if(bone->flag & BONE_CONNECTED)
1558                         pchan->child= get_pose_channel(pose, bone->name);
1559         }
1560         
1561         return counter;
1562 }
1563
1564 /* only after leave editmode, duplicating, validating older files, library syncing */
1565 /* NOTE: pose->flag is set for it */
1566 void armature_rebuild_pose(Object *ob, bArmature *arm)
1567 {
1568         Bone *bone;
1569         bPose *pose;
1570         bPoseChannel *pchan, *next;
1571         int counter=0;
1572                 
1573         /* only done here */
1574         if(ob->pose==NULL) {
1575                 /* create new pose */
1576                 ob->pose= MEM_callocN(sizeof(bPose), "new pose");
1577                 
1578                 /* set default settings for animviz */
1579                 animviz_settings_init(&ob->pose->avs);
1580         }
1581         pose= ob->pose;
1582         
1583         /* clear */
1584         for(pchan= pose->chanbase.first; pchan; pchan= pchan->next) {
1585                 pchan->bone= NULL;
1586                 pchan->child= NULL;
1587         }
1588         
1589         /* first step, check if all channels are there */
1590         for(bone= arm->bonebase.first; bone; bone= bone->next) {
1591                 counter= rebuild_pose_bone(pose, bone, NULL, counter);
1592         }
1593
1594         /* and a check for garbage */
1595         for(pchan= pose->chanbase.first; pchan; pchan= next) {
1596                 next= pchan->next;
1597                 if(pchan->bone==NULL) {
1598                         free_pose_channel(pchan);
1599                         free_pose_channels_hash(pose);
1600                         BLI_freelinkN(&pose->chanbase, pchan);
1601                 }
1602         }
1603         // printf("rebuild pose %s, %d bones\n", ob->id.name, counter);
1604         
1605         /* synchronize protected layers with proxy */
1606         if(ob->proxy) {
1607                 object_copy_proxy_drivers(ob, ob->proxy);
1608                 pose_proxy_synchronize(ob, ob->proxy, arm->layer_protected);
1609         }
1610         
1611         update_pose_constraint_flags(ob->pose); // for IK detection for example
1612         
1613         /* the sorting */
1614         if(counter>1)
1615                 DAG_pose_sort(ob);
1616         
1617         ob->pose->flag &= ~POSE_RECALC;
1618         ob->pose->flag |= POSE_WAS_REBUILT;
1619
1620         make_pose_channels_hash(ob->pose);
1621 }
1622
1623
1624 /* ********************** SPLINE IK SOLVER ******************* */
1625
1626 /* Temporary evaluation tree data used for Spline IK */
1627 typedef struct tSplineIK_Tree {
1628         struct tSplineIK_Tree *next, *prev;
1629         
1630         int     type;                                   /* type of IK that this serves (CONSTRAINT_TYPE_KINEMATIC or ..._SPLINEIK) */
1631         
1632         short free_points;                              /* free the point positions array */
1633         short chainlen;                                 /* number of bones in the chain */
1634         
1635         float *points;                                  /* parametric positions for the joints along the curve */
1636         bPoseChannel **chain;                   /* chain of bones to affect using Spline IK (ordered from the tip) */
1637         
1638         bPoseChannel *root;                             /* bone that is the root node of the chain */
1639         
1640         bConstraint *con;                               /* constraint for this chain */
1641         bSplineIKConstraint *ikData;    /* constraint settings for this chain */
1642 } tSplineIK_Tree;
1643
1644 /* ----------- */
1645
1646 /* Tag the bones in the chain formed by the given bone for IK */
1647 static void splineik_init_tree_from_pchan(Scene *scene, Object *UNUSED(ob), bPoseChannel *pchan_tip)
1648 {
1649         bPoseChannel *pchan, *pchanRoot=NULL;
1650         bPoseChannel *pchanChain[255];
1651         bConstraint *con = NULL;
1652         bSplineIKConstraint *ikData = NULL;
1653         float boneLengths[255], *jointPoints;
1654         float totLength = 0.0f;
1655         short free_joints = 0;
1656         int segcount = 0;
1657         
1658         /* find the SplineIK constraint */
1659         for (con= pchan_tip->constraints.first; con; con= con->next) {
1660                 if (con->type == CONSTRAINT_TYPE_SPLINEIK) {
1661                         ikData= con->data;
1662                         
1663                         /* target can only be curve */
1664                         if ((ikData->tar == NULL) || (ikData->tar->type != OB_CURVE))  
1665                                 continue;
1666                         /* skip if disabled */
1667                         if ( (con->enforce == 0.0f) || (con->flag & (CONSTRAINT_DISABLE|CONSTRAINT_OFF)) )
1668                                 continue;
1669                         
1670                         /* otherwise, constraint is ok... */
1671                         break;
1672                 }
1673         }
1674         if (con == NULL)
1675                 return;
1676                 
1677         /* make sure that the constraint targets are ok 
1678          *      - this is a workaround for a depsgraph bug...
1679          */
1680         if (ikData->tar) {
1681                 Curve *cu= ikData->tar->data;
1682                 
1683                 /* note: when creating constraints that follow path, the curve gets the CU_PATH set now,
1684                  *              currently for paths to work it needs to go through the bevlist/displist system (ton) 
1685                  */
1686                 
1687                 /* only happens on reload file, but violates depsgraph still... fix! */
1688                 if ((cu->path==NULL) || (cu->path->data==NULL))
1689                         makeDispListCurveTypes(scene, ikData->tar, 0);
1690         }
1691         
1692         /* find the root bone and the chain of bones from the root to the tip 
1693          * NOTE: this assumes that the bones are connected, but that may not be true...
1694          */
1695         for (pchan= pchan_tip; pchan && (segcount < ikData->chainlen); pchan= pchan->parent, segcount++) {
1696                 /* store this segment in the chain */
1697                 pchanChain[segcount]= pchan;
1698                 
1699                 /* if performing rebinding, calculate the length of the bone */
1700                 boneLengths[segcount]= pchan->bone->length;
1701                 totLength += boneLengths[segcount];
1702         }
1703         
1704         if (segcount == 0)
1705                 return;
1706         else
1707                 pchanRoot= pchanChain[segcount-1];
1708         
1709         /* perform binding step if required */
1710         if ((ikData->flag & CONSTRAINT_SPLINEIK_BOUND) == 0) {
1711                 float segmentLen= (1.0f / (float)segcount);
1712                 int i;
1713                 
1714                 /* setup new empty array for the points list */
1715                 if (ikData->points) 
1716                         MEM_freeN(ikData->points);
1717                 ikData->numpoints= ikData->chainlen+1; 
1718                 ikData->points= MEM_callocN(sizeof(float)*ikData->numpoints, "Spline IK Binding");
1719                 
1720                 /* perform binding of the joints to parametric positions along the curve based 
1721                  * proportion of the total length that each bone occupies
1722                  */
1723                 for (i = 0; i < segcount; i++) {
1724                         if (i != 0) {
1725                                 /* 'head' joints 
1726                                  *      - 2 methods; the one chosen depends on whether we've got usable lengths
1727                                  */
1728                                 if ((ikData->flag & CONSTRAINT_SPLINEIK_EVENSPLITS) || (totLength == 0.0f)) {
1729                                         /* 1) equi-spaced joints */
1730                                         ikData->points[i]= ikData->points[i-1] - segmentLen;
1731                                 }
1732                                 else {
1733                                          /*     2) to find this point on the curve, we take a step from the previous joint
1734                                           *       a distance given by the proportion that this bone takes
1735                                           */
1736                                         ikData->points[i]= ikData->points[i-1] - (boneLengths[i] / totLength);
1737                                 }
1738                         }
1739                         else {
1740                                 /* 'tip' of chain, special exception for the first joint */
1741                                 ikData->points[0]= 1.0f;
1742                         }
1743                 }
1744                 
1745                 /* spline has now been bound */
1746                 ikData->flag |= CONSTRAINT_SPLINEIK_BOUND;
1747         }
1748         
1749         /* apply corrections for sensitivity to scaling on a copy of the bind points,
1750          * since it's easier to determine the positions of all the joints beforehand this way
1751          */
1752         if ((ikData->flag & CONSTRAINT_SPLINEIK_SCALE_LIMITED) && (totLength != 0.0f)) {
1753                 Curve *cu= (Curve *)ikData->tar->data;
1754                 float splineLen, maxScale;
1755                 int i;
1756                 
1757                 /* make a copy of the points array, that we'll store in the tree 
1758                  *      - although we could just multiply the points on the fly, this approach means that
1759                  *        we can introduce per-segment stretchiness later if it is necessary
1760                  */
1761                 jointPoints= MEM_dupallocN(ikData->points);
1762                 free_joints= 1;
1763                 
1764                 /* get the current length of the curve */
1765                 // NOTE: this is assumed to be correct even after the curve was resized
1766                 splineLen= cu->path->totdist;
1767                 
1768                 /* calculate the scale factor to multiply all the path values by so that the 
1769                  * bone chain retains its current length, such that
1770                  *      maxScale * splineLen = totLength
1771                  */
1772                 maxScale = totLength / splineLen;
1773                 
1774                 /* apply scaling correction to all of the temporary points */
1775                 // TODO: this is really not adequate enough on really short chains
1776                 for (i = 0; i < segcount; i++)
1777                         jointPoints[i] *= maxScale;
1778         }
1779         else {
1780                 /* just use the existing points array */
1781                 jointPoints= ikData->points;
1782                 free_joints= 0;
1783         }
1784         
1785         /* make a new Spline-IK chain, and store it in the IK chains */
1786         // TODO: we should check if there is already an IK chain on this, since that would take presidence...
1787         {
1788                 /* make new tree */
1789                 tSplineIK_Tree *tree= MEM_callocN(sizeof(tSplineIK_Tree), "SplineIK Tree");
1790                 tree->type= CONSTRAINT_TYPE_SPLINEIK;
1791                 
1792                 tree->chainlen= segcount;
1793                 
1794                 /* copy over the array of links to bones in the chain (from tip to root) */
1795                 tree->chain= MEM_callocN(sizeof(bPoseChannel*)*segcount, "SplineIK Chain");
1796                 memcpy(tree->chain, pchanChain, sizeof(bPoseChannel*)*segcount);
1797                 
1798                 /* store reference to joint position array */
1799                 tree->points= jointPoints;
1800                 tree->free_points= free_joints;
1801                 
1802                 /* store references to different parts of the chain */
1803                 tree->root= pchanRoot;
1804                 tree->con= con;
1805                 tree->ikData= ikData;
1806                 
1807                 /* AND! link the tree to the root */
1808                 BLI_addtail(&pchanRoot->iktree, tree);
1809         }
1810         
1811         /* mark root channel having an IK tree */
1812         pchanRoot->flag |= POSE_IKSPLINE;
1813 }
1814
1815 /* Tag which bones are members of Spline IK chains */
1816 static void splineik_init_tree(Scene *scene, Object *ob, float UNUSED(ctime))
1817 {
1818         bPoseChannel *pchan;
1819         
1820         /* find the tips of Spline IK chains, which are simply the bones which have been tagged as such */
1821         for (pchan= ob->pose->chanbase.first; pchan; pchan= pchan->next) {
1822                 if (pchan->constflag & PCHAN_HAS_SPLINEIK)
1823                         splineik_init_tree_from_pchan(scene, ob, pchan);
1824         }
1825 }
1826
1827 /* ----------- */
1828
1829 /* Evaluate spline IK for a given bone */
1830 static void splineik_evaluate_bone(tSplineIK_Tree *tree, Scene *scene, Object *ob, bPoseChannel *pchan, int index, float ctime)
1831 {
1832         bSplineIKConstraint *ikData= tree->ikData;
1833         float poseHead[3], poseTail[3], poseMat[4][4]; 
1834         float splineVec[3], scaleFac, radius=1.0f;
1835         
1836         /* firstly, calculate the bone matrix the standard way, since this is needed for roll control */
1837         where_is_pose_bone(scene, ob, pchan, ctime, 1);
1838         
1839         VECCOPY(poseHead, pchan->pose_head);
1840         VECCOPY(poseTail, pchan->pose_tail);
1841         
1842         /* step 1: determine the positions for the endpoints of the bone */
1843         {
1844                 float vec[4], dir[3], rad;
1845                 float tailBlendFac= 1.0f;
1846                 
1847                 /* determine if the bone should still be affected by SplineIK */
1848                 if (tree->points[index+1] >= 1.0f) {
1849                         /* spline doesn't affect the bone anymore, so done... */
1850                         pchan->flag |= POSE_DONE;
1851                         return;
1852                 }
1853                 else if ((tree->points[index] >= 1.0f) && (tree->points[index+1] < 1.0f)) {
1854                         /* blending factor depends on the amount of the bone still left on the chain */
1855                         tailBlendFac= (1.0f - tree->points[index+1]) / (tree->points[index] - tree->points[index+1]);
1856                 }
1857                 
1858                 /* tail endpoint */
1859                 if ( where_on_path(ikData->tar, tree->points[index], vec, dir, NULL, &rad, NULL) ) {
1860                         /* apply curve's object-mode transforms to the position 
1861                          * unless the option to allow curve to be positioned elsewhere is activated (i.e. no root)
1862                          */
1863                         if ((ikData->flag & CONSTRAINT_SPLINEIK_NO_ROOT) == 0)
1864                                 mul_m4_v3(ikData->tar->obmat, vec);
1865                         
1866                         /* convert the position to pose-space, then store it */
1867                         mul_m4_v3(ob->imat, vec);
1868                         interp_v3_v3v3(poseTail, pchan->pose_tail, vec, tailBlendFac);
1869                         
1870                         /* set the new radius */
1871                         radius= rad;
1872                 }
1873                 
1874                 /* head endpoint */
1875                 if ( where_on_path(ikData->tar, tree->points[index+1], vec, dir, NULL, &rad, NULL) ) {
1876                         /* apply curve's object-mode transforms to the position 
1877                          * unless the option to allow curve to be positioned elsewhere is activated (i.e. no root)
1878                          */
1879                         if ((ikData->flag & CONSTRAINT_SPLINEIK_NO_ROOT) == 0)
1880                                 mul_m4_v3(ikData->tar->obmat, vec);
1881                         
1882                         /* store the position, and convert it to pose space */
1883                         mul_m4_v3(ob->imat, vec);
1884                         VECCOPY(poseHead, vec);
1885                         
1886                         /* set the new radius (it should be the average value) */
1887                         radius = (radius+rad) / 2;
1888                 }
1889         }
1890         
1891         /* step 2: determine the implied transform from these endpoints 
1892          *      - splineVec: the vector direction that the spline applies on the bone
1893          *      - scaleFac: the factor that the bone length is scaled by to get the desired amount
1894          */
1895         sub_v3_v3v3(splineVec, poseTail, poseHead);
1896         scaleFac= len_v3(splineVec) / pchan->bone->length;
1897         
1898         /* step 3: compute the shortest rotation needed to map from the bone rotation to the current axis 
1899          *      - this uses the same method as is used for the Damped Track Constraint (see the code there for details)
1900          */
1901         {
1902                 float dmat[3][3], rmat[3][3], tmat[3][3];
1903                 float raxis[3], rangle;
1904                 
1905                 /* compute the raw rotation matrix from the bone's current matrix by extracting only the
1906                  * orientation-relevant axes, and normalising them
1907                  */
1908                 VECCOPY(rmat[0], pchan->pose_mat[0]);
1909                 VECCOPY(rmat[1], pchan->pose_mat[1]);
1910                 VECCOPY(rmat[2], pchan->pose_mat[2]);
1911                 normalize_m3(rmat);
1912                 
1913                 /* also, normalise the orientation imposed by the bone, now that we've extracted the scale factor */
1914                 normalize_v3(splineVec);
1915                 
1916                 /* calculate smallest axis-angle rotation necessary for getting from the
1917                  * current orientation of the bone, to the spline-imposed direction
1918                  */
1919                 cross_v3_v3v3(raxis, rmat[1], splineVec);
1920                 
1921                 rangle= dot_v3v3(rmat[1], splineVec);
1922                 rangle= acos( MAX2(-1.0f, MIN2(1.0f, rangle)) );
1923                 
1924                 /* multiply the magnitude of the angle by the influence of the constraint to 
1925                  * control the influence of the SplineIK effect 
1926                  */
1927                 rangle *= tree->con->enforce;
1928                 
1929                 /* construct rotation matrix from the axis-angle rotation found above 
1930                  *      - this call takes care to make sure that the axis provided is a unit vector first
1931                  */
1932                 axis_angle_to_mat3(dmat, raxis, rangle);
1933                 
1934                 /* combine these rotations so that the y-axis of the bone is now aligned as the spline dictates,
1935                  * while still maintaining roll control from the existing bone animation
1936                  */
1937                 mul_m3_m3m3(tmat, dmat, rmat); // m1, m3, m2
1938                 normalize_m3(tmat); /* attempt to reduce shearing, though I doubt this'll really help too much now... */
1939                 copy_m4_m3(poseMat, tmat);
1940         }
1941         
1942         /* step 4: set the scaling factors for the axes */
1943         {
1944                 /* only multiply the y-axis by the scaling factor to get nice volume-preservation */
1945                 mul_v3_fl(poseMat[1], scaleFac);
1946                 
1947                 /* set the scaling factors of the x and z axes from... */
1948                 switch (ikData->xzScaleMode) {
1949                         case CONSTRAINT_SPLINEIK_XZS_ORIGINAL:
1950                         {
1951                                 /* original scales get used */
1952                                 float scale;
1953                                 
1954                                 /* x-axis scale */
1955                                 scale= len_v3(pchan->pose_mat[0]);
1956                                 mul_v3_fl(poseMat[0], scale);
1957                                 /* z-axis scale */
1958                                 scale= len_v3(pchan->pose_mat[2]);
1959                                 mul_v3_fl(poseMat[2], scale);
1960                         }
1961                                 break;
1962                         case CONSTRAINT_SPLINEIK_XZS_VOLUMETRIC:
1963                         {
1964                                 /* 'volume preservation' */
1965                                 float scale;
1966                                 
1967                                 /* calculate volume preservation factor which is 
1968                                  * basically the inverse of the y-scaling factor 
1969                                  */
1970                                 if (fabs(scaleFac) != 0.0f) {
1971                                         scale= 1.0 / fabs(scaleFac);
1972                                         
1973                                         /* we need to clamp this within sensible values */
1974                                         // NOTE: these should be fine for now, but should get sanitised in future
1975                                         scale= MIN2(MAX2(scale, 0.0001) , 100000);
1976                                 }
1977                                 else
1978                                         scale= 1.0f;
1979                                 
1980                                 /* apply the scaling */
1981                                 mul_v3_fl(poseMat[0], scale);
1982                                 mul_v3_fl(poseMat[2], scale);
1983                         }
1984                                 break;
1985                 }
1986                 
1987                 /* finally, multiply the x and z scaling by the radius of the curve too, 
1988                  * to allow automatic scales to get tweaked still
1989                  */
1990                 if ((ikData->flag & CONSTRAINT_SPLINEIK_NO_CURVERAD) == 0) {
1991                         mul_v3_fl(poseMat[0], radius);
1992                         mul_v3_fl(poseMat[2], radius);
1993                 }
1994         }
1995         
1996         /* step 5: set the location of the bone in the matrix */
1997         if (ikData->flag & CONSTRAINT_SPLINEIK_NO_ROOT) {
1998                 /* when the 'no-root' option is affected, the chain can retain
1999                  * the shape but be moved elsewhere
2000                  */
2001                 VECCOPY(poseHead, pchan->pose_head);
2002         }
2003         else if (tree->con->enforce < 1.0f) {
2004                 /* when the influence is too low
2005                  *      - blend the positions for the 'root' bone
2006                  *      - stick to the parent for any other
2007                  */
2008                 if (pchan->parent) {
2009                         VECCOPY(poseHead, pchan->pose_head);
2010                 }
2011                 else {
2012                         // FIXME: this introduces popping artifacts when we reach 0.0
2013                         interp_v3_v3v3(poseHead, pchan->pose_head, poseHead, tree->con->enforce);
2014                 }
2015         }
2016         VECCOPY(poseMat[3], poseHead);
2017         
2018         /* finally, store the new transform */
2019         copy_m4_m4(pchan->pose_mat, poseMat);
2020         VECCOPY(pchan->pose_head, poseHead);
2021         VECCOPY(pchan->pose_tail, poseTail);
2022         
2023         /* done! */
2024         pchan->flag |= POSE_DONE;
2025 }
2026
2027 /* Evaluate the chain starting from the nominated bone */
2028 static void splineik_execute_tree(Scene *scene, Object *ob, bPoseChannel *pchan_root, float ctime)
2029 {
2030         tSplineIK_Tree *tree;
2031         
2032         /* for each pose-tree, execute it if it is spline, otherwise just free it */
2033         for (tree= pchan_root->iktree.first; tree; tree= pchan_root->iktree.first) {
2034                 /* only evaluate if tagged for Spline IK */
2035                 if (tree->type == CONSTRAINT_TYPE_SPLINEIK) {
2036                         int i;
2037                         
2038                         /* walk over each bone in the chain, calculating the effects of spline IK
2039                          *      - the chain is traversed in the opposite order to storage order (i.e. parent to children)
2040                          *        so that dependencies are correct
2041                          */
2042                         for (i= tree->chainlen-1; i >= 0; i--) {
2043                                 bPoseChannel *pchan= tree->chain[i];
2044                                 splineik_evaluate_bone(tree, scene, ob, pchan, i, ctime);
2045                         }
2046                         
2047                         /* free the tree info specific to SplineIK trees now */
2048                         if (tree->chain) MEM_freeN(tree->chain);
2049                         if (tree->free_points) MEM_freeN(tree->points);
2050                 }
2051                 
2052                 /* free this tree */
2053                 BLI_freelinkN(&pchan_root->iktree, tree);
2054         }
2055 }
2056
2057 /* ********************** THE POSE SOLVER ******************* */
2058
2059 /* loc/rot/size to given mat4 */
2060 void pchan_to_mat4(bPoseChannel *pchan, float chan_mat[4][4])
2061 {
2062         float smat[3][3];
2063         float rmat[3][3];
2064         float tmat[3][3];
2065         
2066         /* get scaling matrix */
2067         size_to_mat3(smat, pchan->size);
2068         
2069         /* rotations may either be quats, eulers (with various rotation orders), or axis-angle */
2070         if (pchan->rotmode > 0) {
2071                 /* euler rotations (will cause gimble lock, but this can be alleviated a bit with rotation orders) */
2072                 eulO_to_mat3(rmat, pchan->eul, pchan->rotmode);
2073         }
2074         else if (pchan->rotmode == ROT_MODE_AXISANGLE) {
2075                 /* axis-angle - not really that great for 3D-changing orientations */
2076                 axis_angle_to_mat3(rmat, pchan->rotAxis, pchan->rotAngle);
2077         }
2078         else {
2079                 /* quats are normalised before use to eliminate scaling issues */
2080                 float quat[4];
2081                 
2082                 /* NOTE: we now don't normalise the stored values anymore, since this was kindof evil in some cases
2083                  * but if this proves to be too problematic, switch back to the old system of operating directly on 
2084                  * the stored copy
2085                  */
2086                 QUATCOPY(quat, pchan->quat);
2087                 normalize_qt(quat);
2088                 quat_to_mat3(rmat, quat);
2089         }
2090         
2091         /* calculate matrix of bone (as 3x3 matrix, but then copy the 4x4) */
2092         mul_m3_m3m3(tmat, rmat, smat);
2093         copy_m4_m3(chan_mat, tmat);
2094         
2095         /* prevent action channels breaking chains */
2096         /* need to check for bone here, CONSTRAINT_TYPE_ACTION uses this call */
2097         if ((pchan->bone==NULL) || !(pchan->bone->flag & BONE_CONNECTED)) {
2098                 VECCOPY(chan_mat[3], pchan->loc);
2099         }
2100 }
2101
2102 /* loc/rot/size to mat4 */
2103 /* used in constraint.c too */
2104 void chan_calc_mat(bPoseChannel *pchan)
2105 {
2106         /* this is just a wrapper around the copy of this function which calculates the matrix 
2107          * and stores the result in any given channel
2108          */
2109         pchan_to_mat4(pchan, pchan->chan_mat);
2110 }
2111
2112 /* NLA strip modifiers */
2113 static void do_strip_modifiers(Scene *scene, Object *armob, Bone *bone, bPoseChannel *pchan)
2114 {
2115         bActionModifier *amod;
2116         bActionStrip *strip, *strip2;
2117         float scene_cfra= (float)scene->r.cfra;
2118         int do_modif;
2119
2120         for (strip=armob->nlastrips.first; strip; strip=strip->next) {
2121                 do_modif=0;
2122                 
2123                 if (scene_cfra>=strip->start && scene_cfra<=strip->end)
2124                         do_modif=1;
2125                 
2126                 if ((scene_cfra > strip->end) && (strip->flag & ACTSTRIP_HOLDLASTFRAME)) {
2127                         do_modif=1;
2128                         
2129                         /* if there are any other strips active, ignore modifiers for this strip - 
2130                          * 'hold' option should only hold action modifiers if there are 
2131                          * no other active strips */
2132                         for (strip2=strip->next; strip2; strip2=strip2->next) {
2133                                 if (strip2 == strip) continue;
2134                                 
2135                                 if (scene_cfra>=strip2->start && scene_cfra<=strip2->end) {
2136                                         if (!(strip2->flag & ACTSTRIP_MUTE))
2137                                                 do_modif=0;
2138                                 }
2139                         }
2140                         
2141                         /* if there are any later, activated, strips with 'hold' set, they take precedence, 
2142                          * so ignore modifiers for this strip */
2143                         for (strip2=strip->next; strip2; strip2=strip2->next) {
2144                                 if (scene_cfra < strip2->start) continue;
2145                                 if ((strip2->flag & ACTSTRIP_HOLDLASTFRAME) && !(strip2->flag & ACTSTRIP_MUTE)) {
2146                                         do_modif=0;
2147                                 }
2148                         }
2149                 }
2150                 
2151                 if (do_modif) {
2152                         /* temporal solution to prevent 2 strips accumulating */
2153                         if(scene_cfra==strip->end && strip->next && strip->next->start==scene_cfra)
2154                                 continue;
2155                         
2156                         for(amod= strip->modifiers.first; amod; amod= amod->next) {
2157                                 switch (amod->type) {
2158                                 case ACTSTRIP_MOD_DEFORM:
2159                                 {
2160                                         /* validate first */
2161                                         if(amod->ob && amod->ob->type==OB_CURVE && amod->channel[0]) {
2162                                                 
2163                                                 if( strcmp(pchan->name, amod->channel)==0 ) {
2164                                                         float mat4[4][4], mat3[3][3];
2165                                                         
2166                                                         curve_deform_vector(scene, amod->ob, armob, bone->arm_mat[3], pchan->pose_mat[3], mat3, amod->no_rot_axis);
2167                                                         copy_m4_m4(mat4, pchan->pose_mat);
2168                                                         mul_m4_m3m4(pchan->pose_mat, mat3, mat4);
2169                                                         
2170                                                 }
2171                                         }
2172                                 }
2173                                         break;
2174                                 case ACTSTRIP_MOD_NOISE:        
2175                                 {
2176                                         if( strcmp(pchan->name, amod->channel)==0 ) {
2177                                                 float nor[3], loc[3], ofs;
2178                                                 float eul[3], size[3], eulo[3], sizeo[3];
2179                                                 
2180                                                 /* calculate turbulance */
2181                                                 ofs = amod->turbul / 200.0f;
2182                                                 
2183                                                 /* make a copy of starting conditions */
2184                                                 VECCOPY(loc, pchan->pose_mat[3]);
2185                                                 mat4_to_eul( eul,pchan->pose_mat);
2186                                                 mat4_to_size( size,pchan->pose_mat);
2187                                                 VECCOPY(eulo, eul);
2188                                                 VECCOPY(sizeo, size);
2189                                                 
2190                                                 /* apply noise to each set of channels */
2191                                                 if (amod->channels & 4) {
2192                                                         /* for scaling */
2193                                                         nor[0] = BLI_gNoise(amod->noisesize, size[0]+ofs, size[1], size[2], 0, 0) - ofs;
2194                                                         nor[1] = BLI_gNoise(amod->noisesize, size[0], size[1]+ofs, size[2], 0, 0) - ofs;        
2195                                                         nor[2] = BLI_gNoise(amod->noisesize, size[0], size[1], size[2]+ofs, 0, 0) - ofs;
2196                                                         add_v3_v3(size, nor);
2197                                                         
2198                                                         if (sizeo[0] != 0)
2199                                                                 mul_v3_fl(pchan->pose_mat[0], size[0] / sizeo[0]);
2200                                                         if (sizeo[1] != 0)
2201                                                                 mul_v3_fl(pchan->pose_mat[1], size[1] / sizeo[1]);
2202                                                         if (sizeo[2] != 0)
2203                                                                 mul_v3_fl(pchan->pose_mat[2], size[2] / sizeo[2]);
2204                                                 }
2205                                                 if (amod->channels & 2) {
2206                                                         /* for rotation */
2207                                                         nor[0] = BLI_gNoise(amod->noisesize, eul[0]+ofs, eul[1], eul[2], 0, 0) - ofs;
2208                                                         nor[1] = BLI_gNoise(amod->noisesize, eul[0], eul[1]+ofs, eul[2], 0, 0) - ofs;   
2209                                                         nor[2] = BLI_gNoise(amod->noisesize, eul[0], eul[1], eul[2]+ofs, 0, 0) - ofs;
2210                                                         
2211                                                         compatible_eul(nor, eulo);
2212                                                         add_v3_v3(eul, nor);
2213                                                         compatible_eul(eul, eulo);
2214                                                         
2215                                                         loc_eul_size_to_mat4(pchan->pose_mat, loc, eul, size);
2216                                                 }
2217                                                 if (amod->channels & 1) {
2218                                                         /* for location */
2219                                                         nor[0] = BLI_gNoise(amod->noisesize, loc[0]+ofs, loc[1], loc[2], 0, 0) - ofs;
2220                                                         nor[1] = BLI_gNoise(amod->noisesize, loc[0], loc[1]+ofs, loc[2], 0, 0) - ofs;   
2221                                                         nor[2] = BLI_gNoise(amod->noisesize, loc[0], loc[1], loc[2]+ofs, 0, 0) - ofs;
2222                                                         
2223                                                         add_v3_v3v3(pchan->pose_mat[3], loc, nor);
2224                                                 }
2225                                         }
2226                                 }
2227                                         break;
2228                                 }
2229                         }
2230                 }
2231         }
2232 }
2233
2234
2235 /* The main armature solver, does all constraints excluding IK */
2236 /* pchan is validated, as having bone and parent pointer
2237  * 'do_extra': when zero skips loc/size/rot, constraints and strip modifiers.
2238  */
2239 void where_is_pose_bone(Scene *scene, Object *ob, bPoseChannel *pchan, float ctime, int do_extra)
2240 {
2241         Bone *bone, *parbone;
2242         bPoseChannel *parchan;
2243         float vec[3];
2244         
2245         /* set up variables for quicker access below */
2246         bone= pchan->bone;
2247         parbone= bone->parent;
2248         parchan= pchan->parent;
2249         
2250         /* this gives a chan_mat with actions (ipos) results */
2251         if(do_extra)    chan_calc_mat(pchan);
2252         else                    unit_m4(pchan->chan_mat);
2253
2254         /* construct the posemat based on PoseChannels, that we do before applying constraints */
2255         /* pose_mat(b)= pose_mat(b-1) * yoffs(b-1) * d_root(b) * bone_mat(b) * chan_mat(b) */
2256         
2257         if(parchan) {
2258                 float offs_bone[4][4];  // yoffs(b-1) + root(b) + bonemat(b)
2259                 
2260                 /* bone transform itself */
2261                 copy_m4_m3(offs_bone, bone->bone_mat);
2262                 
2263                 /* The bone's root offset (is in the parent's coordinate system) */
2264                 VECCOPY(offs_bone[3], bone->head);
2265                 
2266                 /* Get the length translation of parent (length along y axis) */
2267                 offs_bone[3][1]+= parbone->length;
2268                 
2269                 /* Compose the matrix for this bone  */
2270                 if((bone->flag & BONE_HINGE) && (bone->flag & BONE_NO_SCALE)) { // uses restposition rotation, but actual position
2271                         float tmat[4][4];
2272                         /* the rotation of the parent restposition */
2273                         copy_m4_m4(tmat, parbone->arm_mat);
2274                         mul_serie_m4(pchan->pose_mat, tmat, offs_bone, pchan->chan_mat, NULL, NULL, NULL, NULL, NULL);
2275                 }
2276                 else if(bone->flag & BONE_HINGE) {      // same as above but apply parent scale
2277                         float tmat[4][4];
2278
2279                         /* apply the parent matrix scale */
2280                         float tsmat[4][4], tscale[3];
2281
2282                         /* the rotation of the parent restposition */
2283                         copy_m4_m4(tmat, parbone->arm_mat);
2284
2285                         /* extract the scale of the parent matrix */
2286                         mat4_to_size(tscale, parchan->pose_mat);
2287                         size_to_mat4(tsmat, tscale);
2288                         mul_m4_m4m4(tmat, tmat, tsmat);
2289
2290                         mul_serie_m4(pchan->pose_mat, tmat, offs_bone, pchan->chan_mat, NULL, NULL, NULL, NULL, NULL);
2291                 }
2292                 else if(bone->flag & BONE_NO_SCALE) {
2293                         float orthmat[4][4];
2294                         
2295                         /* do transform, with an ortho-parent matrix */
2296                         copy_m4_m4(orthmat, parchan->pose_mat);
2297                         normalize_m4(orthmat);
2298                         mul_serie_m4(pchan->pose_mat, orthmat, offs_bone, pchan->chan_mat, NULL, NULL, NULL, NULL, NULL);
2299                 }
2300                 else
2301                         mul_serie_m4(pchan->pose_mat, parchan->pose_mat, offs_bone, pchan->chan_mat, NULL, NULL, NULL, NULL, NULL);
2302                 
2303                 /* in these cases we need to compute location separately */
2304                 if(bone->flag & (BONE_HINGE|BONE_NO_SCALE|BONE_NO_LOCAL_LOCATION)) {
2305                         float bone_loc[3], chan_loc[3];
2306
2307                         mul_v3_m4v3(bone_loc, parchan->pose_mat, offs_bone[3]);
2308                         copy_v3_v3(chan_loc, pchan->chan_mat[3]);
2309
2310                         /* no local location is not transformed by bone matrix */
2311                         if(!(bone->flag & BONE_NO_LOCAL_LOCATION))
2312                                 mul_mat3_m4_v3(offs_bone, chan_loc);
2313
2314                         /* for hinge we use armature instead of pose mat */
2315                         if(bone->flag & BONE_HINGE) mul_mat3_m4_v3(parbone->arm_mat, chan_loc);
2316                         else mul_mat3_m4_v3(parchan->pose_mat, chan_loc);
2317
2318                         add_v3_v3v3(pchan->pose_mat[3], bone_loc, chan_loc);
2319                 }
2320         }
2321         else {
2322                 mul_m4_m4m4(pchan->pose_mat, pchan->chan_mat, bone->arm_mat);
2323
2324                 /* optional location without arm_mat rotation */
2325                 if(bone->flag & BONE_NO_LOCAL_LOCATION)
2326                         add_v3_v3v3(pchan->pose_mat[3], bone->arm_mat[3], pchan->chan_mat[3]);
2327                 
2328                 /* only rootbones get the cyclic offset (unless user doesn't want that) */
2329                 if ((bone->flag & BONE_NO_CYCLICOFFSET) == 0)
2330                         add_v3_v3(pchan->pose_mat[3], ob->pose->cyclic_offset);
2331         }
2332         
2333         if(do_extra) {
2334                 /* do NLA strip modifiers - i.e. curve follow */
2335                 do_strip_modifiers(scene, ob, bone, pchan);
2336                 
2337                 /* Do constraints */
2338                 if (pchan->constraints.first) {
2339                         bConstraintOb *cob;
2340
2341                         /* make a copy of location of PoseChannel for later */
2342                         VECCOPY(vec, pchan->pose_mat[3]);
2343
2344                         /* prepare PoseChannel for Constraint solving
2345                          * - makes a copy of matrix, and creates temporary struct to use
2346                          */
2347                         cob= constraints_make_evalob(scene, ob, pchan, CONSTRAINT_OBTYPE_BONE);
2348
2349                         /* Solve PoseChannel's Constraints */
2350                         solve_constraints(&pchan->constraints, cob, ctime);     // ctime doesnt alter objects
2351
2352                         /* cleanup after Constraint Solving
2353                          * - applies matrix back to pchan, and frees temporary struct used
2354                          */
2355                         constraints_clear_evalob(cob);
2356
2357                         /* prevent constraints breaking a chain */
2358                         if(pchan->bone->flag & BONE_CONNECTED) {
2359                                 VECCOPY(pchan->pose_mat[3], vec);
2360                         }
2361                 }
2362         }
2363         
2364         /* calculate head */
2365         VECCOPY(pchan->pose_head, pchan->pose_mat[3]);
2366         /* calculate tail */
2367         VECCOPY(vec, pchan->pose_mat[1]);
2368         mul_v3_fl(vec, bone->length);
2369         add_v3_v3v3(pchan->pose_tail, pchan->pose_head, vec);
2370 }
2371
2372 /* This only reads anim data from channels, and writes to channels */
2373 /* This is the only function adding poses */
2374 void where_is_pose (Scene *scene, Object *ob)
2375 {
2376         bArmature *arm;
2377         Bone *bone;
2378         bPoseChannel *pchan;
2379         float imat[4][4];
2380         float ctime;
2381         
2382         if(ob->type!=OB_ARMATURE) return;
2383         arm = ob->data;
2384         
2385         if(ELEM(NULL, arm, scene)) return;
2386         if((ob->pose==NULL) || (ob->pose->flag & POSE_RECALC)) 
2387            armature_rebuild_pose(ob, arm);
2388            
2389         ctime= bsystem_time(scene, ob, (float)scene->r.cfra, 0.0);      /* not accurate... */
2390         
2391         /* In editmode or restposition we read the data from the bones */
2392         if(arm->edbo || (arm->flag & ARM_RESTPOS)) {
2393                 
2394                 for(pchan= ob->pose->chanbase.first; pchan; pchan= pchan->next) {
2395                         bone= pchan->bone;
2396                         if(bone) {
2397                                 copy_m4_m4(pchan->pose_mat, bone->arm_mat);
2398                                 VECCOPY(pchan->pose_head, bone->arm_head);
2399                                 VECCOPY(pchan->pose_tail, bone->arm_tail);
2400                         }
2401                 }
2402         }
2403         else {
2404                 invert_m4_m4(ob->imat, ob->obmat);      // imat is needed 
2405                 
2406                 /* 1. clear flags */
2407                 for(pchan= ob->pose->chanbase.first; pchan; pchan= pchan->next) {
2408                         pchan->flag &= ~(POSE_DONE|POSE_CHAIN|POSE_IKTREE|POSE_IKSPLINE);
2409                 }
2410                 
2411                 /* 2a. construct the IK tree (standard IK) */
2412                 BIK_initialize_tree(scene, ob, ctime);
2413                 
2414                 /* 2b. construct the Spline IK trees 
2415                  *  - this is not integrated as an IK plugin, since it should be able
2416                  *        to function in conjunction with standard IK
2417                  */
2418                 splineik_init_tree(scene, ob, ctime);
2419                 
2420                 /* 3. the main loop, channels are already hierarchical sorted from root to children */
2421                 for(pchan= ob->pose->chanbase.first; pchan; pchan= pchan->next) {
2422                         /* 4a. if we find an IK root, we handle it separated */
2423                         if(pchan->flag & POSE_IKTREE) {
2424                                 BIK_execute_tree(scene, ob, pchan, ctime);
2425                         }
2426                         /* 4b. if we find a Spline IK root, we handle it separated too */
2427                         else if(pchan->flag & POSE_IKSPLINE) {
2428                                 splineik_execute_tree(scene, ob, pchan, ctime);
2429                         }
2430                         /* 5. otherwise just call the normal solver */
2431                         else if(!(pchan->flag & POSE_DONE)) {
2432                                 where_is_pose_bone(scene, ob, pchan, ctime, 1);
2433                         }
2434                 }
2435                 /* 6. release the IK tree */
2436                 BIK_release_tree(scene, ob, ctime);
2437         }
2438                 
2439         /* calculating deform matrices */
2440         for(pchan= ob->pose->chanbase.first; pchan; pchan= pchan->next) {
2441                 if(pchan->bone) {
2442                         invert_m4_m4(imat, pchan->bone->arm_mat);
2443                         mul_m4_m4m4(pchan->chan_mat, imat, pchan->pose_mat);
2444                 }
2445         }
2446 }