There are several challenges associated with inverse problems in which the unknown parameters can be modeled as piecewise constant functions. We model the unknown parameter using multiple level sets to represent the piecewise constant function. Adopting a Bayesian approach, we impose prior distributions on both the level set functions that determine the piecewise constant regions as well as the parameters that determine their magnitudes. We develop a Gauss-Newton approach with a backtracking line search to efficiently compute the maximum a priori (MAP) estimate as a solution to the inverse problem. We use the Gauss-Newton Laplace approximation to construct a Gaussian approximation of the posterior distribution and use preconditioned Krylov subspace methods to sample from the resulting approximation. To visualize the uncertainty associated with the parameter reconstructions we compute the approximate posterior variance using a new diagonal estimator (LanczosMC) that combines Lanczos and Monte Carlo methods. We illustrate that the incorporation of Lanczos iterations act as a variance reduction technique with a computational example. We will also demonstrate the benefits of our Bayesian level set approach and solvers on synthetic test problems as well as an application to X-ray imaging with real data.