Euler angles (and Tait-Bryan angles) are evil. They can be defined in a number of ways, and even when people remember to specify the axes and the order they are rotated, it all gets too hairy too fast. Instead, I show how to easily deal with this using only elementary vector algebra.
Let's say you have an ellipsoid centered at $\vec{c}$ with semi-principal axes $\vec{r}_1 $, $\vec{r}_2$, and $\vec{r}_3$:
$$\vec{r}_1 = \left [ \begin{array}{c} r_{11} \\ r_{21} \\ r_{31} \end{array} \right ], \qquad \vec{r}_2 = \left [ \begin{array}{c} r_{12} \\ r_{22} \\ r_{32} \end{array} \right ], \qquad \vec{r}_3 = \left [ \begin{array}{c} r_{13} \\ r_{23} \\ r_{33} \end{array} \right ]$$
The surface of the ellipsoid passes through points $\vec{c} \pm \vec{r}_1$, $\vec{c} \pm \vec{r}_2$, and $\vec{c} \pm \vec{r}_3$.
The semi-principal axes are all perpendicular: $\vec{r}_1 \cdot \vec{r}_2 = 0$, $\vec{r}_1 \cdot \vec{r}_3 = 0$, and $\vec{r}_2 \cdot \vec{r}_3 = 0$.
The axis-aligned bounding box for the above ellipsoid is defined by
$$\begin{array}{c}
x = c_x \pm \sqrt{r_{11}^2 + r_{12}^2 + r_{13}^2} \\
y = c_y \pm \sqrt{r_{21}^2 + r_{22}^2 + r_{23}^2} \\
z = c_z \pm \sqrt{r_{31}^2 + r_{32}^2 + r_{33}^2} \end{array}$$
Note that this applies also to the case where the homogenous transformation matrix $\mathbf{M}$ that transforms the unit sphere to the ellipsoid is
$$\mathbf{M} = \left [ \begin{array}{cccc} r_ {11} & r_{12} & r_{13} & c_x \\
r_{21} & r_{22} & r_{23} & c_y \\ r_{31} & r_{32} & r_{33} & c_z \\ 0 & 0 & 0 & 1 \end{array} \right ]$$
We can describe the transformation of the unit sphere to the above ellipsoid using
$$\vec{p}(u, v, w) = \vec{p}(\vec{t}) = \vec{c} + u \vec{r}_1 + v \vec{r}_2 + w \vec{r}_3 = \vec{c} + \mathbf{R} \vec{t} \tag{1}\label{1}$$
where $\vec{t} = ( u , v , w )$ fulfills
$$\left\lVert \vec{t} \right\rVert = 1 \qquad \iff \qquad u^2 + v^2 + w^2 = 1$$
on the surface of the ellipsoid, and
$$\left\lVert \vec{t} \right\rVert \lt 1 \qquad \iff \qquad u^2 + v^2 + w^2 \lt 1$$
within the ellipsoid, and matrix $\mathbf{R}$ has the semi-principal axes as column vectors,
$$\mathbf{R} = \left [ \begin{array}{ccc} \vec{r}_1 & \vec{r}_2 & \vec{r}_3 \end{array} \right ] = \left [ \begin{array}{ccc} r_{11} & r_{12} & r_{13} \\ r_{21} & r_{22} & r_{23} \\ r_{31} & r_{32} & r_{33} \end{array} \right ]$$
Essentially, we scale the unit sphere by the semi-principal axes, rotate it to the new orientation, and translate the center.
The above operation is invertible, if the ellipsoid is not degenerate; that is, if $\left\lVert\vec{r}_1\right\rVert \gt 0$, $\left\lVert\vec{r}_2\right\rVert \gt 0$, and $\left\lVert\vec{r}_3\right\rVert \gt 0$.
Since rotation matrices are orthonormal, its inverse is its transpose. To counteract the scaling, we scale by the inverses of the semi-principal axes.
The scaling vectors are now row vectors of the transformation matrix, with lengths inverse to the corresponding semi-principal axis length; i.e.
$$\vec{s}_n = \left ( \frac{\vec{r}_n}{\left\lVert \vec{r}_n \right\rVert^2} \right)^{T}, \qquad n = 1, 2, 3$$
so that $$\left\lVert \vec{s}_n \right\rVert = \left\lVert \vec{r}_n \right\rVert^{-1}, \qquad n = 1, 2, 3$$
i.e.
$$\vec{s}_1 = \frac{\vec{r}_1}{\left\lVert\vec{r}_1\right\rVert^2} = \left [ \begin{array}{ccc} s_{11} & s_{12} & s_{13} \end{array} \right ]$$
$$\vec{s}_2 = \frac{\vec{r}_2}{\left\lVert\vec{r}_2\right\rVert^2} = \left [ \begin{array}{ccc} s_{21} & s_{22} & s_{23} \end{array} \right ]$$
$$\vec{s}_3 = \frac{\vec{r}_3}{\left\lVert\vec{r}_3\right\rVert^2} = \left [ \begin{array}{ccc} s_{31} & s_{32} & s_{33} \end{array} \right ]$$
$$s_{i j} = \frac{r_{j i}}{ r_{1 i}^2 + r_{2 i}^2 + r_{3 i}^2 }, \qquad i = 1, 2, 3; \; j = 1, 2, 3$$
$$\mathbf{S} = \left [ \begin{array}{ccc} s_{11} & s_{12} & s_{13} \\ s_{21} & s_{22} & s_{23} \\ s_{31} & s_{32} & s_{33} \end{array} \right ]$$
Now,
$$\vec{t}(\vec{p}) = \vec{s}_1 \cdot ( \vec{p} - \vec{c} ) + \vec{s}_2 \cdot ( \vec{p} - \vec{c} ) + \vec{s}_3 \cdot ( \vec{p} - \vec{c} ) = \mathbf{S}\left(\vec{p} - \vec{c}\right) \tag{2}\label{2}$$
Essentially, $\eqref{1}$ and $\eqref{2}$ are the inverse of each other.
If you insist on using Euler or Tait-Bryan angles or other similar angular convention, I recommend you use them to calculate the semi-principal axes $\vec{r}_1$, $\vec{r}_2$, and $\vec{r}_3$ only, and then apply the above vector solution.
