I have three quaternions A, B, C.
I also have a 3D vector P which describes the barycentric coordinate to interpolate.
I want to interpolate the three quaternions with respect to the barycentric coordinate P
This Image shows how to do it would A, B and C be usual vectors.
but how to calculate the quaternion interpolation? Can I use the "slerp" command somehow? How would the glm code look like?
The weighted arithmetic mean $\mu$, of a set of vectors, minimizes squared L2 norm of the error ($v_i$) to all elements of the set.
$$ \begin{equation} \begin{split} v_i &= \mu - x_i \\ \mu &= \underset{\mu \in \mathbb{R^n}}{\text{argmin}} \; \frac{1}{2}\sum\limits_{v_i} w_i ||v_i||^2_2 \\ \end{split} \end{equation} $$
This empowers us to define special means, where we use different notions of distance.
A reasonable notion of distance for SO(3) is the angle separating two rotations. That's $||\log(\mu \ q^{-1})||$ Your geometry library probably implements this $log$, and expressions for it can be found here. It is not the log one is accustomed to on the real numbers.
Intuition: Recall that every sequence of rotations corresponds to a single rotation about a single axis, so the difference between a pair of rotations can always be summarized by a single angle.
$$ \begin{equation} \begin{split} v_i &= log(\mu \ x_i^{-1}) \\ \mu &= \underset{\mu \in \mathbb{SO3}}{\text{argmin}} \sum\limits_{v_i} w_i || v_i||^2_2 \\ \end{split} \end{equation} $$
There's no analytical form for this minimum for more than two rotations. We can, however, treat this as a least squares problem and compute the result by Gauss-Newton.
To apply Gauss-Newton in general, for some error, we need to be able to find $\mathbf{J_i} = \frac{d}{d\mu}v_i(\mu)$, which requires differentiating $v_i$. Then we can take a step as
$$ \begin{equation} \begin{split} \Delta \mu &= \sum\limits_{i}(w_i \mathbf{J_i}^T \mathbf{J_i})^{-1} \sum\limits_{i} w_i \mathbf{J_i}^T v_i \end{split} \end{equation} $$
Our error function cannot readily be differentiated, since its argument is a rotation $\mu$ and not a vector. We invent a parameter $\delta \in \mathbb{R^3}$ which we will $Exp$ into a rotation. This, we can differentiate
$$ \begin{equation} \begin{split} v_i &= log(Exp(\delta)\; \mu \ x_i^{-1}) \\ \mathbf{J_i} &= \frac{d}{d\delta} v_i = \mathbf{I} \\ \end{split} \end{equation} $$
So the update to $\mu$ is
$$ \begin{equation} \begin{split} \delta &= \sum\limits_{i}(w_i \mathbf{J_i}^T \mathbf{J_i})^{-1} \sum\limits_{i} w_i \mathbf{J_i}^T v_i\\ \delta &= \sum\limits_{i} w_i \sum\limits_{i} w_i v_i \\ \delta &= \sum\limits_{i} w_i v_i \end{split} \end{equation} $$
$$ \begin{equation} \begin{split} \delta = \sum\limits_{i} w_i v_i &= \sum\limits_{i}w_i log(\mu_{k}^{\mathstrut} x_{i}^{-1}) \\ \mu_{k + 1} &= Exp(\delta) \mu_{k}\\ \end{split} \end{equation} $$ This will converge quadratically to a mean.
Caveat emptor: for rotations, there may be many "means". By analogy: What is the average point, on Earth's surface, of the north and south poles? Every point on the equator has equal right to be called the average. This is the case for all equispaced rotations.
