We propose a new nonnegative matrix underapproximation model for images. We design an algorithm that runs in linear time in the dimensions of the input matrix. This allows us to extract sequentially localized and spatially coherent features. We illustrate the effectiveness of our approach on a synthetic data set, facial and hyperspectral images. We show that it competes favorably with comparable state-of-the-art techniques.