High-dimensional manifold modeling increases the precision and performance of cortical morphometry analysis by densely sampling on the grey matters. But this also brings redundant information and increased computational burden. Gaussian process regression has been used to tackle this problem by learning a mapping to a low-dimensional subspace. However, current methods may not take relevant morphometric properties, usually measured by geometric features, into account, and as a result, may generate morphometrically insignificant selections. In this paper, we propose a morphometric Gaussian process (M-GP) as a novel Bayesian model on the gray matter tetrahedral meshes. We also implement an M-GP regression landmarking algorithm as a manifold learning method for non-linear dimensionality reduction. The definition of M-GP involves a scale-invariant wave kernel signature distance map measuring the local similarities of geometric features, and a heat flow entropy which implicitly embeds the global curvature flow. With such a design, the prior knowledge fully encodes the geometric information so that a posterior predictive inference is morphometrically significant. In experiments, we use 518 grey matter tetrahedral meshes generated from structural magnetic resonance images of a publicly available Alzheimer's disease imaging cohort to empirically and numerically evaluate our method. The results verify that our method is theoretically and experimentally valid in selecting a representative subset from the original massive data. Our work may benefit any studies involving large-scale or iterative computations on extensive manifold-valued data, including morphometry analyses and general medical data processing.