The non-perturbative nature of QCD at hadronic scales implied the development of phenomenological approaches such as quark models or, more recently, computer-based calculations using Lattice QCD. However, the unique properties of heavy quarkonium systems allow an entire calculation in terms of non-relativistic perturbative QCD. In this work, the bottomonium spectrum, up to n=3, and the ground state charmonium states, are analyzed in the framework of Non-Relativistic Quantum Chromodynamics at N$^3$LO. For bottomonium, finite charm quark mass effects in the QCD potential and the $\overline{MS}$-pole mass relation are incorporated to the highest known order, $\mathcal{O}(\varepsilon^3)$ in the $\Upsilon$-scheme counting. The bottom quark pole mass is expressed in terms of the MSR mass, a low-scale short-distance mass which cancels the u=1/2 renormalon of the static potential. We study the $n_\ell=3$ and $n_\ell=4$ schemes, finding a negligible difference between the two if finite $m_c$ effects are smoothly incorporated in the MSR mass definition. We find that bottomonium n=3 states are not well behaved within perturbative NRQCD. Hence, fitting to the n=1,2 $b\bar b$ states we obtain $\overline{m}_b(\overline{m}_b) = 4.216\pm 0.039$ GeV. Similarly, from the lowest lying charmonium states we find $\overline{m}_c(\overline{m}_c)=1.273 \pm 0.054$ GeV.