Harmonic phase (HARP) MRI is used to measure myocardial motion and strain from tagged MR images. HARP MRI uses limited number of samples from the spectrum of the tagged images to reconstruct motion and strain. The HARP strain maps, however, suffer from artifacts that limit the accuracy of the computations and degrade the appearance of the strain maps. Causes of these, so called 'zebra', artifacts include image noise, Gibbs ringing, and interference from other Fourier spectral peaks. Computing derivatives of the HARP phase, which are needed to estimate strain, further accentuates these artifacts. Previous methods to reduce these artifacts include 1-D and 2-D nonlinear filtering of the HARP derivatives, and a 2-D linear filtering of unwrapped HARP phase. A common drawback among these methods is the lack of proper segmentation of the myocardium from the blood pool. Because of the lack of segmentation, the noisy phase values from the blood pool enter into the computation in the smoothed strain maps, which causes artifacts. In this work, we propose a smoothing method based on anisotropic diffusion that filters the HARP derivatives strictly within the myocardium without the need for prior segmentation. The information about tissue geometry and the strain distribution is used to restrict the smoothing to within the myocardium, thereby ensuring minimum distortion of the final strain map. Preliminary results demonstrate the ability of anisotropic diffusion for better artifact reduction and lesser strain distortion than the existing methods.