Standard 3D imaging systems process only a single return at each pixel from an assumed single opaque surface. However, there are situations when the laser return consists of multiple peaks due to the footprint of the beam impinging on a target with surfaces distributed in depth or with semi-transparent surfaces. If all these returns are processed, a more informative multi-layered 3D image is created. We propose a unified theory of pixel processing for Lidar data using a Bayesian approach that incorporates spatial constraints through a Markov Random Field with a Potts prior model. This allows us to model uncertainty about the underlying spatial process. To palliate some inherent deficiencies of this prior model, we also introduce two proposal distributions, one based on spatial mode jumping, the other on a spatial birth/death process. The different parameters of the several returns are estimated using reversible jump Markov chain Monte Carlo (RJMCMC) techniques in combination with an adaptive strategy of delayed rejection to improve the estimates of the parameters.