The human cortical folding is intriguingly complex in its variability and regularity across individuals. Exploring the principal patterns of cortical folding is of great importance for neuroimaging research. The term-born neonates with minimum exposure to the complicated environments are the ideal candidates to mine the postnatal origins of principal cortical folding patterns. In this work, we propose a novel framework to study the gyral patterns of neonatal cortical folding. Specifically, first, we leverage multi-view curvature-derived features to comprehensively characterize the complex and multi-scale nature of cortical folding. Second, for each feature, we build a dissimilarity matrix for measuring the difference of cortical folding between any pair of subjects. Then, we convert these dissimilarity matrices as similarity matrices, and nonlinearly fuse them into a single matrix via a similarity network fusion method. Finally, we apply a hierarchical affinity propagation clustering approach to group subjects into several clusters based on the fused similarity matrix. The proposed framework is generic and can be applied to any cortical region, or even the whole cortical surface. Experiments are carried out on a large dataset with 600+ term-born neonates to mine the principal folding patterns of three representative gyral regions.