# Scaling Up SLIC Superpixels Using a Tile-Based Approach

Dawa Derksen, Jordi Inglada, Julien Michel

<figure><img src="image_1.png" /></figure>

## To cite this version:

Dawa Derksen, Jordi Inglada, Julien Michel. Scaling Up SLIC Superpixels Using a Tile-Based Approach. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57 (5), pp.3073-3085. 10.1109/TGRS.2018.2880248. [hal-04886083](https://hal-INRAeS.sciencedirect.com/10.1109/TGRS.2018.2880248)

HAL Id: hal-04886083

https://ut3-toulouseinp.hal.science/hal-04886083v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Scaling Up SLIC Superpixels Using a Tile-Based Approach

Dawa Derksen<sup>1,∗</sup>, Jordi Inglada, and Julien Michel

**Abstract**— Image segmentation techniques are challenging to apply to large-size remote sensing imagery. Indeed, if the data to be processed are larger than the computer's available memory, it must be split into smaller pieces. Without precaution, segmentation errors appear along the edges of these pieces. The goal of this paper is to present a tilewise processing method to overcome this issue for superpixel segmentation, applied in particular to the simple linear iterative clustering algorithm. Incidentally, tilewise methods allow for several pieces of the image to be processed simultaneously, which enables the deployment of these methods in a parallel processing environment. Estimations of the speed-up when using multiple processors are provided. Then, it is demonstrated that the result of the tilewise segmentation is equivalent to the segmentation of the complete image, with respect to a number of global unsupervised segmentation criteria. Finally, experimental results on a large-size Sentinel-2 time series validate the method's feasibility.

**Index Terms**— Image processing, image segmentation, scalability, superpixel segmentation.

## I. INTRODUCTION

RECENT satellite-based Earth observation systems such as Pléiades, WorldView, GeoEye, and Sentinel-1 and 2 generate massive amounts of data, which can sometimes be challenging to process. The growth in the data size is mainly due to a gain in spatial, spectral, and temporal resolutions, as satellite systems become more and more sophisticated. For example, the Sentinel-2 optical Earth observation constellation provides images with a ground sampling of 10 m in four channels, 20 m in six channels, and 60 m in three other channels. Combined with a temporal revisit time of 5 days, a time series of one scene (300 km²) over a year can reach up to 1 TB of memory. Complex memory management mechanisms must, therefore, be implemented to enable computation on such large images. In many cases, the image can be split into smaller regions, usually rectangular tiles, so that each region fits in the computer's memory. The final result is a stitching

of the individually processed regions. This mechanism, known as piecewise or tilewise processing, works well for most classic pixelwise or fixed neighborhood algorithms. This paper addresses the issue of tilewise processing applied to superpixel segmentation methods. A superpixel segmentation method is a kind of image segmentation, in which the objective is to isolate small, homogeneous regions of similar size, equally spread across the image [1]. This is different from other image segmentation methods, in which the objective is to extract real objects from an image. A more complete definition of superpixel segmentation will be given in Section III.

Recent studies suggest that superpixel segmentation can be used for a finer analysis of the image than a pixel-based or an object-based analysis, by providing small compact homogeneous regions which can be used at various steps in a processing algorithm. This has proven to be useful in several applications, for instance, to aid in hyperspectral image classification [2], [3], or even in the field of convolutional neural networks [4].

This paper arises from the need to apply superpixel segmentation to very large remote sensing images, which could otherwise not be processed in this way due to limitations in memory. In this paper, the main objective is to demonstrate a tilewise method to apply superpixel segmentation to images of any size, while ensuring the same segmentation quality as if the image had been processed in one piece. Although the segmentation quality will be measured, the aim of this paper is not to improve the segmentation itself. Actually, segmentation quality will only be used as a quantitative measure of similarity between the image segmented with and without tilewise processing.

The main contribution of this paper is a novel tilewise processing framework to apply a superpixel algorithm without memory constraints. In addition, some modifications to the original simple linear iterative clustering (SLIC) algorithm, in order to enable tilewise processing, are also proposed.

Furthermore, these following experimental results are shown.

1) SLIC applied on a small Pléiades scene, proves the method validity.

2) SLIC applied on a large Sentinel-2 time series, to prove the feasibility of large data sets.

The rest of this paper is structured in the following manner. Section II details the previous studies on large-scale segmentation. Next, Section III explains the superpixel paradigm, as well as the SLIC algorithm and its modifications.

Manuscript received September 27, 2017; revised March 5, 2018 and June 11, 2018; accepted August 13, 2018. Date of publication January 9, 2019; date of current version April 22, 2019. The work of D. Derksen was supported in part by the Centre National d'Etudes Spatiales and in part by ATOS under Grant 2714. (Corresponding author: Dawa Derksen.)

D. Derksen is with CESBIO, Centre National d'Etudes Spatiales, 31401 Toulouse, France (e-mail: joederksen1@gmail.com).

J. Inglada is with DCT/SI/CB, Centre National d'Etudes Spatiales, 31401 Toulouse, France.

J. Michel is with DCT/SI/AP, Centre National d'Etudes Spatiales, 31401 Toulouse, France.

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TGRS.2018.2880248
The tilewise processing mechanism applied to this algorithm will be detailed in Section IV. Finally, the experimental results on Pléiades and Sentinel-2 images will be shown in Section V.

## II. PREVIOUS WORK

The scaling up of segmentation algorithms is a challenging task, and one unifying framework suitable for all algorithms has not yet been found. Two main groups of methods to deal with this issue can be distinguished. The first group contains methods based on adding a postprocessing step to fix the errors linked to tiling. Such methods operate on the tile edge regions after the segmentation step is finished. For example, in [5], segment fusion along the tile edges is based on user-defined topological criteria such as the length of the contact surface between the segments to be merged. However, this does not completely eliminate errors on the edges as the edge can still appear in certain parts of the segments. Another idea from the same group of methods is to reprocess the region edge areas after the algorithm is finished using a more sophisticated method than simple region fusing. In [6], the segments on the tile edges are merged by the full lambda schedule algorithm, first, on the vertical borders, and then, on the horizontal ones. However, if this method is applied to a superpixel segmentation, fusing regions will create abnormally large superpixels across the tile edges. A similar idea is developed in [7], in which the segment border areas are processed separately at the end of each iteration of the algorithm. Unfortunately, reading and writing the image into disk memory at each iteration is very costly.

The other group of methods for dealing with the tile edges issue involves acting on the tiling itself before applying the segmentation algorithm. In [8], the image is split into nonrec- tangular tiles. The authors propose a tiling following gradients in the image. However, there is no guarantee that the actual segments would have followed this tiling, if the method had been applied to the whole image. Furthermore, image gradients are difficult to use in the case very high dimensional data, such as a Sentinel-2 time series. Another idea is to pad each tile by a margin on all sides, instead of using nonoverlapping tiles. By doing so, the hope is that segmentation results in the central area of the tile are less affected (even unaffected) by the tile border. A certain rule can then be defined to choose how to stitch the margin areas together. Michel et al. [9] propose a method that guarantees precisely the same results with and without tilewise processing, applied to the mean shift segmentation algorithm. For this, they use the concept of stability. A segmentation algorithm is considered stable if, when applying it to any subregion of an image, it produces the same result as if the image had been first segmented, then cropped to the subregion. Marginal stability has a looser constraint: the results must be identical within a central zone of the subregion. If an algorithm is shown to be marginally stable, and this margin can be calculated, the tiling with margins method (described earlier) is guaranteed to provide exactly the same result as a segmentation without tiling. In a further study [10], these concepts are applied to the region merging family of algorithms. The tiling process developed in

this paper is similar to the methodology developed in [10]. Section IV explains the details of how this can be applied to SLIC and to superpixel methods in general.

## III. SUPERPIXEL SEGMENTATION

Segmentation refers to the general problem of dividing an image into nonoverlapping regions in a way that these regions exhibit certain properties, for example, spectral homogeneity or compacity. Often, it is used as a way of extracting objects from an image. Some view the object extraction via segmentation problem as ill-posed, as the definition of an object in an image can vary from user to user. For instance, in a remote sensing image, one user might want to single out trees, while another might want to isolate forests. This means the segmentation problem does not have a unique solution for each image. Superpixel segmentation is a less ambitious version of the segmentation problem, in which the purpose is to divide the image into small, homogeneous, compact, and similarly sized segments. Superpixels can be seen as conglomerates of several pixels, but that adhere to the natural object borders in the image. As opposed to usual segmentation methods, superpixels almost always oversegment the objects in an image. Figs. 1–3 show an example of a superpixel segmentation of a remote sensing image.

There are different ways of obtaining a superpixel segmentation of an image. With the right parameters, some well-known segmentation methods such as watershed [11], mean shift, or region merging can generate superpixels. However, these provide fewer guarantees and less control on the superpixel parameters than custom-made superpixel methods.

There are two main families of superpixel methods: graph based and gradient-ascent based. Graph-based methods consider each pixel of the image as a node in a graph. Neighboring nodes are connected by a vertex, the weight of which depends on the similarity between the pixels. By minimizing a global objective function over the image, superpixels can be created. Gradient-ascent-based methods, on the other hand, start with a poor initial segmentation, usually a regular square lattice, and iteratively relabel the pixels to optimize an objective function.

An example of gradient-ascent-based methods are the lattice methods, in which one considers the segmentation as a lattice (horizontal and vertical lines) which is then iteratively adapted to the image gradients [12], [13].

The most popular superpixel method, introduced by Achanta et al. [1], is called SLIC and belongs to the group of gradient-ascent-based methods. It provides better segment adherence to boundaries when compared to other state-of-the-art superpixel methods and was, therefore, chosen for this paper. The following paragraph describes the SLIC algorithm as in [1]. This algorithm is referred to in the rest of this paper as the original SLIC algorithm.

### A. SLIC

The following notations will be used throughout this paper. Let $I$ be an image with $N$ pixels of dimensions $X \times Y$, in which each pixel $p_n$ at position $(x_n, y_n)$ contains $D$ features called $f_n^d, n = [1..N], d = [1..D]$. For example, in a remote
sensing image, the spectral bands as well as any number of derived indices normalized differential vegetation index can be used as features.

SLIC is an algorithm similar to the K-means clustering algorithm. Each segment is represented by a centroid $C_k = [X_k, Y_k, F_k^d]$, $k = [1..K]$, $d = [1..D]$. The couple $(X_k, Y_k)$ represents the spatial mean (geometric center of gravity), while $F_k^d$ are the mean features on the segment $k$.

At the initialization step, the image is split into K square segments following a regular grid. The initial grid size (length of the side of the square grid elements) is a parameter of the algorithm, called the spatial width (SpW) parameter. This very important parameter defines the initial size of the superpixels and, therefore, the total number of superpixels. The algorithm then proceeds as follows. During the iterative step, each pixel is associated with its nearest centroid, following a similar principle as in K-means. The distance in this step is calculated as a weighted sum of the spatial and feature Euclidean distances, shown in (1), which is similar to the definition in [14]. The distance weight parameter, DW allows the user to give more importance to either segment compacity or segment homogeneity

$$ \begin{aligned} \text{dist}(p_n, C_k) &= \text{dist}_{\text{feature}} + \text{DW} \times \text{dist}_{\text{spatial}} \\ &= \sqrt{\sum_{d=1}^{D} (F_k^d - f_n^d)^2} \\ &\quad + \text{DW} \sqrt{(X_k - x_n)^2 + (Y_k - y_n)^2}. \end{aligned} \quad (1) $$

In the original algorithm, the candidate centroids are chosen in a limited search radius, to avoid calculating unnecessary distances [1]. Once all the pixels are updated, the centroids are recalculated, and the next iterative step is started. Iterations stop when a convergence criterion is met: either a maximal number of iterations or a total residual (average square distance between the centroids of two consecutive steps) smaller than a given value.

The most significant parameters of the algorithm are defined as follows.

1) SpW is the length of the side of each segment in the initial square grid. This parameter controls the total number of superpixels and their average size.

2) The distance weight, DW, is the weight of the spatial distance with respect to the feature distance. The combined distance is given in (1).

3) The search radius, SR, gives the search radius for candidate centroids.

SLIC, like K-means applied to an image, does not necessarily generate connected segments. This is due to the composite distance criteria: a pixel can be similar in terms of features to a centroid but spatially distant, meaning that other pixels can isolate a pixel from its centroid. Achanta et al. [1] propose a postprocessing step which involves removing the isolated regions by fusing them with a neighboring segment, based on the same composite distance. A region is considered isolated if the centroid of its segment does not belong to the region itself. Simply put, this postprocessing step will keep only the

regions that contain the centroids. Such a postprocessing step has two downsides as follows.

1) If the segment centroid happens to be outside the segment, for instance, in the case of a crescent-shaped segment, it will be incorrectly fused with a neighbor.

2) This postprocessing step is costly, as it requires the isolated regions to be found and distances with their neighbors to be calculated.

## B. SLIC With Simple Points

In exploring the application of SLIC to videos, Chang et al. [15] expose a way to overcome the issue of the postprocessing step. The solution is based on the concept of simple points, introduced in [16] and more recently in [17]. A point is characterized as simple if changing its label cannot alter the topology of its four-neighboring segments. For instance, a pixel in the middle of a segment is not simple, because if its label is changed, a hole will appear in the segment. On the other hand, a pixel on a straight edge between two segments will be simple, because it can be safely changed without altering the segments' topology. Calculating if a point is simple can be done in constant time, as it depends only on the labels in the eight neighborhood, as shown in [17]. In SLIC, the initial segments are fully connected, so if only the simple points are changed, this connectivity will be insured until the end of the algorithm. The other concept introduced by [15] is the log-likelihood of the segmentation. This evaluates the quality of a Gaussian fit on the data, both in the feature and the spatial domains. The log-likelihood is used to find the best neighbor for a label change during the algorithm.

In the algorithm proposed by [15], the iterative step is modified in the following way. For each pixel,

1) check if simple;

2) if simple, find which four-neighbor fusion brings the highest change in log-likelihood;

3) assign label of that neighbor; and

4) in case of change, update centroids of the two segments.

In this paper, following the logic of the original SLIC algorithm from [1], we have chosen to associate the pixel to the nearest centroid in terms of composite distance in the four neighborhood, instead of calculating the rather costly log-likelihood which contains both variance and mean terms.

Therefore, in the modified algorithm proposed by this paper, the iterative step becomes, for each pixel,

1) check if simple;

2) if simple, find which of the four neighbor's centroid is the nearest in terms of composite distance;

3) assign label of that neighbor; and

4) in case of change, update centroids of the two segments.

Compared to the original SLIC algorithm, the modified algorithm has the following advantages.

1) Fewer distances need to be calculated.

a) Distances are only calculated if the pixel is a simple point.

b) The original SLIC algorithm already restricts the nearest centroid search space by only calculating
<figure><img src="image_3.png" /><figcaption>Fig. 1. SLIC applied on a small 4000 × 4000 Pléiades image, with a superpixel width parameter of 100 and a distance weight of 1, with 20 iterations. Segments are equally spread across the image.</figcaption></figure>

distances to centroids in a limited search radius. In the modified algorithm, the search space is reduced even further, to the four neighboring segment centroids.

2) The search radius for candidate centroids parameter, SR, is no longer needed.

3) No postprocessing step is required.

In other words, the algorithm presented in this paper is an updated version of the original SLIC from [1], with a few of the several changes made in [15]. The impact of these changes is not quantified in this paper, as they are all taken from [15], which already validated these changes in terms of segmentation performance. Furthermore, the objective, here, is not to provide a better segmentation than other versions of SLIC but rather to demonstrate the tile-based upscaling of a superpixel segmentation.

Fig. 1 shows the modified version of SLIC applied on a Pléiades image. The superpixels are equally spread across the image and are all quite compact. Superpixel segmentations are known to oversegment image objects. However, Fig. 2 shows a zoom on an urban area, to show that superpixels can also undersegment some parts of the image, when the image objects are smaller than the desired superpixel size. Fig. 3 shows a close-up on an agricultural area. The different fields are relatively well separated, but hedges lying on the edge between fields will not be isolated in the segmentation.

For comparison, Fig. 4 shows the result of a region merging algorithm with Baatz and Schäpe criterion [10] applied to the same image. The segments are larger in homogeneous places (fields) and smaller in textured places (urban area) on the image.

## IV. TILEWISE SEGMENTATION

The main difficulty is finding an efficient way to scale up the segmentation algorithm. Indeed, if the image does not fit in

<figure><img src="image_4.png" /><figcaption>Fig. 2. Zoomed-in view of Fig. 1 over an urban area, the urban elements (houses and streets) are undersegmented and the surrounding fields are oversegmented.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Zoom-in view of Fig. 1, fields are well distinguished but thin hedge-like structures are not isolated.</figcaption></figure>

memory, it must be split into smaller pieces which can be dealt with individually. This mechanism is called tilewise processing. If no precautions are taken, the tile edges themselves will appear as segment edges in the final segmentation result. To counteract this effect, it is possible to initialize SLIC's segmentation on the entire image. The initialization is a square regular grid, and the initial labels can be calculated anywhere, even if the whole image is not loaded in the memory. Then, the algorithm can run as usual. In this case, the tile edges will not directly appear in the final segmentation, but other kinds of anomalies will, as shown in Fig. 9(c).

1) Sharp edges will appear along the tile edge zones.

2) A segment that is on two adjacent tiles will end up being twice as big as the average segment size.

3) Disjoint segments can be found along the edges, as the simple point constraint is no longer present between the two independent tiles.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Region merging with Baatz and Shäpe criterion.</figcaption></figure>

The solution that is developed here is inspired by the work in [9], in which, Michel et al. take margins around each processing tile. Their objective was to create a tilewise processing mechanism that guaranteed precisely the same segmentation result as if the image had been processed with no tiling. This was shown to work for algorithms that were stable. An image processing algorithm is considered stable if when applying it to any region of the image, the same result is obtained as if you had applied the algorithm on the whole image and, then, extracted the result in that region. SLIC is not stable following the notion of stability developed in that study, because the result depends on the order in which the pixels are considered. If applying it to only a subregion of the image, the pixels will not be scanned in the same order as on the entire image. Even the original version of SLIC is not stable, because the postprocessing step also depends on which order the isolated regions are read. In fact, the scanning order is totally arbitrary, [15] scans pixels in a random order, but for processing speed reasons, a lexicographical order was chosen for our application. This scanning order on the full image does not provide better results than scanning each region independently, because both orders are arbitrary. That is why in this paper, the identical results requirement from [9] and [10] was loosened, and replaced by a similarity requirement, which involves two aspects.

1) No segmentation errors should appear along tile edges.

2) The segmentations with and without tilewise processing should show similar overall segmentation characteristics, as defined by the four criteria described in the next paragraph.

A. Segmentation Quality Criteria

In this part, a method to verify if the two segmentations (with and without tilewise processing) have the same overall characteristics is developed. In [1], superpixel quality is given

by two measures, adherence to boundaries and boundary recall. These measures, such as the Hoover metrics [18], are based on dense validation data, in which all the pixels are labeled. The two segmentations could be compared to a reference segmentation, to ensure they have the same quality. Unfortunately, dense validation data are not available in our case. In [10], the Hoover metrics were used with the segmentation on the full image as reference data, because the objective was to achieve an identical segmentation. However, as mentioned previously, SLIC is not stable. Therefore, the differences between the segmentations with and without tilewise processing are not due only to tiling but are also linked to the arbitrary pixel scanning order. The Hoover metrics will be sensitive to these differences. For this reason, unsupervised evaluation metrics have been chosen. Superpixels aim at being compact, homogeneous and similarly sized, so these three criteria will be measured. Compacity shall be calculated two ways. First, based on the perimeter and area, [19], given in (2). This measure is always between 0 and 1, 1 being the value for a circle (maximal compacity)

$$C_p = 4\pi \frac{a^2}{p} \qquad (2)$$

The perimeter compacity measure from (2) is sensitive to segments with rough edges, so another complementary measure is proposed, based on the statistical spread of the position vectors (x, y) of the pixels of a segment. First, the 2 × 2 symmetric covariance matrix of the positions x and y of the pixels in the segment is calculated. Next, this matrix is diagonalized (covariance matrices are always semidefinite positive). These eigenvalues are calculated by the following formulas:

$$S_1 = \sigma_x + \sigma_y + \sqrt{\sigma_x\sigma_y - \sigma_{xy}^2} \qquad (3a)$$

$$S_2 = \sigma_x + \sigma_y - \sqrt{\sigma_x\sigma_y - \sigma_{xy}^2} \qquad (3b)$$

The normalized ratio of the first eigenvalue to the second eigenvalue, (4), gives a measure of the statistical elongation of the segment. This metric is between 0 and 1, with 1 being the maximal compacity (minimal elongation). It is rotation invariant because the eigenvalues are rotation invariant, and scale invariant because it measures a normalized ratio

$$C_s = \frac{2S_1S_2}{S_1^2 + S_2^2} \qquad (4)$$

Spectral homogeneity can be measured in several ways [19]. The selected metric, taken from [12], is the explained variance, as shown in (5). It measures the variance retained in the image after replacing each pixel value by the mean value in the segment. By normalizing this value by the total variance in the image, a value between 0 and 1 is obtained, that is, insensitive to the feature dimension and image content. A value of 1 indicates maximal homogeneity in the segments

$$EV = \frac{1}{D} \sum_{d} \left( \frac{\sigma_{d, \text{seg}}}{\sigma_{d, \text{im}}} \right) \qquad (5)$$

Finally, the similarity of segment sizes is measured by the relative standard deviation (RSD) of the areas of the segments, i.e., the standard deviation divided by the mean segment
<figure><img src="image_7.png" /><figcaption>Fig. 5. Chessboard margin strategy, white tiles take margins on the four surrounding black tiles and on the two white tiles of the above row.</figcaption></figure>

size, (6). The higher this value, the more the segment sizes are spread out

$$RSD_{size} = \frac{\sigma_{size}}{\mu_{size}}. \qquad (6)$$

These measures will be used to compare the characteristics of a segmentation with and without tilewise processing, the results of which are given in Section V.

## B. Tilewise Processing Procedure

The key of this tilewise processing procedure lies in taking margins around certain tiles. A margin is defined as a padding of a given width, applied on the tile in all directions. The idea is that the overlapping parts of two neighboring tiles can be combined to eliminate errors in these areas.

The process starts by splitting the image into $Nx \times Ny$ tiles, along with the x- and y-directions. Margins do not need to be taken around all the tiles. Indeed, a tile edge is shared by two tiles, so only one of the two tiles needs to take a margin on the other, to fix the errors around their common edge. Margins are, therefore, taken only on one in two tiles. Visually, the tiling can be compared to a chessboard with alternating black and white squares. The tiles around which margins are taken will be referred to as white tiles, the others as black tiles. This is illustrated in Fig. 5. White tiles take margins on the four surrounding black tiles and the corners of two white tiles above them.

This whole process is described step-by-step on a 2 × 2 tiling example in Fig. 6.

1) The first step consists in running the segmentation algorithm on all the black tiles. These tiles are treated independently, so they can even be processed in parallel.

2) Once all black tiles areas are segmented, the white tile processing can start. White tiles of the first row take margins on the surrounding black tiles. They load into memory not only the image in this area but also the segmentation result provided by the surrounding tiles.

3) The key is then to freeze all the segments on the margin outer edges, meaning they will not be altered by the segmentation algorithm. Then, the remaining parts of the white tiles are processed.

<figure><img src="image_8.png" /><figcaption>Fig. 6. Step-by-step tilewise segmentation process on a 2 × 2 example. A detailed description is available in Section IV-B. 1: all the black tiles are segmented first. 2: white tiles take margins on surrounding black tiles. 3: freeze the segments on margin edges. 4: no anomalies on the tile edges. 5: white tiles need to take a margin on the white tiles of the row above. 6: segmentation is complete.</figcaption></figure>

4) Because the outer edge segments have been frozen, the segments in that area will be coherent with the ones in the black tile and will not exhibit any edge anomalies.

5) The white tile of the next row is processed. If all the white tiles are processed independently like the black tiles, errors will appear around the corners of the tiles. For this reason, each white tile takes margins not only on the surrounding black tiles but also as on the white tiles of the row above. This means the white tiles can only be processed row by row, and the white tiles in one row can be processed in parallel.

6) The processing is finished, and the frozen segments of the final white tile fit well with the segments of all surrounding tiles.

The condition for this mechanism to function correctly is that the margin must be wider than the maximal width of a segment. If a segment crosses the whole width of the margin, it will be frozen with the tile edge still intact, and a straight line will appear along the edge in the final result. For general segmentation methods, it is difficult to calculate a maximal segment width, as segments can have any size. Fortunately, superpixels are known to be compact, and their average width is an entry parameter of the algorithm. Therefore, it is sufficient to take a margin of several superpixel widths. Experimental results show that for three superpixel widths, no more edge anomalies appear in the final segmentation result.

## C. Parallel Processing

Tiling the image as independent pieces opens up the possibility of processing several pieces at the same time, to speed-up computation. One property of the tiling algorithm as described in Section IV is that all the black tiles can be processed in parallel. Furthermore, all the white tiles in one row can be processed in parallel, however, the white rows must be processed in order. Therefore, the speed-up
will not necessarily be a linear function of the number of processors, as the algorithm has sequential steps. It is possible to calculate the theoretical speed-up as a function of the tiling parameters. For example, say $N_p$ processors have to perform $N$ tasks. We define $N_{eq}$ as the number of times $N_p$ tasks have to be performed in order to finish all $N$ tasks, in other words, the equivalent number of tasks to perform when several processors are involved. $N_{eq}$ can be seen as the total computation time expressed in a unit of time equal to the time needed to accomplish one task. Equation (7) gives the formula for calculating $N_{eq}$ as a function of $N$ and $N_p$. In the following equations, $E(.)$ is the integer floor function:

$$N_{eq} = E \left( \frac{N-1}{N_p} \right) + 1. \quad (7)$$

The achievable speed-up $s$ is $s = (N/N_{eq})$. For example, if five processors wish to perform 12 tasks, they will first perform five in parallel, then the next five in parallel, then two processors will perform the last two tasks, while the others are idle. This could be the case for the black tiles in a 6×4 layout, or the white tiles of one row in a 24 × 24 layout. In this case, $N_{eq} = 3$, and $s = (12/3) = 4$, meaning that using five processors only brings a speed up of 4. Estimating these speed-ups for the particular case of tiling with black and white squares can be done by calculating the total number of tasks for the different steps of the algorithm. In the first step of the algorithm, all the black tiles can be processed in parallel. The total number of black tiles is

$$N_{bl} = E \left( \frac{N_x N_y + 1}{2} \right)$$

meaning $N_{eq}$ for the first phase can be calculated using (7).

To calculate the achievable speed-up for the white tile processing step, it is necessary to first derive the number of white tiles in odd and even rows, the details of this calculation can be found in Appendix A.

The dotted lines in Fig. 7 represent the theoretical speed-up from (9) for a given number of processors, which is an upper bound for the actual speed-up. Memory reduction is the ratio of the size of the image to the size of a tile. Best performances are achieved in the upper-left region of the graph, where the speed-up is high for a low memory reduction. As the number of processors increases, the curves are nearer to this region. Although the tiling mechanism does not allow for all tiles to be processed in parallel, Fig. 7 shows that the maximal achievable speed-up does approach the number of processors when the memory reduction is high, i.e., when the number of tiles becomes large. This graph shows that having more processors does not always mean more speed-up. More processors become beneficial when the need in memory reduction is strong, that is, to say when dealing very large images, or when very little memory is available. The local maxima of the curves are points where the parallelization degree is the highest. For instance, the first local maximum of the five processors curve is placed at a memory reduction of 100, which is a 10×10 tiling, meaning each row has exactly five white tiles, so the processors are never idle. The saw-tooth behavior is due to the fact that some processors will be idle

<figure><img src="image_9.png" /><figcaption>Fig. 7. Memory reduction versus speed-up for different numbers of processors. Memory reduction is the ratio of the size of one tile to the size of the whole image. Solid lines represent the measured speed-ups, while the dotted lines represent the theoretical speed-up as given in (9) in Appendix A. The measured speed-up is lower than the predicted speed-up, but they exhibit the same trends.</figcaption></figure>

when the number of tasks is not a multiple of the number of processors.

The full lines in Fig. 7 show the measured memory reduction for different numbers of processors. The image used for these experiments is a 11000 × 11000 time series of Sentinel-2 images, sampled at 10 m, with 33 dates each containing 10 spectral bands. The total size of this image is 75 Gb, making it impossible to process without tiling. The experiments are run on a shared memory architecture, with Intel Xeon CPU E5-2650 v4 @ 2.20 GHz processors that communicate using the message passing interface protocol. For each tiling layout, the computation time with several processors is divided by the computation time with one processor to obtain the speed-up value. As predicted, the measured speed-ups are always lower than the theoretical speed-ups, but clearly, follow the same trends, and show the same local extrema, meaning the theoretical curves can be used to calculate the optimal memory reduction parameter for a given number of processors. The difference is stronger as the number of processors gets higher, because of the need for synchronization in a parallel computing environment.

## V. EXPERIMENTAL RESULTS

The experimental validation is done in two parts as follows.

1) Validating the tilewise method by comparing the results with and without tilewise processing, on a small 4000 × 4000 Pleiades image, with four spectral bands and a total size of 123 Mb.

2) Validating the scalability by testing the procedure on a large image: a 11000 × 11000 Sentinel-2 time series, with 330 features per pixel, and a total size of 75 Gb.

### A. Tiling Method Validation

For these experiments, a Pleiades image of 4000 × 4000 with four channels (R, G, B, and NIR) was chosen. First, a visual comparison between SLIC with and without margins exhibits the nature of the anomalies on tile intersection areas. Fig. 9 shows the results of the method on three intersection zones,
<figure><img src="image_10.png" /><figcaption>Fig. 8. Different tiling layouts for the tests on the Pleiades image. The tile edges cover a variety of inhomogenous terrains. Some of the edges of the layouts cannot be seen in this figure as they are covered by others. Two regions were selected for the illustrations in Fig. 9.</figcaption></figure>

covering relatively homogeneous natural areas, as well as more heterogeneous urban areas. The layouts and the zoomed regions can be seen in Fig. 8. In a second experiment, many different layouts are tested to show the robustness of the method to a variety of terrains at the tile edges. The different layouts for the second experiment are shown in Fig. 8.

The zoomed regions and the tile intersections are shown in Fig. 9(a). Fig. 9(b) shows the result of SLIC applied to the whole image, with no tilewise processing. The image is small and fits in memory, so this result can be generated for validation purposes. Fig. 9(c) contains the result of SLIC when applied tilewise with no margins. The three types of anomalies cited in Section IV are found in this example.

1) Straight edges along the border zone.

2) Abnormally large superpixels, especially on the intersection of all four tiles.

3) A disjoint segment (marked with a star).

Fig. 9(d) shows the segmentation result of SLIC applied with the tilewise strategy explained in Section IV. None of the previously mentioned errors are seen in this result, and it is very similar to the result without tilewise processing, as shown in Fig. 9(b). In Fig. 9(e), the difference between tilewise processing with no margins and applying the algorithm to the whole image, i.e., between Fig. 9(b) and (c) is shown. A pixel is colored white if the two labels are different. Fig. 9(f) shows the difference between Fig. 9(b) and (d). The segmentation result when using margins is clearly much closer to the original. However, as predicted in Section IV, there are still minor differences. Fig. 9 shows that the method is robust to heterogeneous terrains such as urban areas. In fact, there are generally fewer differences in such areas, because the result of a superpixel segmentation is more sensitive to the initial superpixel grid in homogeneous areas, as can be seen in Fig. 9(f).

To see if these differences have an impact on the overall segmentation statistics, the unsupervised segmentation criteria from Section IV-A have been measured on the three cases mentioned earlier, for different values of SpW, with a distance

weight of 1. In order for the tile edges to cover a variety of terrains, tiling layouts ranging from 2 × 2 to 5 × 5, as well as a rectangular layout, 1 × 8, were tested, as shown in Fig. 8. In these cases, the tiles size range from 2000 × 2000 to 800 × 800. In the 1 × 8 case, tiles are 4000 × 500. These results are shown in Fig. 10. Fig. 10(a), (c), (e), and (g) shows the results when processing without margins. The overall criteria clearly vary according to the tiling layout. Generally speaking, the differences are more important when the tiles are smaller (5 × 5 tiling layout), because the number of segments affected by the tile edges increases. Fig. 10(e) shows that the explained variance is slightly lower when processing without margins, which indicates less homogeneous regions. The differences in RSD of segment size are explained by the abnormally large segments on tile edges. In some cases, the values without margins are equal to the values without tilewise processing, for instance, 2 × 2 tiling layout with an SpW of 50. These points appear when the tile size is a multiple of SpW, meaning there are no segments that are initially on two tiles. The perimeter compacity increase is due to the straight borders on tile edges, as opposed to the rough edges which appear on most segments. Fig. 10(b), (d), (f), and (h) shows that when using the tilewise processing strategy with margins, the overall characteristics are indistinguishable from when the algorithm is applied on the whole image. While the segmentation result is not precisely the same, the margin strategy guarantees that these two cases will be similar in terms of their global characteristics.

## B. Scalability Validation

The definition of scalability has two aspects. First, the method should have linear complexity, meaning that the processing time should be a linear function of the size of input data. This is the case for SLIC, as each pixel is scanned a fixed number of times. The second aspect is that the method should function on an input of any size without degrading the quality of the result. A method is said to be scalable if it can process an image larger than the available memory, in linear time. This is done by tiling the image, i.e., splitting it into smaller regions and processing the regions independently, as explained in Section IV.

1) *Calculating the Optimal Splitting Parameters:* The objective is to provide a fully automatic scalable method that works on any image without user interference. For this, we must calculate the number of necessary regions for the processing to fit in the memory. For this calculation, we will assume the following parameters are known (input parameters).

1) The image *I* is of size *X* × *Y*.

2) The margin value is *m* = 3 × SpW in pixels.

3) The available memory, *T* bits.

4) Input pixel size, *p*<sub>in</sub> bits.

5) Output pixel size (label coded on an integer) *p*<sub>out</sub> bits.

The tiling is entirely defined by two parameters, $N_x$ and $N_y$, which represent the number of divisions along the *x* and *y*-axes of the image. For example, in the illustration in Fig. 5, $N_x$ = 5 and $N_y$ = 5. The problem can be formulated as finding the values of $N_x$ and $N_y$ that satisfy the memory constraint.
<figure><img src="image_11.png" /><figcaption>Fig. 9. (a) Zoom on the intersection regions between tiles. The intersection areas cover a variety of both agricultural and urban covers, to ensure the robustness of the method to heterogeneous terrains. (b) Result of SLIC applied to the entire image with no tilewise processing. (c) SLIC applied on each tile independently with no margins. Linear tile edges and abnormally large superpixels can be seen, as well as a disjointed segment marked with a star. (d) SLIC applied with tilewise procedure, with margins of three superpixel widths. No more anomalies can be seen on the edge areas. (e) Difference between SLIC on entire image and tilewise SLIC with no margins (a white pixel means the segments are different). There are many differences between the two segmentations, especially around the tile edges. (f) Difference between SLIC on entire image and SLIC with tilewise processing and margins of three superpixel widths. There are still minor differences between the two segmentations. These are concentrated in smooth terrains, where small differences in the initial superpixel grid can cause large differences in the segment boundaries, as there are no strong gradients for them to adhere to. The differences are almost all eliminated in rougher terrains, such as built-up urban areas.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 10. Comparison of segmentation statistics of SLIC applied to a Pleiades image, for different tiling layouts (2 × 2 to 5 × 5 and 1 × 8) and different superpixel sizes (SpWs). If no margin is taken (left column), anomalies around tile edges change the statistics for different tiling layouts. When a margin is taken around the tiles (right column), the statistics become identical to when the image is processed without tiling, and no longer depend on the tiling layout. (a) Perimeter compacity, no margins. (b) Perimeter compacity, margin = 3. (c) Variance compacity, no margins. (d) Variance compacity, margin = 3. (e) Explained variance, no margins. (f) Explained variance, margin = 3. (g) RSD of segment area, no margins. (h) RSD of segment area, margin = 3.</figcaption></figure>
<table><caption>TABLE I
MEMORY PROFILING RESULTS ON SENTINEL-2 TILE</caption><thead><tr><th>SpW (px)</th><th>Optimal tiling</th><th>Tested tiling</th><th>Peak mem. (GB)</th></tr></thead><tbody><tr><td>20</td><td>10x10</td><td>10x10</td><td>1.893</td></tr><tr><td>20</td><td>10x10</td><td>9x9</td><td>2.277</td></tr><tr><td>2</td><td>11x11</td><td>11x11</td><td>2.000</td></tr><tr><td>2</td><td>11x11</td><td>10x10</td><td>2.350</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Demonstration of SLIC applied to a 11 000 × 11 000 Sentinel-2 time series of 33 dates zoomed on a small coastal area of 1300 × 1000. The background is the visible bands of the first date of the series.</figcaption></figure>

The finer the tiling, the higher the total margin surface, which means extra memory overhead and in the end, longer processing times. Therefore, the objective is finding the largest tile size that fits in memory. The total margin $m_{tot}$ can be seen as a cost function, and the problem as a constrained integer optimization problem.

1) Variables: $N_x$, $N_y$.

2) Minimize

$$m_{tot} = m(XN_y + YN_x)$$

3) Constraint

$$T > \alpha \frac{XY}{N_x N_y} + \beta \left( \frac{X}{N_x} + \frac{Y}{N_y} \right)$$

where $\alpha$ and $\beta$ depend on the input parameters

$$\alpha = 2p_{in}$$

$$\beta = (2p_{in} + p_{out})2m.$$

The difficulty is finding an integral solution to this problem. Dedicated libraries could be used to find an exact solution, however, in this particular case, as the constraint function is strictly decreasing with respect to the input variables, a valid approximate solution is the ceiling of the real solution. The real solution is simply the solution of a second-degree polynomial. This provides us with a solution that is always valid, and while not necessarily optimal, it is a good approximation of the optimal solution

$$N_x^* = E \left( \frac{X\alpha}{\sqrt{\alpha T + \beta^2} - \beta} \right) + 1$$

$$N_y^* = E \left( \frac{Y\alpha}{\sqrt{\alpha T + \beta^2} - \beta} \right) + 1.$$

2) Validating the Optimal Splitting Parameters: To validate this method, it has been applied to a large Sentinel-2 time series. This image is the stacking of 33 dates of Sentinel-2 acquisitions, in 10 spectral channels, on a 110 × 110 km region. All channels have been sampled at 10 m. The total size of this image is around 75 GB, which clearly exceeds the 16 GB of available memory on the computer used for testing. Table I shows the peak memory measured on a run of the SLIC algorithm. These measurements were made using the Valgrind massif tool [20]. The optimal tiling calculation was set with a maximal memory of 2 GB for validation purposes. The table shows a comparison of the case where the optimal tiling is used, and the case where bigger tiles are used. It shows that calculating the tiling parameters allows the algorithm to best respect the memory constraints.

Fig. 11 shows a small region (1300 × 1000) of the result the application of the tiled algorithm to an entire Sentinel-2 time series. The linear and round elements that can be seen in the water are linked to the temporal interpolation of cloudy images on other dates in the year. On land, the superpixels adhere to the natural boundaries and manage to separate the different agricultural plots from urban cover, and from forests.

## VI. CONCLUSION

In this paper, we formalized a method for the tilewise processing of a superpixel algorithm. The method does not guarantee a result exactly the same as if the image had been processed in one piece. However, experimental tests using a modified version of SLIC on a Pléiades image show that the results are of the same quality in terms of four unsupervised segmentation criteria: perimeter compacity, statistical compacity, explained variance, and RSD of segment sizes. Finally, the scalability was validated by applying this method with SLIC on a real-size Sentinel-2 time series. This method allows us to process large images that were previously impossible to process on machines with limited memory. As the image is processed tilewise, some tiles can be processed in parallel. We have shown that although the speed-up is not necessarily, a linear function of the number of processors, if the tiling is chosen wisely, a maximal speed-up can be obtained. Further studies could develop a method to increase the parallelization degree of the tiling method, by designing a more complex way of organizing the tile processing. The tilewise processing method should also be validated by other superpixel segmentation methods. All of the results were generated using the Orfeo Toolbox 6.0.0 open-source software [21]. The source codes can be found at [22].

## APPENDIX

### A. Speed-Up of White Tiles

The total number of white tiles is

$$N_{wh} = E \left( \frac{N_x N_y}{2} \right).$$
The number of white tiles in an even row is

$$N_{\text{wh,even}} = E \left( \frac{N_{\text{wh}}}{N_{\text{y}}} \right).$$

The number of white tiles in an odd row is

$$N_{\text{wh,odd}} = E \left( \frac{N_{\text{wh}} + N_{\text{y}} - 1}{N_{\text{y}}} \right).$$

The total number of even rows is

$$N_{\text{even}} = E \left( \frac{N_{\text{y}} + 1}{2} \right).$$

The total number of odd rows is

$$N_{\text{odd}} = E \left( \frac{N_{\text{y}}}{2} \right).$$

Reminder: $N_p$ is the number of processors. The $N_{\text{wh,even}}$ tiles of any even row can be processed in parallel. There is a total of $N_{\text{even}}$ even rows, which must be processed sequentially. The equivalent number of tasks for white tiles on even rows is

$$N_{\text{eq,even}} = \left( E \left( \frac{N_{\text{wh,even}} - 1}{N_p} \right) + 1 \right) N_{\text{even}}.$$

The same logic applies to odd rows, the equivalent number of tasks for white tiles on odd rows is

$$N_{\text{eq,odd}} = \left( E \left( \frac{N_{\text{wh,odd}} - 1}{N_p} \right) + 1 \right) N_{\text{odd}}.$$

The total equivalent number of tasks for all white tiles is

$$N_{\text{eq,wh}} = N_{\text{eq,even}} + N_{\text{eq,odd}}.$$

The total equivalent number of tasks for all tiles is

$$N_{\text{eq}} = N_{\text{eq,bl}} + N_{\text{eq,wh}}.$$

Combining all of the previous equations, the equivalent number of tasks for the entire process is

$$\begin{aligned} N_{\text{eq}} = & E \left( \frac{E \left( \frac{N_x N_y + 1}{2} \right) - 1}{N_p} \right) + 1 \\ & + E \left( \frac{E \left( \frac{E \left( \frac{N_x N_y}{N_y} \right)}{N_p} \right) - 1}{N_p} \right) + 1 \right) E \left( \frac{N_y + 1}{2} \right) \\ & + E \left( \frac{E \left( \frac{E \left( \frac{N_x N_y}{N_y} \right) + N_y - 1}{N_p} \right) - 1}{N_p} \right) + 1 \right) E \left( \frac{N_y}{2} \right). \end{aligned} \quad (8)$$

The speed-up shown in the dotted lines of Fig. 7 is based on (8)

$$s = \frac{N_x N_y}{N_{\text{eq}}}. \quad (9)$$

## REFERENCES

[1] R. Achanta, A. Shaji, K. Smith, A. Lucchi, P. Fua, and S. Süsstrunk, "SLIC superpixels compared to state-of-the-art superpixel methods," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 34, no. 11, pp. 2274-2282, Nov. 2012.

[2] T. Lu, S. Li, L. Fang, L. Bruzzone, and J. A. Benediktsson, "Set-to-set distance-based spectral-spatial classification of hyperspectral images," *IEEE Trans. Geosci. Remote Sens.*, vol. 54, no. 12, pp. 7122-7134, Dec. 2016.

[3] X. Jin and Y. Gu, "Superpixel-based intrinsic image decomposition of hyperspectral images," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 8, pp. 4285-4295, Aug. 2017.

[4] M. Volpi and D. Tuia, "Dense semantic labeling of subdecimeter resolution images with convolutional neural networks," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 2, pp. 881-893, Feb. 2016.

[5] J. Michel et al., "Open tools and methods for large scale segmentation of very high resolution satellite images," in *Proc. OGRS*, 2012, pp. 179-184.

[6] Z. Yang, "Tiling and merging framework for segmenting large images," U.S. Patent 8086037, Dec. 27, 2011.

[7] P. N. Happ, R. S. Ferreira, C. Bentes, G. A. O. P. Costa, and R. Q. Feitosa, "Multiresolution segmentation: A parallel approach for high resolution image segmentation in multicore architectures," *Int. Arch. Photogram., Remote Sens. Spatial Inf. Sci.*, vol. 38, no. 4, p. C7, 2010.

[8] T. S. Korting, E. F. Castejon, and L. M. G. Fonseca, "Divide and segment: An alternative for parallel segmentation," in *Proc. GeoInfo*, 2011, pp. 97-104.

[9] J. Michel, D. Youssefi, and M. Grizonnet, "Stable mean-shift algorithm and its application to the segmentation of arbitrarily large remote sensing images," *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 2, pp. 952-964, Feb. 2015.

[10] P. Lassalle, J. Inglada, J. Michel, M. Grizonnet, and J. Malik, "A scalable tile-based framework for region-merging segmentation," *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 10, pp. 5473-5485, Oct. 2015.

[11] L. Vincent and P. Soille, "Watersheds in digital spaces: An efficient algorithm based on immersion simulations," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 13, no. 6, pp. 583-598, Jun. 1991.

[12] A. P. Moore, S. J. D. Prince, J. Warrell, U. Mohammed, and G. Jones, "Superpixel lattices," in *Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR)*, Jun. 2008, pp. 1-8.

[13] D. Chai, Y. Huang, and Y. Bao, "IRSL: Iterative refining superpixel lattice," *IEEE Geosci. Remote Sens. Lett.*, vol. 14, no. 3, pp. 344-348, Mar. 2017.

[14] C.-Y. Hsu and J.-J. Ding, "Efficient image segmentation algorithm using SLIC superpixels and boundary-focused region merging," in *Proc. 9th Int. Conf. Inf., Commun. Signal Process. (ICICS)*, Dec. 2013, pp. 1-5.

[15] J. Chang, D. Wei, and J. W. Fisher, III, "A video representation using temporal superpixels," in *Proc. IEEE Conf. Comput. Vis. Pattern Recognit.*, Jun. 2013, pp. 2051-2058.

[16] G. Bertrand, "Simple points, topological numbers and geodesic neighborhoods in cubic grids," *Pattern Recognit. Lett.*, vol. 15, no. 10, pp. 1003-1011, Oct. 1994.

[17] X. Han, C. Xu, and J. L. Prince, "A topology preserving level set method for geometric deformable models," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 25, no. 6, pp. 755-768, Jun. 2003.

[18] A. Hoover et al., "An experimental comparison of range image segmentation algorithms," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 18, no. 7, pp. 673-689, Jul. 1996.

[19] H. Zhang, J. E. Fritts, and S. A. Goldman, "Image segmentation evaluation: A survey of unsupervised methods," *Comput. Vis. Image Understand.*, vol. 110, no. 2, pp. 260-280, 2008.

[20] N. Nethercote and J. Seward, "Valgrind: A framework for heavyweight dynamic binary instrumentation," *ACM SIGPLAN Notices*, vol. 42, no. 6, pp. 89-100, 2007.

[21] M. Grizonnet, J. Michel, V. Poughon, J. Inglada, M. Savinaud, and R. Cresson, "Orfeo ToolBox: Open source processing of remote sensing images," *Open Geospatial Data, Softw. Standards*, vol. 2, no. 1, p. 15, Jun. 2017, doi: 10.1186/s40965-017-0031-6.

[22] D. Derksen. (Jul. 2017). *Simple Linear Iterative Clustering (SLIC) Superpixel Method Applicable on Large Images*. [Online]. Available: https://doi.org/10.5281/zenodo.830171
<figure><img src="image_14.png" /></figure>

Dawa Derksen received the master's degree in aerospace engineering from the Institut Supérieur de l'aéronautique et de l'Espace (ISAE-Supaéro), Toulouse, France, in 2016. He is currently pursuing the Ph.D. degree with the Centre d'Etudes Spatiales de la Biosphère (CESBIO) Laboratory, Toulouse, under the supervision of the Centre National d'Etudes Spatiales (French Space Agency), Toulouse.

He is currently involved in the operation production of image processing algorithms applied to the large scale classification of Earth observation images for land cover mapping.

<figure><img src="image_15.png" /></figure>

Julien Michel received the master's degree in telecommunications engineering from the École Nationale Supérieure des Télécommunications de Bretagne, Brest, France, in 2006.

From 2006 to 2010, he was with Communications et Systèmes, Toulouse, France, where he was involved in the studies and developments of remote sensing image processing. He is currently with the Centre National d'Etudes Spatiales, Toulouse, where he is involved in the development of image processing algorithms for the exploitation of Earth

observation images.

Mr. Michel is a member of the Orfeo Toolbox Project Steering Committee and a Charter Member of the OSGEO Foundation.

<figure><img src="image_16.png" /></figure>

Jordi Inglada received the master's degree in telecommunications engineering from the Universitat Politècnica de Catalunya, Barcelona, Spain, and the École Nationale Supérieure des Télécommunications de Bretagne, Brest, France, in 1997, and the Ph.D. degree in signal processing and telecommunications from the Université de Rennes 1, Rennes, France, in 2000.

He is currently with the Centre National d'Etudes Spatiales (French Space Agency), Toulouse, France, where he is involved in the field of remote sensing image processing at the Centre d'Etudes Spatiales de la Biosphère (CESBIO) Laboratory. He is involved in the development of image processing algorithms for the operational exploitation of Earth observation images, mainly in the field of multitemporal image analysis for land use and cover change.