Details are presented of an efficient formalism for calculating transmission and reflection matrices from first principles in layered materials. Within the framework of spin density functional theory and using tight-binding muffin-tin orbitals, scattering matrices are determined by matching the wave-functions at the boundaries between leads which support well-defined scattering states and the scattering region. The calculation scales linearly with the number of principal layers N in the scattering region and as the cube of the number of atoms H in the lateral supercell. For metallic systems for which the required Brillouin zone sampling decreases as H increases, the final scaling goes as H^2*N. In practice, the efficient basis set allows scattering regions for which H^{2}*N ~ 10^6 to be handled. The method is illustrated for Co/Cu multilayers and single interfaces using large lateral supercells (up to 20x20) to model interface disorder. Because the scattering states are explicitly found, ``channel decomposition of the interface scattering for clean and disordered interfaces can be performed.