The first extension of the virtual fields method to the reconstruction of heterogeneous stiffness properties from 3D bulk full-field displacement data is presented in this paper. Data are provided by Magnetic Resonance Imaging (MRI). Two main issues are addressed: 1. the identification of the stiffness ratio between two different media in a heterogeneous solid; 2. the reconstruction of stiffness heterogeneities buried in a heterogeneous solid. The approach is based on a finite element discretization of the equilibrium equations. It is tested on experimental full-field data obtained on a phantom with the stimulated echo MRI technique. The phantom is made of a stiff spherical inclusion buried within a lower modulus material. Preliminary independent tests showed that the material of the inclusion was four times stiffer than the surrounding material. This ratio value is correctly identified by our approach directly on the phantom with the MRI data. Moreover, the modulus distribution is promisingly reconstructed across the whole investigated volume. However, the resulting modulus distribution is highly variable. This is explained by the fact that the approach relies on a second order differentiation of the data, which tends to amplify noise. Noise is significantly reduced by using appropriate filtering algorithms.