Brain extraction is a standard preprocessing step for subsequent tasks such as bias field correction, tissue segmentation, and cortical surface reconstruction. This process is challenging to fully automate, due to the anatomical variations of the brain and the imperfections in MR images. Existing approaches could be categorized into boundary-based¹, region-based², atlas-based³, learning-based⁴, and hybrid⁵ methods, all of which have strengths and weaknesses (e.g., sensitivity to noise for boundary/region-based, computational costs for atlas-based, feature/classifier design for learning-based methods, etc.) We herein present a generic supervised learning framework for fast, accurate, and robust brain extraction with extendibility to multi-contrast data.Our framework is based on multi-resolution classification on uni/multimodal MR data. First, standard preprocessing steps are applied to normalize the intensity distributions for the input data and affinely register it to a preselected template. Next, the training stage is triggered to build binary classifiers at four spatial resolution levels (i.e., the original images and the ones downsampled by 2, 4, and 8), which involves data sampling, feature extraction, and random forests learning for each level. Specifically, a sampling region is estimated by morphologically processing the ground truth mask to highlight a narrow band along the boundaries. For each randomly drawn sample from this region, the conventional spatial features, multi-scale intensity contextual features⁶, and spatial prior features (i.e., contextual features computed on the averaged group mask at the coarsest level) are extracted, followed by the use of random forests⁷ as binary classifiers. For an unseen image, starting with the average group mask at the coarsest level, we hierarchically refine the boundary by classifying voxels in a fixed-size narrow band along previous estimations. At the finest level, we couple the voxelwise classification with a narrow band level set approach⁸ using Chan-Vese⁹ region force to dynamically determine the test voxels. As the front propagates, the narrow band shifts accordingly and classification scores are computed only for the newly appeared voxels. This allows the surface to recover from previous mistakes without examining a large search region. We also induce a curvature term in the front propagation as regularization to maintain a smooth closed surface.We evaluate our framework on multiple datasets with different acquisition sequences and scanner types using single or multi-contrast images. We first consider the LONI-LPBA40¹⁰ dataset, which comprises 40 normal subjects each with a 1.5T T1-weighted image and a manually delineated brain mask. We randomly split the dataset into 20 for training and the other 20 for testing. The worst, an average, and the best results are illustrated in Figure 1, together with those using BET, a publicly available tool in FSL. Our segmentations are accurate while BET consistently oversegments the brain near the hippocampus. This observation is reflected in Figure 2 for all test subjects, as our method yields higher Dice, Jaccard, and specificity values, while BET achieves higher sensitivity due to oversegmentation. The same classifiers are then tested on a traumatic brain injury dataset of 10 subjects, each with a 3T T1-weighted image, and acceptable results are obtained (Figure 3). Additional tests are performed on the ATAG¹¹ dataset of 53 normal subjects, each with a 7T T1-weighted image. We randomly select 20 subjects to retrain the classifiers and test on the remaining 33. As there is no manual annotation, we use the brain masks produced by the CBS¹² tools for MIPAV as pseudo ground truth. We again obtain results that are better than the ones produced (whenever reasonable) by the CBS tools (Figure 4). In our final experiment, we consider a subset of the HCP¹³ datasets, which contains 50 normal subjects of 3T T1-weighted and T2-weighted images. To adapt our algorithm for multi-modal images, we extract the intensity contextual features from each co-registered image sequence independently and concatenate them to form the entire feature set. Considering the brain masks computed by the HCP processing pipeline as a pseudo ground truth, we retrain our models using randomly selected 20 cases and evaluate those models on the remaining 30. With a mean Dice coefficient of 0.98 and Jaccard index of 0.96, our multi-contrast framework achieves high accuracy, with a typical case shown in Figure 5.We present a generic framework for MR brain extraction using hierarchical voxelwise classification coupled with narrow band level set. Through a series of experiments, we demonstrate that our method is accurate and robust across various scanners at different strength of magnetic field with natural extendibility for multi-contrast data. It is also very fast, taking only 30 seconds to strip the skull on a standard preprocessed image.