There are several challenges associated with inverse problems in which we\nseek to reconstruct a piecewise constant field, and which we model using\nmultiple level sets. Adopting a Bayesian viewpoint, we impose prior\ndistributions on both the level set functions that determine the piecewise\nconstant regions as well as the parameters that determine their magnitudes. We\ndevelop a Gauss-Newton approach with a backtracking line search to efficiently\ncompute the maximum a priori (MAP) estimate as a solution to the inverse\nproblem. We use the Gauss-Newton Laplace approximation to construct a Gaussian\napproximation of the posterior distribution and use preconditioned Krylov\nsubspace methods to sample from the resulting approximation. To visualize the\nuncertainty associated with the parameter reconstructions we compute the\napproximate posterior variance using a matrix-free Monte Carlo diagonal\nestimator, which we develop in this paper. We will demonstrate the benefits of\nour approach and solvers on synthetic test problems (photoacoustic and\nhydraulic tomography, respectively a linear and nonlinear inverse problem) as\nwell as an application to X-ray imaging with real data.\n