This paper considers the numerical solution of elliptic boundary value problems with a complicated (nonperiodic) diffusion matrix which is smooth but highly oscillating on very different scales. We study the influence of the scales and amplitudes of these oscillations to the regularity of the solution. We introduce weighted Sobolev norms of integer order, where the (p + 1)st seminorm is weighted by properly scaled pth derivative of the diffusion coefficient. The constants in the regularity estimates then turn out to depend only on global lower and upper bounds of the diffusion matrix but not on its derivatives; in particular, the constants are independent of the scales of the oscillations. These regularity results give rise to error estimates for hp-finite element discretizations with scale-adapted distribution of the mesh cells. The adaptation of the mesh is explicit with respect to the local variations of the diffusion coefficient. Numerical results illustrate the efficiency of these oscillation adapted finite elements, in particular, in the preasymptotic regime.