We discuss in detail a method to study transverse momentum dependent parton distribution functions (TMDs) using lattice QCD. To develop the formalism and to obtain first numerical results, we directly implement a bi-local quark-quark operator connected by a straight Wilson line, allowing us to study T-even, process-independent TMDs. Beyond results for x-integrated TMDs and quark densities, we present a study of correlations in x and transverse momentum. Our calculations are based on domain wall valence quark propagators by the LHP collaboration calculated on top of gauge configurations provided by MILC with 2+1 flavors of asqtad-improved staggered sea quarks.