The parabolic leaf spring is widely used in modern vehicle suspension systems because it has many desirable features, such as weak interleaf friction and light weight. In this paper, the parabolic leaf spring is analyzed using automatic geometric data acquisition (AGDA) and the finite element (FE) absolute nodal coordinate formulation (ANCF). In order to account for manufacturing considerations in developing the virtual models, the relaxed uniform cubic B-Spline is used to represent the leaf spring profile curves and geometry. Three-dimensional scanning techniques based on structured light and multiple images are explored in this study to automatically extract the leaf spring complex geometric data. A new procedure is proposed to develop the ANCF/FE mesh from the physical object. Both double-leaf uniform-thickness and double-leaf parabolic spring models are developed and analyzed using different ANCF elements. Using ANCF geometry, piecewise linearly tapered parabolic leaf spring models are constructed, accounting for the leaf pre-stress. The interleaf contact is enforced using a penalty approach and a smoothed Coulomb friction model. It is shown that the fully-parameterized low-order beam and plate elements suffer from locking problems, while the thin plate element can lead to less accurate results. The use of the new strain-split method (SSM) as a locking alleviation technique is also examined in this investigation. It is shown that while the current SSM implementation can be effective in solving the locking problem in the case of symmetric bending-dominant loading, it may not produce accurate results in the case of torsional loading. The comparative study performed demonstrates that the higher order ANCF beam element is more suitable for developing the leaf spring models compared to other ANCF elements considered in this investigation. The numerical results obtained show that the friction effect in parabolic leaf springs is much weaker than that in the leaf spring with uniform-thickness.