This paper introduces a new approach for the joint alignment of a large collection of segmented images into the same system of coordinates while estimating at the same time an optimal common coordinate system. The atlas resulting from our group-wise alignment algorithm is obtained as the hidden variable of an Expectation-Maximization (EM) estimation. This is achieved by identifying the most consistent label across the collection of images at each voxel in the common frame of coordinates. In an iterative process, each subject is iteratively aligned with the current probabilistic atlas until convergence of the estimated atlas is reached. Two different transformation models are successively applied in the alignment process: an affine transformation model and a dense non-rigid deformation field. The metric for both transformation models is the mutual information that is computed between the probabilistic atlas and each subject. This metric is optimized in the affine alignment step using a gradient based stochastic optimization (SPSA) and with a variational approach to estimate the non-rigid atlas to subject transformations. A first advantage of our method is that the computational cost increases linearly with the number of subjects in the database. This method is therefore particularly suited for a large number of subjects. Another advantage is that, when computing the common coordinate system, the estimation algorithm identifies weights for each subject on the basis of the typicality of the segmentation. This makes the common coordinate system robust to outliers in the population. Several experiments are presented in this paper to validate our atlas construction method on a population of 80 brain images segmented into 4 labels (background, white and gray matters and ventricles). First, the 80 subjects were aligned using affine and dense non-rigid deformation models. The results are visually assessed by examining how the population converges closer to a central tendency when the deformation model allows more degrees of freedom (from affine to dense non-rigid field). Second, the stability of the atlas construction procedure for various sizes of population was investigated by starting from a subset of the total population which was incrementally augmented until the total population of 80 subjects was reached. Third, the consistency of our group-wise reference (hidden variable of the EM algorithm) was also compared to the choice of an arbitrary subject for a subset of 10 subjects. According to William’s index, our reference choice performed favorably. Finally, the performance of our algorithm was quantified on a synthetic population of 10 subjects (generated using random B-Spline transformations) using a global overlap measure for each label. We also measured the robustness of this measure to the introduction of noisy subjects in the population.