extfieldsofstudy:
Computer Science
Mathematics

text:
Homotopy Continuation Approaches for Robust SV Classification and Regression

In support vector machine (SVM) applications with unreliable data that contains a portion of outliers, non-robustness of SVMs often causes considerable performance deterioration. Although many approaches for improving the robustness of SVMs have been studied, two major challenges remain in robust SVM learning. First, robust learning algorithms are essentially formulated as non-convex optimization problems. It is thus important to develop a non-convex optimization method for robust SVM that can find a good local optimal solution. The second practical issue is how one can tune the hyperparameter that controls the balance between robustness and efficiency. Unfortunately, due to the non-convexity, robust SVM solutions with slightly different hyper-parameter values can be significantly different, which makes model selection highly unstable. In this paper, we address these two issues simultaneously by introducing a novel homotopy approach to non-convex robust SVM learning. Our basic idea is to introduce parametrized formulations of robust SVM which bridge the standard SVM and fully robust SVM via the parameter that represents the influence of outliers. We characterize the necessary and sufficient conditions of the local optimal solutions of robust SVM, and develop an algorithm that can trace a path of local optimal solutions when the influence of outliers is gradually decreased. An advantage of our homotopy approach is that it can be interpreted as simulated annealing, a common approach for finding a good local optimal solution in non-convex optimization problems. In addition, our homotopy method allows stable and efficient model selection based on the path of local optimal solutions. Empirical performances of the proposed approach are demonstrated through intensive numerical experiments both on robust classification and regression problems.

Introduction
The support vector machine (SVM) has been one of the most successful machine learning algorithms [1,2,3].
However, in recent practical machine learning applications with less reliable data that contains a portion of outliers (e.g., consider situations where the labels are automatically obtained by semi-supervised learning [4,5] or manually annotated in crowdsourcing framework [6,7]

Limitations of Existing Robust Classification and Regression
Although a great deal of efforts have been devoted to improving the robustness of SVM and other similar learning algorithms [8,9,10,11,12,13,14,15,16,17,18], the problem of learning SVM under massive noise still poses two major challenges.
First, robust learning algorithms are essentially formulated as non-convex optimization problems because the loss function must be designed to alleviate the effect of outliers (see the robust loss functions for classification and regression problems in Figure2 and Figure3). Since there could be many local optimal solutions in non-convex optimization problems, it is important to develop non-convex optimization tricks such as simulated annealing [19] that enables us to find good local optimal solution.
The second practiall issue is how to make a balance between robustness and efficiency. Any robust SVM formulations contain an additional hyper-parameter for controling the trade-off between robustness and efficiency. Since such a hyper-parameter gorverns the influence of outliers on the model, it must be carefully tuned based on the property of noise contained in the data set. In practice, we need to empirically select an appropriate value for the hyper-parameter, e.g., by cross-validation, because we usually do not have sufficient knowledge about the noise in advance. Furthermore, due to the non-convexity, robust SVM solutions with slightly different hyper-parameter values can be significantly different, which makes model selection highly unstable.
In this paper, we address these two issues simultaneously by introducing a novel homotopy approach to robust SV classification (SVC) and SV regression (SVR) learnings 1 . Our basic idea is to consider parametrized formulations of robust SVC and SVR which bridge the standard SVM and fully robust SVM via a parameter that gorverns the influence of outliers. We use homotopy methods [20,21,22,23] for tracing a path of solutions when the influence of outliers is gradually decreased. We call the parameter as the robustness parameter and the path of solutions obtained by tracing the robustness parameter as the robustification path. Figure2 and Figure3 illustrate how the robust loss functions for classification and regression problems can be gradually robustified, respectively.

Our contributions
Our first technical contribution is in analyzing the properties of the robustification path for both classification and regression problems. In particular, we derive the necessary and sufficient conditions for SVC and SVR solutions to be locally optimal (note that the well-known Karush-Khun Tucker (KKT) conditions are only necessary, but not sufficient). Interestingly, the analyses indicate that the robustification paths contain a finite number of dicontinuous points. To the best of our knowledge, the above property of robust learning has not been known previously. 1 For regression problems, we study least absolute deviation (LAD) regression. It is straightforward to extend it to original SVR formulation with ε-insensitive loss function. In order to simplify the description, we often call LAD regression as SV regression (SVR). In what follows, we use the term SVM when we describe common properties of SVC and SVR.    Our second contribution is to develop an efficient algorithm for actually computing the robustification path based on the above theoretical investigation of the geometry of robust SVM solutions. Here, we use parametric programming technique [20,21,22,23], which is often used for computing the regularization path in machine learning literature. The main technical challenge here is how to handle the discontinuous points in the robustification path. We develop an algorithm that can precisely detect such discontinuous points, and jump to find a strictly better local optimal solution. Unlike solution path of convex problems [20,21,22,23], our local optimal solution paths is shown to have a finite number of discrete points. We overcome this difficulty by precisely analyzing the necessary and sufficient conditions for the local optimality.
Many existing studies on robust SVM employ Concave Convex Procedure (CCCP) [24] or a variant called Difference of Convex (DC) programming [9,10,11,12,14,15]. In these methods, a non-convex loss function is decomposed into the concave and convex parts. A local optimal solution will be obtained by iteratively solving a sequence of convex optimization problems. Local optimal solutions found by CCCP are sensitive to hyper-parameter values (for controling the robustness and efficiency balance), which makes model selection highly challenging.
Experimental results indicate that our outlier approach can find better robust SVM solutions more efficiently than alternative approaches based on CCCP. We conjecture that there are two reasons why favorable results can be obtained. At first, the robustification path shares similar advantage to simulated annealing [19]. Simulated annealing is known to find better local solutions in many non-convex optimization problems by solving a sequence of solutions along with so-called temprature parameter. If we regard the robustness parameter as the temprature parameter, our robustification path algorithm can be interpreted as simulated annealing with infinitesimal step size. Another possible explanation for favorable performances of our method is the ability of stable and efficient model selection. Since our algorithm provides the path of local solutions, unlike other non-convex optimization algorithms such as CCCP, two solutions with slightly different robustness parameter values tend to be similar, which makes model selection stable. Accoding to our experiments, choice of the robustness parameter is quite sensitive to the generalization performances.
Thus, it is important to finely tune the robustness parameter. Since our algorithm can compute the path of solutions, it is much more computationally efficient than running CCCP many times at different robustness parameter values.

Structure of This Paper
After we formulate robust SVC and SVR as parametrized optimization problems in § 2, we derive in § 3 the necessary and sufficient conditions for a robust SVM solution to be locally optimal, and show that there exist a finite number of discontinuous points in the local solution path. We then propose an efficient algorithm in § 4 that can precisely detect such discontinuous points and jump to find a strictly better local optimal solution. In § 5, we experimentally demonstrate that our proposed method, named the robustification path algorithm, outperforms the existing robust SVM algorithm based on CCCP or DC programming. Finally, we conclude in § 6.
This paper is an extented version of our preliminary conference paper presented at ICML 2014 [25].
In this paper, we have extended our robusitification path framework to the regression problem, and many more experimental evaluations have been conducted. To the best of our knowledge, the homotopy method [20,21,22,23] is first used in our preliminary conference paper in the context of robust learning, So far, homotopy-like methods have been (often implicitly) used for non-convex optimization problems in the context of sparse modeling [26,27,28] and semi-supervised learning [29].
In this section, we first formulate robust SVMs for classification and regression problems, which we denote by robust SVC (SV classification) and robust SVR (SV regression), respectively. Then, we introduce parameterized formulation both for robust SVC and SVR, where the parameter gorverns the influence of outlieres to the model. The problem is reduced to ordinary non-robust SVM at one end of the parameter, while the problem corresponds to fully-robust SVM at the other end of the parameter. In the following sections, we develop algorithms for computing the path of local optimal solutions when the parameter is changed form one end to the other.

Robust SV Classification
Let us consider a binary classification problem with n instances and d features. We denote the training set 1} is the binary class label, and the notation N n := {1, . . . , n} represents the set of natural numbers up to n. We write the decision function as where φ is the feature map implicitly defined by a kernel K, w is a vector in the feature space, and ⊤ denotes the transpose of vectors and matrices.
We introduce the following class of optimization problems parameterized by θ and s: where C > 0 is the regularization parameter which controls the balance between the first regularization term and the second loss term. The loss function ℓ is characterized by a pair of parameters θ ∈ [0, 1] and s ≤ 0 as where [z] + := max{0, z}. We refer to θ and s as the homotopy parameters. Figure 2 shows the loss functions for several θ and s. The first homotopy parameter θ can be interpreted as the weight for an outlier: θ = 1 indicates that the influence of an outlier is the same as an inlier, while θ = 0 indicates that outliers are completely ignored. The second homotopy parameter s ≤ 0 can be interpreted as the threshold for deciding outliers and inliers.
In the following sections, we consider two types of homotopy methods. In the first method, we fix s = 0, and gradually change θ from 1 to 0 (see the top five plots in Figure 2). In the second method, we fix θ = 0 and gradually change s from −∞ to 0 (see the bottom five plots in Figure 2). Note that the loss function is reduced to the hinge loss for the standard (convex) SVC when θ = 1 or s = −∞. Therefore, each of the above two homotopy methods can be interpreted as the process of tracing a sequence of solutions when the optimization problem is gradually modified from convex to non-convex. By doing so, we expect to find good local optimal solutions because such a process can be interpreted as simulated annealing [19]. In addition, we can adaptively control the degree of robustness by selecting the best θ or s based on some model selection scheme.

Robust SV Regression
Let us next consider a regression problem. We denote the training set of the regression problem as where the input x i ∈ X is the input vector as the classification case, while the output y i ∈ R is a real scalar. We consider a regression function f (x) in the form of (1). SV regression is formulated as min w where C > 0 is the regularization parameter, and the loss function ℓ is defined as The loss function in (5) has two parameters θ ∈ [0, 1] and s ∈ [0, ∞) as the classification case. Figure 3 shows the loss functions for several θ and s.

Local Optimality
In order to use the homotopy approach, we need to clarify the continuity of the local solution path. To this end, we investigate several properties of local solutions of robust SVM, and derive the necessary and sufficient conditions. Interestingly, our analysis reveals that the local solution path has a finite number of discontinuous points. The theoretical results presented here form the basis of our novel homotopy algorithm given in the next section that can properly handle the above discontinuity issue. We first discuss the local optimality of robust SVC in detail in § 3.1 and § 3.2, and then present the corresponding result of robust SVR briefly in § 3.3.

Conditionally Optimal Solutions (for Robust SVC)
The basic idea of our theoretical analysis is to reformulate the robust SVC learning problem as a combinatorial optimization problem. We consider a partition of the instances N n := {1, . . . , n} into two disjoint sets I and O. The instances in I and O are defined as Inliers and Outliers, respectively. Here, we restrict that the margin y i f (x i ) of an inlier should be larger than s, while that of an outlier should be smaller than s. We denote the partition as P := {I, O} ∈ 2 Nn , where 2 Nn is the power set 2 of N n . Given a partition P, the above restrictions define the feasible region of the solution f in the form of a convex polytope 3 : Using the notion of the convex polytopes, the optimization problem (2) can be rewritten as where the objective function J P is defined as 4 When the partition P is fixed, it is easy to confirm that the inner minimization problem of (7) is a convex problem.
Definition 1 (Conditionally optimal solutions) Given a partition P, the solution of the following convex problem is said to be the conditionally optimal solution: The formulation in (7) is interpreted as a combinatorial optimization problem of finding the best solution from all the 2 n conditionally optimal solutions f * P corresponding to all possible 2 n partitions 5 . Using the representer theorem or convex optimization theory, we can show that any conditionally optimal solution can be written as where {α * j } j∈Nn are the optimal Lagrange multipliers. The following lemma summarizes the KKT optimality conditions of the conditionally optimal solution f * P .

Lemma 2
The KKT conditions of the convex problem (8) is written as 3 Note that an instance with the margin y i f (x i ) = s can be the member of either I or O. 4 Note that we omitted the constant terms irrelevant to the optimization problem. 5 For some partitions P, the convex problem (8) might not have any feasible solutions.
The proof is omitted because it can be easily derived based on standard convex optimization theory [30].

SVC)
From the definition of conditionally optimal solutions, it is clear that a local optimal solution must be conditionally optimal within the convex polytope pol(P; s). However, the conditional optimality does not necessarily indicate the local optimality as the following theorem suggests.
Theorem 3 For any θ ∈ [0, 1) and s ≤ 0, consider the situation where a conditionally optimal solution f * P is at the boundary of the convex polytope pol(P; s), i.e., there exists at least an instance such that y i f * P (x i ) = s. In this situation, if we define a new partitionP := {Ĩ,Õ} as then the new conditionally optimal solution f * P is strictly better than the original conditionally optimal solution f * P , i.e., JP (f * P ; θ) < J P (f * P ; θ).
The proof is presented in Appendix A. Theorem 3 indicates that if f * P is at the boundary of the convex polytope pol(P; s), i.e., if there is one or more instances such that y i f * P (x i ) = s, then f * P is NOT locally optimal because there is a strictly better solution in the opposite side of the boundary.
The following theorem summarizes the necessary and sufficient conditions for local optimality. Note that, in non-convex optimization problems, the KKT conditions are necessary but not sufficient in general.
Theorem 4 For θ ∈ [0, 1) and s ≤ 0, are necessary and sufficient for f * to be locally optimal.
The proof is presented in Appendix B. The condition (13e) indicates that the solution at the boundary of the convex polytope is not locally optimal. Figure 4 illustrates when a conditionally optimal solution can be locally optimal with a certain θ or s.
Theorem 4 suggests that, whenever the local solution path computed by the homotopy approach encounters a boundary of the current convex polytope at a certain θ or s, the solution is not anymore locally optimal. In such cases, we need to somehow find a new local optimal solution at that θ or s, and restart the local solution path from the new one. In other words, the local solution path has discontinuity at that θ or s. Fortunately, Theorem 3 tells us how to handle such a situation. If the local solution path arrives at the boundary, it can jump to the new conditionally optimal solution f * P which is located on the opposite side of the boundary. This jump operation is justified because the new solution is shown to be strictly better than the previous one. Figure 4 (c) and (d) illustrate such a situation.

Local optimality of SV Regression
In order to derive the necessary and sufficent conditions of the local optimality in robust SVR, with abuse of notation, let us consider a partition of the instances N n into two disjoint sets I and O, which represent inliers and outliers, respectively. In regression problems, an instance (x i , y i ) is regarded as an outlier if the aboslute residual |y i − f (x i )| is sufficiently large. Thus, we define inliers and outliers of regression problem Given a partition P := {I, O} ∈ 2 Nn , the feasible region of the solution f is represented as a convex polytope: Then, as in the classification case, the optimization problem (4) can be rewritten as min P min f ∈pol(P;s) where the objective function J P is defined as Since the inner problem of (15) is a convex problem, any conditionally optimal solution can be written as The KKT conditions of f * P (x) are written as Based on the same discussion as §3.2, the necessary and sufficient conditions for the local optimality of robust SVR are summarized as the following theorem: Theorem 5 For θ ∈ [0, 1) and s ≥ 0, if (y i f (x i ) = s ∀ i ∈ N n ) then

6:
Run C-step. We omit the proof of this theorem because they can be easily derived in the same way as Theorem 4.

Outlier Path Algorithm
Based on the analysis presented in the previous section, we develop a novel homotopy algorithm for robust SVM. We call the proposed method the outlier-path (OP) algorithm. For simplicity, we consider homotopy path computation involving either θ or s, and denote the former as OP-θ and the latter as OP-s. OP-θ computes the local solution path when θ is gradually decreased from 1 to 0 with fixed s = 0, while OP-s computes the local solution path when s is gradually increased from −∞ to 0 with fixed θ = 0.
The local optimality of robust SVM in the previous section shows that the path of local optimal solutions has finite discontinuous points that satisfy (13e) or (18d). Below, we introduce an algorithm that appropriately handles those discontinuous points. In this section, we only describe the algorithm for robust SVC. All the methodologies described in this section can be easily extended to robust SVR counterpart.

Overview
The main flow of the OP algorithm is described in Algorithm 1. The solution f is initialized by solving the standard (convex) SVM, and the partition P := {I, O} is defined to satisfy the constraints in (6). The algorithm mainly switches over the two steps called the continuous step (C-step) and the discontinuous step (D-step).
In the C-step (Algorithm 2), a continuous path of local solutions is computed for a sequence of gradually

Algorithm 2 Continuous
Step (C-step) Solve the sequence of convex problems, min f ∈pol(P;s) for gradually decreasing θ in OP-θ or gradually increasing s in OP-s. The C-step can be implemented by any homotopy algorithms for solving a sequence of quadratic problems (QP). In OP-θ, the local solution path can be exactly computed because the path within a convex polytope can be represented as piecewise-linear functions of the homotopy parameter θ. In OP-s, the C-step is trivial because the optimal solution is shown to be constant within a convex polytope. In § 4.2 and § 4.3, we will describe the details of our implementation of the C-step for OP-θ and OP-s, respectively.
In the D-step, we only need to solve a single quadratic problem (QP). Any QP solver can be used in this step. We note that the warm-start approach [31] is quite helpful in the D-step because the difference between two conditionally optimal solutions in adjacent two convex polytopes is typically very small. In § 4.4, we describe the details of our implementation of the D-step. Figure 5 illustrates an example of the local solution path obtained by OP-θ.

Continuous-Step for OP-θ
In the C-step, the partition P := {I, O} is fixed, and our task is to solve a sequence of convex quadratic problems (QPs) parameterized by θ within the convex polytope pol(P; s). It has been known in optimization literature that a certain class of parametric convex QP can be exactly solved by exploiting the piecewise linearity of the solution path [23]. We can easily show that the local solution path of OP-θ within a convex polytope is also represented as a piecewise-linear function of θ. The algorithm presented here is similar to the regularization path algorithm for SVM given in [32].
Let us consider a partition of the inliers in I into the following three disjoint sets: For a given fixed partition {R, E, L, O}, the KKT conditions of the convex problem (8) indicate that The KKT conditions also imply that the remaining Lagrange multipliers {α i } i∈E must satisfy the following linear system of equations: where Q ∈ R n×n is an n × n matrix whose (i, j) th entry is defined as Q ij := y i y j K(x i , x j ). Here, a notation such as Q EL represents a submatrix of Q having only the rows in the index set E and the columns in the index set L. By solving the linear system of equations (19), the Lagrange multipliers α i , i ∈ N n , can be written as an affine function of θ.
Noting that y i f (x i ) = j∈Nn α j y i y j K(x i , x j ) is also represented as an affine function of θ, any changes of the partition {R, E, L} can be exactly identified when the homotopy parameter θ is continuously decreased.
Since the solution path linearly changes for each partition of {R, E, L}, the entire path is represented as a continuous piecewise-linear function of the homotopy parameter θ. We denote the points in θ ∈ [0, 1) at which members of the sets {R, E, L} change as break-points θ BP .
Using the piecewise-linearity of y i f (x i ), we can also identify when we should switch to the D-step. Once we detect an instance satisfying y i f (x i ) = s, we exit the C-step and enter the D-step.

Continuous-Step for OP-s
Since θ is fixed to 0 in OP-s, the KKT conditions (10) yields This means that outliers have no influence on the solution and thus the conditionally optimal solution f * P does not change with s as long as the partition P is unchanged. The only task in the C-step for OP-s is therefore to find the next s that changes the partition P. Such s can be simply found as

Discontinuous-Step (for both OP-θ and OP-s)
As mentioned before, any convex QP solver can be used for the D-step. When the algorithm enters the D-step, we have the conditionally optimal solution f * P for the partition P := {I, O}. Our task here is to find another conditionally optimal solution f * P forP := {Ĩ,Õ} given by (11). Given that the difference between the two solutions f * P and f * P is typically small, the D-step can be efficiently implemented by a technique used in the context of incremental learning [33].
Let us define Then, we consider the following parameterized problem with parameter µ ∈ [0, 1]: and α (bef) be the corresponding α at the beginning of the D-Step. We can show that fP (x i ; µ) is reduced to f P (x i ) when µ = 1, while it is reduced to fP (x i ) when µ = 0 for all i ∈ N n . By using a similar technique to incremental learning [33], we can efficiently compute the path of solutions when µ is continuously changed from 1 to 0. This algorithm behaves similarly to the C-step in OP-θ. The implementation detail of the D-step is described in Appendix C.

Numerical Experiments
In this section, we compare the proposed outlier-path (OP) algorithm with conventional concave-convex procedure (CCCP) [24] because, in most of the existing robust SVM studies, non-convex optimization for robust SVM training are solved by CCCP or a variant called difference of convex (DC) programming [9,10,11,12,14,15].

Setup
We used several benchmark data sets listed in Tables 1 and 2. We randomly divided data set into training (40%), validation (30%), and test (30%) sets for the purposes of optimization, model selection (including the selection of θ or s), and performance evaluation, respectively. For robust SVC, we randomly flipped 15% of the labels in the training and the validation data sets. For robust SVR, we first preprocess the input and output variables; each input variable was normalized so that the minimum and the maximum values are −1 and +1, respectively, while the output variable was standardized to have mean zero and variance one.
Then, for the 5% of the training and the validation instances, we added an uniform noise U (−2, 2) to input variable, and a Gaussian noise N (0, 10 2 ) to output variable, where U (a, b) denotes the uniform distribution between a and b and N (µ, σ 2 ) denotes the normal distribution with mean µ and variance σ 2 .

Generalization Performance
First, we compared the generalization performance. We used the linear kernel and the radial basis function where γ is a kernel parameter fixed to γ = 1/d with d being the input dimensionality. Model selection was carried out by finding the best hyper-parameter combination that minimizes the validation error. We have a pair of hyper-parameters in each setup. In all the setups, the regularization parameter C was chosen from {0.01, 0.1, 1, 10, 100}, while the candidates of the homotopy parameters were chosen as follows: with f SVC being the ordinary non-robust SVR.
• In CCCP-s, the homotopy parameter s was selected from s ∈ {s C , 0.75s C , 0.5s C , 0.25s C , 0} for robust SVC, while it was selected from s ∈ {s R , 0.8s R , 0.6s R , 0.4s R , 0.2s R } for robust SVR.
Note that both OP and CCCP were initialized by using the solution of standard SVM.
Tables 3-6 represent the average and the standard deviation of the test errors on 10 different random data splits. These results indicate that our proposed OP algorithm tends to find better local solutions and the degree of robustness was appropriately controlled.

Computation Time
Finally, we compared the computational costs of the entire model-building process of each method. The results are shown in Figure 6. Note that the computational cost of the OP algorithm does not depend on

Conclusions
In this paper, we proposed a novel robust SVM learning algorithm based on the homotopy approach that allows efficient computation of the sequence of local optimal solutions when the influence of outliers is gradually emphasized. The algorithm is built on our theoretical findings about the geometric property and the optimality conditions of local solutions of robust SVM. Experimental results indicate that our algorithm tends to find better local solutions possibly due to the simulated annealing-like effect and the stable control of robustness. One of the important future works is to adopt scalable homotopy algorithms [28] or approximate parametric programming algorithms [34] for further improving the computational efficiency.

A Proof of Theorem 3
Although f * P is a feasible solution, it is not a local optimum for θ ∈ [0, 1) and s ≤ 0 because violate the KKT conditions (10) forP. These feasibility and sub-optimality indicates that JP (f * P ; θ) < J P (f * P ; θ), we arrive at (12). Q.E.D.

B Proof of Theorem 4
Sufficiency: If (13e) is true, i.e., if there are NO instances with y i f * P (x i ) = s, then any convex problems defined by different partitionsP = P do not have feasible solutions in the neighborhood of f * P . This means that if f * P is a conditionally optimal solution, then it is locally optimal. (13a)-(13d) are sufficient for f * P to be conditionally optimal for the given partition P. Thus, (13) is sufficient for f * P to be locally optimal. Necessity: From Theorem 3, if there exists an instance such that y i f * P (x i ) = s, then f * P is a feasible but not locally optimal. Then (13e) is necessary for f * P to be locally optimal. In addition, (13a)-(13d) are also necessary for local optimality, because of every local optimal solutions are conditionally optimal for the given partition P. Thus, (13) is necessary for f * P to be locally optimal. Q.E.D.

C Implementation of D-step
In D-step, we work with the following convex problem where,P is updated from P as (11).
If we write the conditionally optimal solution as f * P (x) := j∈Nn α * j y j K(x, x j ), (24) {α * j } j∈Nn must satisfy the following KKT conditions s < y i f * P (x i ) < 1 ⇒ α * i = C (25c) At the beginning of the D-step, f * P (x i ) violates the KKT conditions by where α (bef) is the corresponding α at the beginning of the D-step, while ∆ I→O and ∆ O→I denote the difference inP and P defined as ∆ I→O := {i ∈ I | y i f P (x i ) = s}, Then, we consider the following another parametrized problem with a parameter µ ∈ [0, 1]: fP(x i ; µ) := fP (x i ) + µ∆f i ∀ i ∈ N n .
In order to always satisfy the KKT conditions for fP (x i ; µ), we solve the following linear system where A := {E,Ĩ ′ ,Õ ′ }. This linear system can also be solved by using the piecewise-linear parametric programming while the scalar parameter µ is continuously moved from 1 to 0.
In this parametric problem, we can show that f * P (x i ; µ) = f * P (x i ) if µ = 1 and f * P (x i ; µ) = f * P (x i ) if µ = 0 for all i ∈ N n .
Since the number of elements in ∆ I→O and ∆ O→I are typically small, the D-step can be efficiently implemented by a technique used in the context of incremental learning [33].


Toward Depth Estimation Using Mask-Based Lensless Cameras

Recently, coded masks have been used to demonstrate a thin form-factor lensless camera, FlatCam, in which a mask is placed immediately on top of a bare image sensor. In this paper, we present an imaging model and algorithm to jointly estimate depth and intensity information in the scene from a single or multiple FlatCams. We use a light field representation to model the mapping of 3D scene onto the sensor in which light rays from different depths yield different modulation patterns. We present a greedy depth pursuit algorithm to search the 3D volume and estimate the depth and intensity of each pixel within the camera field-of-view. We present simulation results to analyze the performance of our proposed model and algorithm with different FlatCam settings.

I. INTRODUCTION
Lens-based cameras are standard vision sensors in system that records visual information. However, lens-based cameras are bulky, heavy, and rigid-partly because of the size and material of a lens. The shape of a lens-based camera is also fixed in a cube-like form because of the physical constraints on placing the lens at a certain distance from the sensor. A lensless camera can potentially be very thin and lightweight, can operate over a large spectral range, can provide an extremely wide field of view, and can have curved or flexible shape.
Recently, a new lensless imaging system, called FlatCam, was proposed in [1]. FlatCam consists of a coded binary mask placed at a small distance from a bare sensor. FlatCam can be viewed as an example of a coded aperture system in which the mask was placed extremely close to the sensor [2]. The mask pattern was selected in a way that the image formation model takes a linear separable form. Image reconstruction from the sensor measurements requires solving a linear inverse problem.
One limitation of the imaging model in [1] is that it assumes a 2D scene that consists of a single plane at a fixed distance from the camera. In this paper, we present a new imaging model for FlatCam in which the scene consists of multiple planes at different (unknown) depths. We use a light field representation in which light rays from different depths yield different modulation patterns. We use the lightfield representation to analyze the sensitivity of FlatCam to the sampling pattern and depth mismatch. We present a greedy algorithm that jointly estimates the depths and intensity of each pixel. We present simulation results to demonstrate the performance of our algorithm under different settings. : 3D imaging with a mask-based lensless camera that consists of a bare sensor with a fixed, binary mask on top of it. Every light source from within the camera field-of-view casts a shadow of the mask on sensor, resulting in a multiplexed image on the sensor. The shadow of any light source depends on its 3D location with respect to the mask-sensor assembly. A depth-selective pursuit algorithm reconstructs the 3D image of the scene.

II. BACKGROUND AND RELATED WORK
A pinhole camera is a classical example of a lensless camera in which an opaque mask with a single pinhole is placed in front of a light-sensitive surface. A pinhole camera can potentially take an arbitrary shape. However, a major drawback of a pinhole camera is that it only allows a tiny fraction of the ambient light to pass through the single pinhole; therefore, it typically requires very long exposure times. Coded aperture imaging systems extend the idea of a pinhole camera by using a mask with multiple pinholes [3], [2], [4], [5], [6], [7]. However, the image formed on the sensor is a linear superposition of images from multiple pinholes. We need to solve an inverse problem to recover the underlying scene image from sensor measurements. The primary purpose of the coded aperture is to increase the amount of light recorded at the sensor. A coded aperture system offers another advantage by virtue of encoding light from different directions and depths differently. Note that a bare sensor can provide the intensity of a light source but not its spatial location. A mask in front of the sensor encodes directional information of the source in the sensor measurements. In a coded aperture system, every light source in the scene casts a unique shadow of the mask onto the sensor. Therefore, sensor measurements encode information about locations and intensities of all the light sources in the scene. Consider a single light source with a dark background; the image formed on the sensor will be a shadow of the mask. If we change the angle of the light source, the mask shadow on the sensor will shift. Furthermore, if we change the depth of the light source, the size of the shadow will change (see Figure 2). Thus, we can represent the relationship between all the points in the scene and the sensor measurements as a linear system, which depends on the pattern and the placement of the mask. We can solve this system using an appropriate computational algorithm to recover the image of the scene. The depth-dependent imaging capability in coded aperture systems is known since the pioneering work in this domain [8], [2]. The following excerpt in [2] summarizes it well: "One can reconstruct a particular depth in the object by treating the picture as if it was formed by an aperture scaled to the size of the shadow produced by the depth under consideration." However, the classical methods assume that the scene consists of a single plane at known depth. In this paper, we assume that the depth scene consists of multiple depth planes and the true depth map is unknown at the time of reconstruction.
Coded-aperture cameras have traditionally been used for imaging wavelengths beyond the visible spectrum (e.g., Xray and gamma-ray imaging), for which lenses or mirrors are expensive or infeasible [3], [2], [4], [5], [6], [7]. Maskbased lensless designs have been proposed for flexible field-ofview selection in [9], compressive single-pixel imaging using a transmissive LCD panel [10], and separable coded masks [11]. In recent years, coded masks and light modulators have been added to lens-based cameras in different configurations to build novel imaging devices that can capture image and depth [12] or 4D light field [13], [14] from a single coded image. Light field imaging by moving a lensless cameras has been demonstrated in [15] and 3D imaging with a single-shot diffuser-based lensless camera was recently demonstrated in [16].

III. FLATCAM: REPLACING LENSES WITH CODED MASKS AND COMPUTATIONS
FlatCam is a coded aperture system that consists of a bare, planar sensor and a binary mask [1]. Coded-aperture cameras have traditionally been used for imaging wavelengths beyond the visible spectrum (e.g., X-ray and gamma-ray imaging), for which lenses or mirrors are expensive or infeasible [3], [2], [6]. A bare sensor can provide information about the intensity of a light source but not its spatial location. By adding a mask in front of the sensor, we can encode directional information of the source in the sensor measurements.
The imaging model in [1] assumes that the scene consists of a single plane parallel to mask-sensor planes. Let us consider a 1D imaging system, shown in Fig. 3a, in which a 1D mask is placed at a distance d in front of a 1D sensor array with M pixels, and the scene consists of a single plane at distance D from the sensor with N scene pixels. Let us denote a scene pixel as l(θ), where θ is uniformly distributed along an angular interval [θ − , θ + ]. We can represent the measurement at sensor pixel s as where ϕ(θ, s) denotes the modulation coefficient of the mask for a light ray between scene pixel l(θ) and the sensor pixel at location s. We can write (1) in a compact form as where Φ denotes an M × N system matrix that maps N scene pixels to M sensor pixels. This is a linear system that we can solve using an appropriate computational algorithm to recover the image of the scene (more details can be found in [1]).

A. Imaging model
The model in (1) assumes a 2D scene that consists of a single plane at some fixed depth. The system matrix Φ encodes mapping of scene points from that plane to the sensor pixels. In our new model, we consider a 3D scene that consists of multiple planes, each of which contribute to the sensor measurements. Without loss of generality, we consider a 1D imaging system in Fig. 3a and assume that the sensor plane is centered at the origin and the mask plane is placed in front of it at distance d. The scene consists of K planes at depths [D 1 , . . . , D K ] and the scene pixels are distributed uniformly along angles in interval [θ − , θ + ], as before. We can describe the measurement at sensor pixel s as where ϕ(s, θ, D) denotes modulation coefficient of the mask for a light ray between light source l(θ, D) and the sensor pixel at location s. We can write (3) in a compact form as where each l D denotes intensity of N pixels in a plane at depth D and Φ D is an M × N matrix that represents the mapping of l D onto the sensor. For the case of 3D imaging, let us represent the light distribution as L(θ x , θ y , z) for z > d. Measurements for a sensor pixel located at (s u , s v ) can then be described using (a) Imaging system geometry. Mask and sensor planes are separated by distance d. Each point in the scene (at any angle θ and distance D) contributes to the sensor measurements. Our goal is to jointly estimate depth and intensity of each pixel within the field-of-view.
where mask [u, v] denotes the transparency value of the mask at location (u, v) within the mask plane. If the mask pattern is symmetric and separable in (θ x , θ y ) space, then we can write the 2D measurements in (5) as where Φ D is an M × N matrix and L D denotes N × N light distribution corresponding to plane at depth D. Furthermore, we can include multiple cameras at different locations and orientations with respect to a reference frame. Such a system will provide us multiple sensor measurements of the form where I c denotes sensor measurements at camera c and Φ D,c is matrix that represents mapping of L D onto camera c.

B. Joint image and depth reconstruction
We estimate the depth and intensity of each pixel within the field-of-view of our cameras using a greedy depth-selective algorithm, in which we assume a sparse prior that L(θ x , θ y , D) has nonzero value only for one depth. Our proposed algorithm is inspired by structured sparse recovery algorithms in modelbased compressive sensing [17].
To simplify the presentation, let us represent (6) or (7) as the following general linear system: where L is an N × N × K light distribution with N × N spatial resolution and K depth planes, A denotes the linear measurement operator in (6) or (7), and I denotes the sensor measurements. Suppose our current estimate of L isL with exactly one depth assigned to each pixel. Let us denote the initial depth map as Ω. In our experiments, we initialize the depth estimate with the farthest plane in the scene. Our proposed depthselective pursuit algorithm is an iterative method that performs the following three main steps at every iteration: Compute proxy depth estimate. We first select new candidate depth for each pixel by picking the maximum magnitude in the following proxy map corresponding to each pixel: Let us denote the new depth map as Ω.
Merge depths and estimate image. We first merge the original depth estimate with the proxy depth estimate. Let us denote the merged depth support as T = {Ω ∪Ω}. Then we solve a least-squares problem over the merged depth support asL = arg min L I T − A T L 2 2 . Prune depth and threshold image. We prune the depth estimate at every spatial location by picking the depth corresponding to higher pixel intensity inL. Finally, we threshold L toL that has only one nonzero depth per spatial location.

C. Depth sampling and sensitivity
For a single camera, sensor measurements for a single point source at location (θ, D) can be described as From this lightfield expression we note that the slope of a line corresponding to any light source is inversely proportional to its depth. As a light source moves farther from the mask-sensor assembly, its line would rotate around the center. As a light source moves along a plane at a fixed depth, its line would shift with the same slope. Therefore, in our imaging model, we select depth planes in a given range by sampling lightfield at uniform angles, which results in planes at non-uniform depths.

V. EXPERIMENTAL RESULTS
To validate the performance of our proposed imaging model and reconstruction algorithm, we performed extensive simulations under different settings of FlatCam parameters.
First, we show results of a simple simulation in which our scene consists of three depth planes as shown in Fig. 4. We simulated an imaging system in which a binary mask is placed at 1mm distance from the sensor. We simulated a 3D voxel space with ten depth planes within a depth range of 100mm and 3m. We chose the ten depth planes by uniformly sampling the lightfield representation. We simulated a scene with 128 × 128 spatial resolution, a mask with a binary random sequence, and a sensor with 256 × 256 pixels. To generate sensor measurements, we assumed that the 3D scene consists of three planes chosen at random out of ten fixed planes. We added a small amount of Gaussian noise (a) Test scene with three cards placed at three different depths picked out of K = 10 depth planes at random.
(b) Image reconstructed by assuming that the scene consists of a single plane at a fixed depth. PSNR = 15 dB (c) Image reconstructed using depth-selective algorithm that jointly estimates the depth and intensity of each pixel. PSNR = 33.5 dB (d) Image reconstructed with three cameras via depth-selective algorithm that jointly estimates the depth and intensity of each pixel. PSNR = 39 dB in the true measurements. To reconstruct the 3D scene, we solved the depth-selective pursuit algorithm described in the previous section. The results are presented in Fig. 4, where (a) denotes pixel intensities for three planes in a 3D scene, (b) denotes image reconstructed by assuming that all pixels belong to a plane at fixed depth, and (c) denotes images reconstructed by solving the depth-pursuit algorithm on the same measurements.
Next we discuss an experiment that demonstrates the robustness of our proposed model and method against mismatch in the locations of the original depth planes and those used for reconstruction. In this experiment, we simulated imaging system with K = {5, 10, 15, 20, 25} depth planes chosen uniformly in the lightfield representation. We selected the three depth planes for the scene at random and calculated the peak signal to noise ratio (PSNR) for the recovered images. We present PSNR for each test (averaged over ten instances) in Table I. We see that the quality of reconstruction remains almost the same as we increase the number of depth planes in our model. The computational complexity, however, slightly increases as we increase the number of depth planes.
Finally, we present an experiment in which we simulated a system with three lensless cameras in a convex geometry, where one camera is used as a reference to generate depth planes in the scene. The other two cameras see tilted planes in their field of view. An advantage of such a configuration is that the depth information of the pixels is converted into angular information. However, this configuration also makes a strong assumption that the scene only consists of the finite  number of tilted planes. The simulation results of three camera system are also summarized in Table I. An example of image reconstruction for this case is shown in Fig. 4(d).


Smallest Univalue Segment Assimilating Nucleus approach to Brain MRI Image Segmentation using Fuzzy C-Means and Fuzzy K-Means Algorithms

Image segmentation still remains an important task in image processing and analysis. Sequel to any segmentation process, preprocessing activities carried out on the images have a great effect on the accuracy of the segmentation task. This paper therefore laid emphasis on the preprocessing stage of brain Magnetic Resonance Imaging (MRI) images using Smallest Univalue Segment Assimilating Nucleus (SUSAN) and bias field correction algorithms. Subsequently, brain tissue extraction tool was employed in extracting non-brain tissues from the brain image. Afterwards, Fuzzy K-Means (FKM) and Fuzzy C-Means (FCM) segmentation algorithms were employed for segmenting brain MRI images acquired from four different MRI databases into their White Matter (WM), Gray Matter (GM) and Cerebrospinal Fluid (CSF) constituents. Evaluation metrics such as cluster validity functions using partition coefficients and partition entropy; area error metrics such as false positive, true positive, true negative and false negative (FN); similarity index, sensitivity and specificity were used to evaluate the performance of both techniques. A comparative analysis of the experimental results revealed that in most instances, FKM segmentation technique is preferable to FCM segmentation technique for brain MRI segmentation task.

INTRODUCTION
For accurate information elicitation from medical images, image segmentation is a crucial task that cannot be undermined. It is an essential routine because the eventual outcome of the analysis will determine the success of the pattern image segmentation stage. This involves measuring and visualizing the brain's anatomical structures, analyzing brain changes, delineating pathological regions, surgical planning and image-guided interventions [1]. Generally, image segmentation entails partitioning an image of interest into its constituent regions that are non-overlapping whilst possessing similar attributes that are semantically meaningful. With focus on brain Magnetic Resonance Image (MRI) and iris images, segmentation of MRI image entails classifying the image into its constituent components: the white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF) while segmenting IRIS image entails classifying the Iris image into iris setosa, iris versicolour, and iris virginica. An efficient image segmentation technique must be able to uniquely and correctly classify these dataset into their constituent components to an acceptable degree.
Generally, image segmentation techniques could be classified into layer and block-based segmentation methods [2]. While layered approach defines image shape masks and explain their appearance and depth, the block based approach focuses on extracting features such as colour, textual and spatial information. Block based approach also provides information about the pixels, regions, edges or boundaries of extracted features. This research employs fuzzy variants of C and K means algorithms mostly denoted as Fuzzy C-Means (FCM) and Fuzzy K-Means (FKM) algorithms; a block based approach in segmenting brain MRI and IRIS images. The fuzzy nature of fuzzy algorithms is established in the fact that when an image is broken down into segments, each segment can belong to more than one cluster, though different segments have varying membership function which indicates their strength of association between the element and other clusters. Fuzzy segmentation technique involves assigning membership levels to the segments and then using these levels to assign data elements to one or more clusters. The advantages of fuzzy segmentation approach includes yielding regions more homogeneous than other methods; reducing the spurious blobs; removing noisy spots; reduced sensitivity to noise compared to other techniques [3]. However, fuzzy approach requires prior knowledge about the number of clusters present or possible in a particular dataset; this may be difficult to obtain a priori in a new dataset. In this research work, I S S N 2 2 7 7 -3061 V o l u m e 1 6 N u m b e r 7 I n t e r n a t i o n a l j o u r n a l o f C o m p u t e r s a n d T e c h n o l o g y 7066 | P a g e D e c e m b e r , 2017 http://cirworld.com/ DOI: 10.24297/ijct.v16i7.6170 emphasis will be laid on Fuzzy C-Means (FCM) and Fuzzy K-Means (FKM) clustering algorithms for brain MRI segmentation.

RELATED WORKS
FCM has been widely used for image segmentation in several fields such as feature extraction, pattern and image recognition, fuzzy identification etc. it is widely accepted because it preserves more information from the original image than other segmentation methods [4]. FCM segmentation technique was adopted in [5] for the identification and segmentation of lung nodules in Computed Tomography (CT) images. Their results showed that FCM was able to effectively separate the parenchyma from the mediastinum and the thoracic wall from the lung nodules. A study to determine the most suitable cluster validity index needed to achieve optimal clustering was carried out in [3]. The study employed the topmost 18 validity index in turn to synthetic images corrupted with noise of varying levels, and simulated volumetric MRI datasets. The results obtained showed that the various indexes have different outcomes for various noise levels, while some indexes can help in determining the accurate number of clusters existing in a dataset. As no algorithm is entirely flawless, several improvements of FCM have been proposed by researchers, these have led to the emergence of various versions of FCM. Some of these variations reported in [6] include Bias Corrected A comprehensive performance analysis among these variations was carried and their performance was validated on real MRI images as well as synthetic images. Partition coefficient (Vpc), partition entropy (Vpe), time complexity and segmentation accuracy were used as the performance evaluation metrics. High value of Vpc implies better performance while low value of Vpe implies better clustering. SFCM, TEFCM and MDFCM yielded a high Vpc and low Vpe values which denotes better performance and clustering. When considering feature structure, RFCMK yielded an appreciable performance while WIPFCM, MDFCM and KWFLICM also produce a better performance. Fuzzy c means was modified for segmenting T1-T2 weighted brain MRI images in [4], the modified algorithm called Modified Robust Fuzzy c-Means with weight Bias Estimation (MRFCM-wBE) was employed to solve the intensity inhomogeneity and noise affecting the segmentation process of computer assisted segmentation algorithms across differentiating borders between tissues of medical images. Fuzzy C-means was modified by reducing the number of iterations using dist-max initialization algorithm after which the algorithm was executed iteratively. Running time for obtaining clusters, number of iterations for completion of clustering the datasets, and silhouette width for clustering accuracy were the yardstick used to measure the performance of MRFCM-wBE. Checkerboard and lung cancer datasets were used as test data. The performance of the algorithm was compared with other segmentation algorithms such as Gaussian Kernel based Fuzzy c-Means algorithm (GKFCM), BCFCM, Improved Fuzzy Segmentation (IFS) algorithm and SFCM. The algorithm showed an appreciable improvement in running time over the remaining algorithms. In [7], Micro calcification clusters in digitized mammograms were detected using FCM and Possibilistic FCM segmentation techniques. The approach adopted involved detecting and separating micro calcifications which are small deposits of calcium in breast tissue. The results obtained showed that the algorithms adopted could detect and separate micro calcifications for normal tissues though PFCM performed better since its obtained results depend solely on the threshold value being the most appropriate. FCM has also been used in [8] for the detection of leukemia from blood images. This employs morphological contour segmentation approach to identify the edges of nuclei in white blood cells and further separates them from the microscopic blood image. The approach adopted yielded a better result in terms of accuracy and time consumption when compare to normal hematologist's visual classification. In addition to leukemia detection, FCM segmentation technique could be used in jaw tumor detection as carried out in [9] using panoramic X-ray images as test data. The approach employed involved grouping and clustering jaw's X-ray image pixels into background region and lesion region. Comparing the results obtained with other segmentation techniques, an appreciable result was obtained in terms of specificity, Similarity Index, Area error metrics and Sensitivity value.
Fuzzy K-Means (FKM)a variant of k-means algorithm has also been widely used for image segmentation. In FKM segmentation technique, the data of interest is partitioned into k-groups of disjoint cluster. To begin with, the k-centroid point is calculated after which the cluster with the nearest distance to the centroid point is determined. A research to determine if FKM actually divides a data into k groups and if the k-clusters chosen is accurate was carried out in [10]. Their work also looks into the color space at which KCM segments data better. Experiments were carried out on L*a*b colour space and RGB colour space images using 1, 2 and k-cluster color images (k>3). Silhouette Width which is a function of cluster's tightness and separation likewise Ground truth images were used to evaluate the adopted approach. Results showed that KCM considerably segment images much better in L*a*b* colour space when compared to RGB colour space.
The initial clusters centre in FKM algorithms are generated randomly; this leads to non-repeatable clustering results that may be hard to comprehend. As a result of this, a kernel density based clustering technique where the data is mapped to a high-dimensional space for classification was proposed in [11]. This involved selecting an initial seed value with maximum variance from the data matrix, then employing Gaussian kernel for the estimation of the initial density. Afterwards, subsequent seed values are selected that has densities equivalent to the initial selected points. In this way, k points with similar density with respect to each other are obtained. The proposed technique was able to avoid the dead and trapped centre limitation of clustering based segmentation approach.
The performance of three fuzzy algorithms viz: FCM, Gustafson-Kessel (GK) and FKM algorithms were analyzed based on their clustering output criteria in [12]. Liver disorders and wine datasets were used to test the performance of these algorithms. With FKM and FCM using Euclidian distance measure and GK using fuzzy covariance matrix in their distance metrics, the obtained experimental results showed that GK produced a similar result to FCM but FKM outperforms both FCM and GK algorithms. Thus, the efficiency of FKM was reported to be better than FCM and GK algorithms. Few comparative work between FKM and FCM have been carried out so far, hence, this research work attempts to further test FKM and FCM algorithm on brain MRI images with a view to establish which algorithms perform better.

Materials and Methods
The proposed method consists of pre-processing and the actual segmentation stage. The result of the segmentation stage to a large extent is dependent on the success of the preprocessing stage; hence, efforts were put in place to ensure that the pre-processing stage was meticulously carried out.

Pre-processing stage
The preprocessing stage involves three steps which are: noise reduction, bias field correction and brain extraction.

Noise Reduction
To eliminate the effect of noise on the image and improve local features using local smoothing, Smallest Univalue Segment Assimilating Nucleus (SUSAN) algorithm illustrated in [13] was used. SUSAN algorithm perfectly works for carrying out edge detection, corner detection and structure-preserving image noise reduction. SUSAN algorithm achieves this by placing a predetermined window (mask) at each point of an image and then calculates the brightness threshold of each pixel within the mask at each point. The calculated brightness was compared with the center point or nucleus of the pixel using equation (1): Where is the position of the nucleus in the image and is the position at any other point within the mask. I is the brightness of any pixel with "t" as the brightness difference threshold and "c" is the output of the comparison. This comparison was done for each pixel within the mask using equation (2): Where n is the number of pixels in the USAN. The comparison gives the USAN's area which will be ultimately minimized. The brightness difference threshold "t" helps in knowing the minimum contrast of features which will be detected likewise the maximum amount of noise which will be ignored. The next step of SUSAN algorithm is to carry out a comparison between "n" and a fixed threshold "g". There would have been no need for the geometric threshold if noise were to be absent. Nonetheless, to ensure noise is totally eradicated or reduced to the nearest minimum, "g" is set to 3nmax /4 where nmax is the maximum value of "n". The initial edge response is calculated using equation (3): Depending on the type of edge point being expected (typical straight edge or standard edge points), the direction of the edges associated to an image point must be calculated. If it were to be a typical straight edge, equation (3) can be used to compute the required edge else the direction of the edge must be calculated by first computing the center of gravity followed by finding the longest axis of symmetry. The center of gravity could be computed using equation (4) while the edge direction is projected by finding the sum of equations (5, 6 and 7): If the computed USAN area is smaller than the mask diameter then the edge belongs to the typical edge category however, if the USAN area is larger then, it belongs to the standard edge category [1].

Bias Field Correction
Bias field arises as a result of spatial inhomogeneity of the magnetic field, deviations in the sensitivity of the reception coil, and contact between the magnetic field and the human body [22]. Scanning MRI at 0.5T could make the effect of bias field minimal, hence, its effect could be neglected but when MRI are scanned with MR scanners with magnetic field above 1.5T, the effect of the bias field is always strong and could affect MRI analysis [4], at this instance the effect needs to be corrected. To correct the effect of this signal, the modified fuzzy-c means approach reported in [21] was adopted. This entails the following steps: (a). Choosing an initial cluster prototypes and setting the bias field to a very small positive value (b). Updating partition matrix where using equation (8): Where and the distance expression . is the observed intensity at the pixel, is the bias field at the pixel, is the cluster prototype, is the set of neighbours that exist in a window around the true intensity, is the cardinality of this set of neighbours and the effect of these neighbours is controlled by the parameter .
(c). Obtaining new prototypes of the clusters in the form of weighted averages of the patterns using: (d). Compute the bias term using: (e). Steps (b) -(d) will be repeated till termination criteria in (11) is satisfied (11) Where V is a vector of cluster centers and ϵ is a small number that can be taken during the initialization process.

Removal of Non-brain Tissue
Since emphasis is on brain tissue segmentation, removal of non-brain and extra-cerebral tissues like skull, neck, fat, bones, skin, eyeballs e. t. c. is necessary. These tissues can distort the segmentation process and affect the accuracy of the segmentation result. There are several brain tissue extraction algorithms such as: statistical parametric mapping, brain surface extractor, minneapolis consensus strip, threshold morphologic brain extraction and Brain Extraction Tool (BET). This research work employs BET; a widely known publicly available brain tissue extraction tool in extracting non-brain tissues from the brain image. The "generate binary brain mask image" option of BET was used to effect the brain tissue extraction.

Segmentation Stage
The segmentation task was carried out using FCM and FKM clustering algorithms.

FUZZY C-MEANS Technique
With FCM, a pattern or data sample can belong to more than one cluster while the membership value assigned to each pattern is used to indicate the similarity. In this research, FCM clustering begins by supplying a set of data and a number of clusters into Matlab environment using the Matlab function fcm. The function in turn returns the best cluster centers and membership value for each pattern. FCM clustering is an iterative process, hence, the cluster centers was repeatedly moved to the right location within the dataset by minimizing the objective function. The iteration terminates when the objective function enhancement between two consecutive iterations is less than the minimum amount of improvement specified [16]. FCM algorithm as illustrated in [3,4,17] expects a dataset of interest "x" with n data samples such that to be partitioned into clusters "c". FCM clustering task begins by minimizing the objective function J which is defined as: is the membership value of dataset for each class "I" and pixel "j", is the cluster while is the number of clusters. Since FCM employs an iterative process, this is achieved by with updating the membership value and cluster center repeatedly using equations (13) and (14) Where "·" is a norm metric. The membership value must satisfy the constraint: (15)

FUZZY K-MEANS Technique
Integrating fuzzy logic into K-Means clustering algorithm yields FKM [16]. FKM as illustrated in [13,15] implies that the data set of interest well be partitioned into k-clusters such that each cluster is compact and far from other clusters. This can be achieved by minimizing the distance between the cluster centers and the patterns that belong to that cluster. To establish the fuzziness of the pattern and the cluster center, the degree of the belongingness was represented by a characteristic function . This function must fulfill the constraint: (19) The objective function J must be minimized with respect to and using equation (20), such that  (20) Where the number of patterns is denoted by N, the measured Euclidean distance between and is denoted by while the number of clusters is denoted by k. Typically, the patterns are not expected to overlap but if they do, then each pattern belongs to more than one cluster i.e In this instance, should be referred to as a membership function rather than a characteristic function. As a result of this, equation (20) was modified as  (21) Where m remains the fuzzifier parameter which controls the degree of fuziness. Minimizing the objective function J may produce an insignificant solution, hence, constraints illustrated in equations (22) and (23)

Dataset used
The segmentation approach was validated using four different brain MRI datasets. Real brain MRI datasets were downloaded from [18]. The dataset was acquired using a physical phantom on a 3T MRI scanner with a turbo spin-echo sequence. The data was acquired with a 220 mm x 292 mm field of view on a 256×340 Cartesian sampling grid and was collected with an array of receiver coils, leading to 16 channels of information. The data was provided in Matlab format and were referred to as "D1". Synthetic brain MRI dataset was also downloaded from [23]. Each image in the downloaded dataset (referred to as D2) has a slice thickness of 1mm, resolution of 1mmX1mmX1mm, 20% intensity non-uniformity (RF), noise levels ranging from 0% to 9%, and 181X217X181 matrix voxels. All these parameters are provided before the final dataset is downloaded. Manually labeled brain MRI images [D3] were also downloaded from the Internet Brain Segmentation Repository (IBSR) [24]. The labelled images were not created by an algorithm but by neuroanatomical experts. The images in the datasets have different levels of difficulty which can be used to evaluate the segmentation task in different conditions. Finally, T1 and T2 brain MRI dataset from [25] with 1x1x1 mm 3 resolution and a voxel size of 2x2x2 mm 3 were also used.

Evaluation Metrics
The performance of both FCM and FKM algorithms were measured using various performance evaluation metrics such as cluster validity functions, Area error metrics, similarity index, sensitivity and specificity

Cluster Validity Function
Cluster validity function uses fuzzy partition to determine the efficiency of the segmentation task. Fuzzy partition is a function of the partition coefficient and partition entropy .
Partition Coefficients: Partition coefficients ranges from 0 to 1 with a high value indicating a better performance.

Fig 1: Fuzzy Partition Coefficient using D1
As shown in Figure 1, the results showed that both FKM and FCM could correctly classify the phantom dataset into three clusters but at different partition coefficient. FKM having the highest partition coefficient 0.85 at the third cluster center shows that FKM segments the dataset better at a higher partition coefficient when compared to FCM at 0.65. Partition coefficients of the remaining three datasets are highlighted in table 1. The results show that both techniques have the peak partition coefficients at the 3 rd cluster points across all the datasets (D2, D3 and D4) used. This showed that all the algorithms could at least identify the three clusters correctly.

Partition Entropy
A high partition coefficient is expected for an accurate clustering while low partition entropy is expected. It is computed using equation (27)

Area Error Metrics
Using manually labeled brain MRI images [D3] as the reference image, the area ratio of the False Positive (FP), True Positive (TP), True Negative (TN) and False Negative (FN) as used in [9] was employed to determine how correctly FKM and FCM could classify dataset D3 into its constituent WM, GM and CSF component. The area metrics are given as: The manually measured cluster area is represented by while those computed using FCM and FKM are represented by . The results obtained as illustrated in Table 2 showed that FKM performed preferably better when compared with FCM on dataset "D3". FKM having a high TP at 94% when compared to that of FCM at 90% showed that FKM could identify more area for each clusters than FCM. Also, the FP of FKM being lower than FCM showed that FKM was able to recognize clusters accurately than FCM. The P-value of each metric is less than 0.005 and this statistically showed that the result is highly significant.

Similarity Index, Sensitivity and Specificity
To further determine the effectiveness of the algorithms, Similarity Index (SI), sensitivity and specificity of each technique were computed using the following equations: SI could also be used to know if there is an overlap among the clusters or not. Specificity can be used to determine if an algorithm could classify an image correctly into the expected clusters (WM, GM, CSF in this sense) while sensitivity calculates the number of images that were correctly classified. The mean of the computed metrics on D3 is as described in Table 3. The result showed that FKM segmentation technique is more accurate and reliable in WM, GM and CSF segmentation when compared to FCM, though the results of FCM segmentation is not totally bad.

CONCLUSION
A comparative analysis aimed at asserting a superior segmentation technique between FKM and FCM algorithm has been carried out in this paper. Using partition coefficients as a metric, both techniques have the peak partition coefficients at the 3 rd cluster points across all the datasets. This showed that the techniques could correctly classify the datasets into three clusters, though, FKM records a higher value when compared to FCM across all the datasets used which showed that FKM classifies at a faster rate when compared to FCM. Also, when the partition entropy for each technique was computed, FKM attained lower partition entropy when compared to FCM and this further established the superiority of FKM segmentation technique. With the manually labeled brain MRI images as a reference image, the area error metrics (FP, TN, TP, FN) of the techniques were computed to further know the technique that classifies the datasets better. FKM achieved good results which further establish the fact that it's a better technique when compared to FCM. The P-value of each metric being less than 0.005 statistically showed that the computed result is highly significant. Conclusively, similarity index, sensitivity and specificity of each technique using dataset D3 were computed. The mean of the computed metrics revealed that FKM segmentation technique is more accurate and reliable in WM, GM and CSF segmentation when compared to FCM. Based on the results of the various evaluation metrics used, it could be conjectured that FKM is a good segmentation technique and its preferable to FCM segmentation technique.