Abstract: We extend the four-dimensional unsubtraction method, which is based on the loop-tree duality (LTD), to deal with processes involving heavy particles. The method allows to perform the summation over degenerate IR configurations directly at integrand level in such a way that NLO corrections can be implemented directly in four space-time dimensions. We define a general momentum mapping between the real and virtual kinematics that accounts properly for the quasi-collinear configurations, and leads to an smooth massless limit. We illustrate the method first with a scalar toy example, and then analyse the case of the decay of a scalar or vector boson into a pair of massive quarks. The results presented in this paper are suitable for the application of the method to any multipartonic process.