A method for the reconstruction of a spectrum of mechanical properties from phase-contrast, MR displacement data is developed and implemented via the finite element method. To overcome the significant computational load associated with reconstructive imaging methods of this nature, the process is implemented on reasonably sized subdomains, or “subzones”, of the total region of interest. This reduction in computational size within the algorithm allows full measure to be taken of the data rich environment the MR generated motion field affords, permitting pixel to sub-pixel resolution within the reconstructed elastic property image. The algorithm's development is chronicled from initial 2D formulations to current 3D, multi-parameter elastic and viscoelastic approaches with experimental and early clinical results presented as well. Current results show that the method is capable of accurately deducing a quantitative measure of the elastic property distribution within gel phantoms as well as estimating the elasticity field within soft tissue. While the algorithm retains significant computational intensity, most noticeable in runtime considerations, a first effort at drastically reducing the time to image completion is made by taking advantage of obvious points for “macro-parallelization” within the process, making the effective runtime of the algorithm the quotient of the total processing time needed for image generation by the total number of processors available for computation. Additionally, a separate calculation for determining the covariance matrix for the generated parameter solution is formulated so that measures of image accuracy can be obtained in addition to the elastic property images themselves.