We state vertebral body (VB) segmentation in MRI as a distribution-matching problem, and propose a convex-relaxation solution which is amenable to parallel computations. The proposed algorithm does not require a complex learning from a large manually-built training set, as is the case of the existing methods. From a very simple user input, which amounts to only three points for a whole volume, we compute a multi-dimensional model distribution of features that encode contextual information about the VBs. Then, we optimize a functional containing (1) a feature-based constraint which evaluates a similarity between distributions, and (2) a total-variation constraint which favors smooth surfaces. Our formulation leads to a challenging problem which is not directly amenable to convex-optimization techniques. To obtain a solution efficiently, we split the problem into a sequence of sub-problems, each can be solved exactly and globally via a convex relaxation and the augmented Lagrangian method. Our parallelized implementation on a graphics processing unit (GPU) demonstrates that the proposed solution can bring a substantial speed-up of more than 30 times for a typical 3D spine MRI volume. We report quantitative performance evaluations over 15 subjects, and demonstrate that the results correlate well with independent manual segmentations.