First Thought (probably not the fastest)
Let us assume you have a vector space in $R^{3}$ with a quaternion defined as:
$$
\mathbf{q} = q_{F}^{*} \ q_{B} \\
= a + b \hat{\mathbf{x}} + c \hat{\mathbf{y}} + d \hat{\mathbf{z}}
$$
where $(a, b, c, d)$ are the Euler parameters and $(\hat{\mathbf{x}}, \hat{\mathbf{y}}, \hat{\mathbf{z}})$ defines the reference unit basis set.
If we define the axis of rotation as $\mathbf{n}$ and the angle through which we rotate as $\zeta$, then the Euler parameters are defined as:
$$
a = \cos{\left( \frac{\zeta}{2} \right)} \\
b = n_{x} \ \sin{\left( \frac{\zeta}{2} \right)} \\
c = n_{y} \ \sin{\left( \frac{\zeta}{2} \right)} \\
d = n_{z} \ \sin{\left( \frac{\zeta}{2} \right)} \\
$$
Thus, if you know $\mathbf{q}$, or rather $(a, b, c, d)$, you can find $\mathbf{n}$ and $\zeta$. Once you know the axis of rotation and the angle of rotation, you can determine the Euler angles. First we define the cross product matrix as:
$$
\left[ \mathbf{n} \right]_{x} = \left[
\begin{array}{ c c c }
0 & - n_{z} & n_{y} \\
n_{z} & 0 & - n_{x} \\
- n_{y} & n_{x} & 0
\end{array} \right]
$$
and the outer product of $\mathbf{n}$ with itself given by:
$$
\left[ \mathbf{n} \otimes \mathbf{n} \right] = \left[
\begin{array}{ c c c }
n_{x} \ n_{x} & n_{x} \ n_{y} & n_{x} \ n_{z} \\
n_{y} \ n_{x} & n_{y} \ n_{y} & n_{y} \ n_{z} \\
n_{z} \ n_{x} & n_{z} \ n_{y} & n_{z} \ n_{z}
\end{array} \right]
$$
Then we can define the rotation matrix as:
$$
\overleftrightarrow{\mathbf{R}} = \cos{\zeta} \ \overleftrightarrow{\mathbf{I}} + \sin{\zeta} \ \left[ \mathbf{n} \right]_{x} + \left( 1 - \cos{\zeta} \right) \ \left[ \mathbf{n} \otimes \mathbf{n} \right]
$$
where $\overleftrightarrow{\mathbf{I}}$ is the unit or identity matrix.
Second Thought (probably faster/easier)
An easier method is to follow the procedure given here. Following that procedure, we define:
$$
\alpha = \frac{ 2 \left( a \ b + c \ d \right) }{ 1 - 2 \left( b^{2} + c^{2} \right) } \\
\beta = 2 \left( a \ c - d \ b \right) \\
\gamma = \frac{ 2 \left( a \ d + b \ c \right) }{ 1 - 2 \left( c^{2} + d^{2} \right) }
$$
which gives us the Euler angles:
$$
\phi = \tan^{-1}{ \alpha } \\
\theta = \sin^{-1}{ \beta } \\
\psi = \tan^{-1}{ \gamma }
$$
Since you already have $\mathbf{q}$ and you can numerically/analytically determine $(\phi, \theta, \psi)$, then I would just take the time derivative of each of these angles to find $(\dot{\phi}, \dot{\theta}, \dot{\psi})$ rather than using the angular velocities.
The derivation you quoted assumes the angles $\psi,\theta,\phi$ depend on time, i.e.
$$\begin{cases}\psi=\psi(t)\\\theta=\theta(t)\\\phi=\phi(t)\end{cases}$$
In this case, at the moment $t+\delta t$ you can approximate
$$\begin{cases}\psi(t+\delta t)=\psi(t)+\delta\psi(t)\\\theta(t+\delta t)=\theta(t)+\delta\theta(t)\\\phi(t+\delta t)=\phi(t)+\delta\phi(t)\end{cases}$$
In other words, an infinitesimal increment in time leads to infinitesimal increments of the angles captured by the vector
$$\delta\boldsymbol{E}=\begin{pmatrix}\delta\psi\\\delta\theta\\\delta\phi\end{pmatrix}$$
Best Answer
Using the basis-vectors of the inertial frame, rotational kinetic energy (assuming center of mass is fixed) is
$$ T = \tfrac{1}{2} \boldsymbol{\omega}^\intercal \boldsymbol{L} = \tfrac{1}{2} \boldsymbol{\omega}^\intercal \mathbf{I} \,\boldsymbol{\omega} $$
where $\boldsymbol{\omega}$ is the rotational velocity vector, $\boldsymbol{L}$ the angular velocity vector, and $\mathbf{I}$ the mass moment of inertia tensor. All quantities are summed about the center of mass.
Now if you know the orientation (rotational transformation) of the body $\mathrm{R}(t)$ in terms of body-to-inertial matrix, and the body-fixed mass moment of inertia matrix $\mathbf{I}_{\rm body}$ then the above is
$$ T = \tfrac{1}{2} \boldsymbol{\omega}^\intercal \left( \mathrm{R} \mathbf{I}_{\rm body} \mathrm{R}^\intercal \right) \boldsymbol{\omega} $$
Now if you know that axis of rotation in the inertial frame, then $\boldsymbol{\omega} = \omega \,\boldsymbol{\hat{n}}$
Kinetic energy is thus calculated with
$$ T = \tfrac{1}{2} \omega^2 \left( \boldsymbol{\hat{n}}^\intercal \mathrm{R} \mathbf{I}_{\rm body} \mathrm{R}^\intercal \boldsymbol{\hat{n}} \right) $$
Here the part in the parenthesis changes with time, as the orientation matrix changes. You can use Euler-Angles to encode the rotation, or use quaternions. Either way you need to calculated $\mathrm{R}$ first before doing any calculations on a rigid body.
If you know that axis of rotation in the body frame, then $\boldsymbol{\omega} = \omega \, \mathrm{R}\, \boldsymbol{\hat{n}}_{\rm body}$
Kinetic energy is thus calculated with
$$ T = \tfrac{1}{2} \omega^2 \left( \boldsymbol{\hat{n}}_{\rm body}^\intercal \mathbf{I}_{\rm body} \boldsymbol{\hat{n}}_{\rm body} \right) $$ Note that if the axis of rotation is fixed to the body, then the part in the parenthesis is a constant.
I think you are asking about the case where the axis of rotation is not fixed to the body, and KE is expressed in body-fixed coordinates (like above). In this case $ \boldsymbol{\hat{n}}_{\rm body} = \mathrm{R}^\intercal \boldsymbol{\hat{n}}$
Specifically you are using $\boldsymbol{\hat{n}} = \pmatrix{0 \\ 0 \\ 1}$ and the spherical rotation matrix
$$ \mathrm{R} = \mathrm{rot}_y(-\theta) \mathrm{rot}_z(-\phi) = \begin{pmatrix}\cos\phi\cos\theta & \sin\phi\cos\theta & -\sin\theta\\ -\sin\phi & \cos\phi & 0\\ \cos\phi\sin\theta & \sin\phi\sin\theta & \cos\theta \end{pmatrix} $$
Note that ${\rm R}^\intercal {\rm R} = 1$ and that when the angles all are zero ${\rm R} = 1$.
Such that $$\boldsymbol{\hat{n}}_{\rm body} = {\rm R}^\intercal \boldsymbol{\hat{n}} = \begin{pmatrix}\cos\phi\sin\theta\\ \sin\phi\sin\theta\\ \cos\theta \end{pmatrix}$$
Take the equation above and plug in the above to get
$$ T = \tfrac{1}{2} \omega^2 \left( I_1 \cos^2 \phi \sin^2 \theta + I_2 \sin^2 \phi \sin^2\theta + I_3 \cos^2 \theta \right) $$
where $I_1$, $I_2$, and $I_3$ are the three diagonal terms of the body-fixed MMOI tensor (the principal MMOI values).