We present a comprehensive methodology for the simulation of astronomical images from optical survey telescopes. We use a photon Monte Carlo approach to construct images by sampling photons from models of astronomical source populations, and then simulating those photons through the system as they interact with the atmosphere, telescope, and camera. We demonstrate that all physical effects for optical light that determine the shapes, locations, and brightnesses of individual stars and galaxies can be accurately represented in this formalism. By using large scale grid computing, modern processors, and an efficient implementation that can produce 400,000 photons/second, we demonstrate that even very large optical surveys can be now be simulated. We demonstrate that we are able to: 1) construct kilometer scale phase screens necessary for wide-field telescopes, 2) reproduce atmospheric point-spread-function moments using a fast novel hybrid geometric/Fourier technique for non-diffraction limited telescopes, 3) accurately reproduce the expected spot diagrams for complex aspheric optical designs, and 4) recover system effective area predicted from analytic photometry integrals. This new code, the photon simulator (PhoSim), is publicly available. We have implemented the Large Synoptic Survey Telescope (LSST) design, and it can be extended to other telescopes. We expect that because of the comprehensive physics implemented in PhoSim, it will be used by the community to plan future observations, interpret detailed existing observations, and quantify systematics related to various astronomical measurements. Future development and validation by comparisons with real data will continue to improve the fidelity and usability of the code.