We introduce a novel computational method for a Mumford–Shah functional, which decomposes a given image into smooth regions separated by closed curves. Casting this as a shape optimization problem, we develop a gradient descent approach at the continuous level that yields nonlinear PDE flows. We propose time discretizations that linearize the problem and space discretization by continuous piecewise linear finite elements. The method incorporates topological changes, such as splitting and merging for detection of multiple objects, space–time adaptivity, and a coarse-to-fine approach to process large images efficiently. We present several simulations that illustrate the performance of the method and investigate the model sensitivity to various parameters.