The dipole subtraction method for calculating next-to-leading order corrections in QCD was originally only formulated for massless partons. In this paper we extend its definition to include massive partons, namely quarks, squarks and gluinos. We pay particular attention to the quasi-collinear region, which gives rise to terms that are enhanced by logarithms of the parton masses, M. By ensuring that our subtraction cross section matches the exact real cross section in all quasi-collinear regions we achieve uniform convergence both for hard scales Q ∼ M and Q ≫ M. Moreover, taking the masses to zero, we exactly reproduce the previously-calculated massless results. We give all the analytical formulae necessary to construct a numerical program to evaluate the next-to-leading order QCD corrections to arbitrary observables in an arbitrary process.