Confusion about atomic/molecular orbital terminology
The short answer is that orbitals are defined to within a unitary transformation.



Canonical Hartree-Fock equations

The application of variational principle to the Slater determinant leads to the system of Hartree-Fock equations
$$
    \newcommand{\op}[1]{\hat{#1}}
    \newcommand{\mat}[1]{\mathbf{#1}}
    \op{F} \psi_{i}(1)
    =
    \sum\limits_{j=1}^{n} \lambda_{ij} \psi_{j}(1)
    \quad \text{for} \quad i = 1, 2, \dotsc, n \, ,
$$
which can be simplified by choosing a unitary transformation $\mat{U}$, defined by the way it transforms spin orbitals,
$$
    \psi_{k}' = \sum\limits_{j=1}^{n} u_{kj} \psi_{j} \quad \text{for} \quad k = 1, 2, \dotsc, n \, ,
$$
that makes the matrix of Lagrange multipliers $\lambda_{ij}$ diagonal and leads to the system of the so-called canonical Hartree-Fock equations
$$
    \op{F} \psi_{k}'(1) = \lambda_{kk}' \psi_{k}'(1) \, , \quad k = 1, 2, \dotsc, n \, ,
$$
or, in a much more familiar form using symbol $\varepsilon_k$ for $λ_{kk}'$ and dropping the primes
$$
    \op{F} \psi_{k}(1) = \varepsilon_{k} \psi_{k}(1) \, , \quad k = 1, 2, \dotsc, n \, .
$$
The canonical Hartree-Fock equations can be thought of the eigenvalue ones, for solving which there exist a number of well-developed algorithmic approaches, thus, they are so convenient to work with.

Slater determinant invariance to a linear transformation

It can be shown that Slater determinants composed of original and transformed spin orbitals are related as follows,
$$
    \Phi' = \det(\mat{U}) \Phi \, ,
$$
which by taking the absolute values and squaring both sides leads to
$$
    |\Phi'|^{2} = |\det(\mat{U})|^{2} |\Phi|^{2} \, .
$$
This means that the transformed Slater determinant $\Phi'$ describe exactly the same system as the original Slater determinant $\Phi$ in every observable way. Or, in other words, the Slater determinant is invariant to any linear transformation of spin orbitals (except for normalization). The only advantage of the unitary transformation from that perspective is that $|\det(\mat{U})|^{2} = 1$, so that there is no even need for renormalization.

Transformed canonical orbitals

Now, as it was mentioned earlier, it is convenient to do calculations with canonical orbitals, however, once the Hartree-Fock equations are solved, the resulting canonical spin orbitals can be transformed by any arbitrary transformation making no changes with respect to observable properties of the many electron system in question. And it advantageous to use unitary transformations since there is no need for renormalization of the Slater determinant.

So, other sets of orbitals can be constructed by transforming the canonical orbitals, i.e. by forming linear combinations of them. The way these linear combinations are formed (i.e. the choice of a unitary transformation) is dictated by what is the ultimate goal of the transformation. For instance, canonical orbitals are usually fairly delocalized over a molecule, rather than concentrated in limited spatial regions of it as anticipated from general chemistry models. Thus, of a particular interest for a chemists is a transformation which localizes the orbitals. The way localized orbitals are obtained depends on a particular choice of the unitary transformation and there are in fact few different approaches.

But just localizing the orbitals might not be the end of the story. For instance, one might want to not just localize them, but do so in a way that the maximum-occupancy character in localized 1-center and 2-center regions of the molecule. This procedure results in the most accurate possible natural Lewis structure corresponding to a given Slater determinant and is the basis of the natural bond orbital model.