We introduce a novel approach for segmenting articulated spine shape models from medical images. A nonlinear low-dimensional manifold is created from a training set of mesh models to establish the patterns of global shape variations. Local appearance is captured from neighborhoods in the manifold once the overall representation converges. Inference with respect to the manifold and shape parameters is performed using a higher-order Markov random field (HOMRF). Singleton and pairwise potentials measure the support from the global data and shape coherence in manifold space respectively, while higher-order cliques encode geometrical modes of variation to segment each localized vertebra models. Generic feature functions learned from ground-truth data assigns costs to the higher-order terms. Optimization of the model parameters is achieved using efficient linear programming and duality. The resulting model is geometrically intuitive, captures the statistical distribution of the underlying manifold and respects image support. Clinical experiments demonstrated promising results in terms of spine segmentation. Quantitative comparison to expert identification yields an accuracy of 1.6 ± 0.6 mm for CT imaging and of 2.0 ± 0.8 mm for MR imaging, based on the localization of anatomical landmarks.