"NBO diagrams" versus MO diagrams
A molecular orbital diagram is a schematic representation of how we interpret bonding in certain species. It is as much an accurate representation for a specific bonding situation as a Lewis structure for the the molecular geometry.
When we construct a molecular orbital diagram, usually we resort to shapes of hydrogen like orbitals. We most of the times do not even use mathematical descriptions, accurate energy levels, without changing the descriptive value much. Oftentimes we completely omit core orbitals and don't draw nodal planes/spheres for orbitals with $n=2$.
Such diagrams, however, provide qualitatively good descriptions of what otherwise is almost impossible to fathom. They even provide a solid starting point for quantum chemical calculations.
As such, the choice of atomic orbitals does not matter at all; a picture will never contain enough detail to work out the subtle differences of these choices.

By the use of computational chemistry methodologies we can (1.) produce a molecular wavefunction of a diatomic molecule as a linear combination of basis functions (assume they are centred on each nuclei). This would already satisfy the definition above and generate a suitable MO diagram.

Most basis sets indeed contain basis functions which are centred at the respective nucleus and they resemble atomic orbitals (they don't have to). Through a linear combination of such functions we can obtain suitable molecular orbitals within the respective approximations. The wave function is then obtained as a set of orthogonal molecular orbitals. This in principle is the essence of molecular orbital calculations, which is in turn essentially the same as calculations on the Hartree-Fock level of theory. We know that this level of theory doesn't contain nearly enough precision to accurately describe the most simple compounds because of the neglecting most of the correlation energy. Yet, it is impossible to represent the result of such a calculation as a diagram.

We could express molecular orbitals as a linear combination of optimised atomic orbitals of each atom, which would imply an extra calculation for each element in a molecule. A change of basis for the MOs, from the basis set to the atomic orbitals, does the rest.

Yes you can do that, but apart from adding significant computational cost, this would not really increase the accuracy of the calculation. First and foremost you need to account for a larger portion of the correlation energy, for example with density functional approximations, post-HF methods like perturbations theory MPn, coupled cluster, or configuration interaction. Instead of optimising the basis set for each atom in the molecule, it is much more effective adding more functions to the basis set until you reach the complete basis set limit. There are several types of basis set libraries in use now, which have been designed to effectively handle such calculations. I think it was John Pople who laid out the advantages of this system. And while his basis sets have significant shortcomings, they are still used today.
In any case, what you will not get is a better or more descriptive MO diagram.

The natural bonding orbitals (NBOs) method allows one to produce natural atomic orbitals as well, which could also be used as initial orbitals centred on the atoms.

Indeed natural bond orbital theory uses natural atomic orbitals (and other types), however, these are representations of the atomic orbital within the molecular environment.
The theory itself is a localisation technique of canonical orbitals that are already expressed in the chosen basis set. As a result we get a bonding picture which closely resembles a Lewis structure. While they are different from the actual basis set they still have the same functions and would still lead to the same solution.
It is best to treat MO diagrams as what they are: qualitative tools to describe bonding.