A three-dimensional boundary integral method for deformable drops in viscous flows at low Reynolds numbers is presented. The method is based on a new nonsingular contour-integral representation of the single and double layers of the free-space Green's function. The contour integration overcomes the main difficulty with boundary-integral calculations: the singularities of the kernels. It also improves the accuracy of the calculations as well as the numerical stability. A new element of the presented method is also a higher-order interface approximation, which improves the accuracy of the interface-to-interface distance calculations and in this way makes simulations of polydispersed foam dynamics possible. Moreover, a multiple time-step integration scheme, which improves the numerical stability and thus the performance of the method, is introduced. To demonstrate the advantages of the method presented here, a number of challenging flow problems is considered: drop deformation and breakup at high viscosity ratios for zero and finite surface tension; drop-to-drop interaction in close approach, including film formation and its drainage; and formation of a foam drop and its deformation in simple shear flow, including all structural and dynamic elements of polydispersed foams.