Patterns and nonlinear waves, such as spots, stripes, and rotating spirals, arise prominently in many natural processes and in reaction-diffusion models. Our goal is to compute boundaries between parameter regions with different prevailing patterns and waves. We accomplish this by evolving randomized initial data to full patterns and evaluate feature functions, such as the number of connected components or their area distribution, on their sublevel sets. The resulting probability measure on the feature space, which we refer to as pattern statistics, can then be compared at different parameter values using the Wasserstein distance. We show that arclength predictor-corrector continuation can be used to trace out transition and bifurcation curves in parameter space by maximizing the distance of the pattern statistics. The utility of this approach is demonstrated through a range of examples involving homogeneous states, spots, stripes, and spiral waves.