The cortical folding of the human brain is highly complex and variable across individuals. Mining the major patterns of cortical folding from modern large-scale neuroimaging datasets is of great importance in advancing techniques for neuroimaging analysis and understanding the inter-individual variations of cortical folding and its relationship with cognitive function and disorders. As the primary cortical folding is genetically influenced and has been established at term birth,neonates with the minimal exposure to the complicated postnatal environmental influence are the ideal candidates for understanding the major patterns of cortical folding. In this paper,for the first time,we propose a novel method for discovering the major patterns of cortical folding in a large-scale dataset of neonatal brain MR images (N = 677). In our method,first,cortical folding is characterized by the distribution of sulcal pits,which are the locally deepest points in cortical sulci. Because deep sulcal pits are genetically related,relatively consistent across individuals,and also stable during brain development,they are well suitable for representing and characterizing cortical folding. Then,the similarities between sulcal pit distributions of any two subjects are measured from spatial,geometrical,and topological points of view. Next,these different measurements are adaptively fused together using a similarity network fusion technique,to preserve their common information and also catch their complementary information. Finally,leveraging the fused similarity measurements,a hierarchical affinity propagation algorithm is used to group similar sulcal folding patterns together. The proposed method has been applied to 677 neonatal brains (the largest neonatal dataset to our knowledge) in the central sulcus,superior temporal sulcus,and cingulate sulcus,and revealed multiple distinct and meaningful folding patterns in each region.