Polycrystalline materials, such as metals, are comprised of heterogeneously oriented crystals. Observed crystal orientations are modelled as a sample from an orientation distribution function (ODF), which determines a variety of material properties and is therefore of great interest to practitioners. Observations consist of quaternions, 4dimensional unit vectors reflecting both orientation and rotation of a single crystal. Thus, an ODF must account for known crystal symmetries as well as satisfy the unit length constraint. A popular method for estimating ODFs non-parametrically is symmetrized kernel density estimation. However, disadvantages of this approach include difficulty in interpreting results quantitatively, as well as in quantifying uncertainty in the ODF. We propose to use a mixture of symmetric Bingham distributions as a flexible parametric ODF model, inferring the number of mixture components, the mixture weights, and scale and location parameters based on crystal orientation data. Furthermore, our Bayesian approach allows for structured uncertainty quantification of the parameters of interest. We discuss details of the sampling methodology and conclude with analyses of various orientation datasets, interpretations of parameters of interest, and comparison with kernel density estimation methods.