Vibrational spectroscopy and normal modes are based on the Born-Oppenheimer approximation and the harmonic approximation.
The BO approximation yields the adiabatic electronic states and the corresponding adiabatic energy eigenvalues as functions of the nuclear coordinates. The adiabatic energy eigenvalues of the electronic states are also known as potential energy surfaces $E_\alpha(R)$. Associated with this electronic adiabatic state is a vibronic Hamiltonian for the nuclear wave function,
$$
H_\alpha(R) = \frac{1}{2}\nabla^T_RM^{-1}\nabla_R + E_\alpha(R)
$$
From now on I will fix the electronic state to be the groundstate and drop the index $\alpha=0$, that defines the electronic state. In pure vibrational spectroscopy the electronic state remains the same and is thus of no further interest.
If we knew the exact potential energy surface $E(R)$, then we could try to solve the problem numerically without the harmonic approximation. But doing so is very expensive for a large amount of degrees of freedom.
To avoid this, we can make the harmonic approximation and expand the potential energy surface to second order around a equilibrium structure, i.e. the nuclear coordinates where the potential energy surface has a minimum. Let this structure be $R_{EQ}$,
$$
E(R)\approx E(R_{EQ}) + \sum^{3N}_{i=1}\frac{\partial E} {\partial R_i}|_{R_{EQ}} (R-R_{EQ})_i +\frac{1}{2}\sum^{3N,3N}_{i=1,j=1}(R-R_{EQ})_i \frac{\partial^2 E} {\partial^2 R_i R_j}|_{R_{EQ}} (R-R_{EQ})_j
$$
We can chose our energy axis such that $E(R_{EQ})=0$. We also are at a minimum of the potential energy surface, which means that the gradient vanishes $\frac{\partial E} {\partial R_i}|_{R_{EQ}}=0 \ \forall \ i$. All that remains in the harmonic approximation at the equilbrium geometry is the second order term with the Hessian evaluated at the equilibrium geometry $H_{ij}=\frac{\partial^2 E_\alpha} {\partial^2 R_i R_j}|_{R_{EQ}} $,
$$
E^{Harmonic}(R) = \frac{1}{2}\sum^{3N,3N}_{i=1,j=1}(R-R_{EQ})_i H_{ij} (R-R_{EQ})_j
$$
With this we can define the vibronic Hamiltonian in the harmonic approximation,
$$\begin{aligned}
H^{Harmonic}(R) &= \frac{1}{2}\nabla^T_RM^{-1}\nabla_R + E^{Harmonic}_\alpha(R)\\
H^{Harmonic}(R)&= \frac{1}{2}\nabla^T_RM^{-1}\nabla_R + \frac{1}{2}\sum^{3N,3N}_{i=1,j=1}(R-R_{EQ})_i H_{ij} (R-R_{EQ})_j
\end{aligned}$$
A short remark on what I have not done is also in order. I have not separated the nuclear coordinates into translation, rotation and vibration thus far. Six(five for linear molecules) of the normal coordinates that we will obtain are going to describe translation and rotation. They are typically easy to recognize by their associated angular frequencies that are almost zero.
Now back to the problem. I will now drop the superscript harmonic. Our harmonic vibronic Hamilton operator is not a sum of 1D Harmonic oscillator Hamilton operators, but we can achieve that by doing some coordinate transformations.
First we introduce mass-weighted displacements
$$
x_i=\sqrt m_i (R-R_{EQ})_i \rightarrow \frac{\partial}{\partial x_i} = \frac{1}{\sqrt m_i}\frac{\partial}{\partial R_i}\rightarrow \nabla_x = M^{-1/2}\nabla_R
$$
In these coordinates we obtain
$$
H(x) = \frac{1}{2}\nabla_x^T\nabla_x + \frac{1}{2}\sum^{3N,3N}_{i=1,j=1}x_i \frac{H_{ij}}{\sqrt{m_i m_j}} x_j
$$
We define now the mass weighted Hessian $F_{ij} = \frac{H_{ij}}{\sqrt{m_i m_j}}$
which leads us to
$$
H(x) = \frac{1}{2}\nabla_x^T\nabla_x + \frac{1}{2}\sum^{3N,3N}_{i=1,j=1}x_i F_{ij} x_j
$$
In this equation we still have coupling between different $x_{i/j}$ due to off-diagonal matrix elements in the mass weighted Hessian. But we can eliminate these by transformation to a coordinate system in which the matrix is diagonal. To find these coordinates we solve the eigenvalue problem
$$
FV = V\Lambda
$$
where the matrix $V$ contains the eigenvectors of $F$ as columns. $\Lambda$ is the corresponding diagonal eigenvalue matrix.
The eigenvector matrix allows us to define the transformation to normal coordinates $Q_i$,
$$
\sum^{3N}_{k=1} V_{ik}Q_k = x_i \quad \rightarrow VQ=x \quad \leftrightarrow Q = V^{-1}x
$$
The eigenvectors of the mass weighted Hessematrix are orthogonal and if normalized form a orthonormal set. This means that $VV^T=1=VV^{-1}$. Which means that $F$ can be diagonalized via
$$
\Lambda = V^TFV
$$
Plugging in the transformation into the Hamiltonian yields
$$\begin{aligned}
H(Q) &= \nabla_Q^T({V^{-1}})^TV^{-1}\nabla_Q + \frac{1}{2} Q^TV^TFVQ \\
H(Q) &= \nabla_Q^T\nabla_Q +\frac{1}{2} Q^T \Lambda Q
\end{aligned}$$
Expanding this equation shows that we now have a simple sum of harmonic 1D Hamiltonians since the Eigenvalues matrix is diagonal $\Lambda_{ij}=\delta_{ij}\Lambda_i$ and no longer couples different coordinates.The last step is recognizing that each 1D equation is the same as a harmonic oscillator equation with $\omega_i^2=\Lambda_i.$
$$\begin{aligned}
H(Q) &= \sum^{3N}_{i=1}\frac{1}{2}\frac{\partial^2 }{\partial Q^2_i} + \sum^{3N}_{i=1}\frac{1}{2}Q_i^2\Lambda_i\\
&=\sum_{i=1}^{3N}\frac{1}{2}\frac{\partial^2}{\partial Q^2_i} + \frac{1}{2}Q_i^2\omega_i^2\\
H(Q) &= \sum_{i=1}^{3N} H^{1D}(Q_i)
\end{aligned}$$
Meaning that our large coupled problem finally has taken the shape of a sum of 3N one dimensional harmonic oscillators, which we know how to solve analytically.
At this point I'd like to repeat that we have not separated translation and rotation. Due to this, the first six frequencies should be zero (or almost zero). The six normal modes corresponding to these zero eigenvalues are not oscillations and have to be omitted when calculating vibrational transitions.
It is possible to remove them before transforming to normal coordinates. How this is done is described in the book of Wilson, Decius and Cross, Molecular vibrations but discussing these points goes beyond the scope of my answer.
Let me now address how all of this connects to the figure that you show in your question. In particular how we know how the atoms are displaced relative to each other for a given normal mode.
By inverting the transformations, we see
$$
\Delta R = M^{-1/2}VQ
$$
with $\Delta R = R-R_{EQ}$.
This means that the columns of the matrix
$$\begin{aligned}
K &= M^{-1/2}V
\end{aligned}$$
describe the relative cartesian displacement amplitudes of the normal modes, easily seen by plugging in a unit vector for Q.
If we had calculated benzene then one of the columns of $K$ would contain the relative displacements as shown in your figure. Some quantum chemistry programs also print these column vectors after some sort of normalization, when they output Cartesian normal modes and use them for visualization.
Now let me say a bit about reduced masses and force constants. First of all reduced mass and force constants are no observables. They are helpful quantities that we can define to make the analogy to a 1D system with a mass attached via spring to a wall. Depending on the particular conventions that we use, we will obtain different values. A common convention is to define a reduced mass based on the Cartesian displacement vector norm. The normal mode vectors as I have defined above are assumed to be normalized to one with regard to mass weighted coordinates. A mass weighted normal coordinate vector of norm one corresponds to a displacement of atoms such that all displacements squared, multiplied with their atomic mass then summed and the root taken yield 1. Take as an example one mode with many heavy atoms contributing to the motion and a mode where mostly light atoms contribute to motion. If both modes have now a normal coordinate value of 1, you would observe spatially when looking at the atom displacements, that the mode with the heavy atoms is only a little displaced while the mode with light atoms shows larger displacement.
We can however introduce a new coordinate that is constructed in a way that a value of 1 means always the same net spatial displacement of the atoms contributing to the mode irregardless of the mass of the contributing atoms.
Let us look at the norm of the Cartesian displacement vector associated with a mass weighted normal mode vector. A unit normal mode vector is of the form
$$
(\hat Q_\alpha)_j = \delta_{j\alpha}
$$
The transformation of this vector to Cartesian displacement coordinates yields
$$
R_\alpha = K\hat Q_\alpha
$$
The Cartesian norm squared that describes the actual displacement of the atoms in space is then,
$$
R_\alpha \cdot R_\alpha = K \hat Q_\alpha \cdot K\hat Q_\alpha
= \sum_{nml} K_{nm} (\hat Q_\alpha)_m K_{nl} (\hat Q_\alpha)_l \\
= \sum_{nml} K_{nm} \delta_{\alpha m} K_{nl} \delta_{\alpha l} \\
= \sum_{nml} K_{n\alpha} K_{n\alpha} \\
|R_\alpha|^2 = \sum_{nml} (K_{n\alpha})^2 \\
$$
We define now the reduced mass for mode $\alpha$ as
$$
\mu_\alpha = \frac{1}{|R_\alpha|^2}
$$
and the new set of coordinates
$$
y_\alpha = \frac{1}{\sqrt \mu_\alpha}Q_\alpha
$$
In this coordinates we obtain the Hamiltonians
$$
H(y_\alpha) = \sum_\alpha \frac{1}{2\mu_\alpha}\frac{\partial^2 }{\partial y_\alpha^2} + \frac{1}{2}\mu_\alpha \omega^2_\alpha y_\alpha^2
$$
and can define the force constants for each mode
$$
k_\alpha = \mu_\alpha \omega^2_\alpha.
$$
I am afraid you need a computer to do all of this or at least a very long breath and patience. I think for some systems the diagonalization can be done manually if you use symmetry, breaking down the matrices that have to be diagonalized to managable size, but generally you are going to let a computer do the linear algebra parts. I am pretty sure that is also the reason why almost all texts only discuss the trivial linear case where everything can be done by hand, although that is problematic since the simple case doesn't give you much intuition for the more general case.
All of this is also explained in the white paper https://gaussian.com/wp-content/uploads/dl/vib.pdf
I would not recommend the book from Wilson, Delcius and Cross as first reference. The book is good if you already have some understanding of the matter but as first exposition it might be difficult. It is however a good source if your are interested in internal coordinates and a treatment based on internal coordinates and the problem of separating vibrational and rotational degrees of freedom.
