Accepted Manuscript
Persistent homology for object segmentation in multidimensional grayscale images Rabih Assaf, Alban Goupil, Valeriu Vrabie, Thomas Boudier, Mohammad Kacim PII: DOI: Reference:
S0167-8655(18)30396-9 https://doi.org/10.1016/j.patrec.2018.08.007 PATREC 7271
To appear in:
Pattern Recognition Letters
Received date: Revised date: Accepted date:
7 October 2017 10 July 2018 5 August 2018
Please cite this article as: Rabih Assaf, Alban Goupil, Valeriu Vrabie, Thomas Boudier, Mohammad Kacim, Persistent homology for object segmentation in multidimensional grayscale images, Pattern Recognition Letters (2018), doi: https://doi.org/10.1016/j.patrec.2018.08.007
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Pattern Recognition Letters Authorship Confirmation Please save a copy of this file, complete and upload as the “Confirmation of Authorship” file.
CR IP T
As corresponding author, I, Rabih Assaf, hereby confirm on behalf of all authors that: 1. This manuscript, or a large part of it, has not been published, was not, and is not being submitted to any other journal. 2. If presented at or submitted to or published at a conference(s), the conference(s) is (are) identified and substantial justification for re-publication is presented below. A copy of conference paper(s) is(are) uploaded with the manuscript.
AN US
3. If the manuscript appears as a preprint anywhere on the web, e.g. arXiv, etc., it is identified below. The preprint should include a statement that the paper is under consideration at Pattern Recognition Letters. 4. All text and graphics, except for those marked with sources, are original works of the authors, and all necessary permissions for publication were secured prior to submission of the manuscript. 5. All authors each made a significant contribution to the research reported and have read and approved the submitted manuscript.
M
Signature: Rabih Assaf Date: 04/07/2018
PT
ED
List any pre-prints:
CE
Relevant Conference publication(s) (submitted, accepted, or published):
AC
Justification for re-publication:
ACCEPTED MANUSCRIPT
Graphical Abstract (Optional) To create your abstract, please type over the instructions in the template box below. Fonts or abstract dimensions should not be changed or altered.
AC
CE
PT
ED
M
AN US
CR IP T
Persistent homology for object segmentation in multidimensional grayscale images Rabih Assaf
In this graphical abstract, we explain the steps of our proposed method for object segmentation in multidimensional grayscale images. From the pixels of an image, we start with a combinatorial representation of the pixels, with each pixel in the image being considered as a vertex or a 0-cell. Then, we build 1-cells or edges from each two adjacent pixels considering the four neighbors of a pixel. We then join squares as 2-cells, where we have 4 adjacent pixels. In the next step, we build a filtration scheme to compute the persistent homology on this combinatorial representation. In this scheme, cells will hold the maximum values of their boundaries, and Ki corresponds to the complex that contains cells with values less than or equal to i. In this example, a homology class is born at the fifth step and dies after at the ninth step. As a result of the persistent homology computation, this class is recorded by a persistence diagram that shows its birth and death times. Highlighting this class allow the object to be segmented.
CR IP T
ACCEPTED MANUSCRIPT
Research Highlights (Required)
To create your highlights, please type the highlights against each \item command.
It should be short collection of bullet points that convey the core findings of the article. It should include 3 to 5 bullet points (maximum 85 characters, including spaces, per bullet point.)
AN US
• Demonstration of the potential of the algebraic topology in performing object segmentation tasks. • New combination between the topological constructions and superpixels of images.
• Insensitivity of the proposed method to continuous deformations and stability against perturbations of the input function.
AC
CE
PT
ED
M
• Capability to be applied on multidimensional images and to perform object segmentation without the need of prior parameters.
ACCEPTED MANUSCRIPT 1
Pattern Recognition Letters journal homepage: www.elsevier.com
Persistent homology for object segmentation in multidimensional grayscale images
a CReSTIC,
Université de Reims Champagne-Ardenne, 51100 Reims, France Sorbonne Universités, UPMC Univ Paris 06, UJF, CNRS, IMT, NUS, 138632, Singapore c Mathematics Department, Holy Spirit University of Kaslik, B.P. 446 Jounieh, Lebanon b IPAL,
ABSTRACT
CR IP T
Rabih Assaf1,∗∗, Alban Goupil1 , Valeriu Vrabie1 , Thomas Boudierb , Mohammad Kacim1
ED
M
AN US
In this paper, we develop a methodology originating from algebraic topology, and we demonstrate its capability of performing multidimensional object segmentation without the need of prior parameters. Persistent homology is a method used in algebraic topology to study qualitative features of data that persist across varying scales. The construction of a topological complex on the image is followed by a filtration scheme that consists of composing a nested sequence of cell complexes on which the persistent homology is computed. The most persistent homology classes are extracted by identifying 1D and 2D chains with large lifespans, which allows salient objects in 2D and 3D images to be segmented and detected. A comparison between this method and other segmentation techniques on a synthetic image shows the advantages of the proposed method. The strength of this technique is reflected in its insensitivity to continuous deformations and perturbations of the input function and in its independence of prior parameters. The results obtained on real and biomedical 2D and 3D images also demonstrate the potential of this method. c 2018 Elsevier Ltd. All rights reserved.
1. Introduction
AC
CE
PT
Object segmentation is a challenging technique that has been considered as a key step in image processing and that remains as a long-standing problem in the field, depending highly on the nature of applications and their types (Singh et al., 2017; Bibiloni et al., 2016). The process of object segmentation consists of separating an object in an image from its other components. A large variety of techniques cover this task. The classical techniques are mostly based on mathematical or statistical concepts such as thresholding (Serna et al., 2014) and watershed transformations (Liu et al., 2010) or on the use of soft-computing techniques (Javed et al., 2016). Segmentation is a difficult and challenging task for many reasons. The main reason is the complexity of the employed algorithms that depend on many parameters and metrics that control the segmentation procedure and cause the lack of generic “off-the-shelf” solutions. The rapid advancements in image processing tools and techniques ensure the development of new methods for image anal∗∗ Corresponding
author: Tel.: +0033613970770; fax: +0-000-000-0000; e-mail:
[email protected] (Rabih Assaf)
ysis. In this context, algebraic topology represents an interesting field that produces powerful techniques for analyzing images (Carlsson, 2014). Among the algebraic topology tools, persistent homology is one of the most powerful yet computationally feasible algebraic techniques for measuring the topological features of functions. It is an algebraic invariant that captures topological features at varying spatial resolutions. Specifically, persistent homology probes topological properties of a space from a set of sampled points such as pixels. It can track the appearance and disappearance of topological features due to the changes in the nested space constructed by an operation called filtration, in which a parameter scale, the intensity of pixels in our case, is increased to detect the changes in the studied space over a range of varying scales. Persistent homology and its statistical results were used in many works on images to classify the studied data. In (Adcock et al., 2014), the authors present a methodology for classifying hepatic lesions using persistent homology, the bottleneck distance, and a support vector machine. The authors in (Dunaeva et al., 2016) use the statistical outcomes of persistent homology to classify patterns of stomach images. In this work, the pixels are binarized depending on a threshold to compute persistent homology. A novel tumor detection tool is presented in (Qaiser
ACCEPTED MANUSCRIPT 2
Fig. 1: An example of a cell complex. Chains R, B and G from the text represent the edges of the red, blue and green paths, respectively.
CR IP T
since we do not need the orientation on our complex, we can work on the Z2 field. For example, in Fig. 1, combinations a, a + d and e + f are 0-chains. In a more general manner, each np P p-chain may be expressed uniquely as c p = αi σi , where σi i=1
represents p-cells, n p is the number of p-cells in the space and αi ∈ Z2 . The set of all p-chains together with the addition operation form a vector space C p . The collection of the (p − 1)dimensional faces of a p-cell σ, which is a (p − 1)-chain, is the boundary ∂ p σ of σ. The boundary of the p-chain is the sum of np P the boundaries of the cells in the chain, i.e., ∂ p c p = αi ∂ p (σi ).
AN US
et al., 2016), using the novel idea of persistent homology profiles. At the level of image segmentation, a split-and-merge algorithm is studied in (Letscher and Fritts, 2007), where the accurate split and merge operations are controlled by two parameters that must be correctly chosen to perform the desired segmentation. A method based on integrating topological prior knowledge into random field image models in practical image segmentation tasks is introduced in (Chen et al., 2011). In (Gao et al., 2013), the authors propose an algorithm based on persistent homology to capture fine structures in CT images. A 3D segmentation technique that relies on many preprocessing steps using zeroth homology groups to find clusters of points of point clouds is discussed in (Beksi and Papanikolopoulos, 2016). In this paper, we start by introducing in a pedagogic manner the formalism of topological spaces needed to build our technique. Then, we show the algebraic potential of homology classes in object segmentation compared to other techniques. The various advantages of the proposed technique over existing methods are demonstrated. It is independent of prior parameters, and its strength is reflected in its invariance to continuous deformations and transformations of the pixel values. Moreover, the topological constructions on superpixels show another aspect of the persistent homology computation on images. All these criteria make the proposed technique perform reasonably well compared to existing methods. 2. Methods
2.1. Combinatorial representation
ED
M
The proposed methodology follows a workflow that consists of transforming the input image into a topological complex called a cell complex using the combinatorial representation. Then, the filtration scheme aims to construct a nested sequence of cell complexes where persistent homology is computed at the end.
AC
CE
PT
Spaces given by point clouds with a notion of neighborhoodlike images can be manipulated as combinatorial spaces or topological complexes that permit the computation of homology classes. There are many types of complexes, such as cubical and simplicial, that represent a triangulation of the space. Since topological spaces are not always a cubical or simplicial complex, we can use a more generalized term called cell complex. We will not provide a detailed explanation of cell complexes, and the reader is referred to (Hatcher, 2001; Edelsbrunner and Harer, 2010) for more details. Following the descriptions, disks, squares, and triangles are considered as 2-cells; edges are considered as 1-cells; and vertices are considered as 0-cells. The relations between these patches are given by their boundaries. The set of all patches and the gluing information provided by the boundary will form the cell complex. A cell τ that is on the bounds of the cell σ is called a face of σ, while σ is a coface of τ. The topological structure of a complex is encoded by vector spaces used by homology theory. We can build a linear structure on top of the complex by defining p-chains c p as formal sums of p-cells. We also choose a coefficient group, and
i=1
Relations between chains are ensured by their boundaries. The boundary operator ∂ p is a linear map between chains of different dimensions, ∂ p : C p → C p−1 . The set of vector spaces C p and the boundary operator ∂ p between them are called a chain complex, denoted as: ∂p
∂ p−1
∂2
∂1
∂0
0 → C p −−→ C p−1 −−−→ . . . −→ C1 −→ C0 −→ 0.
(1)
The fundamental relation of a chain complex determines that the boundary of a boundary is void, that is, ∂ p ∂ p+1 = 0 for all p. This relation will play an important role in the definition of homology groups. 2.2. Filtration of the cell complex After describing the combinatorial representation of the studied data in a cell complex and before developing the idea of filtration, we introduce the concept of homology groups, and we present some examples to illustrate this concept. 2.2.1. Homology Homology is a way to uncover p-dimensional holes in a cell complex. The idea is to find chains that surround holes without being a boundary of a cell of higher dimension. Note that for the needs of the example, B, R and G are the three chains of dimension 1 in Fig. 1, which are in strong blue, red and green lines, respectively. First, a chain that surrounds a subspace or a hole is necessarily without boundaries, i.e., null boundaries. Boundaryless p-chains are interesting and form a subgroup of C p that we call the p-th cycle group Z p . The set Z p of all p-cycles is defined as the subspace of C p of chains without boundary: Z p = {x ∈ C p | ∂ p x = 0} = ker ∂ p .
(2)
Among these cycles, we consider the ones that are boundaries of other chains of higher dimension while remaining on
ACCEPTED MANUSCRIPT 3 the surface. They form a subgroup called the p-boundary group Bp: (3)
Fig. 2: Example on a cell complex on the superpixels together with values over cells for filtration scheme.
an inclusion sequence called filtration: ∅ ⊂ K0 ⊂ K1 ⊂ . . . ⊂ Ki ⊂ . . . ⊂ K.
(5)
We can compute the homology H p (Ki ) for all sublevel sets Ki of (5) to depict the evolution of the number of topological features of input data. However, we lose the information concerning the evolution of each particular cycle. Detecting the evolution of these topological features along the filtration procedure described in (5) to obtain nested sequences of cell complexes can provide more important details, and importantly, it imparts our method with independence from intensity values.
ED
M
AN US
For example, the chain G = (cd)+(dn)+(nm)+(mc) displayed in green in Fig. 1 is a 1-cycle. Indeed, ∂G = (c)+(d)+(d)+(n)+ (n) + (m) + (m) + (c) = 0 as each point appears twice, and we work on the binary field Z2 . Chains B and R in Fig. 1 are also 1cycles. Thus, they belong to Z1 . Additionally, G is a boundary of the 2-chain (cdnm), the dark surface in Fig. 1, and it belongs to B1 . Indeed, ∂2 (cdnm) = (cd) + (dn) + (nm) + (mc) = G. The goal of homology is to discard cycles that are also boundaries because they cannot contain voids. For this purpose, we build an equivalence relation on Z p . The equivalence relation ∼ defined above partitions the p-cycles Z p into a union of disjoint subsets, called homology classes. A chain that belongs to Z p and B p at the same time, meaning that it is a cycle and a boundary of a p + 1-chain and thus does not contain a void, is not considered interesting. In contrast, a chain that belongs to Z p but not to B p because it is a cycle that contains a void and is not a boundary of a higher-dimensional chain is thus interesting. The homology group H p keeps the count of essentially different cycles that are interesting by distributing all cycles into equivalent classes. Thus, an element of H p gathers together these equivalent cycles, and H p will be represented by one cycle. Mathematically, two cycles z1 and z2 ∈ Z p are said to be homologous or equivalent, written as z1 ∼ z2 , if they differ by a boundary, i.e., z1 − z2 ∈ B p . We say that z1 and z2 belong to the same class [z]. We let [z] denote the homology class of z ∈ Z p and define the p-th homology of a space as the quotient of the vector space Z p by the subspace B p , which is a set of homology classes:
CR IP T
B p = {x ∈ C p | ∃y ∈ C p+1 , ∂ p+1 y = x} = im ∂ p+1 .
H p = Z p /B p = {[z], z ∈ Z p }.
(4)
AC
CE
PT
For example, the chain G in Fig. 1 is a cycle that does not surround a hole. It belongs to Z1 and B1 at the same time. On the other hand, chains R and B encircle a hole, and thus, they are not boundaries of any 2-chain. Hence, they belong to Z1 but not to B1 . Additionally, the interesting 1-cycles R and B are inherently the same, and they are said to be homologous because they surround the same hole. The difference between both cycles is a boundary of a 2-chain. In fact, R− B = Γ, where Γ ∈ B1 and Γ = (ab) + (bk) + (kl) + (ln) + (nd) + (de) + (e j) + ( ji) + (ih + (ha). 2.2.2. Filtration The topological features of a complex may be due to noise. Computing the homology of a single topological space is not rich enough to describe the significant topological feature that we look for. We will now show how persistent homology solves this problem. Given a cell complex K, let f : K → R be a nondecreasing function, which means that if σ is a face of τ, then f (σ) ≤ f (τ). The level subcomplexes are then defined to be K(a) = f −1 (−∞, a]. Denoting by ai the values of f on cells of K in increasing order, the corresponding level subcomplexes define a filtration scheme of K. Let Ki = K(ai ); then, we have
2.3. From image to cell complexes An image is a set of points, or pixels, arranged on a rectangular grid of 2D. Thus, the cell complex built from the image depends on its pixels, and it is constructed as shown in the associated graphical abstract: 1. Beginning the construction with vertices or 0-cells. Each pixel in the image will be considered as a vertex. 2. Then, we build the 1-cells or edges from each two adjacent pixels considering the four neighbors of a pixel. 3. We then join squares as 2-cells, where we have 4 adjacent pixels.
Note that the connectivity of the image measured by the homology groups will not be strongly affected by the type of pixel neighborhood chosen. For 3D images, we form the 3-cells or cubes from the corresponding 6 adjacent squares. At the end, the cell complex K will be formed by squares and their components, edges and vertices in the 2D case, or by cubes and their parts, squares, edges and vertices in 3D images. Apart from its typical pixel-grid form, an image can be regarded as a set of superpixels. Superpixels are regions or clusters that group “similar” pixels into meaningful atomic regions that can be used to replace the rigid structure of the pixel grid to simplify the space without information loss and while respecting boundaries (Chen et al., 2017), which will lead to another way of topological construction other than the pixel approach. In the construction of the complex based on superpixels, vertices are superpixels themselves. For visualization purposes, we choose centers of superpixels as representatives of these vertices. For edges, we iterate in each superpixel (S i ) over the 8-connected pixels of each pixel. Two vertices associated with two superpixels are connected by an edge if two neighbor pixels exist in each superpixel. In Fig. 3, the blue dots in the superpixel S 1 are neighbors of the ones in S 2 ; thus, we connect the superpixels S 1 and S 2 by an edge. The violet dots in S 2 and
ACCEPTED MANUSCRIPT 4
Fig. 3: Filtration scheme from K1 to K5 with the associated connected sequence of homology groups and the persistence diagram.
CR IP T
As the intensity increases during filtration, new cells are added until we obtain the entire complex K. For every i ≤ j, we obtain an inclusion map from Ki to K j . This continuous map is transformed by homology to a linear one. A topological feature in Ki detected by the p-th homology is transported into K j . Hence, by tracking the topological evolution during filtration, we obtain a sequence of homology groups that are connected by linear maps induced by inclusions: H p (K0 ) → H p (K1 ) → . . . → H p (Ki ) → . . . → H p (K).
CE
PT
ED
M
AN US
S 3 connect the superpixels S 2 and S 3 by an edge, and the same applies for the green dots in S 1 and S 3 . Regarding the 2-cells or the triangles, we will also iterate in each superpixel (S i ) over the 8-connected pixels of each pixel. If two of these neighbors belong to two other different superpixels S j and S k , we have a triangle of vertices S i , S j and S k . In Fig. 3, each red dot has two neighbor pixels that belong to two other superpixels, which lead to the creation of a triangle connecting S 1 , S 2 and S 3 . Note that the creation of a triangle connecting three vertices in this complex means that there is definitely an edge connecting each two vertices and that 3 vertices connected two by two do not form a triangle if another superpixel is in the middle. We assume that the superpixels are convex subsets of the plane in order to guarantee the homology conservation as provided by the Nerve theorem. In practice, this assumption is verified when the superpixels are numerous enough. In this way, we prevent the loss of information in the built complex that will cover the whole image, and we ensure that all homology classes are captured during filtration. In the three-dimensional case, the superpixels will swing over the depth of the image as each superpixel will belong to different depths. This time, we iterate over the 26-connected pixels of each pixel, and we add the 1-cells and the 2-cells in the same way as in the 2D case. Additionally, a 3-cell is added to the complex when three of the 26 neighbors belong to other superpixels to three other superpixels (S j ), (S k ) and (S l ). In this case, we add a 3-cell connecting these superpixels that represents a tetrahedron. 2.4. Computing persistent homology
AC
To compute the persistent homology of a complex K, we must build a filtration on the complex constructed on the image. We begin by assigning a value to each vertex in K. In the case of construction on pixels, intensity values arise as a natural function to develop this strategy. In superpixels, vertices will hold the mean of intensities in its corresponding superpixel S i . The value of a p-cell will be the maximum of values of its (p − 1)cells that represent its boundary. This value assigning will be the key of the filtration built on the image, as described in (5). Spaces Ki of the filtration are the sublevel sets, and they contain all cells with values less than or equal to i. In this way, we will ensure for the two cases, construction on pixels and superpixels, that the subcomplex criterion is respected, Ki ⊂ K j ∀i < j. The graphical abstract shows an example of the assigning of values for the complex built over a 4× 3 toy image.
(6)
Going from K0 to the whole complex K through these maps will cancel the need for a specific intensity value as a parameter. These maps in (6) will perturb the life of homology classes through the filtration. Considering the step from H p (Ki ) to H p (Ki+1 ), changes can occur: new homology classes that represent holes can be created or already existing homology classes can merge or vanish, signifying filling of the hole. Persistent homology provides tools to track the appearance/birth or disappearance/death of classes in this sequence. It detects homology classes that contain voids and that are present through many steps in the filtration. The illustration in Fig. 3 shows an example of a nested sequence of spaces going from K1 to K5 . Tracking the topological evolution of homology classes by (6) allows their persistence ability to be detected. Homology classes in blue and green appear in K1 and K2 , respectively, and persist to K3 and K4 . We can detect homology classes at each step Ki , but it is more informative to detect their evolution along the filtration scheme. This mechanism will allow separating the non sufficient persistent classes from important ones with no need to fix an intensity value. The persistence of a topological feature is defined as the difference (e.g., in function value or ordering index) between its death and its birth in the construction. Following this approach, lifespans of homology classes can be computed during filtration and thus stored to be used in applications afterward. This pattern is supported by the stability of the lifespans of classes with respect to continuous perturbations of the input function (Cohen-Steiner et al., 2007) and by invariance under continuous deformations of the image, such as stretching, rotation and scaling. Computing persistent homology can be performed using the matrix reduction algorithm described in (Bauer et al., 2017). The history of changes in homology classes can be registered by a persistence diagram. It associates with each class a point whose abscissa is the time of its birth and the ordinate is the level of its death. The lifetime of each homology class is given
ACCEPTED MANUSCRIPT 5
3. Results 3.1. Synthetic 2D image In this part, we created a 2D synthetic image, and we applied our method to perform a comparison. Then, we based our technique on topological constructions applied on the superpixels of the same image.
(c) All MSER detected objects
(b) Otsu thresholding
(d) First channel of the image in (c)
(e) Second channel of the image in (c)
(f) Third channel of the image in (c)
(g) Objects in volume range 1001000
(h) Persistence diagram of first dimension homology classes given by our method
(i) Object segmentation by keeping classes above (t1 )
(j) Object segmentation by keeping classes above (t2 )
AC
CE
PT
ED
M
AN US
3.1.1. Object segmentation on pixels and comparison The objective here is to perform a comparison among our method, classical methods that perform global thresholding, and segmentation methods based on tree-like structures. For this purpose, we created a synthetic image with a size of 240×120 pixels, as shown in Fig. 4(a), that contains objects of different sizes, shapes and intensities as well as some noise. Some of these objects are overlapping other objects, and all objects are overlaid on a slightly variable background. First, Fig. 4(b) shows the results using the Otsu thresholding, which is the best known global thresholding method. Note that such methods are not able to segment the small objects that are overlaid in the larger one in the seventh object nor the large objects that are the base of the smallest ones in the bottom right because of the use of a threshold. To overcome the use of a threshold, the maximally stable extremal regions (MSER) method described in (Gul-Mohammed et al., 2014) uses different thresholds to iteratively detect objects and store them in a hierarchy or tree-like structure. Hence, for the lowest threshold, only background is detected. With increasing threshold, some objects are detected, and if their volumes fall into a specific range, these objects are extracted and stored in the tree-like structure. If the same object is detected with different thresholds, this object is stored in a branch of the tree. Then, the best object from this branch is extracted. For MSER, it is the object that has the minimal volume variation between two thresholds. In the case where an object will be partitioned into smaller objects at a higher threshold, a division will be created in the treelike structure. The best detected objects for each branch are computed and displayed in 3 channels colored in red, green and blue in the final resulting image. Using MSER, if we consider the entire tree structure, we obtain the result in Fig. 4(c). This image corresponds to the superposition of the objects detected in the three branches of the tree-like hierarchy that we obtain for the synthetic image. In the first channel shown in Fig. 4(d), as shown, this method identifies some of the objects with different sizes. However, this
(a) Synthetic image
CR IP T
by the vertical distance between its corresponding point and the diagonal. Thus, we can recover the lifespans of homology classes from this persistence diagram. The more these classes are present during filtration, the further they are from the diagonal. Points close to the diagonal represent homology classes that could not resist or persist for a long time against the variations in the space and will thus be considered as noise. This motto depends highly on the application which is segmentation of salient objects in our case. In Fig. 3, the persistence diagram records the homology classes of the associated filtration scheme. Clearly, the classes represented in green and blue appear at levels 1 and 2 and disappear later at levels 4 and 5, respectively.
Fig. 4: Comparison between the Otsu global thresholding technique, the MSER method and our method.
ACCEPTED MANUSCRIPT
Fig. 5: (a) Persistence diagram of homology classes of the first dimension after Gaussian smoothing of Fig. 4(a), (b) object segmentation of the smoothed image after highlighting classes above the threshold (t2 ), (c) persistence diagram of homology classes of the first dimension after rotating the image in Fig. 4(a), (d) highlighting the same classes with the threshold (t2 ), (e) persistence diagram of homology classes of the first dimension after scaling the image in Fig. 4(a), and (f) object segmentation of the scaled image after highlighting classes above the threshold (t2 ).
AN US
method also identifies small dots corresponding to noise and a part of the background on the bottom right. In the second channel shown in Fig. 4(e), which contains a larger object emerging from the merging of smaller objects, two other objects that are the base of the smallest ones are highlighted. In the third channel shown in Fig. 4(f), these two objects and the object with the largest size from the first image are embedded into another part of the background. This method falls to identify only the interesting objects since the noise is also detected. Moreover, this method depends on the background level, and the objects are not extracted correctly. One way to segment overlapping objects with their bases and to avoid noise is to use a specific range of volumes as a parameter. In this manner, it will achieve segmentation depending on the given corresponding information. Only objects that are consistent with the corresponding range volume will be detected. For example, Fig. 4(g) shows the objects that correspond to the volume range of 100-1000. A part of the background is not eliminated, and not all objects can be identified. Moreover, if the volume parameter is lowered to extract the small objects, then noise will also appear. If the volume parameter is sufficiently large, we cannot extract the small objects that are overlaid over other objects.
CR IP T
6
AC
CE
PT
ED
M
Considering our method, a topological complex is first built directly on the pixels of the image in Fig. 4(a), as we explained before. In an analogous approach to the tree, the filtration scheme is constructed. Points, edges and 2-dimensional cells are added to the complex as the intensity increases until we obtain the whole complex. Computing the persistent homology on this complex presents the persistence diagram of the homology classes of dimension 1. These homology classes of dimension 1 consist of chains of 1-cells or edges, which is suitable for object segmentation in 2D cases. This persistence diagram, shown in Fig. 4(h), is the main output of our method since it indicates the births and deaths of homology classes of first dimension. Each point carries all the information necessary for homology classes. The importance of the points and what they represent is related to their distance from the diagonal in our application. Thus, points that are close to the diagonal have a small lifespan, meaning that they correspond to noise, while those that are far from the diagonal are interesting. Note that having the construction of the complex, we can find for each point in the persistence diagram the x, y coordinates of the corresponding class. For example, the homology class of the first dimension represented with a * in Fig. 4(h), which is a chain of edges or 1cells, corresponds to the class that identifies the object shown in blue in Fig. 4(i), the homology class represented with * to the one corresponding to the object in red, and so forth. Highlighting the homology classes in the persistence diagram then allows selecting interesting objects by imposing a parallel threshold to the first diagonal. For example, the homology classes that represent points above the parallel (t1 ) are highlighted in Fig. 4(i) by their corresponding colors. Pulling down the parallel to (t2 ) will allow highlighting other homology classes that are represented by all points above (t2 ) as in Fig. 4(j). Note that this threshold can be selected using the gap in a histogram that represents the lifespans of homology classes. Other ways to optimize this threshold are to use statistical concepts such as en-
tropy (Atienza et al., 2016) or landscapes (Bubenik, 2015), but these concepts are more suitable for Vietoris-Rips complexes. Our method allows computing all the homology classes at the same time and highlighting them by simply changing a posteriori threshold. The method is insensitive to continuous transformations of the images, such as rotation, or scaling, and to continuous perturbations of the pixel values. For demonstration purposes, Fig. 5(a) shows the persistence diagram of homology classes of the first dimension after Gaussian smoothing of the synthetic image in Fig. 4(a) by a standard deviation for a Gaussian kernel equal to 1. We also show in Fig. 5(c) and (e) the persistence diagram of the first dimension after rotating the synthetic image by 30 degrees and scaling the image to 150% and 120% of the initial width and height, respectively. The same salient objects as previously are segmented. Indeed, the corresponding homology classes are colored stars in the smoothed, rotated and scaled images in Fig. 5(b), (d) and (f), respectively. The same important homology classes are obtained after Gaussian smoothing, rotating or scaling the image. Some slight changes have occurred in the lifespans of classes, and some noisy classes of null lifespans are present due to the differences in the intensity values caused by these perturbations. However, these changes do not affect the result of the threshold (t2 ) highlighting the same classes in the original, smoothed, rotated and scaled images. Our method also showed its efficiency in the segmentation of irregular shapes. An example is shown in Fig. 6 on the synthetic image but containing some irregular shapes this time. The stars in the persistence diagram in Fig. 6(a) correspond to the homology classes that are highlighted in Fig. 6(b). Note that, the current method can handle images counted in hundred of million of pixels thanks to the efficiency of the persistent homology (Otter et al., 2017).
ACCEPTED MANUSCRIPT 7
Fig. 7: (a) The persistence diagram of homology classes of the first dimension and (b) object segmentation by keeping the classes represented by * after presegmentation of the image in Fig 4(a) to 2000 superpixels rather than 240×120 pixels.
Fig. 9: (a) Persistence diagram for homology classes after its presegmentation to 5000 superpixels rather than 513x282 pixels. (b) Segmentation results of image of dots.
However, we found the same number of homology classes of the first dimension, represented by *, as well as the same distribution of the topological noise. Imposing the same threshold (t2 ) as previously, the same objects as in Fig. 4(j) are segmented, as shown by the results in Fig. 7(b). We observe that the outlines of the segmented objects are the roughest, which is explained by the fact that the 0-cells are the centers of superpixels that are more spaced out than the raw pixels. However, the computation time is substantially shorter and the memory is saved. For the case of raw pixels, the computation time is 7 ms for 114481 cells, whereas it is 800 µs for 17278 cells in the case of superpixels.
CE
PT
ED
M
AN US
The homology classes are not optimized as MSER does by choosing only objects with a specified volume range. This can be done as a postprocessing by classical methods such as active contour or clustering techniques to obtain a pertinent segmentation for the specific application. Meanwhile, methods that rely on the application of extinction values on hierarchical image representation, described in (Xu et al., 2017) and (Salembier et al., 1998), aim to compute a hierarchy of image segmentation from fine to coarse, unlike our method that provides a single segmentation of the image after the computation of persistent homology. Although the hierarchical image representation depends on the notion of min/max trees that are based on the inclusion relationship of connected components of upper level sets or lower level sets, our method deals with these level sets as a way to construct a nested sequence of topological complexes where cells of different dimensions are added to finally have the entire complex. Moreover, candidate regions in these methods are characterized with attributes that describe them. Extinction values are used to measure the persistence of increasing attributes of these regions, while our notion of persistence is concluded directly from the resistance of homology classes to increasing spaces and does not depend on attributes.
Fig. 8: (a) Persistence diagram of homology classes of first dimension, and (b) object segmentation of the natural complex image.
CR IP T
Fig. 6: (a) Persistence diagram of homology classes of first dimension, and (b) object segmentation of the image containing some irregular shapes.
AC
3.1.2. Object segmentation on superpixels Now, we benefit from the advantages provided by superpixels, particularly in the reduction of the constructed spaces, to perform object segmentation, and we apply our methodology on images presegmented to superpixels using the SLIC technique (Achanta et al., 2012). We apply the method on the synthetic image in Fig. 4(a). As described previously, we build our topological complex on the set of superpixels, and then we apply the steps of computing persistent homology. The births and deaths of these homology classes of the first dimension extracted from the set of superpixels are shown in Fig. 7(a). Some slight differences are observable when comparing them with the homology classes shown in Fig. 4(h). These differences are because the values of 0-cells are not the same since a superpixel realizes an averaging of several raw pixels.
3.2. Real 2D application First, we consider the image of size 500×305 pixels shown in Figure Fig. 8(b). This image is quite complex since the interesting objects, the four fish, are superimposed on a highly variable background, strong noise being also present. Building the topological complex directly on the pixels and computing the persistent homology, we obtain the persistence diagram of homology classes of first dimension shown in Fig. 8(b). Imposing the parallel threshold (t1 ) to the first diagonal in order to keep the four classes represented with red *, we segment only the interesting objects, the fish. Decreasing the threshold to (t2 ) we will allow highlighting other objects, the fish wings. Next, we consider employing our method on superpixels estimated on the image of size 513×282 shown in Fig. 9(b). The dots in this image do not have the same sizes, forms or intensities, and some of them are difficult to distinguish by the naked eye. Additionally, the spaces between these dots contain abundant noise. The persistence diagram in Fig. 9(a) demonstrates 28 homology classes that are far from the diagonal, which represent the dots. We highlight in Fig. 9(b) the homology classes that represent the dots of the image. Depending on the nature of the input image, this methodology can also be applied using the superfiltration procedure by starting from the maximum intensity in a decreasing way and in an analogous way to subfiltration as described in section 2. Note that once we identify
ACCEPTED MANUSCRIPT 8 References
AN US
the homology classes, we can store the information inside them. Applying a simple postprocessing technique to each content of the homology classes will allow providing a more refined segmentation, while these same techniques fail on the entire image.
CR IP T
Fig. 10: Results of segmentation of a 61×249×308 grayscale biomedical image of nucleolus at depths z = 0, 15, 30, 40, 50 and 60.
Achanta, R., Shaji, A., Smith, K., Lucchi, A., Fua, P., Susstrunk, S., 2012. Slic superpixels compared to state-of-the-art superpixel methods. IEEE Trans. Pattern Anal. Mach. Intell. 34, 2274–2282. Adcock, A., Rubin, D., Carlsson, G., 2014. Classification of hepatic lesions using the matching metric. Computer Vision and Image Understanding 121, 36 – 42. Atienza, N., Gonzalez-Diaz, R., Rucco, M., 2016. Separating topological noise from features using persistent entropy, in: Software Technologies: Applications and Foundations, Springer International Publishing, Cham. pp. 3–12. Bauer, U., Kerber, M., Reininghaus, J., Wagner, H., 2017. Phat persistent homology algorithms toolbox. Journal of Symbolic Computation 78, 76 – 90. Algorithms and Software for Computational Topology. Beksi, W.J., Papanikolopoulos, N., 2016. 3D region segmentation using topological persistence, in: Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems (IROS), pp. 1079–1084. Bibiloni, P., Gonzàlez-Hidalgo, M., Massanet, S., 2016. A survey on curvilinear object segmentation in multiple applications. Pattern Recognition 60. Bubenik, P., 2015. Statistical topological data analysis using persistence landscapes. J. Mach. Learn. Res. 16, 77–102. Carlsson, G., 2014. Topological pattern recognition for point cloud data. Acta Numerica 23, 289–368. Chen, C., Freedman, D., Lampert, C.H., 2011. Enforcing topological constraints in random field image segmentation, in: CVPR 2011, pp. 2089– 2096. Chen, J., Li, Z., Huang, B., 2017. Linear spectral clustering superpixel. IEEE Transactions on Image Processing 26, 3317–3330. Cohen-Steiner, D., Edelsbrunner, H., Harer, J., 2007. Stability of persistence diagrams. Discrete Comput. Geom. 37 (2007) . Dunaeva, O., Edelsbrunner, H., Lukyanov, A., Machin, M., Malkova, D., Kuvaev, R., Kashin, S., 2016. The classification of endoscopy images with persistent homology. Pattern Recognition Letters 83, Part 1, 13 – 22. Geometric, topological and harmonic trends to image processing. Edelsbrunner, H., Harer, J.L., 2010. Computational topology:an introduction. American Mathematical Society. . Gao, M., Chen, C., Zhang, S., Qian, Z., Metaxas, D., Axel, L., 2013. Segmenting the papillary muscles and the trabeculae from high resolution cardiac ct through restoration of topological handles, in: Information Processing in Medical Imaging, Springer Berlin Heidelberg. pp. 184–195. Gul-Mohammed, J., Arganda-Carreras, I., Andrey, P., Galy, V., Boudier, T., 2014. A generic classification-based method for segmentation of nuclei in 3d images of early embryos. BMC Bioinformatics 15, 9. Hatcher, A., 2001. Algebraic topology. Cambridge Univ. Press . Javed, U., Riaz, M.M., Ghafoor, A., Cheema, T.A., 2016. SAR image segmentation based on active contours with fuzzy logic. IEEE Transactions on Aerospace and Electronic Systems 52, 181–188. Letscher, D., Fritts, J., 2007. Image segmentation using topological persistence. Computer Analysis of Images and Patterns . Liu, Z., Li, W., Shen, L., Han, Z., Zhang, Z., 2010. Automatic segmentation of focused objects from images with low depth of field. Pattern Recognition Letters 31, 572 – 581. Otter, N., Porter, M.A., Tillmann, U., Grindrod, P., Harrington, H.A., 2017. A roadmap for the computation of persistent homology. EPJ Data Science 6, 17. Qaiser, T., Sirinukunwattana, K., Nakane, K., Tsang, Y.W., Epstein, D., Rajpoot, N., 2016. Persistent homology for fast tumor segmentation in whole slide histology images. Procedia Computer Science 90, 119 – 124. 20th Conference on Medical Image Understanding and Analysis (MIUA 2016). Salembier, P., Oliveras, A., Garrido, L., 1998. Antiextensive connected operators for image and sequence processing. IEEE Transactions on Image Processing 7, 555–570. Serna, A., Marcotegui, B., Decencière, E., Baldeweck, T., Pena, A.M., Brizion, S., 2014. Segmentation of elongated objects using attribute profiles and area stability: Application to melanocyte segmentation in engineered skin. Pattern Recognition Letters 47, 172 – 182. Singh, R.P., Gupta, S., Acharya, U.R., 2017. Segmentation of prostate contours for automated diagnosis using ultrasound images: A survey. Journal of Computational Science 21, 223 – 231. Xu, Y., Carlinet, E., Géraud, T., Najman, L., 2017. Hierarchical segmentation using tree-based shape spaces. IEEE Transactions on Pattern Analysis and Machine Intelligence 39, 457–469.
PT
ED
M
3.3. Real 3D images To demonstrate the multidimensional aspect of our method, we present in Fig. 10 a segmentation of a 3D biomedical grayscale image of a nucleolus using persistent homology computation on 5000 superpixels. The homology classes of the second dimension are chains of 2-cells or triangles, which make them suitable for object segmentation in 3D images. In Fig. 10, we highlight the seven most persistent homology classes formed by 2-chains. We have here four essential objects. The homology classes that represent them and mentioned in increasing order of lifespans are highlighted in red (the bottom right), green, blue, and red (the top left). Our tool is also capable of detecting the interesting objects inside them.
AC
CE
4. Conclusion The purpose of this work is not only to introduce a new way of performing object segmentation but also to explain a methodology arising from algebraic topology tools and persistent homology. The relation between this methodology and image properties allowed us to define generators of homology classes for the object segmentation process. The power of this method is revealed in its independence of prior parameters and its insensitivity to continuous perturbations of the image. A comparison on a synthetic image demonstrates the advantages of our method over some of the existing ones. The construction of a topological complex on images presegmented to superpixels demonstrated the capability of persistent homology and showed its multidimensional efficiency. Moreover, our methodology can be applied on images of reversed intensities and can be extended to n-D dimensional images. For future work, we would like to involve the multidimensional persistence in our method by combining many criteria in the process of filtration to achieve the segmentation of more complex images. Acknowledgment: The research in this paper was financially supported by the CNRS-L /USEK funding program.