AbstractSeveral recent methods for tomographic reconstruction of stress and strain fields from Bragg‐edge neutron strain images have been proposed in the literature. This paper presents an extension of a previously demonstrated approach based on Gaussian process regression that enforces equilibrium in the method. This extension incorporates knowledge of boundary conditions, primarily boundary tractions, into the reconstruction process. This is shown to increase the rate of convergence and is more tolerant of systematic errors that may be present in experimental measurements. An exact expression for a central calculation in this method is also provided which avoids the need for the approximation scheme that was previously used. Convergence of this method for simulated data is compared to existing approaches and a reconstruction from experimental data is provided. Validation of the results to conventional constant wavelength strain measurements and comparison to prior methods show a significant improvement.