We present a new Riemannian framework for comprehensive statistical shape analysis of 3D objects, represented by their boundaries (parameterized surfaces). By comprehensive framework, we mean tools for registration, comparison, averaging, and modeling of observed surfaces. Registration is analogous to removing all shape preserving transformations, which include translation, scale, rotation and re-parameterization. This framework is based on a special representation of surfaces termed square-root normal fields and a closely related elastic metric. The main advantages of this method are: (1) the elastic metric provides a natural interpretation of shape deformations that are being quantified, (2) this metric is invariant to re-parameterizations of surfaces, and (3) under the square-root normal field transformation, the complicated elastic metric becomes the standard L2 metric, simplifying parts of the implementation. We present numerous examples of shape comparisons for various types of surfaces in different application areas. We also compute average shapes, covariances and perform principal component analysis to explore the variability in different shape classes. These quantities are used to define generative shape models and for random sampling. Specifically, we showcase the applicability of the proposed framework in shape analysis of anatomical structures in different medical applications including Attention Deficit Hyperactivity Disorder and endometriosis.