Adaptive Mathematical Morphology on Irregularly Sampled Signals in Two Dimensions

: This paper proposes a way of better approximating continuous, two-dimensional morphology in the discrete domain, by allowing for irregularly sampled input and output signals. We generalize previous work to allow for a greater variety of structuring elements, both flat and non-flat. Experimentally we show improved results over regular, discrete morphology with respect to the approximation of continuous morphology. It is also worth noting that the number of output samples can often be reduced without sacrificing the quality of the approximation, since the morphological operators usually generate output signals with many plateaus, which, intuitively do not need a large number of samples to be correctly represented. Finally, the paper presents some results showing adaptive morphology on irregularly sampled signals.


Introduction
Commonly, morphology is applied to digital, discrete signals. Often, these signals are the result of sampling some continuous phenomenon (e.g., a real-world scene sampled by a digital camera). Notably, a band-limited continuous signal may be represented on a discrete sampling grid, as long as the sampling density is sufficiently high, as shown by Nyquist and Shannon [18,21]. The morphological operators are defined in the continuous as well as the discrete domain. When applying a morphological operator to a continuous signal, the result will, generally, no longer be band-limited. Thus, the result can no longer be represented on the sampling grid. Therefore, if we attempt to approximate the continuous operators in the discrete domain, we will not be able to represent the result within the same domain. Applying the discrete operators does not change the sampling density, and therefore, the correspondence between the continuous and the discrete domain is broken.
Furthermore, since the morphological operators depend on extrema in the signal, the discrete operators will poorly approximate the continuous ones if the sampling misses extreme values (which is likely to happen, even if the sampling rate is sufficient to represent the continuous, band-limited signal). However, since the continuous signal can, theoretically, be reconstructed, the extrema can be found. Using these extrema should allow for a better approximation of the continuous operators, if the extrema are found before applying a discrete morphological operator (it should suffice to find maxima when applying dilations, and minima when applying erosions).
Finally, morphology depends on a so called structuring element (SE), which itself depends on the sampling grid. The structuring element can be considered a geometric object, which is passed over the image when applying an operator. The shape of the structuring element affects the kind of filtering that will take place. One example would be applying openings with a disk-shaped SE. This could be used to filter out objects, whose sizes are smaller than the SE. In this way, one may, for example, compute a size distribution that gives information on the sizes of objects present in an image. However, as mentioned above, this SE depends on the sampling grid. Specifically, this means that a disk-shaped continuous SE is approximated by a discretized version in discrete morphology, which can lead to significant errors [14] (see Fig. 1).
To better approximate the continuous operators, in the discrete domain, we propose to allow for discrete, irregularly sampled signals. The purpose of allowing irregular sampling is threefold: 1. As described, there are issues with the continuous operators destroying the band-limitedness of signals, thereby breaking the correspondence between the continuous and discrete domain. Allowing for irregularly sampling the output would enable us to selectively increase the sampling density of the transformed signal so as to better preserve the correspondence between the continuous and discrete domain.
Since the output can be sampled irregularly, we hope to avoid an explosion in output samples, which would be necessary in the regular case, since the sampling may not have to be dense everywhere. In fact, the morphological operators generally create a large number of plateaus, which-intuitively-one could hope to represent using only a few samples. 2. The issue of missed extrema in the correctly sampled signal could also be alleviated, in a similar manner as the non-band-limitedness issue, since we can selectively increase the sampling density in the input around points of interest to hopefully capture the extreme values. Again, the fact that we allow irregularly sampled signals means that the input does not have to be densely sampled everywhere, unlike in the regular case. 3. Finally, since there is no sampling grid in the irregular case, the issue of discretizing the structuring element is not as problematic. Allowing for irregularly sampling the structuring element, means that we can add points exactly along its edge, which is where the discretization error normally would crop up.
A more thorough discussion on these issues can be found in a previous paper [3]. Although limited to 1D, most of the observations generalize readily to higher dimensions.

Previous Work
In a previous paper [3], the authors proposed a way of applying one-dimensional, irregular morphology. Experiments indicated that this could indeed lead to a better approximation of the corresponding continuous 1D-operators. Another work [2] extends the idea to two dimensions. In this case, too, experiments indicate that it is possible to gain accuracy when approximating continuous 2D-morphology with irregular, discrete morphology. This paper can be seen as an extension of that work along with a more recent paper that takes a similar approach (though avoiding some problems of the former approach) [4]. Some of the concepts are similar, but here the method has been generalized to allow for-more or less-arbitrary structuring elements (including non-flat and adaptive SEs), while also refining the way structuring elements are sampled.
A related approach dealing with questions of continuous morphology is to model the operators as evolutions of partial differential equations [5-7, 20, 25] and representing the result on the sampling grid. This approach is slow [28], and the fact that the sampling grid is still relied upon means that the issues discussed above still apply.
Another approach was taken by Thurley [23]. In this case, the idea is to apply morphology to irregularly sampled data, by moving samples up and down according to a continuous structuring element. In this approach, the output samples are taken at the same positions as the samples in the input. Therefore, the sampling issues discussed above still apply.
Finally, a paper by Calderon and Boubekeur [8] presents an approach for irregularly sampled data. However, they deal with binary morphology on 3D point clouds, not 2D, grayscale morphology. Moreover, they require explicitly estimating the underlying surface. Even so, in this paper we take inspiration from their approach to defining structuring elements using signed distance fields.

Background
In this section we recall some definions and establish the notation used in this paper.

Flat Morphology
Flat morphology, i.e., morphology using flat structuring elements B(x) where B is a set of points describing the shape of the SE, is called flat because the value of B is constant for all x ∈ B. Some examples of SE shapes are given below:

Square
A square structuring element: where S is half the length of a side of the square.

Extending SE shapes by decomposition
Simpler shapes may be used as SEs in repeated morphological operations to compute operations with a more complex SE.

Disk
A disk with radius r is commonly used as a structuring element. Note that the representation of the disk on a regular sampling grid is not isotropic (i.e. uniform in every direction)

Disk using L p -norm
Generalizing the concept of disks using p-norms: gives us a wide range of structuring elements.

Arbitrary shape
Of course, B can be any arbitrary subset of E 2 (where E 2 is the space in which the signals reside), however, for practical reasons, as we shall see, certain shapes are more straightforward to use together with the proposed algorithm (namely those that are easily described as a signed scalar field and for which a sampling strategy is easily developed), nevertheless, in the sequel we will also consider structuring elements of arbitrary shape.

Non-Flat Morphology
In the general case, the structuring element B(x) is some function that is not constantly 0 over B: where h : B →R is some function that gives a height to each point inside the structuring element. A particularly common choice is a parabolic structuring function: where || · || is the Euclidean distance. Parabolic structuring functions are of interest since they are isotropic and separable [24].

Morphological Operators on Irregularly Sampled Data
The main idea of the proposed algorithm is based on placing the SE at each sample and sampling the SE. Importantly, none of these samples can fall underneath a previously placed SE. In other words, if placing the SE at an input sample of a higher intensity than that of the sample at which the SE is currently located would mean that the umbra [22] of the SE contains the sample, then that sample must be suppressed. This is analogous to taking the supremum at the sample location.

The Structuring Element as a Signed Scalar Field
Similar to the idea proposed by Calderon and Boubekeur [8], we model the SE using a signed scalar field, f (see Fig. 2), such that E.g., the p-disk with radius r is modeled by the field The model is such that for a point inside the structuring element, we have that f (x, y) ≤ 0. Similarly, for a point outside the SE, f (x, y) > 0. To place the SE at a given input sample (x ′ , y ′ ), we translate all the input samples by (−x ′ , −y ′ ), meaning the input sample under consideration ends up at the origin.
Using signed scalar fields allows for easily combining simple shapes into more complex ones, e.g., by taking the max or min of two scalar fields one can get intersections and unions of the shapes defined by the fields. Some examples of such combinations are shown in Fig. 3. One may also generate structuring elements by sampling points of a complicated function and letting f be an interpolating function over these sampled points, however, for complicated shapes, the question of how to sample the resulting transformed signal becomes more difficult.
It should be noted that this idea is essentially the same as the concept of level sets [19], where the signed scalar field is equivalent to the level-set function. In particular, signed distance functions (i.e. the special case where the absolute value of f denotes the distance to the boundary of the represented object), which are commonly used in the level set approach [13] are also linked to morphology, since for example using Euclidean distance, the signed distance function gives dilations/erosions (depending on its sign) with ballshaped SEs of a size determined by the absolute value of the function [10] and v.v.

Adaptive (Shape) SEs
Allowing for adapting the shape of the structuring element at different parts of the signal is quite straightforward in this framework. One simply has to vary f based on some desired criterion. For example, one may vary the shape of the SE from a diamond, to a disk, to a square, by changing the p used by the p-norm. Figure 4 shows an example where p depends on the distance from the center of the image (p = 2 in the center, and p → ∞ as the distance increases). Another example is shown in Fig. 11, where the structuring element is an ellipse that is rotated and shaped depending on directional information at the point under consideration.

Sampling SE
For small structuring elements, the number of samples to represent the SE does not need to be particularly high, however, as the size of the SE increases, we need to increase the number of samples taken to keep the sampling density of the SE approximately constant with size. In a previous paper [2], dealing with a similar issue, we considered the gradient magnitude as a way of determining the necessary sampling density, since-intuitively-the output needs more samples near edges to avoid fuzzy borders, but plateau-like features, which do not need as many samples, will have a gradient magnitude close to zero. Here the proposed approach generates samples in a different way, and does not suffer from the same problems of fuzzy borders. Preliminary results do show that taking into account properties of the input (for example the morphological Laplacian [26]) can reduce the number of samples in the output (at least halving the number of samples) without adversely affecting the approximation of continuous morphology, however, for simplicity, we keep the sampling density constant at each point of the signal, meaning the SE is sampled proportional to its size, but not based on any property of the input signal.

Flat Morphology
In order to implement dilations (which will enable erosions, openings, and closings), the following are necessary: -A signed scalar field f (x, y, xc , yc , v, r) that describes a structuring element centered at position (xc , yc) (for example the p-disk described by Eqn. 8 translated) that depends on some value v (usually derived from the input signal), and whose size is determined by r. -A way of sampling the zero-set f = 0, which we denote ∂f (xc , yc , v, r).
Let us consider, for the sake of simplicity, a function f (x, y, xc , yc , r) that is independent of v, and its zero-set ∂f (xc , yc , r) (we will later come back to the case where f is not independent of v).
To perform a dilation of size r on an input signal of N samples, where x i and y i is the sample position, and z i is its value, do the following (see Fig. 5 for a side-view illustration of the steps of the algorithm): 1. Let I 0 δ = ∅. In further iterations we will fill out sets I i δ with samples from the dilation, as we process samples in I. 2. In this step we select an unprocessed sample (xc , yc , zc) ∈ I and mark it as processed. In Fig. 5(a), the sample marked with a blue "x" is currently being processed. 3. Sample ∂f (xc , yc , r). Let us call this set of samples B (shown in blue in Fig. 5(b)). 4. Sample ∂f (xc , yc , r + ε), where ε > 0 is some small real number. We call this set of samples Bε (shown in red in Fig. 5(c)). 5. For each sample (x ′ , y ′ ) ∈ B, and each sample ( (see the blue star and blue "x" in Fig. 5(d)). LetB denote the remaining set of samples (blue circle), i.e., discard samples from the smaller SE that are overlapped. 6. Repeat the procedure in 5. except replacing B with Bε. LetBε denote the remaining samples (red circle in Fig. 5(d)). The red star in Fig. 5(d) shows a discarded sample. Figure 5(e) shows the remaining samples contained inB (blue) andBε (red).
Finally, let Then, for each sample ( . Conceptually, z h is the value of a given sample of the bigger SE (Bε) after being projected downward until hitting a SE at a lower value (or falling through to −∞). This is illustrated in Fig. 5(f) by the red circle. 8. LetB Meaning we add the samples spawned during processing of (xc , yc , zc) to I i+1 δ , where i is the number of previously processed samples. This process is then iterated until all samples in I have been processed. Then I N δ is the dilated signal.
Conceptually, we center a structuring element on top of each input sample (x, y, z) and sample its border (steps 2-3), additionally, we also sample a slightly larger SE (step 4). We then discard all these samples that are subsumed by SEs placed at input samples at a higher value (steps 5-6). The samples from the slightly larger SE are projected down until they hit a SE placed at some sample below the sample value of zc (step 7). These projected samples help define the border between adjacent SEs. The dilated signal is the union of all these sampled structuring elements (step 8). Note that f may be an adaptive structuring element, taking different shapes at different positions. Also note that when implementing the procedure outlined in the steps above, there are a number of optimizations that may be made (we do not, for example, need to actually consider all samples in I \ {xc , yc , zc} in step 5).
Erosions are easily computed using the same procedure by first inverting the signal (by multiplying the values by −1), dilating, and then inverting that result. Openings and closings may then be performed by composing erosions and dilations.

Adaptive Morphology
Coming back to the case where f is not independent of v, we shall see that this case can be used to implement adaptive morphology. In this case we reintroduce the v parameter to all cases where f is applied in steps 1-8. In short, the SE being sampled changes based on the value of v at the given position (xc , yc) and when projecting samples downward in step 7, we need to check for intersections with the adapted structuring elements at the samples (x i , y i ). Later we will show an example of applying this approach where v is the local structure tensor [11,12,15]. In Figure 4, the parameter v is the distance from the center of the image, and the signed scalar field describes a p-disk, where p is 2 when the distance is 0 and increases proportionally with the distance.

Non-Flat Morphology
We can also extend this approach to non-flat morphology. First, let f − denote the set of points where f is negative (i.e., inside the SE), also, let h(x, y, xc , yc) be a function denoting the height of the structuring function centered on (xc , yc) at (x, y). Then we may replace step 3 in the approach above by: , yc , r). Let us call this set of samples B.
In other words, we also need to sample inside the SE since it is no longer flat. In step 5 (and step 6), we also need to take into account the height of the structuring function: 5 ′ . For each sample (x ′ , y ′ ) ∈ B, and each sample (x i , y i , z i ) ∈ I \ {xc , yc , zc}, compute If f i ≤ 0, and h(x ′ , y ′ , xc , yc) + zc ≤ h(x ′ , y ′ , x i , y i ) + z i (16) discard (x ′ , y ′ ). LetB denote the remaining set of samples.
In step 7, we again need to consider the height of the structuring function: where Finally, let Then, for each sample ( . Similarly to the case for a flat SE, we are projecting the samples from the larger SE down until they hit a SE at a lower sample (or fall through). However, in this case we need to consider the height of the non-flat SE by offsetting the values by the height of the relevant SE, meaning we start dropping the samples from a value offset by the height of the SE, and consider the height of SEs when finding intersections.
Finally, in step 8, we need to change the setB c to the following: Meaning the samples inside the SE have to be offset by its height.
In summary: the overall idea is the same as in the flat case, but we now need to offset our samples according to the height of the structuring function. Moreover, when projecting the samples from the slightly larger SE downward, we also need to consider the offset when finding intersection points with other structuring functions. One additional consideration is that the sampling of f − in step 3 ′ is now over a surface instead of a boundary, so there are other considerations to be taken into account when deciding how to generate the samples of B than in the flat case. Figure 6 shows an example of applying a dilation with a parabolic SE to the cameraman image shown in Fig. 4 using the proposed approach. Note however that this is just an example, and that h is not restricted to being a parabolic function.

Sampling Strategy
In these experiments we have chosen a simple sampling strategy where, in the flat case, each SE is sampled along its border by taking ⌊4πR⌋ samples along the border of the SE (where R is the size of the SE) such that the angle between any pair of adjacent samples is equal and the sum of angles is the entire circle. For simplicity we here stick to the naive approach of sampling the SE with equal density everywhere, however preliminary experiments where the sampling density is varied based on some property of the signal shows that the number of samples in the output can be decreased without hurting the accuracy of the proposed approach to any great degree, while yielding faster results. For example, by computing the morphological Laplacian and adjusting the number of samples at each point proportionally to this value, in some experiments the number of samples in the output is more than halved, without significantly affecting performance.
In most of our experiments we consider a p-disk, where p = 2, in which case the number of samples is twice the circumference of a circle of radius R. This value was chosen empirically.
For the non-flat case, where the interior of the SE also needs to be sampled, we choose to sample the SE as if we were sampling several similar SEs of size 1, 2, . . . , R along their edge as described above, where R is the actual size of the non-flat SE.

Inserting Extrema
As mentioned earlier, one advantage of allowing for irregularly sampled input signals is that one may insert extrema that do not fall on the regular sampling grid without increasing the sampling density everywhere. Figure 7 shows a toy example where the continuous signal cos( √︀ x 2 + y 2 ) 9 + 1 2 is dilated by a square structuring element with sides of length 6. The function is sampled from −2π to 2π with samples 0.5 apart. This sampled image is dilated with the corresponding structuring element. The sampling misses the maximum at (0, 0), so the result of the dilation approximates the continuous dilation poorly. Inserting a single sample at (0, 0) improves the approximation significantly (see Fig. 7). This proves the idea, although it is not obvious in general how to find these extrema. In practice, one could consider an instrument sampling additional points of a real world object as deemed necessary. Alternatively, it is possible that cleverly interpolating could find extrema without also adding a large number of non-essential points.

Comparison of Regular, Discrete Dilation and Proposed Algorithm for Regularly Sampled Images
To compare the performance of the proposed dilation to the performance of regular, discrete dilation, we perform experiments on regularly sampled images. The images are first filtered by a Gaussian kernel (σ = 1.0) in order to remove high frequencies, since we will be subsampling these images. The filtered image is dilated using regular dilation with disks of different radii. These dilations will be used as a proxy for the continuous dilation of the image. The filtered image is subsampled by keeping only every third pixel in every third row. This subsampled image is then dilated using both regular, discrete dilation and the proposed irregular dilation. The dilated images are then interpolated onto the regular grid of the original image, using linear interpolation. Finally, we compute the mean square error (MSE) and the structural similarity [27] (SSIM) between the interpolated images and the dilated original filtered image, which acts as a stand-in for the continuous dilation of the image. Figure 8 shows the images used (unfiltered) and the MSE and SSIM for different radii. Figure 9 shows a closeup of part of the mandrill image shown in Fig. 8 dilated by a disk using the proposed approach (interpolated back onto the regular grid) as well as regular, discrete dilation. The figure also shows the regular, discrete dilation of the image before subsampling, acting as a reference. The artifacts created by discretization are clearly visible in the example showing the regular, discrete dilation. We also perform a similar experiment with non-flat morphology (see Fig. 10) where the SE has a disk shaped support of radius 16 and a parabolic height defined by h(x) = −k||x|| 2 , where k scales the height of the SE. The experiments show the MSE and SSIM for a couple of images using varying values of k. When k = 0, the SE is flat. As k increases both variants yield better results. The discretization errors of the disk are less noticeable for the regular, discrete dilation for larger k, since the non-flat SE helps smooth out the disk near the edge of the support.

Adaptive Elliptical Structuring Element
In a paper by Landström and Thurley [17], adaptive morphology using different elliptical structuring elements based on the local structure tensor (LST) [11,12,15] is introduced for regularly sampled signals. In this section we adapt this work to this irregular morphology framework. The main idea is to compute the LST at each sample and apply a morphological operation using an elliptical structuring element that is rotated and shaped based on information extracted from the LST locally. This can be used to preserve strong edges when filtering, since the ellipses near edges are thin and oriented along the edge. The mean square error comparing the dilation of the subsampled image (for the regular discrete dilation and the proposed approach) with the dilation of the original blurred image. The structural similarity index for the same experiment. The regular, discrete dilation has a number of peaks/valleys, which are a result of the SE being restricted by the sampling grid. In general, the proposed dilation with and without preprocessing perform better than the regular, discrete dilation. The jagged nature of the curves for the regular, discrete case is a product of the necessary discretization of the SE.
The eigenvectors e 1 (x) and e 2 (x) associated with λ 1 (x) and λ 2 (x) point in the direction of largest and smallest variation respectively [9].  Fig. 10). The regular discrete dilation (left) is extremely jagged, since the small disk (r = 2 in the sub-sampled case) is poorly approximated on the discrete grid. The proposed approach (middle) yields a much better approximation of continuous morphology. The rightmost image shows the dilation of the blurred image that has not been subsampled (which is a proxy for the continuous dilation of the image). The idea of Landström and Thurley [17] is to use the eigenvectors e 1 and e 2 at each pixel position x of an input image to shape an elliptical structuring element by letting the length of the major axis depend on λ 2 and the minor axis on λ 1 and rotating the ellipse so that the major axis coincides with e 2 . The ellipses are scaled such that M is the maximum size of the semi-major axis for any SE.

Computing the Local Structure Tensor
In the case of regular, discrete morphology estimating the image gradient at each pixel position is relatively straight-forward (as is applying a Gaussian filter). However, in the case of irregularly sampled signals, where sample positions are not constrained by any sampling grid, estimating the gradient is a bit trickier. One might consider normalized differential convolution [16], but in this paper we make use of a kind of morphological gradients [1]. This fits well into the current framework, since all that is required (as we shall see) is dilations and erosions with certain structuring elements on the irregularly sampled signal, all of which is already possible with the approach presented in Section 3. Using this morphological approach also helps illustrate the versatility of the proposed approach to morphology on irregularly sampled signals.
Estimating the Gradient. In a previous paper [1], we propose a morphological approach to estimating the gradient for a regularly sampled image with missing samples and indicate that essentially the same approach can be applied to estimate the gradient for irregularly sampled images as long as we can apply morphology in the irregular domain. The necessary groundwork has already been laid previously in this paper, allowing us to use SEs of appropriate shape. Specifically, we need to define elliptical SEs and half-planes as signed scalar fields, to correspond to the ellipses and half-planes used in the regular case described in the previously mentioned paper [1].
Thus, we aim to construct four scalar fields g + h , g − h , g + v , and g − v , which describe four half-ellipses. First we define six simpler scalar fields, which we then combine two-and-two by taking the point-wise maximum (see Sec. 3.1). The first two scalar fields, g h , gv are ellipses with their major axis oriented horizontally or vertically: where a, b, and r determine the shape and size of the ellipses. The four remaining scalar fields split the plane into two: hr(x, y) = −x + εe (26) hu(x, y) = −y + εe (27) h d (x, y) = y + εe (28) Here εe > 0 gives the offsets of the split from the origin in a cardinal direction.
Combining the scalar fields for the ellipses with the scalar fields for the split plane one may create scalar fields that represent halves of ellipses (in fact the ellipses are cut off slightly-i.e. εe units-past the halfway point). The following four scalar fields will be used to form structuring elements: Applying dilations using these four functions will allow us to estimate the gradient at each sample of an irregularly sampled signal. The idea is to dilate the image by g + h and g − h and then take the difference. This yields an estimate of the horizontal component of the gradient. Similarly, to estimate the vertical component we take the difference between the two dilations using g + v and g − v as structuring elements. To improve the estimate, one may use additional half-ellipses to estimate directional gradients. The information from all estimates may be combined by setting up a system of equations and solving for the vertical and horizontal component of the gradient. Additionally, combining the result of dilations and erosions and using non-flat structuring elements can also improve the estimate, as detailed in the paper [1]. Note that, to compute the difference between the two dilated signals, their values need to be interpolated at each position of the input signal, in this paper we use linear interpolation.
The scalar fields (23)-(28) depend on a few parameters, namely a, b, r and εe. These can be used to alter the shape of the halved ellipses. Parameters a and b alter the general shape of the ellipse (how thin or thick the ellipse is), r alters the size of the ellipse and ε moves the position of the cut that splits the ellipse in two. One could adapt these values based on the underlying signal so that sparser parts of the signal use larger, more circular ellipses, while the SEs are smaller and more narrow as the density increases. However, in our experiments we will fix a = 2, b = 1, r = 1 and εe = 1. See the paper [1] for a more thorough examination of suitable parameter values.
Gaussian Filter. To apply a Gaussian filter to an irregularly sampled signal, we take the weighted sum of the neighbors of a sample, normalizing the weights so that their sum is 1. Now we have everything we need to estimate the local structure tensor (22) at each input sample. In practice an initial smoothing of the input signal with a Gaussian filter using a small σ is applied before estimating the gradient.
In summary: to estimate the gradient, we first apply a small amount of blurring with a Gaussian with σ = 0.3, we then dilate the set of irregularly sampled points using first g + h and g − h and take the difference (after interpolating back to the original sample points using linear interpolation), followed by dilating the blurred input points using g + v and g − v and taking the difference. These differences estimate the horizontal and vertical gradients. For the local structure tensor, the σ of the Gaussian is set to s/ √︀ 2 log 2, as per the paper by Landström and Thurley [17]. Here s is a parameter that controls how the influence of edges will be propagated.
Although the simplest method of using only the horizontal and vertical estimates outlined above has some rotational bias, thie estimation can be markedly improved by performing essentially the same process, except with differently rotated pairs of ellipses and combining the results. The details of this are left to paper [1]. This yields a more rotationally invariant approximation. Figure 11 shows an example of applying a dilation using adaptive elliptical structuring elements to a regularly sampled image with missing samples. Figure 12 shows an irregularly sampled pile of rocks, and the result of applying an erosion with adaptive, elliptical SEs. In this case, the sampling density is higher along the x-axis than along the y-axis, so in this case the half-ellipses g + v and g − v use a = 4 instead of a = 2 to allow for samples farther away to be captured by the morphological gradient. The other parameters are kept as before.

Conclusions and Future Work
In this paper we have extended previous work on morphology on irregularly sampled signals to allow for more general structuring elements, including non-flat morphology and adaptive morphology. We argue that allowing for irregularly sampled signals when dealing with morphology in the discrete domain can lead to a better approximation of continuous morphology for several reasons and show this by comparing the proposed approach and regular, discrete morphology on synthetic as well as natural images. Additionally, we make use of a previously proposed morphological approach [1] to estimating gradients in order to compute the local structure tensor for irregularly sampled signals and use this to implement elliptical, adaptive structuring elements in the irregular case based on the work by Landström and Thurley [17]. Some preliminary results regarding sampling strategies for the proposed approach are reported, indicating that considering the gradient or the Laplacian of the input signal can be useful in order to improve the approximation of continuous morphology while keeping the sampling density of the transformed signal low. Further investigation regarding how to best adaptively change the sampling density of the output based on properties of the input signal (particularly relating to its smoothness, in some sense) and the SE are of interest. For example, the sampling of the SE should probably take into account its shape by increasing the density near sharp locations on its border and decreasing the sampling density where the border is smooth. Additionally, for non-flat SEs, the sampling should also take into account how quickly the height of the SE varies, as well as extrema of its height. Another point of interest to consider, regarding the sampling of the transformed signal, is its suitability for composing morphological operators.
As mentioned in Section 1.1, related work for dealing with considerations of continuous morphology exists in the form of morphology by partial differential equations [5-7, 20, 25]. However, this approach makes use of a sampling grid, which means that they will suffer from the problems outlined in this paper with regard to representing the transformed band-limited signals on a sampling grid. Namely, for some band-limited input function, the transformation of the continuous function may introduce higher frequency content into the signal, which cannot be represented on the sampling grid. In contrast, the approach proposed in this paper does not make use of any sampling grid, and thus is not restricted by this limitation. Nevertheless, another avenue for future work is to experimentally compare the proposed approach with PDE based approaches.
In our experiments we have made use of linear interpolation where necessary, this is motivated by the fact that the morphological operators generate signals where the derivative is discontinuous. A possible improvement would be to consider different kinds of interpolation, perhaps adapting the approach based on properties of the output signal. For example, at smooth parts of the signal, one might use spline-based interpolation, while linear interpolation may be the best choice near sharply varying parts of the output signal. These ideas are also related to the insertion of extrema described in Section 5.2.
The proposed approach may not, strictly speaking, yield morphological operators. One difficulty when dealing with irregularly sampled signals is that there is no straight-forward way of finding the supremum/infimum of two signals (usually defined by taking the pointwise supremum/infimum, which is not possible in the irregular case, since two irregularly sampled signals may not have pairs of corresponding points). One approach is to consider a kind of interpolation of the irregularly sampled signals and instead deal with the corresponding continuous signals. This approach is taken in a previous paper [3] where we deal with irregularly sampled signals in 1D, and show that the proposed approach yields morphological operators under certain conditions.