Euler angles between two given coordinate systems

Hi. I'm trying to do a coordinate system transformation from a global (0,0,0) to a local element coordinate system. The local system's position and unit vectors are known. My problem is that I need three Euler angles for the transformation matrix. I can only calculate the current angles between the global and local axis. But how do i get Euler angles with respect to the right order (e.g. Euler ZYZ rotation)? Is it possible to compute the Euler angles from that information?

The global coord. system is green, the local is blue. The dashed blue system is the local system shifted to the point (0,0,0).

For the calculation of the transformation matrix, I used _localToGlobal3d.m_ from geom3d:
https://www.mathworks.com/matlabcentral/fileexchange/24484-geom3d  
The transform is defined by:
  %   - CENTER: the position of the local origin into the World coordinate
  %       system
  %   - THETA: colatitude, defined as the angle with the Oz axis (between 0
  %       and 180 degrees), positive in the direction of the of Oy axis.
  %   - PHI: azimut, defined as the angle of the normal with the Ox axis,
  %       between 0 and 360 degrees
  %   - PSI: intrinsic rotation, corresponding to the rotation of the object
  %       around the direction vector, between 0 and 360 degrees
  %
  %   The resulting transform is obtained by applying (in that order):
  %   - Rotation by PSI   around he Z-axis
  %   - Rotation by THETA around the Y-axis
  %   - Rotation by PHI   around the Z-axis
  %   - Translation by vector CENTER
  %   This corresponds to Euler ZYZ rotation, using angles PHI, THETA and
  %   PSI.

Thank you for your help!

[James Tursa]

Do you have the unit axes for both coordinate systems? Or the angles between each axis of one system vs each axis of the other system?

[¿ ?]

Yes, I think so. I can calculate the current angles between the local system's axis (e1,e2,e3) and the base axis:

angleX = atan2d(norm(cross(e1,[1 0 0])), dot(e1,[1 0 0]))
angleY = atan2d(norm(cross(e2,[0 1 0])), dot(e2,[0 1 0]))
angleZ = atan2d(norm(cross(e3,[0 0 1])), dot(e3,[0 0 1]))