We introduce a novel numerical approach for the simulation of soft particles interacting via frictional contacts. This approach is based on an implicit formulation of the Material Point Method, allowing for large particle deformations, combined with the Contact Dynamics method for the treatment of unilateral frictional contacts between particles. This approach is both precise due to the treatment of contacts with no regularization and artificial damping parameters, and robust due to implicit time integration of both bulk degrees of freedom and relative contact velocities at the nodes representing the contact points. By construction, our algorithm is capable of handling arbitrary particle shapes and deformations. We illustrate this approach by two simple 2D examples: a Hertz contact and a rolling particle on an inclined plane. We also investigate the compaction of a packing of circular particles up to a solid fraction well above the jamming limit of hard particles. We find that, for the same level of deformation, the solid fraction in a packing of frictional particles is above that of a packing of frictionless particles as a result of larger particle shape change.