Computing image atlases that represent a large data set is an important first step for statistical analysis of images, quantifying shapes, etc. Most current approaches estimate a single atlas to represent the average of a large population of images. However, a single atlas is not sufficiently expressive to capture distributions of images with multiple modes. In this paper, we present a Gaussian mixture model for building diffeomorphic multi-atlases that can represent sub-populations without knowing the category of each observed data point. In our probabilistic model, we treat diffeomorphic image transformations as latent variables, and integrate them out using a Monte Carlo Expectation Maximization (MCEM) algorithm via Hamiltonian Monte Carlo (HMC) sampling. Another benefit of our model is that the Gaussian mixture modeling simultaneously does automatic clustering of the dataset. This work is submitting to MICCAI 2015.