A novel method for inferring the toroidal current distribution from measurements of magnetic field and flux is presented. The method uses a Bayesian approach that gives the full joint probability distribution over all possible current distributions consistent with the measurements and their uncertainties. The plasma current is modelled as a grid of toroidal current-carrying solid beams with rectangular cross sections and the tomographic inversion corresponds to extracting information, such as mean values, standard deviations and marginal distributions from the joint posterior distribution over beam currents, as well as functions of the currents: plasma boundary, internal flux surfaces, X-points, etc. The method's intended primary usage is as a foundation for a principled probabilistic integration of multiple diagnostics, an example of which will be demonstrated. However, with suitable informative prior distributions, the method can be used on its own, giving in addition to best estimates, 'error bars' on plasma boundary, internal flux surfaces, X-point, etc. A fast approximation, suitable for realtime reconstruction of boundary, flux surfaces and q-profiles, will also be presented.