We present a 3D joint inversion method to estimate two physical parameters, density and magnetization of subsurface materials using field intensity measurements. In the method, we introduce the fuzzy c-means (FCM) clustering technique to relate gravity with magnetic data. In the approach, the subsurface structure is discretized to a set of rectangular prisms. For estimating the density and magnetization of each prism, we minimize the quadratic norm of the residuals between the observed data and the calculated responses from inversion model. Two regularization terms, i.e. the roughness and the similarity of the two physical parameters, are introduced in our joint inversion to control the degree of model roughness and similarity. We determine their regularization parameters using the L-curve criterion. We apply our method to a numerical model which represents submarine massive sulphides (SMS). The joint inversion results, which have the advantages of both gravity and magnetic inversion, show better accuracy and resolution than the individual ones.