# Piecewise Linear Patch Reconstruction for Segmentation and Description of Non-smooth Image Structures

In this paper, we propose a unified energy minimization model for the segmentation of non-smooth image structures. The energy of piecewise linear patch reconstruction is considered as an objective measure of the quality of the segmentation of non-smooth structures. The segmentation is achieved by minimizing the single energy without any separate process of feature extraction. We also prove that the error of segmentation is bounded by the proposed energy functional, meaning that minimizing the proposed energy leads to reducing the error of segmentation. As a by-product, our method produces a dictionary of optimized orthonormal descriptors for each segmented region. The unique feature of our method is that it achieves the simultaneous segmentation and description for non-smooth image structures under the same optimization framework. The experiments validate our theoretical claims and show the clear superior performance of our methods over other related methods for segmentation of various image textures. We show that our model can be coupled with the piecewise smooth model to handle both smooth and non-smooth structures, and we demonstrate that the proposed model is capable of coping with multiple different regions through the one-against-all strategy.

## Authors

• 12 publications
• 5 publications
• ### Piecewise Padé-Chebyshev Reconstruction of Bivariate Piecewise Smooth Functions

We extend the idea of approximating piecewise smooth univariate function...
09/23/2021 ∙ by Akansha Singh, et al. ∙ 0

• ### A nonlocal feature-driven exemplar-based approach for image inpainting

We present a nonlocal variational image completion technique which admit...
09/20/2019 ∙ by Viktor Reshniak, et al. ∙ 9

• ### Hierarchical image simplification and segmentation based on Mumford-Shah-salient level line selection

Hierarchies, such as the tree of shapes, are popular representations for...
03/15/2016 ∙ by Yongchao Xu, et al. ∙ 0

• ### Multiphase Segmentation For Simultaneously Homogeneous and Textural Images

Segmentation remains an important problem in image processing. For homog...
06/29/2016 ∙ by Duy Hoang Thai, et al. ∙ 0

• ### Discrete Potts Model for Generating Superpixels on Noisy Images

Many computer vision applications, such as object recognition and segmen...
03/20/2018 ∙ by Ruobing Shen, et al. ∙ 0

• ### An Efficient Algorithm for the Piecewise-Smooth Model with Approximately Explicit Solutions

This paper presents an efficient approach to image segmentation that app...
12/08/2016 ∙ by Huihui Song, et al. ∙ 0

• ### Automatic Pollen Grain and Exine Segmentation from Microscope Images

In this article, we propose an automatic method for the segmentation of ...
03/19/2015 ∙ by François Chung, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

Computer vision problems are often addressed by using mathematical models. The quality of the solutions to the problems are measured objectively in the mathematical models. With the valid mathematical models, we can elucidate the phenomenon of the natural computations, e.g. by human vision, that try to accomplish the tasks, and we can reproduce the result of the natural computations by computerized simulations.

An objective measure of piecewise smooth image segmentation is the Mumford-Shah functional energy [1]. By minimizing the energy, we expect to achieve high quality of segmentation for piecewise smooth images. The problem of minimization of the Mumford-Shah functional energy is the mathematical abstraction, i.e. the model, of piecewise smooth image segmentation. The model is now known as the Mumford-Shah model. This methodology is different from that of those in [2] [3] which target at the objective evaluation of the segmentation based on the ground-truth results from normal subjects.

In Mumford-Shah model, each image region is modeled as a smooth or constant function with parameters. By minimizing the energy of the Mumford-Shah model, we can restore the smooth image in each region, and we can also obtain the partition boundaries of the segmentation located at the discontinuities in the restored image. However, the image values to be restored may not be piecewise smooth or flat as a whole. The images may contain regions of non-smooth structures. Such as the images in Fig. 1(a). Imposing the smoothness on these images leads to the destructive averaging of the image content. This poses problem to the segmentation with the conventional Mumford-Shah model. For example, it is possible that non-smooth visual patterns different in structure may have similar average image values (See Fig. 1(a)). Consequently, the Mumford-Shah model cannot separate such patterns in the image space.

To cope with non-smooth image data, the segmentation has been considered as a framework of two independent processes, i.e. the feature extraction and the segmentation of the feature image. There exist combinations of the Mumford-Shah model or other active contour model with predefined features [4] [5] [6]. The predefined features does not necessarily match the underlying non-smooth image structures. For example, in Fig. 1

(b) the feature values obtained by image filtering with the Gabor filters that well match the texture pattern, or the local averages of the feature values, can have significant spatial variations. Unsupervised feature selection has also been used

[7] [8] [9]. The major problem concerns us is that these frameworks, formed by separate feature selection and segmentation, do not provide a unified optimization model for segmentation. These concerns motivate us to explore a unified and valid mathematical model for segmentation of non-smooth structures.

In this paper, we propose a novel unified energy minimization model of the segmentation of general non-smooth image structures, motivated by the original Mumford-Shah model. We formulate the segmentation of general non-smooth image structures as a single energy minimization problem of piecewise linear patch reconstruction. The formulation is not heuristic since we prove that the error of segmentation is bounded by the energy of piecewise linear patch reconstruction for a fixed number of bases. Thus, minimizing the energy of piecewise linear patch reconstruction can reduce the error of segmentation. Regarding the energy minimization, we prove that the eigen-patches constitute the global optimal solution to the minimization of the error of linear patch reconstruction. We also explore a more efficient alternative to the eigen-decomposition of matrix for computing eigenvectors. We prove that the linear patch reconstruction based on gradient descent converges to the eigen-patches under some mild assumption on the initializations. The assumption can often be met, and the convergence is linear. We therefore propose the gradient descent algorithm for solving the linear patch reconstruction problem even though the problem is not convex. The segmentation algorithm is based on alternating piecewise linear patch reconstruction and curve evolution. The piecewise linear patch reconstruction can naturally be coupled with the conventional Mumford-Shah model for coping with both smooth and non-smooth structures. A unique feature of the method is that it produces both the segmentation and a dictionary of optimized (orthonormal) descriptors for each segmented region upon completion in the same optimization framework.

The rest of the paper is organized as follows. We reviewed the related works in section II. We study the PC/PS model and propose the piecewise linear patch reconstruction model in sections III-IV. We show our error bound of the segmentation in terms of reconstruction error in section IV-B and we present our proved theoretical claim of the global optimality of the gradient descent for the nonconvex linear patch reconstruction problem in section IV-C. The experimental results are presented in section V. We conclude the paper in section VI.

## Ii Background and related works

### Ii-a The Mumford-Shah model

The segmentation of piecewise smooth images has been formulated as an energy minimization problem in the Mumford-Shah model. For an image defined over the image domain , the two-phase Mumford-Shah model can be formulated as follows.

 minϕE(ϕ,g1,g2),E(ϕ,g1,g2)=∫ΩE1(x,y,g1)H(ϕ)dxdy  +∫ΩE2(x,y,g2)(1−H(ϕ))dxdy  +ν∫Ωδ(ϕ(x,y))∥∇ϕ∥dxdy (1)

where is a Dirac delta function, is a Heaviside (step) function, is a penalty coefficient. and

are the reconstruction estimates,

is the signed distance function that partitions the image into two regions, i.e. for one region, for the other. and are the reconstruction errors of the two region models. assigns either of the two models to all the pixels, and the last term tries to minimize the complexity of the labeling. For piecewise smooth (PS) model, we have , , where is the image value at , is a constant penalty coefficient. Note that, , and are all functions defined on . We will omit behind these functions henceforth if there is no risk of confusion. The fundamental optimization technique for the reconstruction is the Green’s functions solution to the Euler-lagrange equation obtained by Calculus of Variations, which was presented in [10]. The numerical schemes for solving the Euler-lagrange equation can also be found in [11] [12]. The smoothness regularization guarantees the global optimality of the solution obtained by any of the optimization methods for a given partition.

If we assume the reconstruction functions , to be constants, i.e. for all in , and , then the gradient terms in the PS model vanish, and the functional becomes the piecewise constant (PC) Mumford-Shah model, which is the prototype of the region competition [13] and the Chan-Vese model [14].

For minimizing the Mumford-Shah functional, the algorithm is often the alternating (or simultaneous) implementation of segmentation and reconstruction. The reconstruction is achieved by image smoothing within each region of segmentation. The level set method is often used in [14] for the segmentation. In recent years, the linear relaxation [15] and the convex relaxation [16] of the Mumford-Shah model have been proposed, which lead to alternatives to the level set method.

### Ii-B Region-based active contours for unsupervised texture segmentation

The original the Mumford-Shah model was formulated for segmentation of piecewise smooth images. Later the model was extended for image texture segmentation. In [17] [7] [8], the Gabor filtering was incorporated in the Mumford-Shah model. In [17], Lee et al. embedded the manually selected 24 Gabor filters into the PS Mumford-Shah model. In [7], Sandberg et al. adopted the Chan-Vese model as well as the maximum difference of feature means as the criterion for Gabor filter selection. Sagiv et al. [8] adopted the framework in [7] with manually selected filters. These frameworks require the textures to be piecewise smooth or flat in the feature space, but the filter selection for unsupervised segmentation cannot ensure this. Although the ad-hoc postprocessing of high dimensional anisotropic diffusion may be applied to ensure the piecewise smoothness, the reasons behind this remain obscure. Kokkinos et al. [9] proposed to apply the Region Competition model to a type of modulation features. This framework assumes that the textures are globally oscillating, as on zebras and tigers, and it still requires filter selection by dominant component analysis (DCA) for parametric texture modeling. However, it is actually unknown what kind of feature or the principle of feature selection can help the segmentation by Mumford-Shah model without supervision. Therefore, the existing combinations of the feature extraction and the Mumford-Shah model is heuristic and they do not provide a unified mathematical model for the segmentation.

There exist other frameworks of texture segmentation based on other similar active contour models with fixed texture feature, e.g. [6] [4] and [5]. For example, in [4] and [5]

the structure tensor has been chosen as the texture feature, and in

[6] the histogram of image values on overlapping patches is chosen as the feature. These frameworks are more likely to form a unified theory of the segmentation. However, the choice of the structure tensor feature or local histogram has not been validated for general textures.

## Iii The principle of segmentation behind the Mumford-Shah model

In this section, we study the segmentation by Mumford-Shah model to understand the rationale of this model for segmentation. We restricted ourselves to the two-phase model. The generalization of the two-phase framework to multi-phase may follow [18], which is out of the scope here.

Let us consider the two-phase Mumford-Shah model in a simplified form as follows.

 (2)

where and are the reconstruction errors corresponding to the subregions and , such that . Since we focus on the region model in this work, we omit discussing about the prior term of arclength for imposing smoothness. Note that, the smoothness term is important for dealing with noisy data. For more detailed discussions on the smoothness term, we refer the readers to [1] [13].

In the PC model, , , where is the image value, are the regional means of the image value. In PS model, the , , where are the smooth restoration of the two image regions (due to the smoothness constraint). The smoothness regularization terms, and are often omitted if and are solved by the normalized Gaussian convolution, such as in [19], which was interpreted as nonparametric regression [12]. Regarding the minimization of this simplified functional (2), we have the following interesting fact upon fixing the error functions and .

###### Proposition III.1
 minH∫ΩE1H+E2(1−H)dxdy=∫ΩminH{(E1−E2)H}dxdy+C (3)

where and is a constant independent of .

This fact tells that optimizing the global assignment of is equivalent to optimizing the assignment locally. The assignment rule for determining the optimal at each pixel location is therefore the following.

 H(x,y)={1,E1(x,y)≤E2(x,y)0,E1(x,y)≥E2(x,y) (4)

The proof is included in the appendix in a separate report. As in the interpretation in [13], Mumford-Shah model tries to achieve clustering based on the measure of the membership defined by . Hence, the proper choice of the measure is essential to the segmentation by Mumford-Shah functional.

## Iv Piecewise linear patch reconstruction

In what follows, we establish the mathematical model and the associated solution for the segmentation of non-smooth image structures.

### Iv-a The model of piecewise linear patch reconstruction

In the original Mumford-Shah model, each region is modeled as a piecewise smooth surface that minimizes the functional energy as follows.

 g∗l=argmingl∫Ωl(I−gl)2+λ∥∇gl∥2dxdy (5)

where ,

is a penalty coefficient. A major reason for this formulation is that the energy can be minimized by solving standard partial differential equations, i.e. the diffusion equation. Besides, the energy is convex. Hence, if the local optimal solution exists it is also global optimal. The problem is that this formulation imposes the smoothness.

To cope with non-smooth image structure, we consider the image patches as linear combination of patch bases, such as in many image appearance models, as follows.

 (6)

with the varying weights for different patches over the image, where is a set of orthonormal bases that fully reconstruct the patch. where and is a square region centered at .

It is obvious that the exact reconstruction can happen regardless of smoothness of the image patch. To find the bases, we may formulate the reconstruction as energy minimization as follows.

 {vlk}∗=argmin{vlk}\bigintsssΩl∥∥ ∥∥p−K∑k=1⟨p,vlk⟩□vlk∥∥ ∥∥2□dxdy (7)

Replacing the image reconstruction error in simplified two-phase Mumford-Shah functional (3) with the patch reconstruction error formulated above, we obtain the formulation of piecewise linear patch reconstruction as follows.

 min{H, {v1k},{v2k}}∫ΩE□1H+∫ΩE□2(1−H)dxdy (8)

where

 (9)

for , where is the size of the patch. From the above, we can note that the error can be computed by convolutions.

From the assignment rule in (4), we know that the patches will be assigned to a region if the set of bases of this region produces the patch reconstruction error smaller than that of the other regions.

### Iv-B The error of segmentation by piecewise linear patch reconstruction

In the above, we have shown that the piecewise linear patch reconstruction is capable of determining the assignment of the patches according to the patch reconstruction error, but we did not answer whether the assignment is correct. In the following, we establish the theoretical foundation of the proposed formulation. Specifically, we show that the error of the assignment, i.e. the segmentation error, is bounded by the patch reconstruction error for a fixed number of, say , bases. Thus, by minimizing the error of reconstruction, we may reduce the error of segmentation.

The assignment rule of (4) enables us to assess the probabilistic analysis of the correctness of the segmentation. We ask whether the segmentation by (4) is consistent with the truth. Specifically, we wish to know if , whether ; and if , whether , where is the true partition. means the set of patches that have their centers in , where . This requires us to analyze the following segmentation error rate.

 εseg=1|Ω1|∑[x′,y′]T∈Ω11[E1(x′,y′)>E2(x′,y′)]+1|Ω2|∑[x′,y′]T∈Ω21[E2(x′,y′)>E1(x′,y′)] (10)

where , are the sizes of the sets and .

The segmentation error rate can be represented by probability, i.e.

for sufficiently large population of and , where is the probability of certain events due to , is the probability of the events due to . We hope the error rate to be small. For our patch reconstruction model, we can deduce the following error bound for segmentation.

###### Proposition IV.1

If the subregion , are sufficiently large, such that , , then the following holds.

 εseg≤∑l=1,2∥□∥∫ΩE□l(x,y)Hldxdy|Ωl|(∥∥p∥∥□−R23−lK2−2q3−l) (11)

where is the patch norm, , , , are constants, and , . Besides, the denominator in the RHS of the bound is positive.

The proof is included in the appendix in a separate report. The above error bound guarantees in theory that the segmentation error rate due to the assignment rule (4) can be decreased by minimizing the reconstruction error with respect to a fixed number of bases. The number of bases will definitely smaller than the number of dimensions of the the patch, since the denominator in the RHS of the bound is positive.

To conclude, we showed that the segmentation error rate is upper bounded by the total patch reconstruction error for a fixed . The minimization of the patch reconstruction error can therefore minimize the segmentation error.

### Iv-C Global optimal linear patch reconstruction

In what follows, we address the optimal linear patch reconstruction.

 {vl∗k}=argmin{vlk}∫ΩEl(x,y,{vli})Hldxdys.t.~{}~{}∀i≠j,⟨vli,vlj⟩□=0,∀k,∥∥vlk∥∥□=1 (12)

where , , and , where is the width of a patch, and we consider square patch in this paper.

A useful identity regarding this formulation is the following.

 argmin{vlk}∫ΩEl(x,y,{vlk})Hldxdy=argmax{vlk}  U({vlk}) (13)

where

 U({vlk})=∫ΩK∑k=1⟨pl,vlk⟩2□Hldxdy=K∑k=1∫Ω∫□|I(x+u,y+v)vlk(u,v)|2dudyHldxdy

according to Eq. (9), and , are orthonormal. This identity can be verified by expanding the squared error.

Before we try to solve the reconstruction problem (12) (or (13) equivalently), we also note that the optimization problem defined in Eq. (12) is a problem of minimizing constrained concave function. The formal statements with their proofs are included in the appendix in a separate report. Such a problem is known to have local optimal solutions [20]. However, we are able to show that the global optimal solution to the linear patch reconstruction problem is attainable. The key result regarding the optimality is the following theorem.

###### Theorem IV.2

Given a set of functions defined as follows,

 (14)

where are all the eigenvectors of , is defined by the following.

 ΛHl(u,v,u′,v′)=∫IHl(x+u,y+v)IHl(x+u′,y+v′)dxdy (15)

where , then the following bound is true.

 U({wk})≤K∑k=1λk=U({ek}) (16)

where is defined in Eq. (13), are the first eigenvalues of .

The proof is included in the appendix in a separate report. The above suggests the matrix eigen-decomposition as our solution to the global optimal reconstruction. However, the eigen-decomposition typically requires converting the image to patches, computing the covariance matrix followed by Singular Value Decomposition (SVD). In our segmentation framework, we require the optimal reconstruction iteratively during the segmentation. Thus, the eigen-decomposition of the matrix obtained from the extracted and labeled patches can be time-storage consuming. Alternatively, we suggest the gradient descent as the solution. The gradient descent equation to minimize the error (

12) is the following.

 (17)

The reasons for this choice are the following.

###### Theorem IV.3

Given the following eigenvalue problem.

 [ΛHl]ek=λkek (18)

where is the eigenvalue, is the eigenvector and , and if there are finitely many, say , eigenvectors of , the following gradient descent procedure for solving (12) (or (13) equivalently) converges to the global optimal solution of (12), with the initial bases for any .

 {while}~{}1≤n≤K {do}  \small{Step} n:⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩% Solve vln by (???)s.t.~{}vln=vln−n−1∑k=1⟨vln,vl∗k⟩□vl∗k,     ∥vln∥□=1{end while} (19)

where for are the solved patch bases.

The proof is included in the appendix in a separate report. Moreover, the convergence of the gradient descent is linear as claimed in the following corollary.

###### Corollary IV.4

If there are eigenvectors of , and the eigenvalues are ordered such that , the rate of convergence for the gradient descent iteration in each step of the greedy procedure defined by Eq. (19) does not exceed . In other words, the procedure converges linearly.

The proof is included in the appendix in a separate report. Besides, the gradient descent does not require computing the matrix or the matrix eigen-decomposition but only the convolutions.

### Iv-D Coupling with piecewise smooth model

It is easy to find two patches, which can be reconstructed equally well (yielding the same residue) by the same set of bases, while having significant difference in intensity. Hence, the proposed linear patch reconstruction model cannot be used for differentiating the image patches that are different only in their illumination and color but similar in their local structure. To cope with both the cases via the same model, we propose to combine our model of linear patch reconstruction with the original piecewise-smooth Mumford-Shah model to form an integrated functional model of segmentation as follows.

 minH∫ΩE⋆1H+E⋆2(1−H)dxdy+ν∫Ω∥∇H∥dxdy (20)

where , for and is a predefined weight. The last term penalizes the complexity of the segmentation.

The numerical solution of to the energy minimization model can be derived from the level set method, in which a signed distance function is used to generate through Heaviside function, i.e. .

 ∂ϕ∂t=−(E⋆1−E⋆2)δ(ϕ)+νdiv(∇ϕ∥∇ϕ∥) (21)

By this equation, it is implied that the reconstruction errors and are fixed when implementing curve evolution. The energy in (20) is minimized by alternatively updating the image partition by solving and updating the reconstruction error in each region by solving and . This energy minimization procedure is in fact the gradient descent method for the separate variables. Hence, the convergence of the procedure is guaranteed.

## V Experiments

### V-a Data preparation and implementation details

Natural textures are examples of the non-smooth structured visual patterns. Evidences show that the texture patches can be reconstructed (modeled or represented) well by patch subspaces, i.e. the eigenfilters, as reported in [21] where the texture classification was investigated more carefully. Therefore, we mainly evaluate our method of segmentation by piecewise linear patch reconstruction on a subset of Brodatz textures [22]. We empirically choose the textures that appear relatively spatially regular with similar size of texture stimuli (textons) for evaluation. Therefore, we obtain a collection of textures: D3, D5-6, D15-22, D24, D34-36, D49, D52-53, D55, D57, D65, D68, D76-77, D79, D81-85, D101-106, which we call the set . Afterwards, we generate two sets of mosaic texture images by paring all different the textures from for evaluation our segmentation method. Each set contains 1260 images. One of the sets is made by the original textures, the other is by the textures with the mean intensity subtracted. The template for paring the textures is shown in Figure 2(a). We also use this template as the ground truth for evaluating the segmentation. Both the sets are challenging for segmentation, especially the second one.

Our segmentation algorithm for minimizing the functional energy (20) can be implemented by alternating an algorithm of image partitioning and the patch reconstruction via either the SVD or the gradient descent procedure in (19) . We adopt the curve evolution governed by Eq. (21) as the image partitioning algorithm in the implementation. All the methods in our experiments are based on the curve evolution for fair comparison. For piecewise smooth image we may choose a small penalty coefficient for the contour length, e.g. in our implementation, for non-smooth image we require a large penalty of the contour length for coping with the randomness, e.g. in our implementation. We use the maximum number of iterations to detect the convergence of the curve evolution algorithms. The maximum iteration number is set to be 600, since it is observed that the curve evolutions in the experiment converge before this iteration number is reached. The convergence of the gradient descent method for patch reconstruction is fast. We set the maximum iteration number to be 5.

### V-B Evaluation of the gradient descent patch reconstruction

In this subsection, we evaluate the proposed gradient descent method for computing the optimal bases. The principle is to compare the error of reconstruction by gradient descent with the error of reconstruction by SVD. We evaluate the gradient descent method for linear patch reconstruction on the Set and Yale face database. We compare the reconstruction errors according to the first 1-20 optimal bases produced by the gradient descent method with the errors according to the bases produced by SVD. The reconstruction errors can be computed by evaluating (12) for the entire image domain. We also present the comparison of the reconstructions by orthogonalized Gabor filter and the SVD solution. We apply the principle of maximum filtering response to select the first Gabor filters from a filter bank of filters. We then orthogonalize the selected filters to compute the energy. We compare the averaged errors but not the total error. The total error is the sum of the reconstruction errors for all the patches in the image. The averaged error is the total error divided by the number of patches.

The results are shown in Figure 3. When visualizing the comparison of the errors, we divide the errors by the value of the maximum averaged errors to form a normalized error. The results by gradient descent are very close to that of the SVD, while the results of Gabor filters do not match the SVD. The initial bases for the gradient descent method are randomly generated.

### V-C Evaluation of the segmentation by linear patch reconstruction

We mainly evaluate the segmentation of our methods, i.e. the SVD and GD based methods, on the two subsets of mosaic textures introduced previously. We compare our methods with the PS Mumford-Shah model [19] [12], the Region-Scalable Fitting (RSF), a.k.a. the Local Binary Fitting (LBF) [23], the Gabor filtering based method [7] which is also adopted in [8] and the local histogram based Chan-Vese model [6], which we denote as HistPC henceforth. The Gabor filtering based method could be viewed as a baseline approach for texture segmentation, and the local histogram-based method is the state-of-the-art approach. We use 8 bases for each region in SVD and GD. We choose the 8 Gabor filters from a bank of 24 filters for each region according to the criterion of maximum filtering response. This filter selection criterion appears like model fitting [24]. The patch size is . The methods are all based on curve evolution. We used a common initial curve for the curve evolutions. The initial curve in the image domain is shown in Figure 2(b). Note that the initial contour crosses the true boundary of the two regions, and the converged contour is expected to outline the region of texture B on the right.

We adopt the pixel-wise segmentation error rate to measure the quality of the segmentation by using the ground truth. The boxplot in Figure 4 shows the segmentation errors for all the methods. We can observe a clear lower error rate by our methods for differentiating the textures that differs only in their structures but not in their intensities. The results by GD and SVD are comparable, but the computational time for the segmentation based on SVD is 0.145 0.048 seconds per iteration for the two datasets, while the computational time for the segmentation based on GD is 0.137 0.006 seconds. This suggests us to use the GD based bases updating scheme for segmentation. The mean error by the local histogram based Chan-Vese active contour for the original textures is small, but the variation of its performance is large. Besides, this local histogram based method is still ineffective for differentiating the different textures having the same mean intensity. We can also observe from the results that Gabor features can deal with textures when the textures differ in their intensities. However, when there is little difference in the intensities of textures, the Gabor features are still powerless. We summarize the quantities in Table I. Besides, we visualize the curve evolutions, as well as the corresponding converged patch bases, for two textured images and a picture from the Berkeley dataset [25] in Figure 5. The picture is composed of different structures, and the piecewise linear patch reconstruction is capable of differentiating the non-smooth structure from the piecewise smooth structure. We may observe that the descriptors are semantic. The results by other models are shown in Fig. 6.

### V-D Segmentation by the coupled model with parameter tuning

To cope with both smooth and non-smooth contents, we may use the coupled model proposed previously. The quality of the result of the segmentation by the proposed coupled model depends on the value of the parameter . It is obvious that the coupled model tends to be the conventional Mumford-Shah model if , and the model tends to be the pure piecewise linear patch reconstruction if . Fig. 8 shows the results of applying the pure piecewise linear patch reconstruction to piecewise smooth images (with noise). We can observe that the optimal bases of the regions are similar. Hence, the model can not deal with such images. We hope to find the such that the coupled model can be used for segmenting piecewise smooth images while its ability for differentiating non-smooth structures is preserved. In our implementation, we run the segmentation with different values on all the mosaic images composed of the original textures (without subtracting the mean off), obtained in the last subsection. We alter the value of from to to obtain the error of segmentation shown in Fig. 7. We observe the steady performance for and , and we observe the clear decay of the performance when changing from to . This means that by selecting , the performance of the coupled model for differentiating non-smooth image structures is almost as good as that of the pure piecewise linear patch reconstruction. We also apply the model with to the piecewise smooth images under noise and the good results are shown in Fig. 9, which indicates that the is already sufficient for segmentation of piecewise smooth images. We therefore chose in the coupled model to optimally balance the two terms in the coupled model.

### V-E One-against-all segmentation

We have so far only considered the presence of two different groups of contents in the image in our formulation, which is known as the two-phase model. Due to the capability of the level set method for handling topology changes, the two-phase model is still capable of partitioning the image into multiple smooth or non-smooth regions of two groups. However, the two-phase model cannot cope with multiple regions of multiple groups. The problem of segmentation of image into multiple different regions can be addressed via the one-against-all strategy. In other words, we may consider a problem of -phase segmentation as subproblems of two-phase segmentation. In each subproblem, the image is to be partitioned into the target region and the background region which is composed of all the other regions.

To evaluate the one-against-all strategy for coping with multiple different groups of regions, we apply this strategy to the mosaic images containing five different textures. The input and output of the segmentation are shown in Fig. 10. The initial rectangular contours are shown in a unique color assigned to a region in one image. The converged contour curves are shown in the same colors in the other image. The optimal bases corresponding to different regions are also visualized. We can observe that these bases well capture the principal structure of the corresponding regions. The error rates of the segmentation are summarized in Table II.

## Vi Discussions and conclusion

### Vi-a Discussions on the limitations

The size of the patches and the number of bases are predetermined, which is empirical. Both of the two properties affect the segmentation. For example, when increasing the patch size, the number of bases needed for small reconstruction error would increase due to the enlarged dimensionality of the patch. Therefore, to achieve a small error of segmentation the number of bases needed has to be large. Larger patch size and more bases would also provide a richer description of the non-smooth structure, which may help segmentation in complex situations. With a small patch size, which means the reconstruction error can be small with few bases, the discrimination might not be clear, since the same set of bases may give similar reconstruction errors to different group of patches in such case.

This work did not really target at natural image segmentation. Specifically, this work explores the mathematical model for addressing an important aspect of the natural image segmentation, i.e. that of coping with non-smooth structures in the segmentation of images which are 2D signals. For natural image segmentation, more sophisticated framework has to be adopted to mimic human vision.

### Vi-B Conclusion

We propose a unified energy minimization model of non-smooth image segmentation without requiring any separate process of feature extraction. The model is the energy minimization of the error of piecewise linear patch reconstruction. The segmentation error rate of the proposed model is proven to be bounded by the patch reconstruction error. The gradient descent method for solving the linear patch reconstruction is proven to be globally optimal under mild conditions on the initialization. The experiments validate our theoretical claim and show the clear supreme performance of our methods over other relevant methods. The linear patch reconstruction in our approach can be viewed as an unsupervised dictionary learning process. Segmentation with features depends largely on the prior knowledge, whilst our approach is mostly data-driven. Our approach is useful when prior knowledge is inapplicable.

## Appendix

• Since and , we have . Thus, for any position we can choose , which is equivalent to at . As a result, for any . Thus, , which completes our proof.

Note that the proof relies on the Axiom of Choice for the continuous domain. The choosing is feasible if we approximate the integral by discretization.

• To prove this proposition, we require the following lemma.

###### Lemma F.1

Given the same condition in Proposition IV.1 then the following holds.

 εseg≤∫ΩE1(x,y)Hdxdy|Ω1|E2(x′,y′)+∫ΩE2(x,y)(1−H)dxdy|Ω2|E1(x′,y′) (A-1)

The proof of this lemma is due to Markov’s inequality. This bound connects the error of segmentation to the error of reconstruction.

Additionally, we require the following fundamental model of functions in the literature of signal analysis, such as wavelets [26] and especially compressive sensing [27] [28].

• Given a (discrete) function , and any subset of fixed system of sorted orthogonal bases such that , then the following bound holds for every ,

 ∥⟨f,vn⟩∥l2≤Rn−q, (A-2)

where , and .

The worst case is when , i.e. .

From Definition of universal energy bound, the linear patch reconstruction error at every pixel could be bounded by using (A-2) as follows.

 (A-3)

where . Substituting (A-3) into the denominator of (A-1), and if the denominator is positive, we complete our proof.

###### Theorem F.2

The optimization problem defined in Eq. (12) is nonconvex.

• First, we shall prove that the objective functional is nonconvex. Taking functional derivatives of the energy twice, we have the following.

 D2Dvli2Ep({vi})=−∫ΩI(x+u′,y+v′)I(x+u,y+v)Hldxdy=−∫ΩI(x+u′,y+v′)HlI(x+u,y+v)Hldxdy=−∫ΩIHl(x+u′,y+v′)IHl(x+u,y+v)dxdy=−ΛHl(u′,v′,u,v) (A-4)

If this functional is convex, will be positive semi-definite. However, we can show that the is actually negative semi-definite as follows.

 −∫Ωh(u′,v′)ΛHl(u′,v′,u,v)h(u,v)dudvdu′dv′=−∫Ω∣∣ ∣∣∫□I(x+u,y+v)h(u,v)dudv∣∣ ∣∣2Hldxdy≤0 (A-5)

Hence, the objective functional is concave. Further, we consider the constraints in terms of orthonormality. Let and both satisfy the constraints. We wonder whether the convex combination also satisfy the constraints. For example, let , where . Then we can verify that , unless , which is generally not true. Therefore, the constraints are also nonconvex. Either of the above two facts can complete our proof.

###### Theorem F.3

In the optimization problem defined in Eq. (12), if is the global optimal solution, then the transformation of

via an orthogonal matrix

is also the global optimal solution.

• First, we define a linearly transformed patch basis

as follows.

 w1(u,v)=[vl1(u,v),vl2(u,v),...,vlK(u,v)]Ta1w2(u,v)=[vl1(u,v),vl2(u,v),...,vlK(u,v)]Ta2⋮wK(u,v)=[vl1(u,v),vl2(u,v),...,vlK(u,v)]TaK (A-6)

where is an mixing matrix. We may simply write the following.

 w1=Va1,w2=Va2,…,wK=VaK (A-7)

where .

Hence, the energy in terms of could be rewritten as follows.

 U({wk})=∫ΩK∑k=1⟨p,wk⟩2□Hldxdy=∫ΩK∑k=1⟨p,Vak⟩2□Hldxdy=∫Ω∥∥AVTp∥∥2Hldxdy=∫Ω(pTV[ATA]VTp)Hldxdy (A-8)

If the mixing matrix is orthogonal, we obtain the following.

 (A-9)

Therefore, if is the global optimal solution, the is also the global optimal solution, which completes the proof.

To prove Theorem IV.3 we require some lemmas.

###### Lemma F.4

The integral operator is symmetric positive semi-definite.

The symmetry is straightforward. The proof of positive semi-definiteness is as in the proof of Theorem F.2.

Hence according to Mercer’s theorem of eigen-decomposition of symmetric nonnegative definite bounded integral operator, we can write

in the form of infinite series of eigenfunctions as follows.

###### Theorem F.5 (Mercer’s representation)

Assuming , then

 ΛHl(u,v,u′,v′)=∞∑i=1λiei(u,v)ei(u′,v′) (A-10)

where are the eigenvalues, is the set of eigenfunctions, and the eigenfunctions form a system of orthogonal basis.

A proof of the above may be found in [29]. Now we are in a position to prove the Theorem IV.3.

• First, we write the gradient descent differential equation for the STEP 1 as follows.

 ∂vl1(u,v,t)∂t=∫ΛHl(u,v,u′,v′)vl1(u′,v′,t)du′dv′ (A-11)

where . The corresponding update equation for the th iteration is the following.

 vl,n+11(u,v)=vl,n1(u,v)  +Δt∫ΛHlvl,n1(u′,v′)du′dv′ (A-12)

Suppose , the update equation for the first iteration could be rewritten as follows.

 vl,11(u,v)=∞∑i=1αiei(u,v)  +Δt∫ΛHl∞∑i=1αiei(u,v)du′dv′=∞∑i=1αiei(u,v)+Δt∞∑i=1λiαiei(u,v)=∞∑i=1(1+Δtλi)αiei(u,v) (A-13)

where we applied Mercer’s representation. Hence, the update equation for the th iteration is the following.

 vl,n1(u,v)=∞∑i=1(1+Δtλi)nαiei(u,v)=(1+Δtλ1)n[α1e1(u,v)    +∞∑i=2(1+Δtλi1+Δtλ1)n≈0αiei(u,v)]≈(1+Δtλ1)nα1e1(u,v) (A-14)

The normalization constraint gives us the desired result.

 limn→∞vl,n1(u,v)=vl,n1(u,v)∥vl1(u,v)∥=±e1(u,v) (A-15)

where we may omit the sign. The proofs for the subsequent steps is hence straightforward, where we only need to replace with

. This setting is due to the Gram-Schmidt process in the constraints. The energy

is therefore the sum of eigenvalues due to Mercer’s representation.

• According to the proof of Theorem IV.3, the th iteration for STEP is the following.

 vl,nk(u,v)=(1+Δtλk)n[αkek(u,v)    +N∑h=k+1(1+Δtλh1+Δtλk)nαheh(u,v)] (A-16)

Let , . Note that the eigenvalues are ordered. Hence, . From these we deduce the following.

 ∥˜vl,n+1k(u,v)−αkek(