Disease heterogeneity is a significant obstacle to understanding pathological processes and delivering precision diagnostics and treatment. Clustering methods have gained popularity in stratifying patients into subpopulations (i.e., subtypes) of brain diseases using imaging data. However, unsupervised clustering approaches are often confounded by anatomical and functional variations not related to a disease or pathology of interest. Semi-supervised clustering techniques have been proposed to overcome this and, therefore, capture disease-specific patterns more effectively. An additional limitation of both unsupervised and semi-supervised conventional machine learning methods is that they typically model, learn and infer from data at a basis of feature sets pre- defined at a fixed scale or scales (e.g, an atlas-based regions of interest). Herein we propose a novel method, Multi-scAle heteroGeneity analysIs and Clustering (MAGIC), to depict the multi-scale presentation of disease heterogeneity, which builds on a previously proposed semi-supervised clustering method, HYDRA. It derives multi-scale and clinically interpretable feature representations and exploits a double-cyclic optimization procedure to drive inter-scale-consistent disease subtypes or neuroanatomical dimensions effectively. More importantly, to fill in the gap of understanding under what conditions the clustering model can estimate true heterogeneity related to diseases, we conducted extensive and systematic semi-simulated experiments to evaluate the proposed method on a sizeable healthy control sample from the UK Biobank (N=4403). We then applied MAGIC to real imaging data of Alzheimers disease (ADNI, N=1728) to demonstrate its potential and challenges in dissecting the neuroanatomical heterogeneity of brain diseases. Taken together, we aim to provide guidelines on when such analyses can succeed or should be taken with caution. The code of the proposed method is publicly available at https://github.com/anbai106/MAGIC.