Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, 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 //XXX #include "nla.h"
38
39 #include "BLI_arithb.h"
40 #include "BLI_blenlib.h"
41
42 #include "DNA_armature_types.h"
43 #include "DNA_action_types.h"
44 #include "DNA_constraint_types.h"
45 #include "DNA_mesh_types.h"
46 #include "DNA_lattice_types.h"
47 #include "DNA_meshdata_types.h"
48 #include "DNA_nla_types.h"
49 #include "DNA_object_types.h"
50 #include "DNA_scene_types.h"
51 #include "DNA_view3d_types.h"
52
53 #include "BKE_armature.h"
54 #include "BKE_action.h"
55 #include "BKE_blender.h"
56 #include "BKE_constraint.h"
57 #include "BKE_curve.h"
58 #include "BKE_deform.h"
59 #include "BKE_depsgraph.h"
60 #include "BKE_DerivedMesh.h"
61 #include "BKE_displist.h"
62 #include "BKE_global.h"
63 #include "BKE_library.h"
64 #include "BKE_lattice.h"
65 #include "BKE_main.h"
66 #include "BKE_object.h"
67 #include "BKE_object.h"
68 #include "BKE_utildefines.h"
69 #include "BIK_api.h"
70 #include "BKE_sketch.h"
71
72 #ifdef HAVE_CONFIG_H
73 #include <config.h>
74 #endif
75
76 /*      **************** Generic Functions, data level *************** */
77
78 bArmature *add_armature(char *name)
79 {
80         bArmature *arm;
81         
82         arm= alloc_libblock (&G.main->armature, ID_AR, name);
83         arm->deformflag = ARM_DEF_VGROUP|ARM_DEF_ENVELOPE;
84         arm->flag = ARM_COL_CUSTOM; /* custom bone-group colors */
85         arm->layer= 1;
86         return arm;
87 }
88
89 bArmature *get_armature(Object *ob)
90 {
91         if(ob->type==OB_ARMATURE)
92                 return (bArmature *)ob->data;
93         return NULL;
94 }
95
96 void free_boneChildren(Bone *bone)
97
98         Bone *child;
99         
100         if (bone) {
101                 
102                 child=bone->childbase.first;
103                 if (child){
104                         while (child){
105                                 free_boneChildren (child);
106                                 child=child->next;
107                         }
108                         BLI_freelistN (&bone->childbase);
109                 }
110         }
111 }
112
113 void free_bones (bArmature *arm)
114 {
115         Bone *bone;
116         /*      Free children (if any)  */
117         bone= arm->bonebase.first;
118         if (bone) {
119                 while (bone){
120                         free_boneChildren (bone);
121                         bone=bone->next;
122                 }
123         }
124         
125         
126         BLI_freelistN(&arm->bonebase);
127 }
128
129 void free_armature(bArmature *arm)
130 {
131         if (arm) {
132                 free_bones(arm);
133                 
134                 /* free editmode data */
135                 if (arm->edbo) {
136                         BLI_freelistN(arm->edbo);
137                         
138                         MEM_freeN(arm->edbo);
139                         arm->edbo= NULL;
140                 }
141
142                 /* free sketch */
143                 if (arm->sketch) {
144                         freeSketch(arm->sketch);
145                         arm->sketch = NULL;
146                 }
147         }
148 }
149
150 void make_local_armature(bArmature *arm)
151 {
152         int local=0, lib=0;
153         Object *ob;
154         bArmature *newArm;
155         
156         if (arm->id.lib==0)
157                 return;
158         if (arm->id.us==1) {
159                 arm->id.lib= 0;
160                 arm->id.flag= LIB_LOCAL;
161                 new_id(0, (ID*)arm, 0);
162                 return;
163         }
164         
165         if(local && lib==0) {
166                 arm->id.lib= 0;
167                 arm->id.flag= LIB_LOCAL;
168                 new_id(0, (ID *)arm, 0);
169         }
170         else if(local && lib) {
171                 newArm= copy_armature(arm);
172                 newArm->id.us= 0;
173                 
174                 ob= G.main->object.first;
175                 while(ob) {
176                         if(ob->data==arm) {
177                                 
178                                 if(ob->id.lib==0) {
179                                         ob->data= newArm;
180                                         newArm->id.us++;
181                                         arm->id.us--;
182                                 }
183                         }
184                         ob= ob->id.next;
185                 }
186         }
187 }
188
189 static void     copy_bonechildren (Bone* newBone, Bone* oldBone)
190 {
191         Bone    *curBone, *newChildBone;
192         
193         /*      Copy this bone's list*/
194         BLI_duplicatelist(&newBone->childbase, &oldBone->childbase);
195         
196         /*      For each child in the list, update it's children*/
197         newChildBone=newBone->childbase.first;
198         for (curBone=oldBone->childbase.first;curBone;curBone=curBone->next){
199                 newChildBone->parent=newBone;
200                 copy_bonechildren(newChildBone,curBone);
201                 newChildBone=newChildBone->next;
202         }
203 }
204
205 bArmature *copy_armature(bArmature *arm)
206 {
207         bArmature *newArm;
208         Bone            *oldBone, *newBone;
209         
210         newArm= copy_libblock (arm);
211         BLI_duplicatelist(&newArm->bonebase, &arm->bonebase);
212         
213         /*      Duplicate the childrens' lists*/
214         newBone=newArm->bonebase.first;
215         for (oldBone=arm->bonebase.first;oldBone;oldBone=oldBone->next){
216                 newBone->parent=NULL;
217                 copy_bonechildren (newBone, oldBone);
218                 newBone=newBone->next;
219         };
220         
221         return newArm;
222 }
223
224 static Bone *get_named_bone_bonechildren (Bone *bone, const char *name)
225 {
226         Bone *curBone, *rbone;
227         
228         if (!strcmp (bone->name, name))
229                 return bone;
230         
231         for (curBone=bone->childbase.first; curBone; curBone=curBone->next){
232                 rbone=get_named_bone_bonechildren (curBone, name);
233                 if (rbone)
234                         return rbone;
235         }
236         
237         return NULL;
238 }
239
240
241 Bone *get_named_bone (bArmature *arm, const char *name)
242 /*
243         Walk the list until the bone is found
244  */
245 {
246         Bone *bone=NULL, *curBone;
247         
248         if (!arm) return NULL;
249         
250         for (curBone=arm->bonebase.first; curBone; curBone=curBone->next){
251                 bone = get_named_bone_bonechildren (curBone, name);
252                 if (bone)
253                         return bone;
254         }
255         
256         return bone;
257 }
258
259
260 #define IS_SEPARATOR(a) (a=='.' || a==' ' || a=='-' || a=='_')
261
262 /* finds the best possible flipped name. For renaming; check for unique names afterwards */
263 /* if strip_number: removes number extensions */
264 void bone_flip_name (char *name, int strip_number)
265 {
266         int             len;
267         char    prefix[128]={""};       /* The part before the facing */
268         char    suffix[128]={""};       /* The part after the facing */
269         char    replace[128]={""};      /* The replacement string */
270         char    number[128]={""};       /* The number extension string */
271         char    *index=NULL;
272
273         len= strlen(name);
274         if(len<3) return;       // we don't do names like .R or .L
275
276         /* We first check the case with a .### extension, let's find the last period */
277         if(isdigit(name[len-1])) {
278                 index= strrchr(name, '.');      // last occurrance
279                 if (index && isdigit(index[1]) ) {              // doesnt handle case bone.1abc2 correct..., whatever!
280                         if(strip_number==0) 
281                                 strcpy(number, index);
282                         *index= 0;
283                         len= strlen(name);
284                 }
285         }
286
287         strcpy (prefix, name);
288
289         /* first case; separator . - _ with extensions r R l L  */
290         if( IS_SEPARATOR(name[len-2]) ) {
291                 switch(name[len-1]) {
292                         case 'l':
293                                 prefix[len-1]= 0;
294                                 strcpy(replace, "r");
295                                 break;
296                         case 'r':
297                                 prefix[len-1]= 0;
298                                 strcpy(replace, "l");
299                                 break;
300                         case 'L':
301                                 prefix[len-1]= 0;
302                                 strcpy(replace, "R");
303                                 break;
304                         case 'R':
305                                 prefix[len-1]= 0;
306                                 strcpy(replace, "L");
307                                 break;
308                 }
309         }
310         /* case; beginning with r R l L , with separator after it */
311         else if( IS_SEPARATOR(name[1]) ) {
312                 switch(name[0]) {
313                         case 'l':
314                                 strcpy(replace, "r");
315                                 strcpy(suffix, name+1);
316                                 prefix[0]= 0;
317                                 break;
318                         case 'r':
319                                 strcpy(replace, "l");
320                                 strcpy(suffix, name+1);
321                                 prefix[0]= 0;
322                                 break;
323                         case 'L':
324                                 strcpy(replace, "R");
325                                 strcpy(suffix, name+1);
326                                 prefix[0]= 0;
327                                 break;
328                         case 'R':
329                                 strcpy(replace, "L");
330                                 strcpy(suffix, name+1);
331                                 prefix[0]= 0;
332                                 break;
333                 }
334         }
335         else if(len > 5) {
336                 /* hrms, why test for a separator? lets do the rule 'ultimate left or right' */
337                 index = BLI_strcasestr(prefix, "right");
338                 if (index==prefix || index==prefix+len-5) {
339                         if(index[0]=='r') 
340                                 strcpy (replace, "left");
341                         else {
342                                 if(index[1]=='I') 
343                                         strcpy (replace, "LEFT");
344                                 else
345                                         strcpy (replace, "Left");
346                         }
347                         *index= 0;
348                         strcpy (suffix, index+5);
349                 }
350                 else {
351                         index = BLI_strcasestr(prefix, "left");
352                         if (index==prefix || index==prefix+len-4) {
353                                 if(index[0]=='l') 
354                                         strcpy (replace, "right");
355                                 else {
356                                         if(index[1]=='E') 
357                                                 strcpy (replace, "RIGHT");
358                                         else
359                                                 strcpy (replace, "Right");
360                                 }
361                                 *index= 0;
362                                 strcpy (suffix, index+4);
363                         }
364                 }               
365         }
366
367         sprintf (name, "%s%s%s%s", prefix, replace, suffix, number);
368 }
369
370 /* Finds the best possible extension to the name on a particular axis. (For renaming, check for unique names afterwards)
371  * This assumes that bone names are at most 32 chars long!
372  *      strip_number: removes number extensions  (TODO: not used)
373  *      axis: the axis to name on
374  *      head/tail: the head/tail co-ordinate of the bone on the specified axis
375  */
376 void bone_autoside_name (char *name, int strip_number, short axis, float head, float tail)
377 {
378         unsigned int len;
379         char    basename[32]={""};
380         char    extension[5]={""};
381
382         len= strlen(name);
383         if (len == 0) return;
384         strcpy(basename, name);
385         
386         /* Figure out extension to append: 
387          *      - The extension to append is based upon the axis that we are working on.
388          *      - If head happens to be on 0, then we must consider the tail position as well to decide
389          *        which side the bone is on
390          *              -> If tail is 0, then it's bone is considered to be on axis, so no extension should be added
391          *              -> Otherwise, extension is added from perspective of object based on which side tail goes to
392          *      - If head is non-zero, extension is added from perspective of object based on side head is on
393          */
394         if (axis == 2) {
395                 /* z-axis - vertical (top/bottom) */
396                 if (IS_EQ(head, 0)) {
397                         if (tail < 0)
398                                 strcpy(extension, "Bot");
399                         else if (tail > 0)
400                                 strcpy(extension, "Top");
401                 }
402                 else {
403                         if (head < 0)
404                                 strcpy(extension, "Bot");
405                         else
406                                 strcpy(extension, "Top");
407                 }
408         }
409         else if (axis == 1) {
410                 /* y-axis - depth (front/back) */
411                 if (IS_EQ(head, 0)) {
412                         if (tail < 0)
413                                 strcpy(extension, "Fr");
414                         else if (tail > 0)
415                                 strcpy(extension, "Bk");
416                 }
417                 else {
418                         if (head < 0)
419                                 strcpy(extension, "Fr");
420                         else
421                                 strcpy(extension, "Bk");
422                 }
423         }
424         else {
425                 /* x-axis - horizontal (left/right) */
426                 if (IS_EQ(head, 0)) {
427                         if (tail < 0)
428                                 strcpy(extension, "R");
429                         else if (tail > 0)
430                                 strcpy(extension, "L");
431                 }
432                 else {
433                         if (head < 0)
434                                 strcpy(extension, "R");
435                         else if (head > 0)
436                                 strcpy(extension, "L");
437                 }
438         }
439
440         /* Simple name truncation 
441          *      - truncate if there is an extension and it wouldn't be able to fit
442          *      - otherwise, just append to end
443          */
444         if (extension[0]) {
445                 int change = 1;
446                 
447                 while (change) { /* remove extensions */
448                         change = 0;
449                         if (len > 2 && basename[len-2]=='.') {
450                                 if (basename[len-1]=='L' || basename[len-1] == 'R' ) { /* L R */
451                                         basename[len-2] = '\0';
452                                         len-=2;
453                                         change= 1;
454                                 }
455                         } else if (len > 3 && basename[len-3]=='.') {
456                                 if (    (basename[len-2]=='F' && basename[len-1] == 'r') ||     /* Fr */
457                                                 (basename[len-2]=='B' && basename[len-1] == 'k')        /* Bk */
458                                 ) {
459                                         basename[len-3] = '\0';
460                                         len-=3;
461                                         change= 1;
462                                 }
463                         } else if (len > 4 && basename[len-4]=='.') {
464                                 if (    (basename[len-3]=='T' && basename[len-2]=='o' && basename[len-1] == 'p') ||     /* Top */
465                                                 (basename[len-3]=='B' && basename[len-2]=='o' && basename[len-1] == 't')        /* Bot */
466                                 ) {
467                                         basename[len-4] = '\0';
468                                         len-=4;
469                                         change= 1;
470                                 }
471                         }
472                 }
473                 
474                 if ((32 - len) < strlen(extension) + 1) { /* add 1 for the '.' */
475                         strncpy(name, basename, len-strlen(extension));
476                 }
477         }
478
479         sprintf(name, "%s.%s", basename, extension);
480 }
481
482 /* ************* B-Bone support ******************* */
483
484 #define MAX_BBONE_SUBDIV        32
485
486 /* data has MAX_BBONE_SUBDIV+1 interpolated points, will become desired amount with equal distances */
487 static void equalize_bezier(float *data, int desired)
488 {
489         float *fp, totdist, ddist, dist, fac1, fac2;
490         float pdist[MAX_BBONE_SUBDIV+1];
491         float temp[MAX_BBONE_SUBDIV+1][4];
492         int a, nr;
493         
494         pdist[0]= 0.0f;
495         for(a=0, fp= data; a<MAX_BBONE_SUBDIV; a++, fp+=4) {
496                 QUATCOPY(temp[a], fp);
497                 pdist[a+1]= pdist[a]+VecLenf(fp, fp+4);
498         }
499         /* do last point */
500         QUATCOPY(temp[a], fp);
501         totdist= pdist[a];
502         
503         /* go over distances and calculate new points */
504         ddist= totdist/((float)desired);
505         nr= 1;
506         for(a=1, fp= data+4; a<desired; a++, fp+=4) {
507                 
508                 dist= ((float)a)*ddist;
509                 
510                 /* we're looking for location (distance) 'dist' in the array */
511                 while((dist>= pdist[nr]) && nr<MAX_BBONE_SUBDIV) {
512                         nr++;
513                 }
514                 
515                 fac1= pdist[nr]- pdist[nr-1];
516                 fac2= pdist[nr]-dist;
517                 fac1= fac2/fac1;
518                 fac2= 1.0f-fac1;
519                 
520                 fp[0]= fac1*temp[nr-1][0]+ fac2*temp[nr][0];
521                 fp[1]= fac1*temp[nr-1][1]+ fac2*temp[nr][1];
522                 fp[2]= fac1*temp[nr-1][2]+ fac2*temp[nr][2];
523                 fp[3]= fac1*temp[nr-1][3]+ fac2*temp[nr][3];
524         }
525         /* set last point, needed for orientation calculus */
526         QUATCOPY(fp, temp[MAX_BBONE_SUBDIV]);
527 }
528
529 /* returns pointer to static array, filled with desired amount of bone->segments elements */
530 /* this calculation is done  within unit bone space */
531 Mat4 *b_bone_spline_setup(bPoseChannel *pchan, int rest)
532 {
533         static Mat4 bbone_array[MAX_BBONE_SUBDIV];
534         static Mat4 bbone_rest_array[MAX_BBONE_SUBDIV];
535         Mat4 *result_array= (rest)? bbone_rest_array: bbone_array;
536         bPoseChannel *next, *prev;
537         Bone *bone= pchan->bone;
538         float h1[3], h2[3], scale[3], length, hlength1, hlength2, roll1=0.0f, roll2;
539         float mat3[3][3], imat[4][4], posemat[4][4], scalemat[4][4], iscalemat[4][4];
540         float data[MAX_BBONE_SUBDIV+1][4], *fp;
541         int a, doscale= 0;
542
543         length= bone->length;
544
545         if(!rest) {
546                 /* check if we need to take non-uniform bone scaling into account */
547                 scale[0]= VecLength(pchan->pose_mat[0]);
548                 scale[1]= VecLength(pchan->pose_mat[1]);
549                 scale[2]= VecLength(pchan->pose_mat[2]);
550
551                 if(fabs(scale[0] - scale[1]) > 1e-6f || fabs(scale[1] - scale[2]) > 1e-6f) {
552                         Mat4One(scalemat);
553                         scalemat[0][0]= scale[0];
554                         scalemat[1][1]= scale[1];
555                         scalemat[2][2]= scale[2];
556                         Mat4Invert(iscalemat, scalemat);
557
558                         length *= scale[1];
559                         doscale = 1;
560                 }
561         }
562         
563         hlength1= bone->ease1*length*0.390464f;         // 0.5*sqrt(2)*kappa, the handle length for near-perfect circles
564         hlength2= bone->ease2*length*0.390464f;
565         
566         /* evaluate next and prev bones */
567         if(bone->flag & BONE_CONNECTED)
568                 prev= pchan->parent;
569         else
570                 prev= NULL;
571         
572         next= pchan->child;
573         
574         /* find the handle points, since this is inside bone space, the 
575                 first point = (0,0,0)
576                 last point =  (0, length, 0) */
577         
578         if(rest) {
579                 Mat4Invert(imat, pchan->bone->arm_mat);
580         }
581         else if(doscale) {
582                 Mat4CpyMat4(posemat, pchan->pose_mat);
583                 Mat4Ortho(posemat);
584                 Mat4Invert(imat, posemat);
585         }
586         else
587                 Mat4Invert(imat, pchan->pose_mat);
588         
589         if(prev) {
590                 float difmat[4][4], result[3][3], imat3[3][3];
591
592                 /* transform previous point inside this bone space */
593                 if(rest)
594                         VECCOPY(h1, prev->bone->arm_head)
595                 else
596                         VECCOPY(h1, prev->pose_head)
597                 Mat4MulVecfl(imat, h1);
598
599                 if(prev->bone->segments>1) {
600                         /* if previous bone is B-bone too, use average handle direction */
601                         h1[1]-= length;
602                         roll1= 0.0f;
603                 }
604
605                 Normalize(h1);
606                 VecMulf(h1, -hlength1);
607
608                 if(prev->bone->segments==1) {
609                         /* find the previous roll to interpolate */
610                         if(rest)
611                                 Mat4MulMat4(difmat, prev->bone->arm_mat, imat);
612                         else
613                                 Mat4MulMat4(difmat, prev->pose_mat, imat);
614                         Mat3CpyMat4(result, difmat);                            // the desired rotation at beginning of next bone
615                         
616                         vec_roll_to_mat3(h1, 0.0f, mat3);                       // the result of vec_roll without roll
617                         
618                         Mat3Inv(imat3, mat3);
619                         Mat3MulMat3(mat3, result, imat3);                       // the matrix transforming vec_roll to desired roll
620                         
621                         roll1= (float)atan2(mat3[2][0], mat3[2][2]);
622                 }
623         }
624         else {
625                 h1[0]= 0.0f; h1[1]= hlength1; h1[2]= 0.0f;
626                 roll1= 0.0f;
627         }
628         if(next) {
629                 float difmat[4][4], result[3][3], imat3[3][3];
630                 
631                 /* transform next point inside this bone space */
632                 if(rest)
633                         VECCOPY(h2, next->bone->arm_tail)
634                 else
635                         VECCOPY(h2, next->pose_tail)
636                 Mat4MulVecfl(imat, h2);
637                 /* if next bone is B-bone too, use average handle direction */
638                 if(next->bone->segments>1);
639                 else h2[1]-= length;
640                 Normalize(h2);
641                 
642                 /* find the next roll to interpolate as well */
643                 if(rest)
644                         Mat4MulMat4(difmat, next->bone->arm_mat, imat);
645                 else
646                         Mat4MulMat4(difmat, next->pose_mat, imat);
647                 Mat3CpyMat4(result, difmat);                            // the desired rotation at beginning of next bone
648                 
649                 vec_roll_to_mat3(h2, 0.0f, mat3);                       // the result of vec_roll without roll
650                 
651                 Mat3Inv(imat3, mat3);
652                 Mat3MulMat3(mat3, imat3, result);                       // the matrix transforming vec_roll to desired roll
653                 
654                 roll2= (float)atan2(mat3[2][0], mat3[2][2]);
655                 
656                 /* and only now negate handle */
657                 VecMulf(h2, -hlength2);
658         }
659         else {
660                 h2[0]= 0.0f; h2[1]= -hlength2; h2[2]= 0.0f;
661                 roll2= 0.0;
662         }
663
664         /* make curve */
665         if(bone->segments > MAX_BBONE_SUBDIV)
666                 bone->segments= MAX_BBONE_SUBDIV;
667         
668         forward_diff_bezier(0.0, h1[0],         h2[0],                  0.0,            data[0],        MAX_BBONE_SUBDIV, 4*sizeof(float));
669         forward_diff_bezier(0.0, h1[1],         length + h2[1], length,         data[0]+1,      MAX_BBONE_SUBDIV, 4*sizeof(float));
670         forward_diff_bezier(0.0, h1[2],         h2[2],                  0.0,            data[0]+2,      MAX_BBONE_SUBDIV, 4*sizeof(float));
671         forward_diff_bezier(roll1, roll1 + 0.390464f*(roll2-roll1), roll2 - 0.390464f*(roll2-roll1),    roll2,  data[0]+3,      MAX_BBONE_SUBDIV, 4*sizeof(float));
672         
673         equalize_bezier(data[0], bone->segments);       // note: does stride 4!
674         
675         /* make transformation matrices for the segments for drawing */
676         for(a=0, fp= data[0]; a<bone->segments; a++, fp+=4) {
677                 VecSubf(h1, fp+4, fp);
678                 vec_roll_to_mat3(h1, fp[3], mat3);              // fp[3] is roll
679
680                 Mat4CpyMat3(result_array[a].mat, mat3);
681                 VECCOPY(result_array[a].mat[3], fp);
682
683                 if(doscale) {
684                         /* correct for scaling when this matrix is used in scaled space */
685                         Mat4MulSerie(result_array[a].mat, iscalemat, result_array[a].mat,
686                                 scalemat, NULL, NULL, NULL, NULL, NULL);
687                 }
688         }
689         
690         return result_array;
691 }
692
693 /* ************ Armature Deform ******************* */
694
695 static void pchan_b_bone_defmats(bPoseChannel *pchan, int use_quaternion, int rest_def)
696 {
697         Bone *bone= pchan->bone;
698         Mat4 *b_bone= b_bone_spline_setup(pchan, 0);
699         Mat4 *b_bone_rest= (rest_def)? NULL: b_bone_spline_setup(pchan, 1);
700         Mat4 *b_bone_mats;
701         DualQuat *b_bone_dual_quats= NULL;
702         float tmat[4][4];
703         int a;
704         
705         /* allocate b_bone matrices and dual quats */
706         b_bone_mats= MEM_mallocN((1+bone->segments)*sizeof(Mat4), "BBone defmats");
707         pchan->b_bone_mats= b_bone_mats;
708
709         if(use_quaternion) {
710                 b_bone_dual_quats= MEM_mallocN((bone->segments)*sizeof(DualQuat), "BBone dqs");
711                 pchan->b_bone_dual_quats= b_bone_dual_quats;
712         }
713         
714         /* first matrix is the inverse arm_mat, to bring points in local bone space
715            for finding out which segment it belongs to */
716         Mat4Invert(b_bone_mats[0].mat, bone->arm_mat);
717
718         /* then we make the b_bone_mats:
719             - first transform to local bone space
720                 - translate over the curve to the bbone mat space
721                 - transform with b_bone matrix
722                 - transform back into global space */
723         Mat4One(tmat);
724
725         for(a=0; a<bone->segments; a++) {
726                 if(b_bone_rest)
727                         Mat4Invert(tmat, b_bone_rest[a].mat);
728                 else
729                         tmat[3][1] = -a*(bone->length/(float)bone->segments);
730
731                 Mat4MulSerie(b_bone_mats[a+1].mat, pchan->chan_mat, bone->arm_mat,
732                         b_bone[a].mat, tmat, b_bone_mats[0].mat, NULL, NULL, NULL);
733
734                 if(use_quaternion)
735                         Mat4ToDQuat(bone->arm_mat, b_bone_mats[a+1].mat, &b_bone_dual_quats[a]);
736         }
737 }
738
739 static void b_bone_deform(bPoseChannel *pchan, Bone *bone, float *co, DualQuat *dq, float defmat[][3])
740 {
741         Mat4 *b_bone= pchan->b_bone_mats;
742         float (*mat)[4]= b_bone[0].mat;
743         float segment, y;
744         int a;
745         
746         /* need to transform co back to bonespace, only need y */
747         y= mat[0][1]*co[0] + mat[1][1]*co[1] + mat[2][1]*co[2] + mat[3][1];
748         
749         /* now calculate which of the b_bones are deforming this */
750         segment= bone->length/((float)bone->segments);
751         a= (int)(y/segment);
752         
753         /* note; by clamping it extends deform at endpoints, goes best with
754            straight joints in restpos. */
755         CLAMP(a, 0, bone->segments-1);
756
757         if(dq) {
758                 DQuatCpyDQuat(dq, &((DualQuat*)pchan->b_bone_dual_quats)[a]);
759         }
760         else {
761                 Mat4MulVecfl(b_bone[a+1].mat, co);
762
763                 if(defmat)
764                         Mat3CpyMat4(defmat, b_bone[a+1].mat);
765         }
766 }
767
768 /* using vec with dist to bone b1 - b2 */
769 float distfactor_to_bone (float vec[3], float b1[3], float b2[3], float rad1, float rad2, float rdist)
770 {
771         float dist=0.0f; 
772         float bdelta[3];
773         float pdelta[3];
774         float hsqr, a, l, rad;
775         
776         VecSubf (bdelta, b2, b1);
777         l = Normalize (bdelta);
778         
779         VecSubf (pdelta, vec, b1);
780         
781         a = bdelta[0]*pdelta[0] + bdelta[1]*pdelta[1] + bdelta[2]*pdelta[2];
782         hsqr = ((pdelta[0]*pdelta[0]) + (pdelta[1]*pdelta[1]) + (pdelta[2]*pdelta[2]));
783         
784         if (a < 0.0F){
785                 /* If we're past the end of the bone, do a spherical field attenuation thing */
786                 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])) ;
787                 rad= rad1;
788         }
789         else if (a > l){
790                 /* If we're past the end of the bone, do a spherical field attenuation thing */
791                 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])) ;
792                 rad= rad2;
793         }
794         else {
795                 dist= (hsqr - (a*a));
796                 
797                 if(l!=0.0f) {
798                         rad= a/l;
799                         rad= rad*rad2 + (1.0f-rad)*rad1;
800                 }
801                 else rad= rad1;
802         }
803         
804         a= rad*rad;
805         if(dist < a) 
806                 return 1.0f;
807         else {
808                 l= rad+rdist;
809                 l*= l;
810                 if(rdist==0.0f || dist >= l) 
811                         return 0.0f;
812                 else {
813                         a= (float)sqrt(dist)-rad;
814                         return 1.0f-( a*a )/( rdist*rdist );
815                 }
816         }
817 }
818
819 static void pchan_deform_mat_add(bPoseChannel *pchan, float weight, float bbonemat[][3], float mat[][3])
820 {
821         float wmat[3][3];
822
823         if(pchan->bone->segments>1)
824                 Mat3CpyMat3(wmat, bbonemat);
825         else
826                 Mat3CpyMat4(wmat, pchan->chan_mat);
827
828         Mat3MulFloat((float*)wmat, weight);
829         Mat3AddMat3(mat, mat, wmat);
830 }
831
832 static float dist_bone_deform(bPoseChannel *pchan, float *vec, DualQuat *dq, float mat[][3], float *co)
833 {
834         Bone *bone= pchan->bone;
835         float fac, contrib=0.0;
836         float cop[3], bbonemat[3][3];
837         DualQuat bbonedq;
838
839         if(bone==NULL) return 0.0f;
840         
841         VECCOPY (cop, co);
842
843         fac= distfactor_to_bone(cop, bone->arm_head, bone->arm_tail, bone->rad_head, bone->rad_tail, bone->dist);
844         
845         if (fac>0.0) {
846                 
847                 fac*=bone->weight;
848                 contrib= fac;
849                 if(contrib>0.0) {
850                         if(vec) {
851                                 if(bone->segments>1)
852                                         // applies on cop and bbonemat
853                                         b_bone_deform(pchan, bone, cop, NULL, (mat)?bbonemat:NULL);
854                                 else
855                                         Mat4MulVecfl(pchan->chan_mat, cop);
856
857                                 //      Make this a delta from the base position
858                                 VecSubf (cop, cop, co);
859                                 cop[0]*=fac; cop[1]*=fac; cop[2]*=fac;
860                                 VecAddf (vec, vec, cop);
861
862                                 if(mat)
863                                         pchan_deform_mat_add(pchan, fac, bbonemat, mat);
864                         }
865                         else {
866                                 if(bone->segments>1) {
867                                         b_bone_deform(pchan, bone, cop, &bbonedq, NULL);
868                                         DQuatAddWeighted(dq, &bbonedq, fac);
869                                 }
870                                 else
871                                         DQuatAddWeighted(dq, pchan->dual_quat, fac);
872                         }
873                 }
874         }
875         
876         return contrib;
877 }
878
879 static void pchan_bone_deform(bPoseChannel *pchan, float weight, float *vec, DualQuat *dq, float mat[][3], float *co, float *contrib)
880 {
881         float cop[3], bbonemat[3][3];
882         DualQuat bbonedq;
883
884         if (!weight)
885                 return;
886
887         VECCOPY(cop, co);
888
889         if(vec) {
890                 if(pchan->bone->segments>1)
891                         // applies on cop and bbonemat
892                         b_bone_deform(pchan, pchan->bone, cop, NULL, (mat)?bbonemat:NULL);
893                 else
894                         Mat4MulVecfl(pchan->chan_mat, cop);
895                 
896                 vec[0]+=(cop[0]-co[0])*weight;
897                 vec[1]+=(cop[1]-co[1])*weight;
898                 vec[2]+=(cop[2]-co[2])*weight;
899
900                 if(mat)
901                         pchan_deform_mat_add(pchan, weight, bbonemat, mat);
902         }
903         else {
904                 if(pchan->bone->segments>1) {
905                         b_bone_deform(pchan, pchan->bone, cop, &bbonedq, NULL);
906                         DQuatAddWeighted(dq, &bbonedq, weight);
907                 }
908                 else
909                         DQuatAddWeighted(dq, pchan->dual_quat, weight);
910         }
911
912         (*contrib)+=weight;
913 }
914
915 void armature_deform_verts(Object *armOb, Object *target, DerivedMesh *dm,
916                            float (*vertexCos)[3], float (*defMats)[3][3],
917                                                    int numVerts, int deformflag, 
918                                                    float (*prevCos)[3], const char *defgrp_name)
919 {
920         bArmature *arm= armOb->data;
921         bPoseChannel *pchan, **defnrToPC = NULL;
922         MDeformVert *dverts = NULL;
923         bDeformGroup *dg;
924         DualQuat *dualquats= NULL;
925         float obinv[4][4], premat[4][4], postmat[4][4];
926         int use_envelope = deformflag & ARM_DEF_ENVELOPE;
927         int use_quaternion = deformflag & ARM_DEF_QUATERNION;
928         int bbone_rest_def = deformflag & ARM_DEF_B_BONE_REST;
929         int invert_vgroup= deformflag & ARM_DEF_INVERT_VGROUP;
930         int numGroups = 0;              /* safety for vertexgroup index overflow */
931         int i, target_totvert = 0;      /* safety for vertexgroup overflow */
932         int use_dverts = 0;
933         int armature_def_nr = -1;
934         int totchan;
935
936         if(arm->edbo) return;
937         
938         Mat4Invert(obinv, target->obmat);
939         Mat4CpyMat4(premat, target->obmat);
940         Mat4MulMat4(postmat, armOb->obmat, obinv);
941         Mat4Invert(premat, postmat);
942
943         /* bone defmats are already in the channels, chan_mat */
944         
945         /* initialize B_bone matrices and dual quaternions */
946         if(use_quaternion) {
947                 totchan= BLI_countlist(&armOb->pose->chanbase);
948                 dualquats= MEM_callocN(sizeof(DualQuat)*totchan, "dualquats");
949         }
950
951         totchan= 0;
952         for(pchan = armOb->pose->chanbase.first; pchan; pchan = pchan->next) {
953                 if(!(pchan->bone->flag & BONE_NO_DEFORM)) {
954                         if(pchan->bone->segments > 1)
955                                 pchan_b_bone_defmats(pchan, use_quaternion, bbone_rest_def);
956
957                         if(use_quaternion) {
958                                 pchan->dual_quat= &dualquats[totchan++];
959                                 Mat4ToDQuat(pchan->bone->arm_mat, pchan->chan_mat, pchan->dual_quat);
960                         }
961                 }
962         }
963
964         /* get the def_nr for the overall armature vertex group if present */
965         for(i = 0, dg = target->defbase.first; dg; i++, dg = dg->next)
966                 if(defgrp_name && strcmp(defgrp_name, dg->name) == 0)
967                         armature_def_nr = i;
968
969         /* get a vertex-deform-index to posechannel array */
970         if(deformflag & ARM_DEF_VGROUP) {
971                 if(ELEM(target->type, OB_MESH, OB_LATTICE)) {
972                         numGroups = BLI_countlist(&target->defbase);
973                         
974                         if(target->type==OB_MESH) {
975                                 Mesh *me= target->data;
976                                 dverts = me->dvert;
977                                 target_totvert = me->totvert;
978                         }
979                         else {
980                                 Lattice *lt= target->data;
981                                 dverts = lt->dvert;
982                                 if(dverts)
983                                         target_totvert = lt->pntsu*lt->pntsv*lt->pntsw;
984                         }
985                         /* if we have a DerivedMesh, only use dverts if it has them */
986                         if(dm)
987                                 if(dm->getVertData(dm, 0, CD_MDEFORMVERT))
988                                         use_dverts = 1;
989                                 else use_dverts = 0;
990                         else if(dverts) use_dverts = 1;
991
992                         if(use_dverts) {
993                                 defnrToPC = MEM_callocN(sizeof(*defnrToPC) * numGroups,
994                                                         "defnrToBone");
995                                 for(i = 0, dg = target->defbase.first; dg;
996                                     i++, dg = dg->next) {
997                                         defnrToPC[i] = get_pose_channel(armOb->pose, dg->name);
998                                         /* exclude non-deforming bones */
999                                         if(defnrToPC[i]) {
1000                                                 if(defnrToPC[i]->bone->flag & BONE_NO_DEFORM)
1001                                                         defnrToPC[i]= NULL;
1002                                         }
1003                                 }
1004                         }
1005                 }
1006         }
1007
1008         for(i = 0; i < numVerts; i++) {
1009                 MDeformVert *dvert;
1010                 DualQuat sumdq, *dq = NULL;
1011                 float *co, dco[3];
1012                 float sumvec[3], summat[3][3];
1013                 float *vec = NULL, (*smat)[3] = NULL;
1014                 float contrib = 0.0f;
1015                 float armature_weight = 1.0f;   /* default to 1 if no overall def group */
1016                 float prevco_weight = 1.0f;             /* weight for optional cached vertexcos */
1017                 int       j;
1018
1019                 if(use_quaternion) {
1020                         memset(&sumdq, 0, sizeof(DualQuat));
1021                         dq= &sumdq;
1022                 }
1023                 else {
1024                         sumvec[0] = sumvec[1] = sumvec[2] = 0.0f;
1025                         vec= sumvec;
1026
1027                         if(defMats) {
1028                                 Mat3Clr((float*)summat);
1029                                 smat = summat;
1030                         }
1031                 }
1032
1033                 if(use_dverts || armature_def_nr >= 0) {
1034                         if(dm) dvert = dm->getVertData(dm, i, CD_MDEFORMVERT);
1035                         else if(dverts && i < target_totvert) dvert = dverts + i;
1036                         else dvert = NULL;
1037                 } else
1038                         dvert = NULL;
1039
1040                 if(armature_def_nr >= 0 && dvert) {
1041                         armature_weight = 0.0f; /* a def group was given, so default to 0 */
1042                         for(j = 0; j < dvert->totweight; j++) {
1043                                 if(dvert->dw[j].def_nr == armature_def_nr) {
1044                                         armature_weight = dvert->dw[j].weight;
1045                                         break;
1046                                 }
1047                         }
1048                         /* hackish: the blending factor can be used for blending with prevCos too */
1049                         if(prevCos) {
1050                                 if(invert_vgroup)
1051                                         prevco_weight= 1.0f-armature_weight;
1052                                 else
1053                                         prevco_weight= armature_weight;
1054                                 armature_weight= 1.0f;
1055                         }
1056                 }
1057
1058                 /* check if there's any  point in calculating for this vert */
1059                 if(armature_weight == 0.0f) continue;
1060                 
1061                 /* get the coord we work on */
1062                 co= prevCos?prevCos[i]:vertexCos[i];
1063                 
1064                 /* Apply the object's matrix */
1065                 Mat4MulVecfl(premat, co);
1066                 
1067                 if(use_dverts && dvert && dvert->totweight) { // use weight groups ?
1068                         int deformed = 0;
1069                         
1070                         for(j = 0; j < dvert->totweight; j++){
1071                                 int index = dvert->dw[j].def_nr;
1072                                 pchan = index < numGroups?defnrToPC[index]:NULL;
1073                                 if(pchan) {
1074                                         float weight = dvert->dw[j].weight;
1075                                         Bone *bone = pchan->bone;
1076
1077                                         deformed = 1;
1078                                         
1079                                         if(bone && bone->flag & BONE_MULT_VG_ENV) {
1080                                                 weight *= distfactor_to_bone(co, bone->arm_head,
1081                                                                              bone->arm_tail,
1082                                                                              bone->rad_head,
1083                                                                              bone->rad_tail,
1084                                                                              bone->dist);
1085                                         }
1086                                         pchan_bone_deform(pchan, weight, vec, dq, smat, co, &contrib);
1087                                 }
1088                         }
1089                         /* if there are vertexgroups but not groups with bones
1090                          * (like for softbody groups)
1091                          */
1092                         if(deformed == 0 && use_envelope) {
1093                                 for(pchan = armOb->pose->chanbase.first; pchan;
1094                                     pchan = pchan->next) {
1095                                         if(!(pchan->bone->flag & BONE_NO_DEFORM))
1096                                                 contrib += dist_bone_deform(pchan, vec, dq, smat, co);
1097                                 }
1098                         }
1099                 }
1100                 else if(use_envelope) {
1101                         for(pchan = armOb->pose->chanbase.first; pchan;
1102                             pchan = pchan->next) {
1103                                 if(!(pchan->bone->flag & BONE_NO_DEFORM))
1104                                         contrib += dist_bone_deform(pchan, vec, dq, smat, co);
1105                         }
1106                 }
1107
1108                 /* actually should be EPSILON? weight values and contrib can be like 10e-39 small */
1109                 if(contrib > 0.0001f) {
1110                         if(use_quaternion) {
1111                                 DQuatNormalize(dq, contrib);
1112
1113                                 if(armature_weight != 1.0f) {
1114                                         VECCOPY(dco, co);
1115                                         DQuatMulVecfl(dq, dco, (defMats)? summat: NULL);
1116                                         VecSubf(dco, dco, co);
1117                                         VecMulf(dco, armature_weight);
1118                                         VecAddf(co, co, dco);
1119                                 }
1120                                 else
1121                                         DQuatMulVecfl(dq, co, (defMats)? summat: NULL);
1122
1123                                 smat = summat;
1124                         }
1125                         else {
1126                                 VecMulf(vec, armature_weight/contrib);
1127                                 VecAddf(co, vec, co);
1128                         }
1129
1130                         if(defMats) {
1131                                 float pre[3][3], post[3][3], tmpmat[3][3];
1132
1133                                 Mat3CpyMat4(pre, premat);
1134                                 Mat3CpyMat4(post, postmat);
1135                                 Mat3CpyMat3(tmpmat, defMats[i]);
1136
1137                                 if(!use_quaternion) /* quaternion already is scale corrected */
1138                                         Mat3MulFloat((float*)smat, armature_weight/contrib);
1139
1140                                 Mat3MulSerie(defMats[i], tmpmat, pre, smat, post,
1141                                         NULL, NULL, NULL, NULL);
1142                         }
1143                 }
1144                 
1145                 /* always, check above code */
1146                 Mat4MulVecfl(postmat, co);
1147                 
1148                 
1149                 /* interpolate with previous modifier position using weight group */
1150                 if(prevCos) {
1151                         float mw= 1.0f - prevco_weight;
1152                         vertexCos[i][0]= prevco_weight*vertexCos[i][0] + mw*co[0];
1153                         vertexCos[i][1]= prevco_weight*vertexCos[i][1] + mw*co[1];
1154                         vertexCos[i][2]= prevco_weight*vertexCos[i][2] + mw*co[2];
1155                 }
1156         }
1157
1158         if(dualquats) MEM_freeN(dualquats);
1159         if(defnrToPC) MEM_freeN(defnrToPC);
1160         
1161         /* free B_bone matrices */
1162         for(pchan = armOb->pose->chanbase.first; pchan; pchan = pchan->next) {
1163                 if(pchan->b_bone_mats) {
1164                         MEM_freeN(pchan->b_bone_mats);
1165                         pchan->b_bone_mats = NULL;
1166                 }
1167                 if(pchan->b_bone_dual_quats) {
1168                         MEM_freeN(pchan->b_bone_dual_quats);
1169                         pchan->b_bone_dual_quats = NULL;
1170                 }
1171
1172                 pchan->dual_quat = NULL;
1173         }
1174 }
1175
1176 /* ************ END Armature Deform ******************* */
1177
1178 void get_objectspace_bone_matrix (struct Bone* bone, float M_accumulatedMatrix[][4], int root, int posed)
1179 {
1180         Mat4CpyMat4(M_accumulatedMatrix, bone->arm_mat);
1181 }
1182
1183 /* **************** Space to Space API ****************** */
1184
1185 /* Convert World-Space Matrix to Pose-Space Matrix */
1186 void armature_mat_world_to_pose(Object *ob, float inmat[][4], float outmat[][4]) 
1187 {
1188         float obmat[4][4];
1189         
1190         /* prevent crashes */
1191         if (ob==NULL) return;
1192         
1193         /* get inverse of (armature) object's matrix  */
1194         Mat4Invert(obmat, ob->obmat);
1195         
1196         /* multiply given matrix by object's-inverse to find pose-space matrix */
1197         Mat4MulMat4(outmat, obmat, inmat);
1198 }
1199
1200 /* Convert Wolrd-Space Location to Pose-Space Location
1201  * NOTE: this cannot be used to convert to pose-space location of the supplied
1202  *              pose-channel into its local space (i.e. 'visual'-keyframing) 
1203  */
1204 void armature_loc_world_to_pose(Object *ob, float *inloc, float *outloc) 
1205 {
1206         float xLocMat[4][4];
1207         float nLocMat[4][4];
1208         
1209         /* build matrix for location */
1210         Mat4One(xLocMat);
1211         VECCOPY(xLocMat[3], inloc);
1212
1213         /* get bone-space cursor matrix and extract location */
1214         armature_mat_world_to_pose(ob, xLocMat, nLocMat);
1215         VECCOPY(outloc, nLocMat[3]);
1216 }
1217
1218 /* Convert Pose-Space Matrix to Bone-Space Matrix 
1219  * NOTE: this cannot be used to convert to pose-space transforms of the supplied
1220  *              pose-channel into its local space (i.e. 'visual'-keyframing)
1221  */
1222 void armature_mat_pose_to_bone(bPoseChannel *pchan, float inmat[][4], float outmat[][4])
1223 {
1224         float pc_trans[4][4], inv_trans[4][4];
1225         float pc_posemat[4][4], inv_posemat[4][4];
1226         
1227         /* paranoia: prevent crashes with no pose-channel supplied */
1228         if (pchan==NULL) return;
1229         
1230         /* get the inverse matrix of the pchan's transforms */
1231         if (pchan->rotmode)
1232                 LocEulSizeToMat4(pc_trans, pchan->loc, pchan->eul, pchan->size);
1233         else
1234                 LocQuatSizeToMat4(pc_trans, pchan->loc, pchan->quat, pchan->size);
1235         Mat4Invert(inv_trans, pc_trans);
1236         
1237         /* Remove the pchan's transforms from it's pose_mat.
1238          * This should leave behind the effects of restpose + 
1239          * parenting + constraints
1240          */
1241         Mat4MulMat4(pc_posemat, inv_trans, pchan->pose_mat);
1242         
1243         /* get the inverse of the leftovers so that we can remove 
1244          * that component from the supplied matrix
1245          */
1246         Mat4Invert(inv_posemat, pc_posemat);
1247         
1248         /* get the new matrix */
1249         Mat4MulMat4(outmat, inmat, inv_posemat);
1250 }
1251
1252 /* Convert Pose-Space Location to Bone-Space Location
1253  * NOTE: this cannot be used to convert to pose-space location of the supplied
1254  *              pose-channel into its local space (i.e. 'visual'-keyframing) 
1255  */
1256 void armature_loc_pose_to_bone(bPoseChannel *pchan, float *inloc, float *outloc) 
1257 {
1258         float xLocMat[4][4];
1259         float nLocMat[4][4];
1260         
1261         /* build matrix for location */
1262         Mat4One(xLocMat);
1263         VECCOPY(xLocMat[3], inloc);
1264
1265         /* get bone-space cursor matrix and extract location */
1266         armature_mat_pose_to_bone(pchan, xLocMat, nLocMat);
1267         VECCOPY(outloc, nLocMat[3]);
1268 }
1269
1270 /* Remove rest-position effects from pose-transform for obtaining
1271  * 'visual' transformation of pose-channel. 
1272  * (used by the Visual-Keyframing stuff)
1273  */
1274 void armature_mat_pose_to_delta(float delta_mat[][4], float pose_mat[][4], float arm_mat[][4])
1275 {
1276         float imat[4][4];
1277  
1278         Mat4Invert(imat, arm_mat);
1279         Mat4MulMat4(delta_mat, pose_mat, imat);
1280 }
1281
1282
1283 /* **************** The new & simple (but OK!) armature evaluation ********* */ 
1284
1285 /*  ****************** And how it works! ****************************************
1286
1287   This is the bone transformation trick; they're hierarchical so each bone(b)
1288   is in the coord system of bone(b-1):
1289
1290   arm_mat(b)= arm_mat(b-1) * yoffs(b-1) * d_root(b) * bone_mat(b) 
1291   
1292   -> yoffs is just the y axis translation in parent's coord system
1293   -> d_root is the translation of the bone root, also in parent's coord system
1294
1295   pose_mat(b)= pose_mat(b-1) * yoffs(b-1) * d_root(b) * bone_mat(b) * chan_mat(b)
1296
1297   we then - in init deform - store the deform in chan_mat, such that:
1298
1299   pose_mat(b)= arm_mat(b) * chan_mat(b)
1300   
1301   *************************************************************************** */
1302 /*  Computes vector and roll based on a rotation. "mat" must
1303      contain only a rotation, and no scaling. */ 
1304 void mat3_to_vec_roll(float mat[][3], float *vec, float *roll) 
1305 {
1306     if (vec)
1307         VecCopyf(vec, mat[1]);
1308
1309     if (roll) {
1310         float vecmat[3][3], vecmatinv[3][3], rollmat[3][3];
1311
1312         vec_roll_to_mat3(mat[1], 0.0f, vecmat);
1313         Mat3Inv(vecmatinv, vecmat);
1314         Mat3MulMat3(rollmat, vecmatinv, mat);
1315
1316         *roll= (float)atan2(rollmat[2][0], rollmat[2][2]);
1317     }
1318 }
1319
1320 /*      Calculates the rest matrix of a bone based
1321         On its vector and a roll around that vector */
1322 void vec_roll_to_mat3(float *vec, float roll, float mat[][3])
1323 {
1324         float   nor[3], axis[3], target[3]={0,1,0};
1325         float   theta;
1326         float   rMatrix[3][3], bMatrix[3][3];
1327         
1328         VECCOPY (nor, vec);
1329         Normalize (nor);
1330         
1331         /*      Find Axis & Amount for bone matrix*/
1332         Crossf (axis,target,nor);
1333
1334         if (Inpf(axis,axis) > 0.0000000000001) {
1335                 /* if nor is *not* a multiple of target ... */
1336                 Normalize (axis);
1337                 
1338                 theta= NormalizedVecAngle2(target, nor);
1339                 
1340                 /*      Make Bone matrix*/
1341                 VecRotToMat3(axis, theta, bMatrix);
1342         }
1343         else {
1344                 /* if nor is a multiple of target ... */
1345                 float updown;
1346                 
1347                 /* point same direction, or opposite? */
1348                 updown = ( Inpf (target,nor) > 0 ) ? 1.0f : -1.0f;
1349                 
1350                 /* I think this should work ... */
1351                 bMatrix[0][0]=updown; bMatrix[0][1]=0.0;    bMatrix[0][2]=0.0;
1352                 bMatrix[1][0]=0.0;    bMatrix[1][1]=updown; bMatrix[1][2]=0.0;
1353                 bMatrix[2][0]=0.0;    bMatrix[2][1]=0.0;    bMatrix[2][2]=1.0;
1354         }
1355         
1356         /*      Make Roll matrix*/
1357         VecRotToMat3(nor, roll, rMatrix);
1358         
1359         /*      Combine and output result*/
1360         Mat3MulMat3 (mat, rMatrix, bMatrix);
1361 }
1362
1363
1364 /* recursive part, calculates restposition of entire tree of children */
1365 /* used by exiting editmode too */
1366 void where_is_armature_bone(Bone *bone, Bone *prevbone)
1367 {
1368         float vec[3];
1369         
1370         /* Bone Space */
1371         VecSubf (vec, bone->tail, bone->head);
1372         vec_roll_to_mat3(vec, bone->roll, bone->bone_mat);
1373
1374         bone->length= VecLenf(bone->head, bone->tail);
1375         
1376         /* this is called on old file reading too... */
1377         if(bone->xwidth==0.0) {
1378                 bone->xwidth= 0.1f;
1379                 bone->zwidth= 0.1f;
1380                 bone->segments= 1;
1381         }
1382         
1383         if(prevbone) {
1384                 float offs_bone[4][4];  // yoffs(b-1) + root(b) + bonemat(b)
1385                 
1386                 /* bone transform itself */
1387                 Mat4CpyMat3(offs_bone, bone->bone_mat);
1388                                 
1389                 /* The bone's root offset (is in the parent's coordinate system) */
1390                 VECCOPY(offs_bone[3], bone->head);
1391
1392                 /* Get the length translation of parent (length along y axis) */
1393                 offs_bone[3][1]+= prevbone->length;
1394                 
1395                 /* Compose the matrix for this bone  */
1396                 Mat4MulMat4(bone->arm_mat, offs_bone, prevbone->arm_mat);
1397         }
1398         else {
1399                 Mat4CpyMat3(bone->arm_mat, bone->bone_mat);
1400                 VECCOPY(bone->arm_mat[3], bone->head);
1401         }
1402         
1403         /* head */
1404         VECCOPY(bone->arm_head, bone->arm_mat[3]);
1405         /* tail is in current local coord system */
1406         VECCOPY(vec, bone->arm_mat[1]);
1407         VecMulf(vec, bone->length);
1408         VecAddf(bone->arm_tail, bone->arm_head, vec);
1409         
1410         /* and the kiddies */
1411         prevbone= bone;
1412         for(bone= bone->childbase.first; bone; bone= bone->next) {
1413                 where_is_armature_bone(bone, prevbone);
1414         }
1415 }
1416
1417 /* updates vectors and matrices on rest-position level, only needed 
1418    after editing armature itself, now only on reading file */
1419 void where_is_armature (bArmature *arm)
1420 {
1421         Bone *bone;
1422         
1423         /* hierarchical from root to children */
1424         for(bone= arm->bonebase.first; bone; bone= bone->next) {
1425                 where_is_armature_bone(bone, NULL);
1426         }
1427 }
1428
1429 /* if bone layer is protected, copy the data from from->pose */
1430 static void pose_proxy_synchronize(Object *ob, Object *from, int layer_protected)
1431 {
1432         bPose *pose= ob->pose, *frompose= from->pose;
1433         bPoseChannel *pchan, *pchanp, pchanw;
1434         bConstraint *con;
1435         
1436         if (frompose==NULL) return;
1437         
1438         /* exception, armature local layer should be proxied too */
1439         if (pose->proxy_layer)
1440                 ((bArmature *)ob->data)->layer= pose->proxy_layer;
1441         
1442         /* clear all transformation values from library */
1443         rest_pose(frompose);
1444         
1445         /* copy over all of the proxy's bone groups */
1446                 /* TODO for later - implement 'local' bone groups as for constraints
1447                  *      Note: this isn't trivial, as bones reference groups by index not by pointer, 
1448                  *               so syncing things correctly needs careful attention
1449                  */
1450         BLI_freelistN(&pose->agroups);
1451         BLI_duplicatelist(&pose->agroups, &frompose->agroups);
1452         pose->active_group= frompose->active_group;
1453         
1454         for (pchan= pose->chanbase.first; pchan; pchan= pchan->next) {
1455                 if (pchan->bone->layer & layer_protected) {
1456                         ListBase proxylocal_constraints = {NULL, NULL};
1457                         pchanp= get_pose_channel(frompose, pchan->name);
1458                         
1459                         /* copy posechannel to temp, but restore important pointers */
1460                         pchanw= *pchanp;
1461                         pchanw.prev= pchan->prev;
1462                         pchanw.next= pchan->next;
1463                         pchanw.parent= pchan->parent;
1464                         pchanw.child= pchan->child;
1465                         pchanw.path= NULL;
1466                         
1467                         /* constraints - proxy constraints are flushed... local ones are added after 
1468                          *      1. extract constraints not from proxy (CONSTRAINT_PROXY_LOCAL) from pchan's constraints
1469                          *      2. copy proxy-pchan's constraints on-to new
1470                          *      3. add extracted local constraints back on top 
1471                          */
1472                         extract_proxylocal_constraints(&proxylocal_constraints, &pchan->constraints);
1473                         copy_constraints(&pchanw.constraints, &pchanp->constraints);
1474                         addlisttolist(&pchanw.constraints, &proxylocal_constraints);
1475                         
1476                         /* constraints - set target ob pointer to own object */
1477                         for (con= pchanw.constraints.first; con; con= con->next) {
1478                                 bConstraintTypeInfo *cti= constraint_get_typeinfo(con);
1479                                 ListBase targets = {NULL, NULL};
1480                                 bConstraintTarget *ct;
1481                                 
1482                                 if (cti && cti->get_constraint_targets) {
1483                                         cti->get_constraint_targets(con, &targets);
1484                                         
1485                                         for (ct= targets.first; ct; ct= ct->next) {
1486                                                 if (ct->tar == from)
1487                                                         ct->tar = ob;
1488                                         }
1489                                         
1490                                         if (cti->flush_constraint_targets)
1491                                                 cti->flush_constraint_targets(con, &targets, 0);
1492                                 }
1493                         }
1494                         
1495                         /* free stuff from current channel */
1496                         if (pchan->path) MEM_freeN(pchan->path);
1497                         free_constraints(&pchan->constraints);
1498                         
1499                         /* the final copy */
1500                         *pchan= pchanw;
1501                 }
1502         }
1503 }
1504
1505 static int rebuild_pose_bone(bPose *pose, Bone *bone, bPoseChannel *parchan, int counter)
1506 {
1507         bPoseChannel *pchan = verify_pose_channel (pose, bone->name);   // verify checks and/or adds
1508
1509         pchan->bone= bone;
1510         pchan->parent= parchan;
1511         
1512         counter++;
1513         
1514         for(bone= bone->childbase.first; bone; bone= bone->next) {
1515                 counter= rebuild_pose_bone(pose, bone, pchan, counter);
1516                 /* for quick detecting of next bone in chain, only b-bone uses it now */
1517                 if(bone->flag & BONE_CONNECTED)
1518                         pchan->child= get_pose_channel(pose, bone->name);
1519         }
1520         
1521         return counter;
1522 }
1523
1524 /* only after leave editmode, duplicating, validating older files, library syncing */
1525 /* NOTE: pose->flag is set for it */
1526 void armature_rebuild_pose(Object *ob, bArmature *arm)
1527 {
1528         Bone *bone;
1529         bPose *pose;
1530         bPoseChannel *pchan, *next;
1531         int counter=0;
1532                 
1533         /* only done here */
1534         if(ob->pose==NULL) ob->pose= MEM_callocN(sizeof(bPose), "new pose");
1535         pose= ob->pose;
1536         
1537         /* clear */
1538         for(pchan= pose->chanbase.first; pchan; pchan= pchan->next) {
1539                 pchan->bone= NULL;
1540                 pchan->child= NULL;
1541         }
1542         
1543         /* first step, check if all channels are there */
1544         for(bone= arm->bonebase.first; bone; bone= bone->next) {
1545                 counter= rebuild_pose_bone(pose, bone, NULL, counter);
1546         }
1547
1548         /* and a check for garbage */
1549         for(pchan= pose->chanbase.first; pchan; pchan= next) {
1550                 next= pchan->next;
1551                 if(pchan->bone==NULL) {
1552                         if(pchan->path)
1553                                 MEM_freeN(pchan->path);
1554                         free_constraints(&pchan->constraints);
1555                         BLI_freelinkN(&pose->chanbase, pchan);
1556                 }
1557         }
1558         // printf("rebuild pose %s, %d bones\n", ob->id.name, counter);
1559         
1560         /* synchronize protected layers with proxy */
1561         if(ob->proxy)
1562                 pose_proxy_synchronize(ob, ob->proxy, arm->layer_protected);
1563         
1564         update_pose_constraint_flags(ob->pose); // for IK detection for example
1565         
1566         /* the sorting */
1567         if(counter>1)
1568                 DAG_pose_sort(ob);
1569         
1570         ob->pose->flag &= ~POSE_RECALC;
1571         ob->pose->flag |= POSE_WAS_REBUILT;
1572 }
1573
1574
1575 /* ********************** THE POSE SOLVER ******************* */
1576
1577
1578 /* loc/rot/size to mat4 */
1579 /* used in constraint.c too */
1580 void chan_calc_mat(bPoseChannel *chan)
1581 {
1582         float smat[3][3];
1583         float rmat[3][3];
1584         float tmat[3][3];
1585         
1586         /* get scaling matrix */
1587         SizeToMat3(chan->size, smat);
1588         
1589         /* rotations may either be quats, eulers (with various rotation orders), or axis-angle */
1590         if (chan->rotmode > 0) {
1591                 /* euler rotations (will cause gimble lock, but this can be alleviated a bit with rotation orders) */
1592                 EulOToMat3(chan->eul, chan->rotmode, rmat);
1593         }
1594         else if (chan->rotmode == PCHAN_ROT_AXISANGLE) {
1595                 /* axis-angle - stored in quaternion data, but not really that great for 3D-changing orientations */
1596                 AxisAngleToMat3(&chan->quat[1], chan->quat[0], rmat);
1597         }
1598         else {
1599                 /* quats are normalised before use to eliminate scaling issues */
1600                 NormalQuat(chan->quat);
1601                 QuatToMat3(chan->quat, rmat);
1602         }
1603         
1604         /* calculate matrix of bone (as 3x3 matrix, but then copy the 4x4) */
1605         Mat3MulMat3(tmat, rmat, smat);
1606         Mat4CpyMat3(chan->chan_mat, tmat);
1607         
1608         /* prevent action channels breaking chains */
1609         /* need to check for bone here, CONSTRAINT_TYPE_ACTION uses this call */
1610         if ((chan->bone==NULL) || !(chan->bone->flag & BONE_CONNECTED)) {
1611                 VECCOPY(chan->chan_mat[3], chan->loc);
1612         }
1613 }
1614
1615 /* NLA strip modifiers */
1616 static void do_strip_modifiers(Scene *scene, Object *armob, Bone *bone, bPoseChannel *pchan)
1617 {
1618         bActionModifier *amod;
1619         bActionStrip *strip, *strip2;
1620         float scene_cfra= (float)scene->r.cfra;
1621         int do_modif;
1622
1623         for (strip=armob->nlastrips.first; strip; strip=strip->next) {
1624                 do_modif=0;
1625                 
1626                 if (scene_cfra>=strip->start && scene_cfra<=strip->end)
1627                         do_modif=1;
1628                 
1629                 if ((scene_cfra > strip->end) && (strip->flag & ACTSTRIP_HOLDLASTFRAME)) {
1630                         do_modif=1;
1631                         
1632                         /* if there are any other strips active, ignore modifiers for this strip - 
1633                          * 'hold' option should only hold action modifiers if there are 
1634                          * no other active strips */
1635                         for (strip2=strip->next; strip2; strip2=strip2->next) {
1636                                 if (strip2 == strip) continue;
1637                                 
1638                                 if (scene_cfra>=strip2->start && scene_cfra<=strip2->end) {
1639                                         if (!(strip2->flag & ACTSTRIP_MUTE))
1640                                                 do_modif=0;
1641                                 }
1642                         }
1643                         
1644                         /* if there are any later, activated, strips with 'hold' set, they take precedence, 
1645                          * so ignore modifiers for this strip */
1646                         for (strip2=strip->next; strip2; strip2=strip2->next) {
1647                                 if (scene_cfra < strip2->start) continue;
1648                                 if ((strip2->flag & ACTSTRIP_HOLDLASTFRAME) && !(strip2->flag & ACTSTRIP_MUTE)) {
1649                                         do_modif=0;
1650                                 }
1651                         }
1652                 }
1653                 
1654                 if (do_modif) {
1655                         /* temporal solution to prevent 2 strips accumulating */
1656                         if(scene_cfra==strip->end && strip->next && strip->next->start==scene_cfra)
1657                                 continue;
1658                         
1659                         for(amod= strip->modifiers.first; amod; amod= amod->next) {
1660                                 switch (amod->type) {
1661                                 case ACTSTRIP_MOD_DEFORM:
1662                                 {
1663                                         /* validate first */
1664                                         if(amod->ob && amod->ob->type==OB_CURVE && amod->channel[0]) {
1665                                                 
1666                                                 if( strcmp(pchan->name, amod->channel)==0 ) {
1667                                                         float mat4[4][4], mat3[3][3];
1668                                                         
1669                                                         curve_deform_vector(scene, amod->ob, armob, bone->arm_mat[3], pchan->pose_mat[3], mat3, amod->no_rot_axis);
1670                                                         Mat4CpyMat4(mat4, pchan->pose_mat);
1671                                                         Mat4MulMat34(pchan->pose_mat, mat3, mat4);
1672                                                         
1673                                                 }
1674                                         }
1675                                 }
1676                                         break;
1677                                 case ACTSTRIP_MOD_NOISE:        
1678                                 {
1679                                         if( strcmp(pchan->name, amod->channel)==0 ) {
1680                                                 float nor[3], loc[3], ofs;
1681                                                 float eul[3], size[3], eulo[3], sizeo[3];
1682                                                 
1683                                                 /* calculate turbulance */
1684                                                 ofs = amod->turbul / 200.0f;
1685                                                 
1686                                                 /* make a copy of starting conditions */
1687                                                 VECCOPY(loc, pchan->pose_mat[3]);
1688                                                 Mat4ToEul(pchan->pose_mat, eul);
1689                                                 Mat4ToSize(pchan->pose_mat, size);
1690                                                 VECCOPY(eulo, eul);
1691                                                 VECCOPY(sizeo, size);
1692                                                 
1693                                                 /* apply noise to each set of channels */
1694                                                 if (amod->channels & 4) {
1695                                                         /* for scaling */
1696                                                         nor[0] = BLI_gNoise(amod->noisesize, size[0]+ofs, size[1], size[2], 0, 0) - ofs;
1697                                                         nor[1] = BLI_gNoise(amod->noisesize, size[0], size[1]+ofs, size[2], 0, 0) - ofs;        
1698                                                         nor[2] = BLI_gNoise(amod->noisesize, size[0], size[1], size[2]+ofs, 0, 0) - ofs;
1699                                                         VecAddf(size, size, nor);
1700                                                         
1701                                                         if (sizeo[0] != 0)
1702                                                                 VecMulf(pchan->pose_mat[0], size[0] / sizeo[0]);
1703                                                         if (sizeo[1] != 0)
1704                                                                 VecMulf(pchan->pose_mat[1], size[1] / sizeo[1]);
1705                                                         if (sizeo[2] != 0)
1706                                                                 VecMulf(pchan->pose_mat[2], size[2] / sizeo[2]);
1707                                                 }
1708                                                 if (amod->channels & 2) {
1709                                                         /* for rotation */
1710                                                         nor[0] = BLI_gNoise(amod->noisesize, eul[0]+ofs, eul[1], eul[2], 0, 0) - ofs;
1711                                                         nor[1] = BLI_gNoise(amod->noisesize, eul[0], eul[1]+ofs, eul[2], 0, 0) - ofs;   
1712                                                         nor[2] = BLI_gNoise(amod->noisesize, eul[0], eul[1], eul[2]+ofs, 0, 0) - ofs;
1713                                                         
1714                                                         compatible_eul(nor, eulo);
1715                                                         VecAddf(eul, eul, nor);
1716                                                         compatible_eul(eul, eulo);
1717                                                         
1718                                                         LocEulSizeToMat4(pchan->pose_mat, loc, eul, size);
1719                                                 }
1720                                                 if (amod->channels & 1) {
1721                                                         /* for location */
1722                                                         nor[0] = BLI_gNoise(amod->noisesize, loc[0]+ofs, loc[1], loc[2], 0, 0) - ofs;
1723                                                         nor[1] = BLI_gNoise(amod->noisesize, loc[0], loc[1]+ofs, loc[2], 0, 0) - ofs;   
1724                                                         nor[2] = BLI_gNoise(amod->noisesize, loc[0], loc[1], loc[2]+ofs, 0, 0) - ofs;
1725                                                         
1726                                                         VecAddf(pchan->pose_mat[3], loc, nor);
1727                                                 }
1728                                         }
1729                                 }
1730                                         break;
1731                                 }
1732                         }
1733                 }
1734         }
1735 }
1736
1737
1738 /* The main armature solver, does all constraints excluding IK */
1739 /* pchan is validated, as having bone and parent pointer */
1740 void where_is_pose_bone(Scene *scene, Object *ob, bPoseChannel *pchan, float ctime)
1741 {
1742         Bone *bone, *parbone;
1743         bPoseChannel *parchan;
1744         float vec[3];
1745         
1746         /* set up variables for quicker access below */
1747         bone= pchan->bone;
1748         parbone= bone->parent;
1749         parchan= pchan->parent;
1750         
1751         /* this gives a chan_mat with actions (ipos) results */
1752         chan_calc_mat(pchan);
1753         
1754         /* construct the posemat based on PoseChannels, that we do before applying constraints */
1755         /* pose_mat(b)= pose_mat(b-1) * yoffs(b-1) * d_root(b) * bone_mat(b) * chan_mat(b) */
1756         
1757         if(parchan) {
1758                 float offs_bone[4][4];  // yoffs(b-1) + root(b) + bonemat(b)
1759                 
1760                 /* bone transform itself */
1761                 Mat4CpyMat3(offs_bone, bone->bone_mat);
1762                 
1763                 /* The bone's root offset (is in the parent's coordinate system) */
1764                 VECCOPY(offs_bone[3], bone->head);
1765                 
1766                 /* Get the length translation of parent (length along y axis) */
1767                 offs_bone[3][1]+= parbone->length;
1768                 
1769                 /* Compose the matrix for this bone  */
1770                 if(bone->flag & BONE_HINGE) {   // uses restposition rotation, but actual position
1771                         float tmat[4][4];
1772                         
1773                         /* the rotation of the parent restposition */
1774                         Mat4CpyMat4(tmat, parbone->arm_mat);
1775                         
1776                         /* the location of actual parent transform */
1777                         VECCOPY(tmat[3], offs_bone[3]);
1778                         offs_bone[3][0]= offs_bone[3][1]= offs_bone[3][2]= 0.0f;
1779                         Mat4MulVecfl(parchan->pose_mat, tmat[3]);
1780                         
1781                         Mat4MulSerie(pchan->pose_mat, tmat, offs_bone, pchan->chan_mat, NULL, NULL, NULL, NULL, NULL);
1782                 }
1783                 else if(bone->flag & BONE_NO_SCALE) {
1784                         float orthmat[4][4];
1785                         
1786                         /* get the official transform, but we only use the vector from it (optimize...) */
1787                         Mat4MulSerie(pchan->pose_mat, parchan->pose_mat, offs_bone, pchan->chan_mat, NULL, NULL, NULL, NULL, NULL);
1788                         VECCOPY(vec, pchan->pose_mat[3]);
1789                         
1790                         /* do this again, but with an ortho-parent matrix */
1791                         Mat4CpyMat4(orthmat, parchan->pose_mat);
1792                         Mat4Ortho(orthmat);
1793                         Mat4MulSerie(pchan->pose_mat, orthmat, offs_bone, pchan->chan_mat, NULL, NULL, NULL, NULL, NULL);
1794                         
1795                         /* copy correct transform */
1796                         VECCOPY(pchan->pose_mat[3], vec);
1797                 }
1798                 else 
1799                         Mat4MulSerie(pchan->pose_mat, parchan->pose_mat, offs_bone, pchan->chan_mat, NULL, NULL, NULL, NULL, NULL);
1800         }
1801         else {
1802                 Mat4MulMat4(pchan->pose_mat, pchan->chan_mat, bone->arm_mat);
1803                 
1804                 /* only rootbones get the cyclic offset (unless user doesn't want that) */
1805                 if ((bone->flag & BONE_NO_CYCLICOFFSET) == 0)
1806                         VecAddf(pchan->pose_mat[3], pchan->pose_mat[3], ob->pose->cyclic_offset);
1807         }
1808         
1809         /* do NLA strip modifiers - i.e. curve follow */
1810         do_strip_modifiers(scene, ob, bone, pchan);
1811         
1812         /* Do constraints */
1813         if (pchan->constraints.first) {
1814                 bConstraintOb *cob;
1815                 
1816                 /* make a copy of location of PoseChannel for later */
1817                 VECCOPY(vec, pchan->pose_mat[3]);
1818                 
1819                 /* prepare PoseChannel for Constraint solving 
1820                  * - makes a copy of matrix, and creates temporary struct to use 
1821                  */
1822                 cob= constraints_make_evalob(scene, ob, pchan, CONSTRAINT_OBTYPE_BONE);
1823                 
1824                 /* Solve PoseChannel's Constraints */
1825                 solve_constraints(&pchan->constraints, cob, ctime);     // ctime doesnt alter objects
1826                 
1827                 /* cleanup after Constraint Solving 
1828                  * - applies matrix back to pchan, and frees temporary struct used
1829                  */
1830                 constraints_clear_evalob(cob);
1831                 
1832                 /* prevent constraints breaking a chain */
1833                 if(pchan->bone->flag & BONE_CONNECTED) {
1834                         VECCOPY(pchan->pose_mat[3], vec);
1835                 }
1836         }
1837         
1838         /* calculate head */
1839         VECCOPY(pchan->pose_head, pchan->pose_mat[3]);
1840         /* calculate tail */
1841         VECCOPY(vec, pchan->pose_mat[1]);
1842         VecMulf(vec, bone->length);
1843         VecAddf(pchan->pose_tail, pchan->pose_head, vec);
1844 }
1845
1846 /* This only reads anim data from channels, and writes to channels */
1847 /* This is the only function adding poses */
1848 void where_is_pose (Scene *scene, Object *ob)
1849 {
1850         bArmature *arm;
1851         Bone *bone;
1852         bPoseChannel *pchan;
1853         float imat[4][4];
1854         float ctime;
1855         
1856         if(ob->type!=OB_ARMATURE) return;
1857         arm = ob->data;
1858         
1859         if(ELEM(NULL, arm, scene)) return;
1860         if((ob->pose==NULL) || (ob->pose->flag & POSE_RECALC)) 
1861            armature_rebuild_pose(ob, arm);
1862            
1863         ctime= bsystem_time(scene, ob, (float)scene->r.cfra, 0.0);      /* not accurate... */
1864         
1865         /* In editmode or restposition we read the data from the bones */
1866         if(arm->edbo || (arm->flag & ARM_RESTPOS)) {
1867                 
1868                 for(pchan= ob->pose->chanbase.first; pchan; pchan= pchan->next) {
1869                         bone= pchan->bone;
1870                         if(bone) {
1871                                 Mat4CpyMat4(pchan->pose_mat, bone->arm_mat);
1872                                 VECCOPY(pchan->pose_head, bone->arm_head);
1873                                 VECCOPY(pchan->pose_tail, bone->arm_tail);
1874                         }
1875                 }
1876         }
1877         else {
1878                 Mat4Invert(ob->imat, ob->obmat);        // imat is needed 
1879
1880                 /* 1. clear flags */
1881                 for(pchan= ob->pose->chanbase.first; pchan; pchan= pchan->next) {
1882                         pchan->flag &= ~(POSE_DONE|POSE_CHAIN|POSE_IKTREE);
1883                 }
1884                 /* 2. construct the IK tree */
1885                 BIK_initialize_tree(scene, ob, ctime);
1886
1887                 /* 3. the main loop, channels are already hierarchical sorted from root to children */
1888                 for(pchan= ob->pose->chanbase.first; pchan; pchan= pchan->next) {
1889                         
1890                         /* 4. if we find an IK root, we handle it separated */
1891                         if(pchan->flag & POSE_IKTREE) {
1892                                 BIK_execute_tree(scene, ob, pchan, ctime);
1893                         }
1894                         /* 5. otherwise just call the normal solver */
1895                         else if(!(pchan->flag & POSE_DONE)) {
1896                                 where_is_pose_bone(scene, ob, pchan, ctime);
1897                         }
1898                 }
1899                 /* 6. release the IK tree */
1900                 BIK_release_tree(scene, ob, ctime);
1901         }
1902                 
1903         /* calculating deform matrices */
1904         for(pchan= ob->pose->chanbase.first; pchan; pchan= pchan->next) {
1905                 if(pchan->bone) {
1906                         Mat4Invert(imat, pchan->bone->arm_mat);
1907                         Mat4MulMat4(pchan->chan_mat, imat, pchan->pose_mat);
1908                 }
1909         }
1910 }