A maximum a posteriori (MAP) technique is developed to identify solar features in cotemporal and cospatial images of line-of-sight magnetic flux, continuum intensity, and equivalent width observed with the NASA/National Solar Observatory (NSO) Spectromagnetograph (SPM). The technique facilitates human understanding of patterns in large data sets and enables systematic studies of feature characteristics for comparison with models and observations of long-term solar activity and variability. The method uses Bayes' rule to compute the posterior probability of any feature segmentation of a trio of observed images from per-pixel, class-conditional probabilities derived from independently-segmented training images. Simulated annealing is used to find the most likely segmentation. New algorithms for computing class-conditional probabilities from three-dimensional Gaussian mixture models and interpolated histogram densities are described and compared. A new extension to the spatial smoothing in the Bayesian prior model is introduced, which can incorporate a spatial dependence such as center-to-limb variation. How the spatial scale of training segmentations affects the results is discussed, and a new method for statistical separation of quiet Sun and quiet network is presented.