Radius-optimized efficient template matching for lesion detection from brain images

Computer-aided detection of brain lesions from volumetric magnetic resonance imaging (MRI) is in demand for fast and automatic diagnosis of neural diseases. The template-matching technique can provide satisfactory outcome for automatic localization of brain lesions; however, finding the optimal template size that maximizes similarity of the template and the lesion remains challenging. This increases the complexity of the algorithm and the requirement for computational resources, while processing large MRI volumes with three-dimensional (3D) templates. Hence, reducing the computational complexity of template matching is needed. In this paper, we first propose a mathematical framework for computing the normalized cross-correlation coefficient (NCCC) as the similarity measure between the MRI volume and approximated 3D Gaussian template with linear time complexity, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{\mathcal{O}}}\left( {{\varvec{a}}_{{{\varvec{max}}}} {\varvec{N}}} \right)$$\end{document}OamaxN, as opposed to the conventional fast Fourier transform (FFT) based approach with the complexity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{\mathcal{O}}}\left( {{\varvec{a}}_{{{\varvec{max}}}} {\varvec{N}}\log {\varvec{N}}} \right)$$\end{document}OamaxNlogN, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{N}}$$\end{document}N is the number of voxels in the image and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{a}}_{{{\varvec{max}}}}$$\end{document}amax is the number of tried template radii. We then propose a mathematical formulation to analytically estimate the optimal template radius for each voxel in the image and compute the NCCC with the location-dependent optimal radius, reducing the complexity to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{\mathcal{O}}}\left( {\varvec{N}} \right)$$\end{document}ON. We test our methods on one synthetic and two real multiple-sclerosis databases, and compare their performances in lesion detection with FFT and a state-of-the-art lesion prediction algorithm. We demonstrate through our experiments the efficiency of the proposed methods for brain lesion detection and their comparable performance with existing techniques.

A lesion in the brain represents a diseased area or one with abnormal tissue due to many factors such as infection, traumatic brain injury, tumors, multiple sclerosis (MS), stroke, etc 1 . The complications resulting from brain lesions are numerous including neurodegeneration, memory or vision loss, and behavioral changes 2 . The presence of a lesion in the brain often demands prompt clinical investigation to initiate the treatment process. Magnetic resonance imaging (MRI) contributes to the diagnosis significantly by visualizing the neuroanatomy with high soft-tissue contrast in vivo. Detection of the exact location of lesion is essential for clinical implications and thus assists in diagnosis and treatment planning. Manual assessment (i.e. visual screening) of individual two-dimensional (2D) slices from the whole MRI volume is the conventional practice of identifying the presence of lesions. This is a cumbersome task, especially in the presence of multiple lesions (e.g., MS, metastasis, etc.) or when performing it on a large dataset. Inter-and intra-rater variabilities are other caveats of manual lesion detection. The limitations of manual screening and volumetric analysis by trained experts have led to the development of semi-automatic and automatic segmentation algorithms in the last two decades, which can assist the clinician in identifying the region of interest (ROI) from the magnetic resonance (MR) image. The study conducted by Egger et al. 3 on fluid attenuation inversion recovery (FLAIR) MRI of MS lesions over 50 cases, reported that the segmentation methods such as lesion growth 4 and lesion prediction algorithm (LPA) 5 underperformed for quantifying the number of lesions on the subjects with high lesion load. According to them, the difficulties in the evaluation of the confluent lesions in periventricular region might be the cause of poor performance, which emphasizes the importance of prior detection of lesions in MRI volume, especially when multiple lesions are present in brain.
(5) f * h  Next, we rewrite the first factor of the denominator of Eq. (1) using the common expansion of the variance , where f p ⇀ X , p = 1, 2 , is calculated similarly to Eq. (8) (by fixing n = 1 ) as: Lastly, we calculate the second factor in the denominator of Eq. (1), again using the common expansion of the variance, as: To that end, we calculate the Fourier transform of the box function as F {h D } = D j=1 sinc aω j , where sincθ = (sin θ )/θ , which leads to the Fourier transform of the kernel F h (n) D = F {h D } n = D j=1 sinc n aω j via the convolution theorem. The integral of the template then is = D j=1 sinc n (0) = 1 , from which the template mean is computed as h (n) D = 1/(2b) D . As for the mean of the square of the template, we use Parseval's theorem as follows: We now define and compute 25,26 : Therefore, h S n,D f := S (n) n,D f , (1) will result in the following formula for the proposed fast approach of NCCC computation: NCCC values range from −1 to 1 . In this subsection, we compute the NCCC for several different values of a . However, note that S n,D f p are volumes independent of a and can be pre-computed, along with the second term in the numerator and the first factor in the denominator of Eq. (15). Furthermore, the second factor in the denominator is a scalar that is computed fast in O(1) for each a . Thus, in the proposed approach to estimate NCCC with the Gaussian template, the bulk of the computational cost is only due to the first term in the numerator of Eq. (15), which can be computed in O(N) for each a , with N the number of voxels in the image. On the contrary, the computational complexity of FFT is O N log N for each a 27 , making the overall computational cost of lesion detection noticeably higher for the FFT-based algorithm (template has to be zero-padded to the size of the image) than the proposed approach. This difference in computational cost is particularly amplified given that the NCCC needs to be repeatedly computed for many values of a = 1, 2, . . . , a max . O(N). Defining NCCC . We intend to optimize the NCCC with respect to the template radius; however, analytical maximization of Eq. (15) with respect to the template radius is difficult. Therefore, we aimed at redefining NCCC, which would be suitable for voxel-wise maximization of template radius. A flowchart of the proposed method is illustrated in Fig. 1. Let us assume the template,

Template matching in
, and the ROI are both Gaussian, where � • � denotes the Euclidean norm. As in the previous section, we solve the NCCC by considering all terms of numerator and denominator, individually. We begin by focusing on the first term of the numerator of Eq. (1). For a Gaussian-shaped lesion, this term can be solved as: Figure 1. A flowchart depicting the proposed O(N) method. We first calculate a closed-form formula for the optimal radius by Gaussian modeling of the lesion, which we then apply to the local image statistics (estimated from a smoothed version of the MRI volume) to compute the optimal radius at each voxel. We then use an approximated Gaussian template (by consecutive convolutions with the boxcar kernel) to compute the NCCC at each voxel, from which we keep a pre-determined number of top detected lesions. www.nature.com/scientificreports/ where, ⇀ µ( ⇀ X) and σ ( ⇀ X) are the mean (with respect to ( ⇀ X) ) and standard deviation (SD) of the ROI, i.e. the location and the radius of the lesion, respectively. We now simplify Eq. (16) by reducing the power of exponential via completing the square as, Now we can rewrite Eq. (16) as, After the change of variables, and assuming the radius of engulfing template to be sufficiently large, i.e. b → ∞ , Eq. (18) can be reduced by using the properties of the error function ( erf −v e −t 2 dt; erf (∞) = 1 ) to the following, which completes the calculation of the first term of the numerator: We now focus on the first term of denominator, i.e. image variance inside the template. We use the expansion , and compute f 2 ( ⇀ X) and f ( ⇀ X) 2 as follows.
is computed in the same way as in Eq. (20), and for a large b , one can see that, Subtracting Eq. (22) from Eq. (21), we get an approximation for the first term of denominator of NCCC for a large b as, Next, we took on the second term of denominator, i.e. the template variance. We followed the expansion solution to initiate its computation. First, the mean of the square of the template can be computed as, Like previous steps, for a large b → ∞ , Eq. (24) can be simplified as, (2b) D g D 2 = √ πa D . Similarly to above, the mean of the template, (2b) D g D 2 , is zero for b → ∞ , and we obtain the approximate solution (for large b ) of the second term of denominator of NCCC as, Next, we compute the second term of numerator of NCCC. We have, g D = 0 ; therefore, Substituting the obtained solutions of all terms associated with NCCC in Eq. (1), and simplifying it further, results in the following formula as the final solution for computing NCCC, when both image and template are Gaussian.
Optimization of template radius. The present form of NCCC in Eq. (26) is suitable for analytical maximization with respect to the template radius, a . The voxel-wise maximizer, a , would allow the computation of NCCC in a single step, thus reducing the computational complexity to the linear time, given that the NCCC no longer needs to be computed for a set of template radii. Let us define α = a/σ ⇀ X . With this notation, Eq. (26) can be simplified as, (20) 2α(α 2 +1) , to zero, resulting in the following analytically optimized voxel-wise template radius, Local statistics of image volume. The local statistics, mean ⇀ µ ( ⇀ X) and SD σ ( ⇀ X) above, are the key parameters to compute the voxel-wise optimized radius. We measured these two parameters from Gaussian smoothed image to counter the small intensity variations in MRI. Gaussian smoothing was performed on the MRI volume with a 3D Gaussian filter of size W × W × W , where W = 2 b ′ + 1 and b ′ < b , and an SD of 2 voxels. We then computed the background minimum image, f BG , via moving minimum filtering over the smoothed image with a kernel of the same size to record the local minimum value inside the sliding window, which we then used in computing local statistics. The time complexity of this step is still kept O(N) by sequentially applying 1D filters in each dimension in both filtering steps.
Mean. We compute the local (moving) mean at each voxel in the smoothed image. The moving minimum is locally subtracted from the local image, in order to remove the effect of the background in the local neighborhood. We define the D-dimensional mean as described below.
where f s ⇀ X is the Gaussian-smoothed image. The second term of the numerator is zero due to symmetry. To simplify the calculation of the above expression, we make use of a linear change of variables. Assume, Standard deviation. SD is computed for each voxel by using a W × W × W kernel in Gaussian-smoothed background-subtracted image. From the definition of variance, the SD ( σ ) can be written as, The covariance of a vector is defined as a D × D matrix. However, given that the template is a radially symmetric Gaussian, we consider a scalar variance here as 1/D of the trace of the covariance matrix. We proceed by computing ⇀ X 2 as, , by a change of variables, Similarly to the case of the mean, we apply our approach to computation of convolution in the above equation (with n = 1) as described below, We compute the second term in the numerator of Eq. (33) the same way as in Eq. (35), by replacing f s Template matching with maximized radius. In this step, the mean and SD are substituted in Eq. (28) to obtain the " a"-maximized volume. Further, the optimized a * volume should be thresholded as a * ← min(a * , b/n − 1) before being substituted in the proposed solution of NCCC as described in Eq. (15). This is how the NCCC computational complexity can be reduced form O(a max N) to O(N) , i.e. we no longer need to use a set of templates with varying radii. In contrast, the FFT approach cannot make use of an optimal radius map to reduce its complexity any lower than O a max NlogN , given that FFT-based convolution cannot be performed with a locally varying kernel.

Results
In this section, we present experimental results comparing the proposed methods and other existing techniques, on 3D synthetic and real MRI data ( D = 3 ). We first describe the FFT-based NCCC computation, to which we will compare our methods next.
FFT-based template matching. We describe here how to use FFT to calculate the NCCC via Eq. (1). We define a truncated Gaussian kernel ( n → ∞ ) as, h (∞) The denominators can be similarly computed. In the FFT-based approach, NCCC needs to be computed for a set of template radii, resulting in the computational complexity of O a max NlogN . www.nature.com/scientificreports/ Experiment on synthetic data. We first evaluated the proposed methods on synthetic data and compared their performances with that of FFT. We compared the detection accuracy and runtime of the three methods. We created twenty artificial volumes of size 513 × 513 × 513 voxels, where each volume contained an enhancing sphere located at a random location with a random radius (varying in between 8 to 21 pixels) to generate the synthetic data. In our first experiment, we aimed to determine how accurately the enhancing spheres could be detected with their approximated radii from the respective volumes. During the experiments, we used the same input parameters for all methods such as b = 50 and a max = 24 (except for the O(N) method, instead b ′ = 24 ).
Note that n is ∞ for FFT and 2 for our methods. a max was calculated as a max = b/n − 1 . During this experiment, we observed that both the FFT-based and our O(a max N) template matching accurately identified the centers of the enhancing spheres in the respective volumes. For the proposed O(N) method, we calculated the analytically optimized a * for each voxel, using which we computed the NCCC in linear time over the entire image. We then simply looked for the voxel with maximal NCCC value to identify the center of the detected enhancing sphere, and the a * at that location was considered to be the detected radius. Since the SD of the Gaussian is proportional-but not necessarily equal-to the radius of the detected binary sphere, we conducted a linear regression analysis to reveal the ratio of the gold-standard binary radius to the detected Gaussian SD. From the experiment, we found their ratio to be 1.61 ± 0.05 and 1.74 ± 0.06 for the O(a max N) and the FFT methods, respectively, both values significant ( p = 9 × 10 −18 ) 12 .
As for the O(N) method, we first obtained the abovementioned ratio as 1.96 ± 0.11 ( p = 5 × 10 −16 ). The centers of the enhancing spheres, however, were not detected by the O(N) method with pinpoint accuracy; the mean distance between the gold-standard and detected centers was 1.6 ± 1.8 voxels. Consequently, we looked up the values of a * at the known centers of the spheres, and recomputed the gold-standard-to-detected ratio as 2.12 ± 0.11. Next, we repeated the O(N) experiment with the corrected a * ← 2.12 × a * . This helped us to reach absolute localization accuracy, with the new gold-standard-to-detected ratio of 1.01 ± 0.03 ( p = 8 × 10 −21 ). www.nature.com/scientificreports/ We performed another experiment on the synthetic data to analyze the runtime. The data were generated by creating 20 different volumes with linearly increasing number of voxels from N = 10 6 (100 × 100 × 100) to N = 1.25 × 10 8 (500 × 500 × 500) , each of which containing a homogeneous bright sphere (all ones) of radius 17 pixels at the center. We ran the scripts of the algorithms 10 times for each synthetic volume for sphere detection with the same parameters as above. We measured the mean and SD of the 10 runtimes and compared them among the different methods, as depicted in Fig. 2. We ran the experiments in MATLAB on a Linux compute cluster with each node including 8 processors with 56 GB shared virtual memory (cluster allots 1 CPU with 7 GB of virtual memory to a single submitted job). To run the scripts on single thread for an unbiased comparison, we used the singleCompThread option, but we found that MATLAB still sometimes multithreaded the codes on multiple cores. Therefore, we report both the wall time (the real-world time elapsed between the start and end of the code) and the CPU time (the total time spent by all cores to execute the script). The O(N) method (green curve in Fig. 2) was respectively 3.3 ± 0.7 and 1.4 ± 0.1 times faster than the FFT and O(a max N) methods in terms of CPU time, and 2.2 ± 0.6 and 1.2 ± 0.2 times faster in terms of wall time. Additionally, the O(a max N) approach was 2.5 ± 0.4 and 1.9 ± 0.4 times faster than the FFT method in terms of CPU and wall times, respectively. The non-monotonicity seen in the FFT runtime is partially because the radix-2 implementation is faster when each dimension is a power of 2. We also observed that the O(a max N) method required 1.7 ± 0.3 and 2.5 ± 0.3 times less resident memory then the FFT and O(N) methods, respectively, and that the FFT-based template matching used 1.5 ± 0.3 times less memory than the O(N) method.
Experiments on clinical MRI. We tested all three methods on in vivo MRI data of human brain for automatic detection of MS lesions. MS is a chronic demyelinating disease that generally affects the central nervous system, with lesions appearing in various regions of the brain, predominantly the white matter (WM) 4,28 . Given that the cerebrospinal fluid (CSF) is attenuated in FLAIR [though sometime lesions are iso-intense to the hyperintense gray matter (GM)], FLAIR MRI is suitable for understanding the lesion burden and measuring the amount of damage due to demyelination 3,6 . We tested the methods over the FLAIR sequence of two different MS lesion MRI databases; the first one is the training data of 'MS lesion segmentation challenge 2008' (MSLS 2008), which was conducted in conjunction with MICCAI 2008 (http:// www. ia. unc. edu/ MSseg), and the other dataset is the training data of the 'MS segmentation challenge using a data management and processing infrastructure' (MSSEG 2016), which was conducted in MICCAI 2016 (https:// portal. fli-iam. irisa. fr/ msseg-chall enge/ overv iew).

MRI scans of MS subjects.
We used 20 cases of training data of MSLS 2008 to evaluate the methods using the gold-standard annotation from the Boston Children's Hospital (BCH) rater as the reference. Isotropic-voxel (0.5 × 0.5 × 0.5 mm 3 ) FLAIR images of size 512 × 512 × 512 29 were preprocessed through N4 bias field correction (https:// www. slicer. org), brain extraction from T1 images using the volBrain tool (http:// volbr ain. upv. es), and then the resulting brain masks were multiplied with bias-corrected FLAIR images to remove the skulls.
The 15 cases of preprocessed (nonlocal means (NLM) filtering, rigid registration, skull stripping, and N4 bias field correction) MSSEG 2016 training data, as provided in the challenge website, were used in this study 6 . We sampled the 3D FLAIR images in the slice direction of each volume (with linear interpolation for FLAIR and nearest-neighbor for the gold-standard lesion and brain masks) to obtain isotropic-voxel (1 × 1 × 1 mm 3 ) data. After sampling, images were transformed into the new sizes, including 315 × 512 × 512, 154 × 224 × 256, and 245 × 336 × 336. The consensus segmentation obtained from all seven delineations (by seven expert raters) was used as the gold-standard for the performance evaluation of all methods implemented in this study 6 .
We min-max normalized the images of both databases before all experiments, and also kept the same input parameters during the experiments, such as b = 18 and a max = 8 for the O(a max N) and FFT-based methods, whereas the O(N) method dealt with all similar parameters, except a max ; instead we had b ′ = 8 . We chose n = ∞ for the FFT methods and n = 2 for the other two methods, but eventually following a similar approach to find a max in the FFT and O(a max N) algorithms.
Performance evaluation on clinical MRI data. The performance of all methods applied for lesion detection were quantified with respect to the expert-made manual annotation of MS lesions in MRI volumes by evaluating the popular statistical metrics such as the true positive fraction (TPF, or sensitivity) and false positive fraction (FPF, or false discovery rate) 9 as defined below, where L A and L M are the label volumes corresponding to automatic detection and manual annotation, respectively, and ¬ denotes negation. A TPF value close to 100% signifies that all voxels representing lesions in the annotated gold-standard are correctly detected, while the 0 value indicates absolutely no detection. Conversely, an FPF value of 100% means no detection, and a value of 0 means that all detected voxels were correctly inside the label. We now discuss how we used the above metrics to evaluate the performance of the O(a max N) , O(N) , and FFT approaches. We first computed the a-maximized NCCC from the FLAIR volume using each method. Many methods in the literature threshold the NCCC to keep putative lesion voxels, on which various performance  10,17,20 . In this work, instead of selecting a threshold, we considered N Top detected voxels with the highest NCCC value. We found the voxel with maximal NCCC from the resulting a-maximized NCCC volume, masked a sphere around it with a radius twice the optimal value of a , updated the NCCC volume by filling the sphere with the value of -∞, then repeated this process to find all N Top maxima. Next, as described in Algorithm 1, we created a sphere of all ones with the radius 1.61 times the optimally detected a , centered at the location of top detection, inside an all-zero volume, M , with the same size as the input image. We compared M with the gold-standard label to calculate the TPF and FPF for the case of N Top = 1 . We then added an enhancing sphere in M for the next detection to compute the TPF and FPF for N Top = 2 and so on.

Algorithm 1: Performance evaluation of algorithms for automatic lesion detection from brain images
Input: -Sets of coordinates, , with corresponding detected radii, , , having maximum c = 1,…,5000 NCCC recorded from top 5000 detection.
Create an enhancing sphere ( ) at with the radius c 1.61

3.
Measure the overlap between the automatic and manual label via the logical 'AND' operation; . ∧

4.
Calculate TPF and FPF for . = end for The mean TPF vs. FPF graphs over the 20 subjects of MSLS 2008 and the 15 subjects of MSSEG 2016 are presented in Fig. 3 for the three methods. We also recorded the computation times and resident memory usage of the three algorithms while running them over the real data as summarized in Table 1. The runtimes and required memory while executing the entire MATLAB scripts (termed as Algorithm) and those of the script till getting the a-maximized NCCC volume (termed as NCCC ), are shown in Table 1    From the results of time complexity experiment over synthetic data (see Fig. 2a,b), it is clear that O(N) method is the fastest of the three template-matching algorithms for sphere detection. We derived a new analytical solution to obtain the voxel-wise optimal template radius from local image statistics. The proposed analytical solution is a computational alternative to exhaustive search for the template size. The O(N) method could detect the bright spheres at their exact centers in the synthetic data (after using the appropriate multiplier). On the other hand, the FFT and O(a max N) based template matching with exhaustive search possess similar accuracy for detection of the bright spheres as expected. The only small deviation persisting between them for the prediction of the radius of the sphere is due to employment of approximated Gaussian template in the O(a max N) technique. The results on synthetic data validate the theoretical prediction of achieving the linear complexity in template matching by our methods (Fig. 2a,b). Regarding the resident memory usage, the O(a max N) method appears to be quite efficient holding less memory than the other two methods. The O(N) method, however, uses high memory. The analytically optimized template radius can be computed from mean and SD, which are computed from the local 3D neighborhood. The O(N) method, therefore, requires storing auxiliary variables (such as mean and variance volumes) in the memory, thereby using more memory than the other two methods. In the exhaustive-search cases, we performed loop operation inside the MATLAB script to execute the similarity measure computation step for a set of template radii. We reduced the memory requirement by updating NCCC with the voxel-wise maximized values, thereby avoiding the storage of a separate NCCC volume for each tried radius. Figure 2c shows that that the memory requirement increases linearly with the number of voxels in the volumes. In case of real data, the results follow a similar trend. The FFT method requires more memory than the O(a max N) method does, because during template matching with the set of varying radii (loop operation) in the FFT method, we need to create a convolution kernel (template) that is as large as the input image, for every tried radius (so they can be multiplied elementwise in the frequency domain).
Regarding the runtimes of the methods on real MR images, Table 1 shows that the CPU time of the O(N) algorithm is higher than that of the O(a max N) method, and close (till NCCC computation) to that of the FFT-based   N) is the fastest of the three on the MSLS 2008 data. These results, which contradict those on the synthetic data, are because NCCC was computed for a set of only 8 different radii for both the O(a max N) and FFT methods and this loop-based operation was proven to be faster than the computation of 10 more convolutions for optimizing the template radius in the O(N) algorithm. In contrast, for a large set of template radii as considered in the experiments on synthetic data, the O(N) method outperformed the others. The CPU times of the MSSEG 2016 experiments were very close for the three algorithms, since for small images, the gap between NlogN and N becomes narrow. Overall, the O(N) method is expected to run faster for large MRI volumes, or for the detection of larger lesions. Ex vivo imaging of human brain, which has attracted researchers to explore microstructural neuroanatomy, can produce brain images with (0.1 mm) 3 isometric voxel resolution 30 . As the resolution increases, the ROI volume size increases proportionally, necessitating the reduction in computational complexity of the algorithms for processing very large MRI volumes to automatically detect lesions. Note that in real clinical scenario, we do not need to test the detection for top 5000 lesions, instead it would be for fewer numbers depending on the characteristics of the lesions and the clinician's needs. Also note that the computation of the a-maximized NCCC can be seen to take up only a small fraction of the entire runtime (see Table 1), and has more variability.
Prior studies on MS lesions, as discussed in the next section, however, generally segment the lesions and report the detection accuracy from their automatic delineation. Therefore, our approach of lesion detection is different from most available tools. During the performance evaluation, we observed the monotonically increasing characteristics of both TPF and FPF with N Top for all three methods (except LPA). Thus, a high sensitivity can be achieved with large N Top , which can be interpreted from Fig. 3. The mean TPF vs. FPF plots (Fig. 3) show very similar performance by FFT and our O(a max N) methods over both MRI databases, similar to synthetic data. We can also infer from Fig. 3 that the performance graph of the O(N) method is slightly better than FFT and O(a max N) , making this method a coherent assistive tool. In addition, the practicality of linear-time template matching can be well appreciated in processing of large 3D volumes with a high b (e.g., b = 50 showed the efficiency of O(N) in synthetic data). Throughout the experiment with the O(a max N) and O(N) methods, we used an approximated 3D Gaussian template, which is almost radially symmetric, and expected to work robustly with spherical or close-to-spherical lesions. On the contrary, the MS lesions are arbitrary in shape and size (largely varying) 31 . Therefore, the algorithms might not detect the exact centers of such arbitrarily shaped lesions in FLAIR images. This reduced the overlap between the detected region and the gold standard, even when the lesions were correctly localized. This may be one of the reasons behind the rise of false positives. In the context of having large false positives, we would like to mention the work of Cabezas et al. 9 . Their unsupervised methodology deeply depends on tissue priors to estimate the mean and SD of the GM distribution in FLAIR images, for defining a threshold in order to segment the MS lesions. They reported a mean FPF higher than 90% (with close to 40% of mean TPF), as accuracy of lesion detection, without using any false detection refinement step. So, the presence of high FPF seems to be common among the performance of the algorithms, especially targeted to lesion detection from MS subjects, which we will discuss more in next section. Contrary to the unsupervised approaches, the prior knowledge learnt by LPA during the training phase seems to have been effective in improving the performance of this supervised-learning method. This analysis helps in tallying the performance of supervised and unsupervised techniques (Fig. 3). In template-matching techniques, a general trend is to fix a threshold over the spectrum of NCCC values so as to evaluate the performance. However, the optimal NCCC threshold may depend on acquisition parameters, and fixing it might lead to the detection performance strongly depending on the threshold. In contrast, our strategy of using maximum NCCC from N Top detection is not threshold-dependent.

Visualization of the detected lesions.
We display the results for top 30 out of 5000 detections in Figs. 4 and 5. Blue, red, and green circles represent the FFT, O(a max N) and O(N) methods, respectively, on the 2D axial slices of the FLAIR volumes. We present the slices in order of lesion detection. (There may be slice repetition due to multiple detections in the same slice.) Thick circles indicate that the centers of the detected spheres are located in the shown slice. Thin circles represent the intersections between the shown slice and the spheres representing top-30 detected lesions centered in nearby slices. Therefore, the lack of a circle around a lesion means that it has not been in the top 30 (but the algorithm may still have detected it). A limitation that we observed for all techniques is the multiple detection of parts of a large lesion, such as those scattered around the lateral ventricles or on the corpus callosum, which are common in MS 31 . When lesions of various sizes and shapes appeared mainly in the WM 31 , multiple lesions were detected over the same large lesion because of the small b that we used. This issue may not be expected to arise in metastatic-lesion detection. The selection of a large b, on the other hand, may cause the algorithm to miss smaller lesions. We can observe the presence of false lesions (cyan circles in Figs. 4 and 5) included in the top 30 detection for all methods. False positives are mostly in the GM, which is hyper-intense in FLAIR images and sometimes confused with lesion by the algorithms. This is the main cause of high FPF in template-matching algorithms. We also found false positives in other brain regions such as the splenium of corpus callosum, hippocampus, and medulla oblongata. The estimated lesion radius in the O(N) method is comparatively less than FFT and O(a max N) methods, which is expected from the gold-standard-todetected ratio radius (1.96 ± 0.11) obtained on synthetic data experiment.
Comparison with existing literature. In this section, we first compare the performance among the existing state-of-the-art methods dedicated for MS lesion segmentation with MSLS 2008 data. We have considered some of the well performing techniques of the last decade (summarized in Tables 2 and 3). Challenge organizers evaluated true positive rate (TPR) and false positive rate (FPR) (with the same definition as TPF and FPF used here) to determine the lesion detection ability 29 Table 3    www.nature.com/scientificreports/ In this study, we intended to develop an algorithm for fast and automatic detection of brain lesions, as opposed to determining lesion contours (but still estimating the radius of the detected lesions). In most of the existing works, multi-contrast MR images were considered to delineate lesions 6 , with the exception of some methods that used only FLAIR images (as we did). Using multiple MRI sequences is attractive as it provides a larger feature space than using a single sequence, therefore multi-channel MR images are widely used in machine leaning and deep neural networks. In this context, our template matching algorithms, focusing on estimating lesion center with an approximated radius, are fundamentally different from the segmentation approaches. During literature review, we observed that Weiss et al. and Manjon et al. avoided computing the whole image to possibly reduce the computational burden, instead they downsampled the image 38 or used candidate ROIs 36 . We found that MS lesions were detected as regions with outlier intensity values compared to the healthy brain tissue in many works 40,41 . Some research groups designed two-stage segmentation tasks; they first determine the putative lesion voxels with the help of various techniques such as unsupervised mixture models 42 , estimation of lesion and tissue label using Markov random field (MRF) 43 , cascaded convolutional neural network (CNN) 44 , etc., and then segment the lesions by classifying the voxels. In another two-step approach, brain tissues, mainly WM, GM, and CSF, were segmented using different techniques such as atlas-based EM 37 , Gaussian mixture model (GMM) 45 , SPM8/12 segmentation algorithm 46 , etc., before delineating the lesion while making use of the segmented tissues. These multistage approaches generally require multi-contrast MRI. Our proposed techniques compute the image at the original resolution and no spatial prior information on the tissue or lesion or candidate ROIs needs to be extracted; instead, it can directly localize the lesion region based on the NCCC measures, making it a single-stage method. Another key point we observed in the existing literature is the inclusion of post-processing operations, which is a common practice to refine the results and reduce false positives. The presence of flow artifacts in the FLAIR sequence creates hyperintensity in the CSF space, thus sometimes mimicking the lesion 53,56 . In addition, the hyperintense border voxels at the GM-WM junction may result in some false positives 53,57 . As a remedy, various post-processing operators have been adopted, such as morphological operations 37 , setting a minimum lesion volume 46,53,58 , connectivity rules, neighborhood similarity 46 , etc., to reduce false detections and improve FPR and PPV. In this study, we did not perform any post-processing, and only evaluated the raw outcome of the proposed algorithms. The performance of various supervised [34][35][36] (Sup. Seg.) and unsupervised 37-39 (UnSup. Seg.) segmentation techniques, evaluated on MSLS 2008 training data, are shown in Table 2. We also report the results of the proposed methods and the FFT baseline approach for four different values of TPF. One can see that O(a max N) method's performance is similar to the FFT approach, consistent with the synthetic-data experiments. The O(N) method seems to perform slightly better than the other two template matching methods. But all three methods yield poor PPV (due to high FPF). Without the post-processing step for MS lesion detection, unsupervised detection has previously shown to result in high FPF 9 . Besides the fact that a single image (as opposed to three varying-contrast images) was used here, overlapping lesion intensity, flow artifacts, and typical shape and size of MS lesions are likely causes of high FPF, when lesions are detected from FLAIR images without refinement strategies. Having no user tunable learning parameters, as opposed to (supervised or unsupervised) segmentation approaches, the tested unsupervised detection approaches are not meant to be directly compared to the segmentation approaches presented in Tables 2 and 4 (discussed below), especially given that a detection method is designed to locate-but not delineate-a lesion. Nonetheless, to paint a broad picture of the performances and applicability of existing works, we present all methods in the tables side by side. We also observe from Table 2 that both TPR and PPV are low for all methods, which emphasizes the presence of high false positive with poor lesion detection ability. Table 3 Tables 2 and 3 were retrieved from the literature and leaderboard results (http:// www. ia. unc. edu/ MSseg/ re-sults_ table. php), and show that high FPR is very common among the methods even when the computation is done with multi-contrast MRI. Ghribi et al. 45 pointed out the poor quality of MR images in MSLS 2008 as the reason behind the high FPR of their method. As expected, supervised learning methods performed well in both tables compared to unsupervised techniques.
The performance summary of some of the existing works on 15 training images of MSSEG 2016 is presented in Table 4. The table includes both supervised 44,59-69 and unsupervised 58,70,71 methods. In the MSSEG 2016 challenge, the lesion detection capability was evaluated using various performance metrics, among which we have chosen to compare the Dice coefficient (DICE), TPR, and PPV 6 in Table 4. From the table, we can see that different groups explored both multi-channel and single FLAIR MRI for processing. Beaumont et al. 58 did not consider lesions appearing at the edges of the brain mask, not dominantly located in the WM, etc., for final refinement of lesion boundaries. In their following work, Beaumont et al. 70 reported augmented accuracy after implementing postprocessing technique. Knight and Khademi 71 showed 3D Gaussian filtering with 0.5 SD to be efficient for noise removal in FLAIR. In this work, we also used the Gaussian filter, however with a higher SD (2.0) in computing local image statistics. They also used three false-positive reduction strategies based on lesion volume, distance from the edge of the brain, and distance from midline of the brain. In other works 59,60 , handcrafted features and post-processing were used for MS lesion segmentation. Most groups favored working with multi-contrast MRI 6 for either designing machine learning algorithms or obtaining tissue prior, and followed with rule-based or morphological-operator-driven post-processing. Table 4 also includes the performances of the proposed and the baseline FFT detection approaches, similarly to MSLS 2008. The results on MSSEG 2016 show the same trend, i.e. the performance of the O(a max N) and FFT methods are very close, and the O(N) method achieves slightly better accuracy. In general, the sensitivity (TPR) of lesion detection seems moderate with a high false detection in both databases. The lower DICE in the three detection methods that we tested (compared to segmentation methods in the literature) is mainly attributed to the facts that: (1) we used a single image (as opposed to multiple images with different contrasts), (2) we performed no post-processing, (3) we employed no supervision, but most importantly (4) our goal was to detect (find the rough location of) as opposed to segment (delineate) the lesion. www.nature.com/scientificreports/ As such, the segmentation performance metrics used for the three detection methods should be compared among themselves, but not directly to those of the segmentation methods. The organizers of MSSEG 2016, Commowick et al. 6 , summarized that the deep-learning based approaches showed promise in narrowing the gap between manual and automatic segmentation of MS lesions. As our proposed method is a fully automatic unsupervised detection technique, the comparison with machine learning (deep neural networks) is not so informative due to the fundamental differences with supervised methods. In summary, the proposed template matching algorithms do not require tissue prior or any lesion-specific information. They are fully automatic and unsupervised models that can detect lesions from the FLAIR sequence with a lower computational complexity than the conventional FFT-based template matching. We also did not apply any post-refinement approaches to remove the false positives; however, false lesions can be ruled out www.nature.com/scientificreports/ during the review by the clinician. In lesion detection, we give a higher priority to avoiding false negatives than false positives.

Conclusions
We have proposed a new mathematical framework to compute a radius-optimized normalized cross-correlation coefficient (NCCC) similarity measure for 3D template matching, with application for automatic lesion detection. The key feature of our method is its low linear-time computational complexity, allowing for fast detection.
Our alternative way to compute the NCCC exploits the approximation of Gaussian smoothing with multiple convolutions with a box kernel. We further proposed an analytical solution for the template radius that replaces the costly exhaustive search, thereby achieving the computational complexity of O(N) , which is not possible with the FFT-based method. The proposed technique provides a new approach to implementing fast template matching for object detection from images. We tested our lesion-detection method on synthetic images as well as two real MRI databases of multiple sclerosis lesions, and compared its performance with existing state-of-theart methods and the conventional FFT-based template matching algorithm. Future work includes extending our approach for anisotropic-voxel images, and also imaging modalities beyond MRI.

Data availability
The proposed methods were implemented in the MATLAB platform.