This paper considers a numerical method for reconstructing piecewise constant Robin coefficients from boundary measurements. A novel variational formulation of the inverse problem based on the celebrated Modica-Mortola functional is proposed, and its Gamma-convergence to the formulation with total variation regularization is shown. The differentiability of the functional is investigated, based on which the concave convex procedure is justified for its numerical solution. The finite element method is utilized for discretizing the functional, and the convergence of the finite element approximation is established. Three sets of numerical experiments are presented to illustrate its distinct features.