This paper presents a robust unsupervised harmonic co-clustering method for profiling arbor morphologies for ensembles of reconstructed brain cells (e.g., neurons, microglia) based on quantitative measurements of the cellular arbors. Specifically, this method can identify groups and sub-groups of cells with similar arbor morphologies, and simultaneously identify the hierarchical grouping patterns among the quantitative arbor measurements. The robustness of the proposed algorithm derives from use of the diffusion distance measure for comparing multivariate data points, harmonic analysis theory, and a Haar-like wavelet basis for multivariate data smoothing. This algorithm is designed to be practically usable, and is embedded into the actively linked three-dimensional (3-D) visualization and analytics system in the free and open source FARSIGHT image analysis toolkit for interactive exploratory population-scale neuroanatomic studies. Studies on synthetic datasets demonstrate its superiority in clustering data matrices compared to recent hierarchical clustering algorithms. Studies on heterogeneous ensembles of real neuronal 3-D reconstructions drawn from the NeuroMorpho database show that the proposed method identifies meaningful grouping patterns among neurons based on arbor morphology, and revealing the underlying morphological differences.