Hypersharpening consists in generating an unobservable high-spatial-resolution hyperspectral image by fusing an observed low-spatial-resolution hyperspectral image with an observed high-spatial-resolution panchromatic or multispectral one. The obtained image preserves the high spectral resolution of the first image and the high spatial resolution of the second one. Unlike standard hypersharpening methods that do not consider the spectral variability phenomenon, in this letter, a new approach, which addresses this phenomenon, is proposed for fusing hyperspectral and multispectral remote sensing images. This approach, linked to linear spectral unmixing methods, is based on an extension of nonnegative matrix factorization (NMF), namely the inertia-constrained pixel-by-pixel NMF (IP-NMF) algorithm. The developed fusion algorithm, called hyperspectral and multispectral data fusion based on IP-NMF (HMF-IPNMF), is applied to synthetic and real data sets. Experimental results clearly show that the developed fusion method yields sharpened hyperspectral images with higher spectral and spatial fidelities when compared to those provided by tested state-of-the-art methods that do not take spectral variability into account.