The O vacancy (Ov) formation energy, EOv, is an important property of a metal-oxide, governing its performance in applications such as fuel cells or heterogeneous catalysis. These defects are routinely studied with density functional theory (DFT). However, it is well-recognized that standard DFT formulations (e.g., the generalized gradient approximation) are insufficient for modeling the Ov, requiring higher levels of theory. The embedded cluster method offers a promising approach to compute EOv accurately, giving access to all electronic structure methods. Central to this approach is the construction of quantum(-mechanically treated) clusters placed within suitable embedding environments. Unfortunately, current approaches to constructing the quantum clusters either require large system sizes, preventing application of high-level methods, or require significant manual input, preventing investigations of multiple systems simultaneously. In this work, we present a systematic and general quantum cluster design protocol that can determine small converged quantum clusters for studying the Ov in metal-oxides with accurate methods, such as local coupled cluster with single, double, and perturbative triple excitations. We apply this protocol to study the Ov in the bulk and surface planes of rutile TiO2 and rock salt MgO, producing the first accurate and well-converged determinations of EOv with this method. These reference values are used to benchmark exchange-correlation functionals in DFT, and we find that all the studied functionals underestimate EOv, with the average error decreasing along the rungs of Jacob's ladder. This protocol is automatable for high-throughput calculations and can be generalized to study other point defects or adsorbates.