Modeling and simulation of rare events is a problem that has been largely ignored by the image processing community, but is of great interest in many areas of science. Rare events in physical systems are responsible for many failure modes, and as such must be precisely understood. We propose a novel method for simulating and modeling rare images using asymptotically efficient importance sampling, and apply it to binary images of interest in statistical mechanics and materials science. These rare images correspond to the occurrence of rare events in systems modeled by Gibbs distributions, more commonly known as Markov random fields (MRFs) in image processing. We will first give a precise definition of a rare event, in terms of a rare event statistic and a rare event region. The rare event statistic of interest here will be the per-site magnetization of a ferromagnet under the Ising model, but this could be replaced by many other statistics for other problems, such as boundary length in two-phase material microstructures, for example. For the given rare event statistic, we estimate the asymptotically efficient importance sampling (AEIS) distribution, which is based on a large deviation principle (LDP); draw samples of rare binary images from this AEIS distribution; and estimate rare event probabilities for several different rare event regions. Theoretically, the AEIS sampling distribution gives an unbiased estimator with the lowest variance asymptotically from a class of importance sampling distributions that are practically feasible for Monte Carlo Markov chain simulation. Finally, we fit large deviations rate functions from simulations using several rare event regions. This allows us to compute probability estimates associated with a given rare event statistic for any rare event region of interest without requiring further simulations.