In this dissertation, we propose new approaches to the surface registration problem by coupling continuous geometry-based methods and combinatorial graph-based methods. On the one hand, geometry-based methods explore the intrinsic properties of the surfaces to simplify the search of correspondences among surfaces undergoing very large deformations. However, these methods are usually based on certain ideal assumptions on the qualities of the input, such as noise-free or no occlusions. Hence they are usually sensitive to uncertainties from the input that are common in real-world data. On the other hand, graph-based methods are better at dealing with uncertainties due to their statistics nature. Nevertheless, without exploring the geometric properties of the surface, discrete graph-based methods usually suffer from discretization error and high computational complexity. Thus, by exploring the relation between the two approaches, we show that our new approaches deal with surface registration problems under very challenging situations. To this purpose, three distinct approaches are explored in this thesis that achieves dense surface registration in different scenarios. In the first approach, we cast the surface registration into a high-order graph matching problem, through the minimization of an energy function based on multiple measurements of geometric/appearance similarities and deformation priors. Our method takes advantage of conformal mapping based method which derives a closed-form solution to dense surface matching. To this end, we design an efficient way to select a finite number of matching candidates for each point of the source surface based on the a sparse set of correspondences, which naturally induces an efficient two-stage optimization approach for the dense surface registration problem. In the sparse matching stage, the high-order interactions among a sparse set of feature points on the surfaces are used to encode the isometric deformation error using conformal mapping. In the dense registration stage, the high-order interactions of a dense set of sampling points are considered to encode the isometric deformation error as well as orientation consistency. Meanwhile, we also propose the first solution to the high-order graph matching problem that solves partial matching. Our method is validated through a series of experiments demonstrating its accuracy and efficiency, notably in challenging cases of large and/or non-isometric deformations, or meshes that are partially occluded. In the second approach, we propose a graph-based formulation for tracking surfaces in a sequence. In order to deal with noises in the input, we propose a robust metric for the cost of matching arbitrary correspondences, which is defined as the lowest feature differences across this set of matchings that cause the particular correspondence to match. We show that for surface tracking applications, the matching cost can be efficiently computed in the conformal mapping domain. Such a matching cost is then integrated into a complete probabilistic tracking framework that enforces spatial and temporal motion consistencies, as well as error drifts and occlusions. Compared to previous 3D surface tracking approaches that either assume isometric deformations or consistent features, our method achieves dense, accurate tracking results, which we demonstrate through a series of dense, anisometric 3D surface tracking experiments. In the third approach, we accurately characterize arbitrary deformations between two surfaces and propose a high-order graphical model for the surface registration problem. From Riemannian geometry, the local deformation at each point of a surface can be characterized by the eigenvalues of a special transformation matrix between two canonically parameterized domains. This local transformation is able to characterize all the deformations (i.e., diffeomorphisms) between surfaces while being independent of both intrinsic (parametrization) and extrinsic (embedding) representations. In particular, we show that existing deformation representations (e.g., isometry or conformality) can be viewed as special cases of the proposed local deformation model. Furthermore, a computationally efficient, closed-form solution is derived in the discrete setting via finite element discretization. Based on the proposed deformation model, the shape registration problem is formulated as a high-order Markov Random Field (MRF) defined on the simplicial complex (e.g., planar or tetrahedral mesh). An efficient high-order MRF optimization algorithm is designed for such a special structured MAP-MRF problem, which can be implemented in a distributed fashion and requires minimal memory. Finally, we demonstrate the speed and accuracy of the proposed approach in the applications of shape registration and tracking.