With the advances in cardiovascular imaging technologies in recent years, 4D (3D+time) patient-specific modeling of the heart has attracted many research interests. Computational modeling approaches such as Computational Fluid Dynamics (CFD) and Finite Element (FE) have been increasingly used to quantitatively diagnose and predict cardiovascular diseases. In these methods, the geometrical reconstruction of the heart anatomy is usually an indispensable step. This work presents a robust method for reconstructing time-varying subdivision surfaces from the segmentation masks of cardiac images. We first reconstruct a 3D mesh for the first time step by iteratively fitting an initial mesh based on error and tension terms. Each subsequent time step uses the model from its previous time step as the control mesh for subdivision surface fitting. This method preserves the 1-to-1 vertex correspondence between meshes in different time steps and allows us to control the mesh quality (i.e. resolution, smoothness, and accuracy). Furthermore, in contrast to contour-based algorithms, our method can handle non-trivial topological changes such as holes and tunnels. The method has been tested on 3D and 4D datasets of different modalities (i.e. CT and MRI), resolutions, and chambers. For creating visually appealing results, we show that synthetic textures can be mapped to the 4D reconstruction due to the vertex correspondence. We also quantitatively evaluate the reconstructed 3D meshes in terms of mesh quality and conformity to the data.