We propose a new multiscale image decomposition which offers a hierarchical, adaptive representation of general images. The starting point is a variational decomposition of an image, f = u0+v0, where [u0,v0] is the minimizer of a J-functional, J(f,c0; X,Y)=inf{u+v=f} {||u||X + c0 ||v||Yp}. Such minimizers are standard tools for denoising, deblurring, compression, ... of images, e.g., [Mumford-Shah] and [Rudin-Osher-Fatemi]. Here, u0 should capture `essential features' of f, to be separated from the spurious components in v0, and c0 is a fixed threshold which dictates separation of scales. To proceed, we iterate the refinement step [uj+1,vj+1] = arginf J(vj,c02j), leading to the hierarchical decomposition, f = ∑j=0k uj + vk. We focus our attention on the particular case of (X,Y)=(BV,L2) decomposition. The resulting hierarchical decomposition, f ~ ∑j uj, is essentially nonlinear. The questions of convergence, energy decomposition, localization and adaptivity are discussed. The decomposition is constructed by numerical solution of successive Euler-Lagrange equations. Numerical results illustrate applications to synthetic and real images (both grayscale and colored images).