Most modern texts spend some time deriving the LSZ reduction formula that connects S matrix elements to time ordered field correlation functions. It seems essential, and really helps clear up what you are calculating. Yet some earlier texts and even some modern texts (e.g. "Student Friendly Quantum Field Theory" by R. Klauber) seem to skip right past this, working everything out in the "interaction" picture. It seems there must be something going wrong with this latter procedure, but I am not quite able to put it together.
The latter route works as long as perturbative expansion works, i.e. states spectrum does not change. This can be checked by hand using said LSZ formula: you will end up calculating correlators $<\Phi(x_1)...>$ of free theory, which can be readily written in terms of creation-annihilation operators, nicely reproducing canonical perturbation series.
Textbooks that adopt canonical approach usually make handwavy arguments about this. An honest approach would probably be using one-particle eigenstates of full Hamiltonian, and then showing that in every order of perturbation theory this gives the same $S$-matrix as free eigenstates.
Obviously, if the spectrum of the interacting theory differs from its free counterpart, as in QCD, canonical expression for $S$-matrix makes no sense. LSZ formula, on the other hand, does.