How to calculate instantaneous rotation axis of rigid body from point velocity vectors? When we have the angular velocity vector and point positions, calculation velocity vectors for point is not a problem. We use formula
$$ \vec v = \vec \Omega \times \vec r $$
When I started writing algorithms for rigid body rotation I encountered the problem, how to calculate angular velocity from point velocity vectors using formulas, without the use of geometric relationships.
Flat example
We have two points of rigid body $m_1(x_1,0,0) ; m_2 (0,y_2,0) $ and velocity vectors for these points $ V_1(0,0,v_z1) ; V_2 (0,0,v_z2)$. How to calculate the instantaneous angular velocity for this example?
A three-dimensional example
We have three points of rigid body $m_1(x_1,0,0) ; m_2 (0,y_2,0); m_3 (0,0,z_3) $ and angular velocity vectors $\Omega (\omega_x, \omega_y, \omega_z) $. So point velocity vectors are
$ V_1(0,\omega_z x_1,-\omega_y x_1) ; V_2 (v_x2,0,v_z2) ; V_3 (v_x3,v_y3,0)$.
How to calculate the instantaneous angular velocity for this example using formulas not geometry?
 A: Since you have been working with rigid body rotations you are probably familiar with some of the unpleasant properties of those rotations, such as the fact that they are not commutative (the order in which you do the rotations matters).
Fortunately, angular velocity is not like that. It adds like an ordinary vector.
(See here or here.)
Note, however, that for vector addition of angular velocities to be physically meaningful, the angular velocities must be physically composed in some way.
In one of the examples in the link there is a platform rotating around one axis with a pendulum swinging around another axis that rotates with the platform.
These angular velocities add as vectors.
You don't need a physical linkage like that; it is enough to identify coordinate frames rotating relative to the inertial frame or to each other.
To reconstruct an angular velocity from components, however, the components must actually be components of the angular velocity relative to a single set of vectors, properly computed for that geometry.
We face the same problem with linear velocity:
if you want to decompose a linear velocity in a plane into one component parallel to the $x$ axis and one component parallel to the line $y=x,$
the calculation that gives the correct $x$ component for the usual orthogonal decomposition will generally give wrong answers for this non-orthogonal decomposition. It can be tricky to get this right, and it is no easier for decomposition of angular velocity than for decomposition of linear velocity.
But since you are working with orthogonal vectors, the decomposition is relatively straightforward.
For your three-dimensional example, first we can establish that the linear velocity of a point at $(0,0,0)$ on this rigid object would be zero,
so any rotation occurs around an axis through that point.
No angular velocity applied to either the $y$ or $z$ axes can give the linear velocity of the point $(0,y_2,0)$ a non-zero $z$ component,
so any $z$ component of that point must be due entirely to the component of angular velocity around the $x$ axis.
Likewise the $x$ component of the linear velocity of $(0,y_2,0)$ must be due entirely to the component of angular velocity around the $z$ axis,
and the $z$ component of the linear velocity of $(x_1,0,0)$ must be due entirely to the component of angular velocity around the $y$ axis.
That is, 
\begin{align}
v_{z2} &= y_2 \omega_x,\\
v_{x2} &= -y_2 \omega_z,\\
v_{z1} &= -x_1 \omega_y,\\
\end{align}
and the vectors $(\omega_x, 0, 0),$ $(0, \omega_y, 0),$
and $(0, 0, \omega_z)$ are components of the angular velocity of the body.
The total angular velocity is $(\omega_x, \omega_y, \omega_z),$
that is, an angular velocity of $\sqrt{\omega_x^2 + \omega_y^2 + \omega_z^2}$
around the axis parallel to the vector $(\omega_x, \omega_y, \omega_z).$
A: It is not true that you never can commutative operation vector product. For the equation we are interested 
$$\vec v = \vec \Omega \times \vec r   \tag 1$$
this can be easily done for a single point rotation (where vectors $\vec v, \vec \Omega, \vec r$ are perpendicular) and the value of the position vector is one $| \vec r|=1$. Then
$$\vec \Omega = \vec r \times \vec v   \tag 2$$
However, if the position vector is different from one $| \vec r|\neq1$, we must upgrades the pattern
$$\vec \Omega = \frac {\vec r \times \vec v}{| \vec r|^2}   \tag 3$$
When I understood what is temporary instantaneous angular velocity then I found a way to commutative the pattern (1) based on formula (3).
I will now use a controversial record a inverse vector. I have already asked about the correctness of such a record here Can I use a fraction of vector to reverse cross product? and because so far no one has proved that such a record is wrong I will still use it.
I write down as a fragment of formula (3) as follows
$$ \frac {\vec r} {|\vec r|^2}=( \frac {r_x} {|\vec r|^2},\frac {r_y} {|\vec r|^2},\frac {r_z} {|\vec r|^2})= \vec {\frac 1  r}  \tag 4$$
And the current pattern (3) record is as follows
$$\vec \Omega = \vec{\frac 1 r} \times \vec v   \tag 5$$
Now we can proceed to the solution of the first example that was in the question. To be correct, the conditions must be met, we must have two perpendicular point location vectors $ \vec r_1 (r_x,0,0), \vec r_y (0,r_y,0)$. And the angular velocity must be on surface of these vectors $\vec \Omega (\omega_x, \omega_y,0)$. We know from properties the rigid body that all points must move according to the formula (1) hence we learn what point velocity vectors are
$$\vec v_1= \vec \Omega \times \vec r_1 = (0,0,- \omega_y r_x) \tag 6 $$
$$\vec v_2= (0,0, \omega_x r_y) $$
The points rotate only at an angular velocity perpendicular to the position vector. So we count vectors of angular velocity of points from the formula (5)
$$\vec \Omega_1 = \vec{\frac 1 r_1} \times \vec v_1 = (0, \frac 1 r_x r_x \omega_y,0)    \tag 7$$
$$\vec \Omega_2 = \vec{\frac 1 r_2} \times \vec v_2 = ( \frac 1 r_y r_y \omega_x,0,0)    $$
So we see that
$$\vec \Omega = \vec \Omega_1 + \vec \Omega_2=  \vec{\frac 1 r_1} \times \vec v_1 +  \vec{\frac 1 r_2} \times \vec v_2 \tag 8 $$
This is the answer to the first example. Everything works in a flat arrangement. To calculate any instantaneous angular velocity in three-dimensional space need some more calculations. 
Answer to the second example
Important! The points must be perpendicular to each other.
Again, we are counting velocity vectors for points
$$\vec v_1= \vec \Omega \times \vec r_1 = (0,\omega_z r_x,- \omega_y r_x) \tag 9 $$
$$\vec v_2= (-\omega_z r_y,0, \omega_x r_y) $$
$$\vec v_3= (\omega_y r_z, -\omega_x r_z,0) $$
In a freely rotating rigid body in the reference system where the center of mass is stationary point velocity vectors are perpendicular to their position vectors.
We now calculate the angular velocities of the points from the formula (5)
$$\vec \Omega_1 = \vec{\frac 1 r_1} \times \vec v_1 = (0, \omega_y, \omega_z)    \tag {10}$$
$$\vec \Omega_2 = \vec{\frac 1 r_2} \times \vec v_2 = (\omega_x,0,\omega_z) $$
$$\vec \Omega_3 = \vec{\frac 1 r_3} \times \vec v_3 = (\omega_x,\omega_y,0) $$
When these vectors are added together, the angular velocity values are doubled, it means that you have to divide the result by two
$$\vec \Omega = \frac {\vec{\frac 1 r_1} \times \vec v_1 + \vec{\frac 1 r_2} \times \vec v_2 + \vec{\frac 1 r_3} \times \vec v_3  } 2 \tag {11} $$
But we don't have to use such a long formula, observing the results (10) we see the information on the coordinates of the instantaneous angular velocity can be drawn using the formula (5) any two points whose position vectors are perpendicular.
If someone claims that it is impossible to commutative operation the vector product, It seems like I've done something impossible:)
