Multi-parent constraint

Ok what is the deal? For example when you animate foot sometimes you want to rotate it about forefoot and sometimes about hindfoot. Or you want to animate things like this http://www.youtube.com/watch?v=lOfaFJq5Wqk . I have simple solution.

I assume that reader is familiar with quaternions and their aplication to rotation if not than please see http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation

So let's get mathematical. Beware everything is quaternion(I will try to stay consistent in notaion with wiki page on quaternion, like conjugation is denoted with star)!
$$\begin{align} t & \dots \text{time} \\ p_0(t) & \dots \text{position of the foot} \\ q_0(t) & \dots \text{quaternion which specifies rotation of foot around}  p_0 \\ & \text{now for  } \, i>0 \\ p_i(t) & \dots \text{positions of points you want to rotate around} \\ q_i(t) & \dots \text{quaterion which specifies rotations around corresponding point} \end{align}$$

Suppose you know $p_0,q_0,p_i$ at some time $t_0$ and $q_i(t)$ for all $t>t_0$. Than you want to find $p_0(t),q_0(t),p_i(t)$. But what should they satisfy??

First denote with $v_i$ relative position of $p_i$ to the $p_0$ at time $t$. So:
$$v_i(t) = p_i(t)-p_0(t)$$
Points $p_i$ should stay fixed relative to foot. Therefore they have to satisfy:
$$p_i(t) = p_0(t) + R_0(t)v_i(t_0)R_0^*(t)$$

where $R_i(t) = q_i(t)q_i^{-1}(t_0)$ which represents change in rotation from time $t_0$ to time $t$. .

Now imagine situation when you rotate only about point $p_1$. So $p_i,q_i$  for $i>1$, $p_1$ are constant in time and only $q_1$ changes in time. What do we get? We rotate point $p_0$ around point $p_1$:

$$p_0(t) = p_1(t_0) - R_1(t) v_1(t_0) R_1^*(t)$$
$$q_0(t) = R_1(t)q_0(t_0)$$

Now replace $p_1(t_0)$ with $p_i(t_0) = p_0(t_0) + v_i(t_0)$ we get:

$$p_0(t)-p_0(t_0) =  v_1(t_0) - R_1(t) v_1(t_0) R_1^*(t)$$

So point $p_0$ changes about $v_i(t_0) - R_1(t) v_1(t_0) R_1^*(t)$ thanks to rotation about point $p_1$. We assumed that $p_1$. When $p_1$ varies in time than previous equation is "sort of OK"(it does not represents accurately the rotation of p_0 around point p_1(which varies over time)) only for small times $\Delta t = t-t_0$.

Now suppose that $p_i,q_i$ all varies in time and $\Delta t$ is small. Than we can generalize our update equation:

$$p_0(t_0+\Delta t) - p_0(t_0) = \sum_i v_i(t_0) - R_i(t_0,\Delta t) v_i(t_0) R_i^*(t,\Delta t)$$

where $R(t,\Delta t) = q_i(t+\Delta t)q_i^{-1}(t)$. But how do we update $q_0$ ?? We could do $q_0(t+\Delta t) = R_1(t,\Delta t)\dots R_n(t,\Delta t)q_0(t)$ But it depends on order of the $R_i$ which is undesirable. So we do better with:
$$q_0(t+\Delta t) = \{R_i(t,\Delta t)\}_i q_0(t)$$

where$\{R_i(t,\Delta t)\}$is what I call normalized anticommutator defined by:

$a_1,\dots,a_n$ are any quaternions than

$$\{a_i\}_i =\frac{\sum_{ \sigma \in \Pi_n } a_{\sigma(1)} \dots a_{\sigma(n)}}{||\sum_{ \sigma \in \Pi_n } a_{\sigma(1)} \dots a_{\sigma(n)}||}$$

where $\Pi_n$ is set of all permutations of size $n$. I think that normalized anticommutator is the most convenient way how to combine $n$ rotations and get again rotation. And quaternions are the good way to do it. Problem with matrices is that when you add two rotational matrices than you hardly get rotation matrix and finding closest rotation matrix to that sum is just pain. With quaternions it is easy, just add them up and normalize.

We are almost finished. We just polish those equations a little bit.

denote $t_n = t_0 + n \Delta t$ than our update equations looks like this:

$$\begin{align*} v_i(t_{n+1})& = p_i(t_{n+1})-p_0(t_{n+1}) = R_0(t_n,\Delta t) v_i(t_n) R_0^*(t,\Delta t)\\ p_0(t_{n+1}) &= p_0(t_n) + \sum_i v_i(t_n) - R_i(t_n,\Delta t) v_i(t_n) R_i^*(t_n,\Delta t) \\ q_0(t_{n+1}) &= \{R_i(t_n,\Delta t)\}_i q_0(t_n)  \\ p_i(t_n) &= p_0(t_n) + R_0(t_0,t_n-t_0)v_i(t_0)R_0^*(t_0,t_n-t_0) \end{align*}$$

Horay! finished.

I implemented it to Autodesk Maya. You can download the plugin here:
http://uloz.to/xgp17H2G/multiparent-zip
or
http://www.4shared.com/zip/MqwiOuUt/multiParent_1.html

Watch video how to use it:
http://www.youtube.com/watch?v=5Oip-YiKuik

I hope you like it! If you find any mistakes please let me know ;)

Further investigation:

Send $\Delta t$ to zero and obtain differential equations from those update equations and analyze their behavior!