The study of myocardial motion is essential for understanding the normal heart function and developing new treatments for cardiovascular diseases. The goals of my PhD research is to develop new methods for cardiac deformation recovery from 3D magnetic resonance (MR) images. The main contribution of my work is that the proposed methods are guaranteed to generate exactly or nearly incompressible deformations. This is a desirable property since the myocardium has been shown to be close to incompressible. From the recovered deformation, one can directly compute a number of clinically useful parameters, including strains. The first method for 3D deformation recovery of the left ventricular wall (LV) from anatomical cine MRI is based on a deformable model that is incompressible. This method is not suitable for the deformation recovery of the biventricular wall. The second method is based on a 3D deformable model that is nearly incompressible. The model uses a matrix-valued radial basis function to represent divergence free displacement fields, which is a first order approximation of incompressibility. This representation allows for deformation modeling of arbitrary topologies with a relatively small number of parameters, which is suitable for representing the motion of the multi-chamber structure of the heart. The two methods have similar performance. A method to obtain a smooth and accurate surface of the myocardium wall is needed to illustrate the cardiac deformation recovery. I present a novel method for the generation of endocardial and epicardial surface meshes. The same algorithm is independently used to generate the surface meshes of the epicardium and endocardium of the four cardiac chambers. It provides smooth meshes despite the strong voxel anisotropy, which is not the case for the marching cubes algorithm. Phase velocity MRI is an acquisition technique that contains more information about the myocardial motion than cine MRI. I present a method to interpolate the velocity vector field in a phase velocity MRI sequence. The method uses an interpolation model that provides a continuous divergence free velocity vector field, which means that the corresponding deformation is incompressible.