A multi-parameter persistence framework for mathematical morphology

The field of mathematical morphology offers well-studied techniques for image processing and is applicable for studies ranging from materials science to ecological pattern formation. In this work, we view morphological operations through the lens of persistent homology, a tool at the heart of the field of topological data analysis. We demonstrate that morphological operations naturally form a multiparameter filtration and that persistent homology can then be used to extract information about both topology and geometry in the images as well as to automate methods for optimizing the study and rendering of structure in images. For illustration, we develop an automated approach that utilizes this framework to denoise binary, grayscale, and color images with salt and pepper and larger spatial scale noise. We measure our example unsupervised denoising approach to state-of-the-art supervised, deep learning methods to show that our results are comparable.


Introduction
Computational topology and the field of topological data analysis offer powerful tools for analyzing structure in data [1,2,3,4,5].Persistent homology, in particular, has offered a means of measuring topological features across a filtration, or sequence of structures built from the data.A one-filtration is a collection of nested sets where the inclusion map enables the tracking of topological information from one set to the next.A multifiltration extends this notion to an indexed collection of sets satisfying an inclusion relation between a pair of sets whenever their indices are related under a specified partial order.This allows for the construction of many structures related to image or point cloud data, where appropriate inclusion relationships allow for the tracking of topological information across the structures.
To date, most studies using multiparameter persistence, the extension of persistent homology to multifiltrations, have focused on point cloud data rather than the cubical/digital image data we study here [6,7,8,9,10,11].A standard filtration for grayscale digital images is the sublevel set filtration obtained via thresholding.In this work, we use arXiv:2103.13013v1[cs.CG] 24 Mar 2021 morphological operations to construct new filtrations, yielding erosion (13), dilation (14), opening (15), and closing filtrations (16) and their variants (17), (18), (19).We then show that under mild assumptions, combinations of these operations form multiparameter filtrations.This establishes a multifiltration framework for the analysis of digital images where features appear at different spatial scales, including noisy images in which the noise is smaller in spatial scale than the underlying structure we wish to uncover.We then demonstrate that it is possible to use this framework and persistent homology to extract information about underlying structure in the images as well as to automate the production of a denoised image.
Digital images may naturally be thought of as functions on Z m , with sets of pixels for regions of interest in the image given as subsets of Z m .In [12], the authors provide an interesting method for smoothing shapes of objects in Z m , with m = 2, 3, which preserves homotopy structure.To achieve this, the authors give the definition of homotopic equivalence of discrete images and construct homotopic thinning/thickening operations for shape smoothing.However, in image denoising tasks, one often aims to remove small scale features that arise due to noise in the image, thus sometimes dramatically changing the topology of the image.Therefore, in our approach, we do not impose homotopic equivalence and instead adopt a goal of intentionally changing the homological type of rendered structures in order to optimize topological and geometric accuracy by removing features most likely due to noise.
Filtrations and persistence lend themselves well to automation.In [13], Chung and Day used persistent homology to track structure in the sublevel set filtrations, developing an automated method for extracting topological measurements and thresholding grayscale images.In this work, we focus on building an algebraic topological framework for the application of the morphological image processing operations of erosion, dilation, opening, and closing.These are well-developed operations for cleaning images by removing small scale features while keeping the remainder of the image relatively constant [14,15,16].When combining morphological operations, the dimension of the constructed multiparameter filtration grows rapidly in the numbers of operations and utilized structuring elements (see Section 2).Furthermore, thresholding may be combined with opening and closing to form a yet larger multifiltration for studying grayscale and, by extension, color images.As illustration, in Section 5, we use opening, closing, and thresholding to construct a multifiltration that we use for denoising images with salt and pepper noise, providing sample results for binary, grayscale, and color images.
In what follows, we introduce necessary definitions and properties for morphological operations (Section 2) and persistent homology (Section 3).We then use morphological operations and thresholding to build multifiltrations in Section 4, presenting our main results in Theorem 1 and Corollary 1.This presents the necessary framework to combine morphological operations and, if appropriate, thresholding, in a single multiparameter filtration for analyzing digital data.We show some illustrative uses of these filtrations in Section 5.

Background on Mathematical Morphology
Mathematical morphology is a field that provides theoretical and practical techniques for processing digital images [17,18,19,20,21].In later sections, we will focus on using morphological operations to measure and track topological features in images.Here, we focus on establishing properties of the operations that are necessary to this topological approach.Two morphological operations that will be of particular interest in what follows are dilation and erosion.In a binary image, dilation enlarges features in the pixel subset for a specified value (e.g. 1, 255, or 'black') while erosion erases small, isolated features in the pixel subset (see [14,15,16] and references therein).
We denote integers, natural numbers, and real numbers by the standard notation Z, N, and R, respectively.We use Z ≥0 to denote N ∪ {0}, the set of all non-negative integers.The symbol R ≥0 represents the set of all non-negative real numbers.Elements in Z m are denoted by boldface letters e.g.u ∈ Z m to distinguish vectors and scalars.We use this notation to build towards formalizing operations on digital images, which we define as follows.Definition (Section 1.1.2.1.[15], p. 6).Let m ∈ N be a positive number, an m-dimensional (digital) image on pixel/voxel set P ⊆ Z m is a non-negative function g : P → R ≥0 .If the range of the function g is {0, 1}, then g is called a binary image.Otherwise, we refer to g as a grayscale image.The set of all images on P , denoted I P , is defined as In practice, one usually considers a rectangular image whose domain can be expressed as The following discussion focuses primarily on binary images g : P → {0, 1}, where a value of 1 means the pixel is white and a value of 0 means the pixel is black, and 8-bit grayscale images g : P → {0, 1, ..., 255} where 225 is white and 0 is black.When feasible, we consider general images g : P → R ≥0 so that properties and theorems stated in the paper hold in this general case.This setting is also convenient when considering re-scaling of pixel values of images with different range sets.
We now establish the following partial order on the space of images and define image and preimage subsets.Given two images f, g : P → R, we say that f ≤ g if and only if f (p) ≤ g(p) for all p ∈ P.
For functions f : S → T , A ⊆ S and B ⊆ T , the sets f (A) := {f (a) | a ∈ A} and f −1 (B) := {s ∈ S | f (s) ∈ B} are the image of A and the preimage of B under f .As an abbreviation, if B = {b} is a singleton set, we write f −1 (b) instead of f −1 ({b}) to denote f −1 (B).Definition (Section 3.1 [14] p. 64).For m ∈ N, a structuring element is a specified finite set B satisfying 0 Remark.In mathematical morphology, structuring elements defined here are called flat structuring elements [14].In certain applications, a non-flat structuring element B is defined as a function from a finite subset B of Z m to R ≥0 which records weighted values for elements in B. In this case, a flat structuring element can be viewed as a characteristic function χ B on a finite set B ⊆ Z m .In this paper, all structuring elements we consider are flat.
Structuring elements will be used to define local windows over which pixel values are considered during processing operations.This requires the Minkowski sum and difference for subsets A, B ⊆ Z m , defined respectively as Since 0 ∈ B, one may think of x + B as the B-neighborhood of x in Z m .If g ∈ I P is a binary image and g(x + B) = {1}, then the B-neighborhood of x is contained in the white region of the image.On the other hand, if g(x + B) = {0, 1}, then the B-neighborhood of x intersects both the white and the black sets in the image.
In this work, we consider four fundamental morphological operations: erosion, dilation, opening, and closing.We next review their formal definitions.Definition (Equations (1.6) and (1.7) in [15] p. 10).For g ∈ I P and structuring element B ⊆ Z m , the erosion of g via B is an image B (g) ∈ I P defined by Similarly, the dilation of g via B is an image δ B (g) ∈ I P defined by Since 0 ∈ B and B is finite, (x + B) ∩ P and (x − B) ∩ P are non-empty, finite sets whenever x ∈ P .Therefore, B (g) and δ B (g) are well-defined.
Erosion and dilation may now be composed to define the operations of opening and closing.Definition (Section 1.2.1 [15], p. 12).Let P, B ⊆ Z m be a pixel set and structuring element respectively.Then opening and closing operations via B, denoted by O B and C B respectively, are functions O B , C B : I P → I P defined as The opening and closing operators may be used to remove structure that is smaller than the scale prescribed by B while minimizing distortion of larger scale features [22,23,24,14].It is clear that if B = {0}, then B = δ B = id I P where id I P : I P → I P denotes the identity function and O B = C B = id I P .
We conclude this section by reviewing some basic properties of these morphological operations.We will use these properties to establish our main result.
Proposition 1 (Properties 3.4 [14], p. 71).Let f, g ∈ I P be images and B ⊆ Z m be a structuring element.If f ≤ g, then the following inequalities hold Proposition 1 states the increasing property, that is, for a fixed structuring element, the basic morphological operations preserve the ordering relation on images.
In Definition 2, we define images as functions.Our main focus in this work is to construct a filtration or collection of sets ordered by set inclusion.We do this for image sublevel sets, i.e. subsets of the pixel set P corresponding to pixels with image values at or below a prescribed threshold value.When we consider binary images, the filtration property of sublevel sets is naturally related to the increasing property as shown in the following proposition.Proposition 2 (Principle 11.1.1[14], p. 318).Let f, g ∈ I P be images.If f ≤ g, then g −1 (0) ⊆ f −1 (0).In addition, if f, g : P → {0, 1} are binary images, then f ≤ g if and only if g −1 (0) ⊆ f −1 (0) and, similarly, f ≤ g if and only if There are many ways to produce a binary image from a grayscale image.Global thresholding of grayscale image g : P → R ≥0 via threshold value t produces the binary image Note that g −1 t (0) = {x ∈ P : g(x) ≤ t}.This set, g −1 t (0) is the t-sublevel set of g.In general, the operations of erosion, dilation, opening, and closing do not commute.However, these four operations do commute with the operation of global thresholding as follows.Proposition 3 (Proposition 1, [25]).For m ∈ N and pixel set P ⊆ Z m , consider the image g ∈ I P .For each threshold t ∈ R ≥0 we define τ t : I P → I P by g → g t i.e., For any structuring element B ⊆ Z m , the following diagrams are commutative: Finally we note the relationship between the partial order on grayscale images and the partial order on the corresponding thresholded images.Lemma 1 (Lemma 1, [25]).For images f, g ∈ I P , f ≤ g if and only if f t ≤ g t for every t ∈ R ≥0 .

One-Parameter Filtrations and Persistent Homology
In this section, we show that the partial order results for morphological operations presented in Section 2 naturally yield the structure necessary for computing persistent homology.
Persistent homology, a foundational tool in the field of topological data analysis (TDA), measures and tracks topological features.It relies on having a one-parameter filtration, a sequence of nested sets.The goal of this section is to introduce two new filtrations based on morphological operations and define and illustrate the meaning of persistence diagrams based on these filtrations.Topological features of interest include connected components (or individual connected pieces of the set), onedimensional holes (holes in 2d or tunnels in 3d) and two-dimensional holes (cavities in 3d).Higher dimensional holes may appear in higher dimensional data, but for illustration, we will focus on two and three dimensional data sets here.The framework we present throughout this work applies to higher dimensional data and holes as well.For the data sets we study, cubical homology may be summarized using Betti numbers.Betti numbers, β k , count holes of various dimensions.More specifically, given a binary image f , if we consider the set of black pixels, X := f −1 (0), then β 0 (X) is the number of connected components, β 1 (X) is the number of 1-dimensional holes, or tunnels, β 2 (X) is the number of 2-dimensional holes, or cavities, etc.They are computed using algebraic structure defined by the cubical structure of X and there are now efficient software packages for performing these calculations.See, for example, [26] and references therein for a discussion of the mathematical theory behind the definition and computation of Betti numbers as well as their interpretation as direct counts of topological features.
Persistent homology extends the topological measurement offered by Betti numbers across a filtration.A one-parameter filtration is a sequence of sets {X i } i∈A , with indexing set A ⊂ Z, satisfying For ease of notation, we will often write {X i } when the indexing set has already been specified.For a one-parameter filtration, persistent homology records the birth and death coordinates at which a given topological feature first appears and first disappears respectively in cubical sets.That is, given a one-parameter filtration {X i } i∈A of cubical sets, X i , a feature with birth/death coordinates (b, d), does not exist in the sets X i with i < b, appears first in X b and persists through all sets X j with b ≤ j < d and disappears in X d .Like Betti numbers, birth/death coordinates are computed using algebraic structure attached to cubical sets, in this case using the inclusion operator to match some features in X n to their preimages in X n−1 .
The collection of birth/death pairs for all topological features, labeled by dimension of the feature, is called a persistence diagram.For a given one-parameter filtration {X i }, the full persistence diagram is the collection of all birth/death pairs and is denoted by P({X i }), with the k-th persistence diagram, P k ({X i }), denoting the subset of pairs measuring kdimensional holes.By this convention, P( For ease of notation, we represent these persistence diagrams by P (respectively P k ) when the filtration {X i } is understood.Betti numbers may be extracted from persistence diagrams as In other words, β k (X m ) counts the number of points in the k-th persistence diagram whose birth/death coordinates indicate that they are present in set X m .By extension, we also write P(m) := ∪ k P k (m).Furthermore, a given feature's lifespan, l = d − b, measures the length of the interval of set indices over which the feature persists.Birth/death coordinates and corresponding lifespans allow us to study the robustness of the feature with respect to changes in the index m.
For more information about persistent homology, see e.g.[2,4], and references therein.In summary, given a oneparameter filtration of cubical sets, persistence diagrams are efficient to compute.The reference [27] provides an overview of current TDA software.In this work, we use Perseus [28] and DIPHA [29] for cubical persistent homology computations.
Researchers have developed several methods for creating one-parameter filtrations.The most fundamental and commonly-used filtration for grayscale digital images is the sublevel set filtration (see, e.g.[30,31]).Using the sublevel set and thresholding operations (5), for thresholds Setting X i = g −1 ti (0) yields a sublevel set filtration {X i }.For binary images, there are techniques for constructing related grayscale images that would then lead to sublevel set filtrations.These include, for example, using a signed distance function or density estimator to define grayscale values [32,33].While [34] considers the changes of size functions (the 0-th persistence diagram of the sublevel set filtration) under the skeleton operation which combines certain morphological operations, our goal here is build a general filtration framework using erosion, dilation, opening, and closing.To the best of our knowledge, the proposed work is the first to use morphological operations and thresholding to construct a filtration directly.
We now use morphological operations to form new filtrations for binary and grayscale images.For ease of discussion, throughout the article we use a sequence of structuring elements, {B i } n i=0 , where each B i is a (i + 1) × (i + 1) square given by where e 1 = (1, 0) and e 2 = (0, 1).These may be depicted as where Note that since B 0 = {(0, 0)}, the erosion/dilation, and opening/closing operations with respect to B 0 are the identity map.Other sequences of structuring elements will also give rise to filtrations.In particular, there is a notion of shift inclusion that may be used to designate a large class of sequences of structuring elements that may be used to form filtrations.That topic is studied in detail in [25].
The first new filtrations we propose are for binary images and use the erosion and dilation operations.We consider a sequence of erosion and dilation operations with respect to {B i } n i=0 , i.e. for each i, consider δ Bi (f ) and Bi (f ) for a given binary image f .Similar to the sublevel set filtration in (10), the desired property is that if i ≤ j (B i ⊆ B j ), then Bi (f ) −1 (0) ⊆ Bj (f ) −1 (0).Thanks to Proposition 2 and Proposition 4, it is straightforward to verify that This shows that for any sequence of nested structural elements, erosion and dilation form filtrations.We call {X i } n i=0 , where X i = Bi (f ) −1 (0), the erosion filtration, and {X δ j } 0 j=−n , where , we may form one extended filtration by taking { Xi } n i=−n , where for i < 0, Xi = X δ i , and for i > 0, Xi = X i .The second new filtration we propose is related to the opening and closing operations.Since opening and closing are compositions of erosion and dilation operations, one may expect that Proposition 4 would extend to the case of opening or closing.However, it is not true in general.We refer readers to [25] for a counter example and more discussion.Essentially, the sequence of structuring elements cannot be arbitrary and requires additional assumptions.[25] presents a sufficient condition called shift inclusion that guarantees the structure necessary for opening and closing to result in appropriately nested sets.
Since our chosen square structuring elements, B i , satisfy shift inclusion [25], O Bi and, separately, C Bi , also form filtrations.
Similar to erosion and dilation filtration, we call , where , we may form one extended filtration by taking { Xi } n i=−n , where for i < 0, Xi = X C i , and for i > 0, Xi = X O i .As a byproduct, applications of ( 15) and ( 16) lead to three additional filtrations based on the commonly used top-hat transformation: the white top hat [35,16,14].More precisely, one has BTH STH For illustration, we now present a relatively simple opening filtration.Consider the modified Kanji image shown in Figure 1.The original binary image is shown in Figure 1(a) and denoted by f .Sample sets from the opening filtration on f , X O i , i = 0, . . ., 18, are shown in the top two rows of Figure 2 and the corresponding 1st level persistence diagram is shown in Figure 2(k).By (8), we know that β 1 (X O 0 ) = #P 1 (0).In particular, P 1 (0) consists of the points on the vertical axis of Figure 1(k).
By construction of the opening filtration, features in the original image 0-level set, X O 0 , (having birth coordinate b = 0), that are small in spatial scale relative to the structuring elements, have a short lifespan (small death coordinate d) whereas larger scale features have a longer lifespan (large d).The separation in scale between spatially small and large features is evident on the left vertical axis (b = 0) in the persistence diagram in Figure 1(f).

Features in X O
0 that disappear quickly under small amounts of opening (that is features with b = 0 and small d), have a small geometric/spatial scale, features in X 0 that persist under a lot of opening (b = 0 and large d) are more robust, and features that only appear after a lot of opening (large b) are most likely spurious.Choosing m = 4 in this example allows us to separate these three groups by drawing horizontal and vertical lines at d = 4 and b = 4 respectively.Correspondingly, X O 4 contains only so-called robust features and β 1 (X O 4 ) is a count of these features.One of the captured features, however, has a birth coordinate b > 0 indicating that it was not present in the original image.
We now return to the more general erosion/dilation and opening/closing one-parameter filtrations presented earlier and extend these to form filtrations on grayscale images.Combining either of these filtrations with the sublevel set filtration in (10), one may obtain a two-parameter filtration, or bi-filtration.We take the opening filtration as an illustration.Given a grayscale image g, by (10), we have g −1 ti (0) ⊆ g −1 tj (0) for any t i ≤ t j .Since for each t, g t is a binary image, by (15) we have that O Bi (g t ) −1 (0) ⊆ O Bj (g t ) −1 (0).By combining both (10) and ( 15) we obtain This is a 2-filtration, an example of a multi-filtration defined in Definition 4. In the next section, we formalize and extend the class of multi-filtrations constructed from opening and closing operations on binary images and opening, closing, and thresholding operations for grayscale images.

Multi-parameter Filtrations
At this point, we have seen that erosion, dilation, opening, and closing each form one-parameter filtrations for binary images and that combining one of these operations with thresholding forms a 2-filtration (see Definition 4 below).As we show in the next example, opening and closing operations may also be combined to form a 2-filtration for a binary image.In fact, this process may be continued to define k-parameter filtrations, the overall goal of this section.
Definition ( [6]).For k ∈ N and u, v ∈ Z k we say that u ≤ v if and only if u i ≤ v i for all i.Given this partial order on Z k , a family of sets {S i } i∈A with indexing set A ⊆ Z k is a multifiltration (or k-parameter filtration) if for any u, v ∈ A with u ≤ v, S u ⊆ S v .
Combining an opening filtration and a closing filtration and invoking Proposition 1 yields the following 2-filtration.
Figure 2 shows sample images from this opening/closing bifiltration as applied to the Kanji example with additive noise also shown in Figure 1(a).As seen in Figure 1, opening operations alone will not allow us to remove the small scale features due to additive noise while preserving the topology of the larger scale features.By visual inspection of the bifiltration depicted in Figure 2, X O,C (−2,3) appears to be the most accurate rendering of the underlying Kanji image.Since X O,C (0,0) ⊆ X O,C (−2,3) , there is no way to compare X O,C (−2,3) directly to the original image X O,C (0,0) using a one-filtration.However, multiple one-filtrations within the 2-filtration may be used to "connect" the two sets.
We discuss an approach for using persistent homology information to search for optimal renderings within multifiltrations, along with extensions of multifiltrations from Section 3 for binary images to a larger multifiltration that handles grayscale images, in Section 5.
We now present a general multiparameter persistence framework using morphological operations.Consider a sequence of operations E i : I P → I P (e.g.erosion) and a sequence of operations, D i : I P → I P (e.g.dilation) satisfying the following: for any f, g ∈ I P and i, j ∈ {0, 1, 2, . . ., n}, When E i = O Bi and D i = C Bi , assumptions (A1), (A2), and (A3) are similar to the sieving axioms in granulometry: anti-extensivity, increasingness, and the absorption property ( [14,16]).Assumption (A1) is the increasing property seen also in Proposition 1. Assumption (A2) is the absorption property ([15] Sec.1.2.6, p.20).Finally, combining assumptions (A2) and (A3) would lead to the anti-extensive or extensive property.
In what follows, let E i and D i , i ∈ {0, 1, 2, . . ., n} be sequences satisfying (A1), (A2), and (A3).Consider i ∈ {0, ±1, ..., ±n} and define a function M E,D i : and denote the 0-level set by When the context is understood, we sometimes abbreviate M E,D i as M i , and X E,D i (g) as X i .The notation ( 22) unifies the operators E and D in the following way.Lemma 2. Let i, j ∈ {0, ±1, ..., ±n} and g ∈ I P be a binary image.Suppose E i and D i satisfy (A1), (A2), and (A3).
If, on the other hand, i < 0, then there are two cases.In the case when j > 0, then M E,D j (g) = E j (g) ≤ g ≤ D |i| (g) = M E,D i (g).In the second case where j ≤ 0, then since i ≤ j ≤ 0, |j| ≤ |i| and we obtain M E,D The essential step in obtaining a multi-parameter filtration is to apply M E,D inductively.This requires us to extend the notation of ( 22) and ( 23) to a multi-index i ∈ Z k .Definition.Let {E i } and {D i } be sequences of morphological operations that satisfy (A1), (A2), and (A3).For k, n ∈ N and i = (i 1 , i 2 , . . ., i k ) ∈ {0, ±1, ..., ±n} k , we define M E,D i : and Similarly, we abbreviate the notation M E,D i as M i and X E,D i (g) as X i if operations E, D and image g are specified.
For example, for i = (−1, 1), M i (g) = M −1 (M 1 (g)) means that the image g is filtered by Motivated by (21), we consider the sets X i formed by the application of alternating E and D operations.Using (22), we see that alternating these operations corresponds to a multi-index consisting of an alternating sequence of integers.
We are now ready to present our main theorem: alternating the operations E and D leads to a multi-parameter filtration.
Proof.Let u = (u 1 , ..., u n ), v = (v 1 , ..., v n ) ∈ A and u ≤ v.By Definition 4, we need to verify that M u (g By Lemma 2, since u n ≤ v n we have that M vn (g) ≤ M un (g).Applying (A1), we see that Since u n−1 ≤ v n−1 by Lemma 2 again, we have Therefore, by combining ( 26) and ( 27), we prove that Finally, by applying the argument inductively one may conclude that By Proposition 2, we conclude that M u (g Remark.For purposes of exposition and to align with common practices in using morphological operations in image smoothing, we focused the composition of operations on alternating sequences (see [16,15,14]).This is inherent in (A2) as well as the stipulation that the indexing set A in Theorem 1 consists of alternating sequences.We note here, however, that the theorem holds true even if A contains sequences that are not alternating.
We now discuss examples to illustrate the framework given in Theorem 1.As a first example, consider erosion and dilation given as E i := Bi and D i := δ Bi .For this pair of operations, (A1) follows from Proposition 1, (A2) follows from Proposition 4, and (A3) is clear.
Therefore, by Theorem 1, X ,δ i i∈A forms a multi-parameter filtration, where A is any set of alternating sequences.
As a second example, consider the opening and closing operations and let forms a multi-parameter filtration.In fact, erosion/closing, and opening/dilation would also lead to multiparameter filtrations.It is important to note that while erosion/dilation and opening/closing lead naturally to multiparameter filtrations, the top-hat transformations do not.As mentioned in Section 3, these transformations do not satisfy (A1) and (A3) in general and, therefore, do not satisfy the hypotheses of Theorem 1.
At this point, g is assumed to be a binary image.If g is a grayscale image, one may combine the sublevel set filtration with the multiparameter filtration described in Theorem 1 to obtain another multiparameter filtration.In the rest of this section, we will formulate this concept as an extension of Theorem 1.
The following is an example of the framework in Theorem 1, While different methods, including the rank invariant function [7] and sheaf theory [36,37], have been developed to study multi-parameter persistence, for purposes of illustration we will focus on computing persistent homology along nondecreasing paths in the constructed multi-filtration.Definition.Define a nondecreasing path in indexing set A as a sequence u 0 , u 1 , . . ., u l ∈ A such that u i ≤ u i+1 for all i = 0, . . ., l.Then for a multifiltration {X u } u∈A and nondecreasing path u 0 , u 1 , . . ., u l in A, {X ui } i is a one-parameter filtration.
As we outline in the following section, this structure allows us to systematically extract information about geometric scale and optimize for certain topological features.Following multiple or successive nondecreasing paths allows for greater exploration of the multifiltration.

Application: A Denoising Algorithm for Salt and Pepper noise
We now use the multiparameter filtration to construct a denoising algorithm aimed at removing salt and pepper (small spatial scale, high amplitude) noise.See, e. g., Figure 1 where one goal is removing small scale white regions from the images in order to focus on the larger scale features.Given a binary image f , we wish to apply alternating opening/closing operations to f .Traditionally, this requires visual inspection to tune the size of the utilized structural elements as well as the number of operations performed.We now seek to automate this process by more fully utilizing the full multiparameter persistence framework.In this section, we will describe details of our proposed algorithm, demonstrate it on synthetic images, and extend it and apply it to grayscale and color images.
The proposed algorithm is iterative.In each iteration, we will use persistence diagrams computed along a nondecreasing path in the multiparameter filtration to guide the choice of a structuring element used for opening or closing.To get started, consider a binary image f contaminated by salt and pepper noise, that is, certain pixels have been switched to either white (salt) or black (pepper).See Figure 4 (b) and (e) for examples of images contaminated by the salt and pepper noise.Each iteration consists two steps: closing operation followed by opening.It is also possible to perform opening operation followed by closing for instead.
We first consider the closing filtration of f , {X C j } 0 j=−n , as shown in (16), where in this case is the largest set in the closing filtration and the persistence diagram can be decomposed into . Similar to the discussion for opening filtration in Section 3, P C 0 contains features (black regions) that are present in the original image X C 0 = f −1 (0), and |b| where b ∈ P C 0 indicates the size of the feature by giving the amount of closing required to remove it from the image.Since salt and pepper noise creates features that are small in spatial scale, we take a conservative route by choosing i c = min{|b| | (b, 0) ∈ P C 0 }.The binary image of the first step is then X ic := C Bi c (f ) −1 (0).To generalize this approach, we note that a gap in the death coordinate values in P C 0 can be used to detect a separation in spatial scales for features in the original image.Using the new binary image C Bi c (f ), we now consider the opening filtration of C Bi c (f ), {X O i } n i=0 , as shown in (15), where in this case ) reveals size information of the white regions.Specifically, P 1 ({X O i } i ) can be decomposed as As demonstrated in Section 3, P O 1 contains features that are present in the original binary image, X O 0 , and d ∈ P O 1 indicates the spatial size of the feature, that is, the amount of opening required to remove the feature.Similar to our approach in the first step using opening, we choose the size of the structuring element for closing to be i o = min{d | (0, d) ∈ P O 1 }.The binary image following the second step is now We repeat this alternating process.The stopping criterion is when the selected structuring element size exceeds a preset maximum, SizeTol.This could be given by the size of the image, or given as an upper bound on the spatial size of noisy features or features we wish to remove.The algorithm is summarized in Algorithm 1. Figure 3 illustrates Algorithm 1 in a schematic way in the multiparameter space.
1 Input: A binary image f ∈ I P and its set companion parameter SizeTol, and maximum number of iterations MaxIter. 2 Output: X O,C i =: X i , where i ∈ {0, ±1, . . ., ±n} l for some l ∈ N.
3 Denote X = f −1 (0). 4 Set g ← f 5 for j = 1, 2, ..., MaxIter do Return X and i To test the proposed algorithm, we again use the 190 × 190 binary image shown in Figure 4(a) as the ground truth.We add salt and pepper noise to the ground truth with various levels of noise densities.The noise density parameter gives the portion, in probability, of pixels whose values have been changed from their original values (to either white or black).We use the Matlab built-in function imnoise along with a specified noise density parameter.For instance, noisy images with noise density 0.4 and 0.7 can be found in Figure 4 We conduct an experiment to further test the proposed algorithm.Again, we take the ground truth image Figure 4(a) and add salt and pepper noise to it with noise densities from 0.1 to 1.0.For each noise density, we construct 1000 noisy images with the prescribed noise density, and for each noisy image, we apply Algorithm 1 with MaxIter=10 and Sizetol=5.The metric we use to compare the ground truth with the computed one is the intersection over union (IOU) defined as |S1∩S2| |S1∪S2| for sets S 1 and S 2 .In this case, we use S 1 = f −1 (0) and S 2 = f −1 (0), the black sets for the ground truth image and the output, denoised image, respectively.Under this metric, high IOU scores close to 1 measure good agreement/high overlap between the sets while numbers closer to 0 indicate that the sets are very different.The numerical results are shown in Table 1(a).When noise density is less than 0.3, the IOU scores are above 0.9, indicating high agreement between the image produced by Algorithm 1 and the ground truth image.For noise density 0.4 ∼ 0.6,  compare our method which is unsupervised with the commercial denoising software: Topaz AI [38] which uses deep learning.The results are shown in Figure 5(h), (i) and (j). Figure 5(k) shows the average extended IOU scores, where the extended score for each grayscale image is the average of the IOU scores over each threshold.For a color image, the extended IOU score is given for each of the three color channels.We observe that Algorithms 2 and 3 outperform Topaz AI on these examples.
Our focus in this section has been to demonstrate that the multiparameter filtration contains useful information and that automation may be used to extract it.Our proposed denoising algorithm works well on binary, grayscale, and color images with salt and pepper noise where the separation in spatial scale between the noise and true features may be used to effectively remove noise.Recently, others have developed salt and pepper denoising algorithms using deep learning e.g.[39,40], where a training process is required.Our unsupervised approach does not require training and it would be interesting to investigate whether combining the two approaches could lead to even better results.4(a) with added salt and pepper noise.For each prescribed noise density, 1000 images were formed, Algorithm 1 with MaxIter=10 and SizeTol=5 was applied to each, as was denoiseImage, a built-in Matlab function for denoising the image by using the deep neural network.Finally, calculated the IOU scores with respect to the original image (Figure 4

Conclusion
In this work, we establish that, under mild conditions, the morphological operations of erosion, dilation, opening, and closing may be combined to form multiparameter filtrations useful for studying binary images.These operations may also be combined with thresholding to form yet larger multiparameter filtrations useful for studying grayscale, and, by extension, color images.The dimension of the filtration grows with the number of operations and structuring elements, forming a potentially high dimensional framework in which to explore image structure and features.As demonstrated in Sections 5, this framework can be used to create automated approaches to image analysis and processing, in our example application leading to methods for removing salt and pepper noise from images.
There is a much broader class of methods for extracting information from multiparameter filtrations than just the approach of calculating persistence along nondecreasing paths that we describe in Definition 4 and use in Section 5. Persistent homology may be generalized as a cellular sheaf defined on a partially ordered set (P, ≤), that is, a functor from P to the category of vector spaces [41,42,43,44].Cellular sheaves were originally developed for studying nerve theory in topology [45] and have recently been used for describing the persistence of objects in applied topology.Because the order in Definition 4 is also a partial order on Z k , persistent homology defined on a multifiltration has a natural cellular sheaf structure.The persistence of the structure is much more complicated since the totally ordered property fails on the new order.However, we do see a variety of approaches for analyzing topological features in this setting, such as sheaf cohomology [42,46,47], zig-zag homology [48], multi-graded Betti numbers [49], and rank invariants [7].The multiparameter filtration we create here offers a constructive class of examples on which to explore these methods.
If either A = {a} or B = {b} are sets of singleton points, we would simply use a + B, a − B, A + b or A − b rather than {a} + B, {a} − B, A + {b} or A − {b} to denote the Minkowski sum or difference of A and Bs.

X and i end Algorithm 1 :
The proposed denoising algorithm .
(b)  and Figure4(d), respectively.The denoised images by Algorithm 1 can be found in Figure4(c) and Figure4(f), respectively; the denoised images by imnoise can be found in Figure4(d) and Figure4(g), respectively Visually, the denoised images by Algorithm 1 are close to the ground truth.Even in the case when noise density is 0.7, the denoised image (Figure4(f)) still recovers much of the core structure of the ground truth Figure4(a).On the other hand, denoised images by imnoise are still pixelated.
(a)) was calculated.IOU scores are listed in the columns below.All scores are recorded by mean ± standard deviation for the 1000 trials.(b) Betti numbers by mean ± standard deviation for all trials.The Betti pair (β 0 , β 1 ) of the original image (Figure4(a)) is(6,5).

Table 1 :
Performance of Algorithm 1 for denoising 2D binary images formed from Figure