Thinking about it further, the accelerometers values are the components of the vector I was looking for... You only have to express the shortest arc between a reference vector and this vector as a quaternion.
For example, if accel
is a (normalized) vector containing the accelerometer values:
reference = Vector3(0, 0, 1)
axis = crossp(accel, reference)
angle = acos(dotp(accel, reference))
q = Quaternion::from_xyzr(axis, angle)
rotation_matrix = q.matrix()
I chose {0, 0, 1} as the reference vector, because it is the values of the accelerometers at "rest position" (Wiimote on a table, pointing towards you).
This gives the same kind of movements the X/Z accelerometers pitch/roll conversion does, but without gimbal lock at vertical positions.
The only problem is you don't get information about rotations made on the earth gravity axis... I guess this is what the MotionPlus is made for.