Brain network analysis is a promising tool in studies of the human brain’s functional organization and neuropsychiatric disorders. For neuroimaging data based brain network analysis, network nodes are typically defined as distinct greymatter regions delineated by anatomical brain atlases, random parcellation of the brain space, or image voxels, resulting in brain network nodes with different spatial scales. As precise functional organization of the brain remains unclear, it is challenging to determine a proper spatial scale in practice. The brain network nodes defined anatomically or randomly do not necessarily possess desired properties of ideal brain network nodes, i.e., functional homogeneity within each individual node, functional distinctiveness across different nodes, and functional consistence of the same node across different subjects. To obtain a definition of brain network nodes with the desired properties, a brain parcellation method based on functional information is proposed to achieve a brain parcellation consistent across subjects and highly in agreement with the functional organization of the brain. Particularly, spatially contiguous voxel-wise functional information of the brain fMRI data recursively aggregate according to inter-voxel/region functional affinity from voxel level to coarser scales, resulting in a brain parcellation with a multi-level hierarchy. A trade-off between functional homogeneity and distinctiveness is determined by identifying the hierarchy level with network measures highly consistent across subjects. The proposed method has been validated on resting-state fMRI datasets for functional network analysis, and the results demonstrate that brain networks constructed with 200~500 nodes could achieve the highest inter-subject consistence.