The fractional Laplacian operator (−∆)s on a bounded domain Ω can be realized as a Dirichlet-to-Neumann map for a degenerate elliptic equation posed in the semi-infinite cylinder Ω × (0,∞). In fact, the Neumann trace on Ω involves a Muckenhoupt weight that, according to the fractional exponent s, either vanishes (s < 1/2) or blows up (s > 1/2). On the other hand, the normal trace of the solution has the reverse behavior, thus making the Neumann trace analytically well-defined. Nevertheless, the solution develops an increasingly sharp boundary layer in the vicinity of Ω as s decreases. In this work, we extend the technology of automatic hp-adaptivity, originally developed for standard elliptic equations, to the energy setting of a Sobolev space with a Muckenhoupt weight, in order to accommodate for the problem of interest. The numerical evidence confirms that the method maintain exponential convergence. Finally, we discuss image denoising via the fractional Laplacian. In the image processing community, the standard way to apply the fractional Laplacian to a corrupted image is as a filter in Fourier space. This construction is inherently affected by the Gibbs phenomenon, which prevents the direct application to “spliced” images. Since our numerical approximation relies instead on the extension problem, it allows for processing different portions of a noisy image independently and combine them, without complications induced by the Gibbs phenomenon.