We have developed an inversion method for recovery of the geometry of an arbitrary number of geologic units using a regularized least-squares framework. The method addresses cases in which each geologic unit can be modeled using a constant physical property. Each geologic unit or group assigned the same physical property value is modeled using the signed distance to its interface with other units. We invert for this quantity and recover the location of interfaces between units using the level-set method. We formulate and solve the inverse problem in a least-squares sense by inverting for such signed distances. The sensitivity matrix to perturbations of the interfaces is obtained using the chain rule, and model mapping from the signed distance is used to recover the physical properties. Exploiting the flexibility of the framework that we develop allows any number of rock units to be considered. In addition, it allows the design and use of regularization incorporating prior information to encourage specific features in the inverted model. We apply this general inversion approach to gravity data favoring minimum adjustments of the interfaces between rock units to fit the data. The method is first tested using noisy synthetic data generated for a model compoed of six distinct units, and several scenarios are investigated. It is then applied to field data from the Yerrida Basin (Australia) where we investigate the geometry of a prospective greenstone belt. The synthetic example demonstrates the proof of concept of the proposed methodology, whereas the field application provides insights into, and potential reinterpretation of, the tectonic setting of the area.