We describe a novel variational segmentation algorithm designed to split an image in two regions based on their intensity distributions. A functional is proposed to integrate the probability density functions of both regions within the optimization process. The method simultaneously performs segmentation and non-parametric density estimation. It does not make any assumption on the underlying distributions, hence it is flexible and can be applied to a wide range of applications. Although a boundary evolution scheme may be used to minimize the functional, we choose to consider an alternative formulation with membership functions. The latter has the advantage of being convex in each of the variables, so that the minimization is faster and less sensitive to initial conditions. Finally, to improve the accuracy and the robustness to low-frequency artifacts, we present an extension for the more general case of local probability densities, allowed to vary in space. The approach readily extends to multi-channel images and 3D volumes, and we show several results on synthetic and photographic images, as well as on 3D medical data.