Simulations of three dimensional ultrasound propagation in heterogeneous media are computationally intensive due to the constraints arising from the large size of the domain, which is on the order of hundreds of wavelengths, and the small size of scatterers, which are much smaller than a wavelength. Consequently, three dimensional ultrasound imaging simulations are currently based on models that simplify the propagation physics... Here the full three dimensional wave physics is simulated with finite differences to generate ultrasound images of the human body based directly on the first principles of propagation and backscattering. The Visible Human project, a 3D data set of the human body that was generated with photographs of 0.33 mm cryosections, is converted into 3D acoustical maps. A full-wave nonlinear acoustic simulation tool is used to propagate ultrasound into the liver with a 2D intercostal ultrasound array in a $93 \times 39 \times 22$ mm domain with $6\times10^8$ points. Imaging metrics, based on the beamplots, root-mean-square phase aberration, spatial coherence lengths, and contrast-to-noise ratio are used to characterize the image quality. It is shown that the harmonic image quality is better than the fundamental image quality due, in part, to a narrower beam profile. The root-mean-square estimate of aberration after propagation through the simulated body wall is shown to be low (23.4 ns), consistently with previous reports of aberration measured experimentally in a human body wall. The spatial coherence measured at the transducer surface indicates that a transducer array element size of $<0.81 \lambda$ would be required to fully sample the acoustic field. These simulated three dimensional ultrasound images based directly on propagation physics provide a platform to investigate the sources of image degradation in three dimensions (included in Part II). (read more)