Medical Image Segmentation Using Multifractal Analysis

Image segmentation plays a crucial role in image analysis processes. The operations performed on a segmented image tend to affect it differently than if they were performed on the original image; therefore, segmenting an image can show radically different results from the original image and successfully do so can yield features and other important information about the image. Proper image analysis is of high importance to the medical community as accurately classifying different conditions and diseases can be facilitated with excellent patient imaging. Multifractal analysis can be leveraged for performing texture classification and image segmentation. In this paper, we propose fusion-based algorithms utilizing multifractal analysis for medical image segmentation. We use two specific multifractal masks: square and quincunx. Our techniques show new insights by using methods such as histogram decomposition in conjunction with new techniques, such as fusion. By fusing different slope images, we can extract more features, thus making our proposed algorithms more robust and accurate than traditional multifractal analysis techniques. These methods are further capable of reliably segmenting medical images by implementing multifractal analysis techniques in coordination with methods such as gaussian blurring and morphological operations. Medical professionals can easily analyze the resulting image for diagnosing medical conditions. The outcomes show that the proposed algorithms extract dominant features that are more encompassing and powerful than classical techniques. Keywords— image segmentation; multifractal; MRI; fusion; quincunx; square mask.


I. INTRODUCTION
The human nervous system is extremely well-adapted to processing visual information quickly and effectively, but there are limitations that technology can address. Image analysis is often employed in situations where human analysis is inhibited, either due to the quantity of incoming information, such as a data-rich environment, or the quality of the information in scenarios where human senses cannot observe the subject matter, such as low-light or radiation. Computers have the advantage of exponentially faster processing power, the ability to add new sensors, and the ability to be placed in dangerous situations unfit for humans. However, they have the disadvantage that they are not naturally adapted to processing image data and must be manually taught how to do so, a task addressed by the field of image analysis. Despite this, after being trained properly, the primary advantages include detecting features normally indistinguishable to the human eye as well as analyzing images at a rate infeasible to humans. In most cases, the aim is either to let the computer process data, which are too numerous to be analyzed by human beings or to be able to detect features in images that are not easily found by the human eye [1].
Magnetic Resonance Imaging, commonly referred to as MRI, is a non-invasive procedure for scanning and imaging the human body. Compared to X-ray Computerized Axial Topography (CAT) images, or Computerized Topography (CT) for short, MRI images provide more detail without the use of potentially dangerous radioactive waves [2,3]. An MRI scan provides high contrast, three-dimensional images that are invaluable for medical diagnosis and treatment [4][5][6]. MRI is based on many scientists' research, including but not limited to E. Purcell, F. Bloch, P. Mansfield, and P. Grannell's work with nuclear magnetic resonance (NMR) [7][8][9]. The concept behind NMR is that materials containing nuclei composed of a non-even number of neutrons, protons, or a mixture of the two have both a nuclear "spin" and a "magnetic moment." Many materials, including biological tissue, have nuclei with the properties mentioned above, and as a result, it is possible to create images of these materials by means of NMR methods [10]. The atoms of these materials emit radio signals when exposed to Radio Frequency (RF) fields. MRI takes advantage of this fact by using an RF field to stimulate these atoms, causing them to emit signals. These signals generate an electric current that is detected by a pickup coil that surrounds the patient or subject to be imaged.
A computer processes the observed signal to create a sampled grayscale pattern, which constitutes an MRI image. Work done by R. Damadian and P. Lauterbur during the 1970s displayed that NMR methods can have applications in medical diagnosis [11,12]. Damadian found that the relaxation time, the period of time an atom continues to emit signals after exposure to an RF field ends, differs between tissue types. The variation in relaxation time is what allows tissues to be distinguished from surrounding tissues in an MRI image. Detectable relaxation times are divided into two types: T1 and T2. A signal received from a particular tissue is typically composed of both T1 and T2 values. This combination determines the characteristics of the signal. For images that are T1 based, tissues or elements of the generated image that have low T1 values will appear as bright spots in the image while those having high T1 values are shown to be darker. Conversely, T2 based images show high valued T2 tissues as being dark and low valued T2 tissues as bright. Traditional techniques cannot effectively process the volume of data that MRIs are capable of producing. Therefore, medical image analysis methods and techniques for visualization are of great value in the medical imaging field. There are three major topics of research in this area: cross-registration, intuitive visualization, and image segmentation. The goal of the latter most of these research topics is to efficiently identify important structural information concerning the subject's pathology and anatomy. Image segmentation is often done manually and thus creates a bottleneck in clinical applications [13,14]. This is standard but unacceptable in situations where it is crucial to identify many organs within the radiological data sets, such as computer assisted neurosurgery, where using CT or MRI images is the current standard [15]. Other situations require the identification of tissue boundaries, especially those in which the relationship between therapeutic actions and morphological changes must be evaluated and understood. Obtaining statistically significant results demands that many data sets must be segmented.
Image segmentation is typically categorized into three general categories: region extraction, labeled characteristic thresholding, and edge detection [16]. Further research was published in 1985 by R. Haralic and L. Shaprio, in which they characterized image segmentation into the following groups: measurement space guided spatial clustering, single linkage region growing scheme, hybrid linkage region growing scheme, centroid linkage region growing scheme, and split and merge schemes. They showed that clustering and segmentation were different because clustering consisted of grouping in measurement space, whereas segmentation consisted of grouping in the spatial domain of the image [17]. In 1988, P. Sahoo, Soltani, and A. Wong updated the work of Fu and Mui by presenting a survey that showcased the performance of segmentation algorithms using criterion functions such as uniformity and shape measures [16,18]. While the aforementioned research has offered great insights, it did not consider range or MRI image segmentation [19]. M. Vannier et al. helped to address the lack of research on MRI tissue segmentation with techniques utilizing statistical classification applied to the intensity of the signal [20,21]. Today, numerous segmentation algorithms are available, yet, despite the many options, there is not one specific technique that is regarded as the best for every image and not every segmentation process can be successfully applied across fields. If image segmentation is performed correctly, the new image will be less complicated and easier to analyze. This is because the objective of image segmentation, in relation to computer vision, is to compartmentalize various subsections of an image. Segmentation facilitates examining an image for essential details, such as tissue identification [19,22,23]. Brian tumor segmentation has been an intense subfield of research. The fusion of multifractional dimension (multi-FD) and fractal and intensity features can greatly increase the accuracy of brain tumor segmentation [24,25]. Over the past few decades, multifractal analysis has been increasingly used in medical signal analysis [26,27]. Multifractal analysis has proven effective for a wide range of multi-modal medical images such as electrocardiogram signals and mammography imaging [28,29]. This analysis is quite effective for image segmentation, which characterizes a given region of the image. One example involves characterizing the pixels distribution heterogeneity of a region of the image. Multifractal models allow us to describe the scale-to-scale propagation of this distribution [30,31]. Previous research has shown promising results for the analysis of mammography images using this technique. Fractal geometry is additionally applicable due to the nature of micro-calcifications to look like light clusters of spots that vary in size and shape. These clusters are embedded within inhomogeneous tissue in the background. The inhomogeneous background also shows the self-similarity of fractal images in that a region containing the microcalcifications clusters can be viewed [34,35]. Multifractal analysis with fusion has shown to increase accuracy and decrease false positives in image analysis [36].
The technique of expanding clinical relevance for diagnosis of medial illnesses by blending numerous images from single or multiple sources is called medical image fusion. A. James and B. Dasarathy recently classified fusion research into three important areas: modality of images, image fusion, and organ imaging. They believe that fusing medical images has large potential to grow and has already proved important in various fields [37]. Fusing images tends to be required in fields such as medical imaging and computer vision where multiple instruments or devices are used to capture data [38]. These fused images are predominantly intended to facilitate human interpretation and observation [39]. We present a multifractal-fusion-based image segmentation algorithm that extracts important features. An overview of this algorithm is shown in Figure 1. Image fusion is a technique of taking an arbitrary amount of images and merging the pertinent information of each image into a single image. The final image ideally contains the most significant features of the individual images; therefore, the single image will be more complete and helpful to scientists or doctors [38]. There are multiple methods for performing image fusion ranging from relatively simple techniques such as high-pass filtering to more complex fusion methods such as Laplacian or Gradient pyramid [40]. While there could be numerous goals in image fusion, there tend to be two important objectives to accomplish: the reduction of noise and decreasing entropy. These objectives are also essential to fusing methods. These techniques can be accomplished by not introducing abnormalities that could divert the attention of the human analyzer to inessential parts of the image. New deformities could lead to inaccurate interpretations of data ranging from test-tube babies to classifying blood cell neutrophils to diagnosing brain tumors [41][42][43]. The final image should have as much of the pertinent features of the original images as possible [38,44,45].
This remainder of the paper is organized as follows: Section 2 presents an overview of multifractal analysis with a derivation of the Hausdorff spectrum for signals and images, including overviews of thresholding and morphological operations. Section 3 provides the empirical framework for our analysis and the underlying methodology, which includes mathematical, visual, and textual information about the quincunx and square masks. That section also details our general algorithm for segmenting an image. In Section 4, our results are displayed and discussed, including sample slope images and selected images that underwent our entire methodology. We compare and analyze each mask to each other and discuss their strengths and weaknesses. The implications and possible applications of our work are then discussed, specifically in the field of medical imaging as well as future avenues of research.

II. MATERIAL AND METHOD
The utilized material begins with a concise overview of multifractals that dives into Choquet's capacity, Hölders exponent, and the Hausdorff spectrum in section A. Section B: Thresholding details the cumulative probability function and various probability mass functions. Section C: Morphological Operations defines equations for erosion, dilation, opening, and closing. Finally, Section D: Rescaled Range Analysis details the separation of significant categories of images identifies features that are not by other features.

A. Multifractal Analysis
The primary differences between multifractal techniques and classical techniques such as morphological approaches or Canny edge detection is in how local irregularities are managed. Vehel et al.'s [4][5][6] definitions provide the background for the multifractal techniques used in this study. Definition: -Let be a set. Paving on E is a set of subset of containing the empty set and stable under finite union and finite intersection. The pair (E, ε) is called a paved space. Let P(E) denote the power set of E.
which is defined when c (I (x))μ(I (x)) ≠ 0 and when this limit exists. This quantity is known as the pointwise Hölders exponent of c at point x with respect toμ, however the standard definition regarding the limit centered at K, is also given below.
Let μ be a Borel measure defined on a compact setP. For each pointx in P, define the local singularity coefficients as: where B M (x) is an open-ball of diameter N centered at the pointKand when the limit exists. Often α(x) is called the Hölders coefficient. It reflects the local behavior of the measure μ around x. Points bearing the same coefficients can be grouped into sets, named Iso-local singularity sets, define as follows: We can define above sets with threshold value as follows: To characterize the above sets, now we will define set dimensions known as the Hausdorff dimension: Where {E c } +c+∞ is a δ-cov e r of E Such that Finally define f(α) = dim l E (α) The description (α, f m (α)) is called the local singularity spectrum, also known as Hölders or Hausdorff spectrum of the multifractal measureμ. Approaching multifractal analysis from a computational standpoint requires the use of discrete forms of the aforementioned measures and capacities [46][47][48]. Because the Hausdorff is not computable, we use the computed box dimension instead. Additionally, we replace the lim α (x) by the slope of the linear regression of loglog of d against μ [49]. Then the abstract descriptions of multi-fractal analysis are useful for creating feasible computational tools for applications in images and signals. In equation 1, the points correspond with pixels that are from the actual image, openballs are associated with blocks (windows) centered about each pixel, and measures are related to functions of the intensities of gray levels. Wherep is the sum of the pixel intensities(i, j) inside of a region which is centered about pixel(x, y). Then h(x, y) represents the amount of gray at (x, y), μ is defined by the equation: The above equation 11 represents an image's intensity properties such as sharpness, spatial distribution, etc. We denote this as a sum measure of a given image.

B. Thresholding
The primary purpose of thresholding images is to extract significant features from background noise. Several methods have been detailed by Sezgin & Sankur that have been applied to image segmentation, such as histogram shapebased, convex hull, or peak and valley thresholding techniques. The definitions in this section are from Sezgin, Sankur, and D. Carabias [50][51]. The cumulative probability function will be defined as follows: where b is the threshold, the probability of any pixel with a gray level p(i) and |(}) is the probability of gray level less than or equal to the threshold. Equations

C. Morphological Operations
In this section, we describe important morphological operations, which are described as the mechanisms for selecting relevant skeletons, boundaries, convex hulls, or other pertinent and important regions shape characteristics of an image. Morphological algorithms for pre-and postprocessing also are essential to image analysis. Erosion and dilation are two elementary but crucial methods of morphological operations. Erosion lessens an image's attribute's while dilation bolsters them. These are shown in equations 19 and 20 where A, B ⊆ Z ‚ with A ⊗ B and A ⊕ B denoting erosion and dilation, respectively. The following definitions in this section are from R. Gonzalez, R. Woods, and J. Angulo, J. Serra [2,52].
Equation 19 and 20 form the basis for many more complex operations, such as opening and closing. The opening usually smooths out the contours of features, including removing small connections between objects and deleting short projections. The closing also smooths out the contours of features, but instead of eliminating small connections, closing strengthens them along with filling in small breaks or cracks. Opening and closing are defined in equations 21 and 22, respectively.

D. Rescaled Range Analysis
An approach to the characterization of texture data that has seen use in a variety of fields is through Fractal dimension (FD) [54][55][56][57]. This allows for the separation of significant image classes and attempts to characterize information that has not already been characterized by the other features. There have been a number of techniques proposed to calculate the FD of images, referred to as characterization, including Cover Dimension (CD), Box-Counting Dimension (BCD), Hausdorff-Besicovitch Dimension (HD), Wavelet-based FD, etc. We will refer to these procedures as FD-estimators. This section explores using the Rescaled Range analysis to calculate the Hurst exponent [58].
Various signals such as seismic activity, stock prices, respiratory rates, etc. can have the appearance of being random, but can be influenced externally. The nature of the external influences, however, can be random and appear as noise. Hurst defined an empirical descriptor of temporal signals describing natural phenomena. The Hurst exponent OEis a measurement of how smooth a fractal object is, with 0 < H < 1: low H values indicate higher degrees of roughness, almost to the point of filling the next-highest dimension, whereas higher values of H suggests a smoothness such that the next-higher dimension is minimally intruded on. Generally, the relationship is defined as H = E + 1-FD, where E = 0 for a point, 1 for a line, and 2 for a surface. Rescaled Range analysis [59,60] is a simple process that is highly data-intensive. Sequentially, the steps are as follows: 1. Start with the whole observed data set that covers and calculate the mean Ā = (1/N) ∑ a c 2. Next, sum the differences from the mean to get the cumulative total at each time point X 5y from the beginning of the period up to any time: X 5t = ∑ a c -Ā , k = 1,2,3, … n 3. Calculate the range R (τ) = max (X 5y )-min (X 5y ) for k = 1, 2, 3 … n

Calculate R/S = R (τ)/S (τ)
6. For the next stage, partition the time interval in to two blocks of size N/2 = τ and repeat the entire procedure , steps 1-5, and determined R/S for each segment of the data set of length N/2, then take averaged value. Repeat, using successively shorter τ's at each stage dividing the data set into nonoverlapping segments and finding the mean R/S of these segments Plot the log-log plot, that is fit Linear Regression Y on X where Y = log (R/S)and X = log N. The exponent H is the slope of the regression line.

E. Methodology
Below, we detail the hierarchy of our multifractal image analysis. The local Hurst exponents for the image are first calculated, which is then used to compute the local FDs. The smaller the FD, the finer the texture, and the larger the FD, the coarser the texture. Depending on the texture of the image, the image can be thresholded in order to detect edges and segment features in the image. One main problem of multifractal analysis is that it can be difficult to find the local FD, and in order to use multifractal techniques, we must first select a mask. This is also not a straightforward task. We describe the two masks we chose later, but for any mask we segment the image in a discrete number of parts, that use several FD parameters such as Hausdorff, box counting, and Hurst exponent. To find the local FD we use a mask at each pixel. We have created two algorithms using Matlab that apply different filters to medical images. Each filter has seven different mathematical processes of creating the new image that are maximum, mean, median, minimum, range and standard deviation [53].
In the first filter named quincunx, we have 50 images outputted for each mathematical process for a total of 350 images per input image. We represent nodes by a distancevector whose length is equal to the number of different distances from its origin. Figure 2 is a pictorial representation of a size 4 quincunx neighborhood with one decimal place approximation and Figure 3 shows a neighborhood of size 5 [53]. As shown in Figures 4 and 5, the square-mask filter has eight images outputted for each mathematical process for a total of 56 images per input image. The number of output pictures indicates the depth level of the mask for each mathematical process. For example, the quincunx filter has a mask depth of 50, and the square mask has a depth of 8. Our method also incorporates a simple form of R/S analysis. Despite R/S analysis typically being defined for one-dimensional time series, it can be adapted for twodimensional images. Our initial findings found that adaptive smoothing improves the quality of the resulting fused images, with the results of the smoothing techniques shown in Figure  6. After computing the slope images, we fused many of the images together, utilizing Matlab. Next, we convert the image into a binary image of white or black pixels. Then we use the morphological operation opening, as described in Section 2.3, to remove extraneous small clusters of pixels from each image. After that, we blur the image and use histogram-shape based thresholding to smooth out the ragged edges and curves. Finally, we calculate the entropy of the selected images that went through this methodology. Figure 7 provides a flow chart for the steps to this methodology. In Figures 8-12, take, for example, Max8Mean8. The "8" refers to the mask depth, and besides max and mean are the two methods fused that produced the respective image. We found that shallower mask depths did not provide enough detail to be useful for segmentation. Lack of detail can even be seen on layer 8, which is the deepest layer of the square mask in images such as Mean8Min8 in Figure 8 and the Median8Min8 of the quincunx mask shown in Figure 9.

Max8Mean8 Mean8Median8
Median8Min8 Min8Range8 Max8Median8 Mean8Min8 Median8Range8 Min8Stdv8 Many fused masks combine important features from each method, such as isolated edges and defined curves. Visual observation of the results in Figures 9 and 10 showcases the results of fusing two different methods in our two different masks. It makes logical sense that fusing maximum and mean would show darker pixels and more well-defined edges than the fusion of median and minimum. It is interesting to note that fusing the standard deviation in different masks and different methods can produce radically altered results. Take, for example, the Min8Stdv8 of Figure  8 that was produced with the Square mask. This figure shows nearly no discernable details besides the outside edges of the skull, and even that is just two different shades of gray. Then compare that to Range8Stdv8 of Figure 9 produced from the Quincunx mask. This image showcases intense detail of all the edges and brings out the minute structures compared to the rest of the fused methods. Although this is visually appealing, it does not help to reduce the clutter of extraneous details like possible tiny speckles of noise in between the fluid and tissue in the brain. We're striving to focus on major features such as the outline of the skull and brain. Because of this, we utilized the algorithm described in Section 3 and applied it to some selected images as detailed after Figure 10.

Max8Mean8 Mean8Median8
Median8Min8 Min8Range8 Range8Stdv8 Max8Median8 Mean8Min8 Min8Stdv8 The major differences we found were in what details our methods extracted. Even similar-looking fused images such as Max8Median8 and Max8Range8 in Figures 10 and 11 have different final results using our same methodology. To complete our methods, we compute the entropy of all three images in Figure 13: These graphs show that our algorithms perform similarly across multiple masks spanning multiple fused methods. The spike of the blurred image would occur due to increasing the amount of noise present in the image, whereas the rest of our approaches reduce entropy as expected. Figure 13 proves that our techniques outlined in this paper reduce the amount of excess information and noise present.

IV. CONCLUSION
Our methods show strong abilities to segment and reduce the level of excessive information in an MRI image. The fusion of many different methods has shown to emphasize features from the two different sources, making combining images ideal for better segmentation. Our method reduces noise, which is ideal for facilitating interpretation by medical professionals. Both masks extracted important information about the image; however, the quincunx mask extracted more than the square mask, which can be seen in the larger number of dark pixels in Figure 10. Although this is an important distinction, the goal of image segmentation is to extract the object from the image and reduce the excess information in an image, which we accomplished. Further research can be done on different kinds of medical images, and more information should be collected on other biological systems such as muscular or skeletal so that more exact segmentation can be conducted. Different kinds of masks, depths, and fusion of any combination thereof can also be researched in addition to applying our techniques to other fields of study such as object detection or biometric scanning such as facial recognition.