Let me explain the intuition briefly. The authors of the paper you've cited assume that the signal $x(t)$ can be written in the form
\begin{align*}
x(t) &= \sum_{k=1}^K a_k(t) \exp(2\pi\mathrm{i} \phi_k(t)),
\end{align*}
where the $a_k$ denote instantaneous amplitudes, the $\phi_k$ denote instantaneous phases and the $\phi'_k$ denote instantaneous frequencies. The continuous wavelet transform enables the visualization of the $\phi'_k$ in the time-scale plane, but the visualization is necessarily blurry due to the Fourier uncertainty principle. The idea is that the peaks of the continuous wavelet transform can give you an idea of what the $\phi'_k$ values are, but the energy is spread out in the time-scale plane, so the exact values of $\phi'_k$ may not be easy to obtain. The goal of the synchrosqueezing transform (SST) is to partially undo the blurring by a frequency reassignment, which can be explained as follows.
First of all, let's note that each scale $s$ corresponds to a natural frequency $\xi$ satisfying the relation $s = c/\xi$, where $c$ is the center frequency of the mother wavelet $\psi$. Now, suppose the time $u$ is fixed. If $\xi = c/s$ is close but not exactly equal to an instantaneous frequency $\phi'_k(u)$, then the coefficient $W(u,s)$ will have some nonzero energy (i.e., $|W(u,s)|^2 > 0$). The idea of synchrosqueezing is to move all this energy away from the frequency $\xi$, and reassign its frequency location closer to the instantaneous frequency $\phi'_k(u)$. Then, we can arrive at a time-frequency representation where the energy is more closely concentrated around the instantaneous frequency curves.
So, for each scale $s$, we compute the frequency reassignment $\omega(u,s)$ of the wavelet transform coefficient $W(u,s)$ by the formula you cited above. As Daubechies, Wu, and Lu discovered in the paper you've linked, it turns out that $\omega(u,s)$ is often a good approximation to the instantaneous frequency curve $\phi'_k(u)$ for values of $(u,s)$ where $\xi = c/s$ is sufficiently close to $\phi'_k(u)$. (In fact, when $\phi'_k(u)$ is a constant function, $\omega(u,s)$ exactly equals $\phi'_k(u)$ for values of $s$ where $\xi = c/s$ is sufficiently close to $\phi'_k(u)$!)
Then, the computation of the SST in the continuous setting is as follows. First, fix a time of interest $u$. Next, compute the frequency reassignment $\omega(u,s)$ for all scale values $s$. Then, for each frequency of interest $\eta$, we compute the SST $T(u,\eta)$ by adding up all values $W(u,s)$ where the reassigned frequency $\omega(u,s)$ equals $\eta$. This means
\begin{align*}
T(u,\eta) &= \int_\mathbb{R} W(u,s) \delta(\eta - \omega(u,s)) s^{-3/2} ds,
\end{align*}
where $\delta$ is the Dirac delta. (To ensure the convergence result for SST, the authors use a formulation for $T$ which includes an approximation to $\delta$ rather than $\delta$ itself, but in computational practice I've never found this approximation to be necessary.) The multiplication by $s^{-3/2}$ is necessary for reconstruction purposes (which you can read more about in the paper).
In practice, one is limited to the discrete setting where we have only finitely many possible frequency bins $\eta_\ell$. In this case, for a fixed frequency $\eta_{\ell_0}$, one finds all the values of $\omega(u, s)$ which are closer to $\eta_{\ell_0}$ than to any other frequency bin $\eta_\ell$. We can more explicitly formulate the SST in the discrete setting as follows.
For simplicity's sake let's assume the possible frequency values $\eta_\ell$ are uniformly spaced by a distance $\Delta \omega$. So, in the discrete setting, the formula above becomes
\begin{align*}
T(u, \eta_\ell) &= \sum_{s: |\omega(u,s) - \eta_\ell| < \Delta\omega/2} W(u,s) s^{-3/2},
\end{align*}
where I haven't bothered to give an explicit notation that shows that the values $u$ and $s$ are discrete.
By doing this reassignment, you end up with a representation $T$ which is sparser than $W$, and hopefully very sharply concentrated about the curves $\phi'_k$. (In practice, for visualization purposes, it's better to reassign not the original coefficient $W(u,s)$ but the magnitude-squared $|W(u,s)|^2$, or even just the constant number $1$. This is because the coefficients $W(u,s)$ are complex-valued, and summing these coefficients may not actually leave you with a very large energy. However, if you additionally want to use the reconstruction formula given in the paper to compute the $x_k$ from the SST, then you need to sum the $W(u,s)$ in order for that formula to work.)
By the way, I am leaving out the fact that we generally throw out coefficients $W(u,s)$ that do not fall above a certain threshold. This is because the computation of $\omega(u,s)$ is not necessarily accurate for small $W(u,s)$.
