In this paper we present an automatic statistical intensity based-approach for extracting the 3D cerebrovascular system from time-of-flight (TOF) magnetic resonance angiography (MRA) data. The voxels of the dataset are classified as either background tissues, which are modeled by a finite mixture of one Rayleigh and two normal distributions, or blood vessels, which are modeled by one normal distribution. We show that the proposed models fit the clinical data properly and result in fewer misclassified vessel voxels. We estimated the parameters of each distribution using the expectation maximization (EM) algorithm. Since the convergence of the EM is sensitive to the initial estimate of the parameters, a novel method for parameter initialization, based on histogram analysis, is provided. A new geometrical phantom motivated by a statistical analysis was designed to validate the accuracy of our method. The algorithm was also tested on 20 in-vivo datasets. The results showed that the proposed approach provides accurate segmentation, especially those blood vessels of small sizes.