We present a new method for the analysis of images, a fundamental task in observational astronomy. It is based on the linear decomposition of each object in the image into a series of localized basis functions of different shapes, which we call ‘shapelets’. A particularly useful set of complete and orthonormal shapelets is that consisting of weighted Hermite polynomials, which correspond to perturbations around a circular Gaussian. They are also the eigenstates of the two dimensional quantum harmonic oscillator, and thus allow us to use the powerful formalism developed for this problem. One of their special properties is their invariance under Fourier transforms (up to a rescaling), leading to an analytic form for convolutions. The generator of linear transformations such as translations, rotations, shears and dilatations can be written as simple combinations of raising and lowering operators. We derive analytic expressions for practical quantities, such as the centroid (astrometry), flux (photometry) and radius of the object, in terms of its shapelet coefficients. We also construct polar basis functions which are eigenstates of the angular momentum operator, and thus have simple properties under rotations. As an example, we apply the method to Hubble Space Telescope images, and show that the small number of shapelet coefficients required to represent galaxy images lead to compression factors of about 40 to 90. We discuss applications of shapelets for the archival of large photometric surveys, for weak and strong gravitational lensing and for image deprojection.