Multi-scale segmentation algorithm for pattern-based partitioning of large categorical rasters

Analyzing large Earth Observation (EO) data on the broad spatial scales frequently involves regionalization of patterns. To automate this process we present a segmentation algorithm designed speciﬁcally to delineate segments containing quasi-stationary patterns. The algorithm is designed to work with patterns of a categorical variable. This makes it possible to analyze very large spatial datasets (for example, a global land cover) in their entirety. An input categorical raster is ﬁrst tessellated into small square tiles to form a new, coarser, grid of tiles. A mosaic of categories within each tile forms a local pattern, and the segmentation algorithm partitions the grid of tiles while maintaining the cohesion of pattern in each segment. The algorithm is based on the principle of seeded region growing (SRG) but it also includes segment merging and other enhancements to segmentation quality. Our key contribution is an extension of the concept of segmentation to grids in which each cell has a non-negligible size and contains a complex data structure (histograms of pattern features). Speciﬁc modiﬁcation of a standard SRG algorithm include: working in a distance space with complex data objects, introducing six-connected “brick wall” topology of the grid to decrease artifacts associated with tessellation of geographical space, constructing the SRG priority queue of seeds on the basis of local homogeneity of patterns, and using a content-dependent value of segment-growing threshold. The detailed description of the algorithm is given followed by an assessment of its performance on test datasets representing three pertinent themes of land cover, topography, and a high-resolution image. Pattern-based segmentation algorithm will ﬁnd application in ecology, forestry, geomorphology, land management, and agriculture. The algorithm is implemented as a module of GeoPAT – an already existing, open source toolbox for performing pattern-based analysis of categorical rasters.


Introduction
The goals and the means of analyzing Earth Observation (EO) data depend on the spatial scale of the data. Most analyses are performed on the scale of a single image which, in geographical terms, depicts a city or a non-urban region of an equivalent spatial extent. At such scale, the goal of the analysis is either to identify individual objects (like, for example, buildings) if the resolution of an image is sufficiently high, or to generalize an image into a map showing land use/land cover classes (LULC). The means for achieving these goals are either pixel-based classification (Li et al., 2014) or the Object-Based Image Analysis (OBIA) (Lang, 2008;Blaschke, 2010), which combines pixel-based segmentation with segments classification.
However, neither individual objects nor LULC classes are useful when analyzing EO data on broad spatial scales (province, country, continent, the entire Earth surface). This point is illustrated in Fig.1(A) which shows a fragment of the National Land Cover Dataset (NLCD) ) (a sixteen-classes LULC map obtained by classification of Landsat-7 images (Jin et al., 2013)) covering a large portion of the U.S. state of Wisconsin. At this scale, the meaningful analysis of the data is to identify and delineate spatial units (hereafter also referred to as segments) containing unique quasistationary (hereafter referred to also as homogeneous) patterns of LULC classes. Arguably, several of such A B 0 50km developed high intensity (24) developed med. intensity (23) developed low intensity (22) developed open space (21) cultivated crops (82) woody wetlands (90) open water (11) deciduous forest (41) pasture/hay (81) evergreen forest (42) ice/snow (12) barren land (31) mixed forest (43) shrub/scrub (52) grassland ( units could be identified visually in Fig.1(A), but to delineate an exhaustive set of such units an algorithmic approach is necessary. The purpose of this paper is to describe an algorithm that starting with a categorical raster as an input (like the NLCD seen in Fig.1(A)), delineates segments characterized by homogeneous patterns of raster categories (like those seen in Fig.1(B)). Note that restricting an input to categorical rasters is a simplifying assumption, but many pertinent datasets are indeed categorical (for example, a land cover), and other pertinent datasets are categorizable (see below).
The proposed segmentation algorithm is intended for application to very large LULC datasets including 30 m U.S.-wide NLCD, 30 m global GLC30 (http://www.globallandcover.com/), 300 m global GlobCover (http://due.esrin.esa.int/globcover/), 300 m global CCI-LC (http://maps.elie.ucl.ac.be/CCI/viewer/), 100 m Europe-wide CORINE, and the Earth Observation for Sustainable Development of Forests (EOSD) (http://tree.pfc.forestry.ca/) which is a map of forest cover in the entire Canada and has the resolution of 25 m. The spatial arrangement (pattern) of land cover affects many environmental processes including ecological processes, biophysical phenomena, and wildlife habitat quality. In addition, mapping of land cover patterns is needed for resource management in order to achieve a balance between ecological goods and services. However, manual delineation of land cover patterns (Wickham and Norton, 1994) is impractical except for a very limited geographical extent. Thus, delineation of land cover pattern units through algorithmic segmentation of LULC raster is a sought-after capability.
In addition, our algorithm is also intended to be applied to non-LULC categorical datasets. For example, it could be applied to topographic datasets (DEMs), such as the U.S. National Elevation Dataset (NED) (http://nationalmap.gov/elevation.html) or the world-wide Shuttle Radar Topography Mission (SRTM) (http://www2.jpl.nasa.gov/srtm/), after they are categorized into landform types using a technique, such as, for example, the geomorphons (Jasiewicz and Stepinski, 2013). It could also be also applied to high-resolution EO RGB images, like the 1 m RGB mosaics of the U.S. (http://viewer.nationalmap.gov/launch/), after they are categorized by quantizing colors to several representative classes (Rubner et al., 2000). Segmenting a high-resolution image covering an entire city provides means for automated delineation of urban structure types (USTs).
Previous research on algorithmic identification and delineation of homogeneous patterns in large categorical rasters is limited and focuses on clustering of local patterns rather than on segmentation (Cardille et al., 2009;Partington and Cardille, 2013). For large rasters, a segmentation has an advantage over clustering because it can be applied to larger datasets and because a degree of pattern homogeneity of the resulting units is controlled. Thus, segmentation is a preferred technique to regionalize patterns in EO datasets. Most of an extensive literature on segmentation pertains to natural (RGB) images (Li et al., 2014) and remotely sensed images (Dey et al., 2010) with focus on achieving homogeneity of segments with respect to color and/or texture. Note that whereas texture is a structure of pixels arranged quasi-randomly and lacking geometric quality, a pattern is a perceptual structure, placement, or arrangement of objects having a geometric quality. Thus, existing, image-oriented segmentation solutions (for example, the popular JSEG algorithm (Deng and Manjunath, 2001)) are not designed to delineate segments with the cohesive pattern even if they are capable of delineating segments with homogeneous texture.
A pattern-oriented segmentation algorithm capable of being applied to large categorical rasters must be computationally efficient. Thus, an algorithm cannot directly use every cell in the raster, instead, it works on local blocks of cells using only histograms of constituent pattern features. Although there are many possible distance measures between histograms (Cha, 2007), experiments on an agreement between those measures and visual perception of similarity between patterns (Rubner et al., 2001) suggests that the Jensen-Shannon divergence (JSD) measure (Lin, 1991) is preferable over the Euclidean distance. Using JSD means that our algorithm must work in a distance space (Ganti et al., 1999) where calculating a distance between two patterns is the only allowable operation. For this reason, we needed to re-design concepts of dissimilarity measure, linkage, and inhomogeneity (all fundamental to a segmentation algorithm) so they can be calculated exclusively from distances between patterns.
The proposed algorithm is based on the principle of seeded region growing (SRG) (Adams and Bischof, 1994) but introduces a number of novel solutions in order to work with patterns. In particular, (1) Algorithm's smallest processing element is not a raster cell but a data object called motifel  -a square block of cells which encapsulates a local pattern of categories. The raster is transformed into a grid of motifels resulting in a reduction of dimensionality by orders of magnitude from the original raster thus making the algorithm computationally efficient.
(2) Because motifels may have a significant spatial extent, we introduce a "brick wall" topology for the spatial organization of the grid of motifels in order to decrease artifacts associated with tessellation of geographical space.
(3) We introduce a novel method of constructing a priority queue for potential seeds (a necessary feature of the SRG algorithm) and an adoption of locally determined growth-stopping criteria for regions growing from these seeds.
The algorithm is implemented as a module within the Geospatial Pattern Analysis Toolbox (GeoPAT)  -an open source collection of GRASS GIS (Neteler et al., 2012) modules for carrying out pattern-based analysis of categorical rasters. The GeoPAT software, including the segmentation module described in this paper, can be downloaded from http://sil.uc.edu/.

Basic concepts
Because the presented algorithm is implemented as a module of the GeoPAT software  it is based on the same core concept of pattern-based analysis of categorical rasters as the rest of GeoPAT modules.
In this concept  the raster is first divided into a regular grid of motifels. While a motifel is geometrically simple (it is a square array of k × k cells) it holds a complex content -a pattern (or motif) corresponding to composition and spatial configuration of k 2 category-labeled cells within its interior. Motifel's linear size, k, sets the scale over which the pattern is sampled. For example, if k = 500 NLCD cells, segmentation algorithm identifies units exhibiting homogeneous patterns having a characteristic scale of up to ∼225 km 2 . Regardless of the value of k, the information from all cells in the raster is used to construct the grid of motifels, and, ultimately, to segment the raster. The number of segments, their character, and their spatial extents change with motifel's size. Performing segmentation using series of different motifel sizes yields a series of non-hierarchical segmentations accounting for multiple scales in the landscape.

Brick wall grid of motifels
Unlike cells, motifels may have relatively large spatial extent so they visibly tessellate geographical space and a topology of their grid matters for the appearance of the resulting map of pattern units. The standard rectangular topology of a grid of motifels has significant shortcomings in application to segmentation. Rectangular grid topology with 4connectivity ( Fig.2(A)) is not isotropic. Consequently, segments can only grow horizontally or vertically resulting in a poor capture of diagonally oriented objects (see Fig.2(E)). Rectangular grid topology with 8-connectivity ( Fig.2(B)) has two types of adjacency, edge-to-edge and vertex-to-vertex. The existence of vertex-to-vertex adjacencies allows for the possibility of forming segments that permeate each other (an extreme example would be a chessboard-like permeation of a segment consisting of "white" motifels with a segment consisting of "black" motifels). Because of this undesirable property 8-connectivity is rarely used in segmentation algorithms. The grid of hexagonal motifels with 6-connectivity ( Fig.2(D)) is isotropic and has no vertexto-vertex adjacencies but a hexagonal motifel cannot be built from an integer number of raster cells.
Our solution is to use the 6-connected brick wall topology grid of square motifels ( Fig.2(C)). We call this topology "brick wall" because square motifels are "laid" in alternative layers with each layer shifted a half motifel length with respect to the previous one like in masonry. This grid is a closer approximation of an isotropic hexagonal grid than a 4-connected rectangular gird is, does not have vertex-to-vertex adjacencies like the 8-connected rectangular grid, and, unlike in a hexagonal grid, its motifels contain an integer number of raster cells. Using brick wall topology our algorithm can form diagonally oriented segments (see Fig.2(F)).

Pattern signature
Numerical representation (signature) of the pattern within a motifel must convey well the salient features of types of patterns present in the raster. Different datasets may require different choices of signature, the only requirement of GeoPAT is that the signature is expressed in the form of a histogram. Thus, after a raster is partitioned into a grid of motifels the next step is to convert this grid into a grid of histogram data.
Our segmentation algorithm is implemented with two histogram-yielding signatures, one based on pattern features derived from a category co-occurrence matrix, and another based on pattern decomposition features. Both signatures are described in Jasiewicz et al. (2015). Briefly, in the category co-occurrence signature (Barnsley and Barr, 1996;Chang and Krumm, 1999) pattern features are the pairs of raster categories assigned to two neighboring cells. Histogram counts and bins the features from eight co-occurrence matrices calculated for eight different displacement vectors along principal directions (see Niesterowicz et al. (2016) for an illustrative example). The result is a histogram with (C 2 + C/2 bins, where C is the number of raster categories; for example, the co-occurrence histogram describing patterns in the NLCD has 136 bins because C = 16. In our evaluations the co-occurrence signature worked well when used in segmentation of LULC and categorized RGB image datasets, but it did not worked equally well when used in segmentation of categorized topography datasets due to a qualitatively different character of topographic patterns. In the decomposition signature (Remmel and Csillag, 2006) a motifel is assumed to have a linear size k = 2 L cells. A motifel is subdivided by a series of partitions into squares with linear sizes of w = 2 i cells where i = 1, . . . , L. Decompositions levels are labeled 0, 1, . . . , L − 1, where 0 corresponds to the entire motifel, and L − 1 corresponds to partitions into squares having size of 2 × 2 cells (there are L 2 /4 such squares in a motifel). Pattern features are category cells, but they are histogrammed L times, once for every level of decomposition. At each level two-dimensional histogram (category, abundance type) is calculated taking into account all partitions. There are only three abundance types, small -less than 1/4 of all cells in the partition, medium -between 1/4 and 1/2 of all cells in the partition, and large -more than 1/2 of all cells in the partition. The total length of histogram is L × C × 3, so, for example, the decomposition histogram describing patterns in the NLCD using motifels having linear size of 256 = 2 8 cells (∼8 km) has 8×16×3 = 384 bins. The decomposition signature informs about the composition of categories at various scales within a motifel. In our evaluations the decomposition signature worked well with all datasets we have tested. Fig.3 illustrates construction of decomposition signature using as an example two NLCD motifels having linear sizes of 2 5 = 32 cells; the first four decomposition levels are shown, partitions at the highest shown level have sizes of 4 × 4 cells. The middle row in Fig.3 shows the parts (48 bins each) of the histograms corresponding to a given level. Consider the brown section of the histogram (NLCD category -cultivated crops) in the level 3. Each bin contains a total number of cells in each abundance type. This category occupies large sections of both motifels, consequently most cells are assigned to "large" type, and only few are assigned to the other two types. However, for the medium green section of the histogram (NLCD category -deciduous forest) motifel 1 is dominated by "large" type but motifel 2 is dominated by "small" and "medium" types. The entire histogram (over all levels of decomposition) representing each motifel is a concatenation of its parts representing different levels of decomposition.

Dissimilarity measure
We use the Jensen-Shannon Divergence (JSD) (Lin, 1991) as a measure of dissimilarity between two motifels represented by corresponding normalized histograms M 1 and M 2 . Normalized histogram means that the sum of all its bins equals to 1. The JSD expresses the informational distance between the two histograms as a deviation between Shannon's entropy of the conjugate of the two histograms (M 1 + M 2 )/2 and the mean entropy of individual histograms M 1 and M 2 . The value of JSD, denoted by d(M 1 , M 2 ), is given by the following formula: where m i is the value of ith bin in the histogram M and |M| is the number of bins (the same for both histograms). For normalized histograms, the JSD dissimilarity always takes values from 0 to 1 with the value of 0 indicating that two motifels are identical, and the value of 1 indicating maximum dissimilarity (none of the classes existing in one motifel can be found in the other).

Linkage and inhomogeneity
The segmentation algorithm not only requires calculating a value of dissimilarity between two motifels but also a value of dissimilarity between two segments (sets of motifels), which we refer to as a linkage. Consider two segments, S 1 = {M 1,1 , . . . , M 1,k1 } consisting of k1 motifels and S 2 = {M 2,1 , . . . , M 2,k2 } consisting of k2 motifels. To measure a dissimilarity between these two segments we use the so-called average linkage or Unweighted Pair Group Method with Arithmetic Mean (UPGNA) (Sokal and Michener, 1958) given by where function d(x, y) is given by eq.(1). The value of D(S 1 , S 2 ) has a range between 0 and 1 because the values of d are restricted to this range.
Inhomogeneity is a property of a single segment; it measures a degree of mutual dissimilarity between all motifels within the segment. As a measure of inhomogeneity we use an average distance between all distinct pairs of motifels in a segment. For a segment S = {M 1 , . . . , M k1 } the inhomogeneity is given as: as there is k1(k1 − 1) distinct pairs of motifels in the segment S . The value of δ has a range between 0 and 1 because values of d are restricted to this range. The small value of δ indicates that all motifels in the segment represent consistent patterns so the segment is patternhomogeneous. Note that segment is considered homogeneous even if its constituent motifels represent complex patterns of categories as long as the pattern of this complexity is approximately the same among all motifels within a segment.

Segmentation algorithm
An overall framework of our pattern-based segmentation algorithm adheres to an existing framework for modern SRG algorithms as applied to segmentation of either natural (RGB) images (for example, see Shih and Cheng (2005)) or remotely sensed images (for example, see Wang and Boesch (2007)) with the following modifications: (a) input data is restricted to categorical rasters, (b) its smallest processing unit is a complex and spatially extent motifel, not a simple cell; this requires using different dissimilarity measures, (c) motifels are organized into a brick wall grid whereas cells are organized into a rectangular grid, and (d) it uses custom method for construction of a priority queue for potential seeds and for calculating seed-specific criteria for stopping the growth of segments. Fig.4 shows a diagram outlining consecutive steps of our algorithm.

Automated ranking of motifels as seeds
The quality of any SRG algorithm depends on the selection of the seeds and on the order in which these seeds are grown into segments. In an unsupervised method, a queue of seeds must be determined automatically from the data (Shih and Cheng, 2005). In our algorithm seeds are individual motifels. The first task of our algorithm is to organize all motifels into a priority queue from which subsequent segments are grown. This part of an algorithm not only outputs the queue of potential seeds but also assigns to each seed a unique threshold  value for stopping the growth of a segment started from this seed.
The idea is that seeds should be queued in the increased order of the inhomogeneity of their immediate neighborhoods. Thus, a priority to grow into segments should be given to seeds located in the parts of the raster exhibiting high homogeneity of a pattern. Also, the segment's growth stopping criterion should be proportional to the inhomogeneity of the neighborhood of the seed from which the segment is grown. This is to say that when a segment starts to grow from a motifel locally surrounded by highly similar motifels it should only expand by accretion of additional highly similar motifels -the growth stopping threshold assigned to such seed should be low. However, when a segment starts to grow from a motifel locally surrounded by less similar motifels it could expand by accretion of additional motifels of lesser similarity -the growth stopping threshold assigned to such seed should be higher.
We assume that an immediate neighborhood, denoted by Ω n , of the nth motifel consists of 18 other motifels arranged around the focus motifel as shown in Fig.5 (2-motifel radius neighborhood in the brick wall topology). The size of the neighborhood is a parameter, we have selected 18 motifels to go one step above the six immediate neighbors. Selecting a larger neighborhood would result in a longer calculation time. The neighborhood may consist of a peer group, Ω n,1 , of motifels highly similar to the focus motifel, and another group, Ω n,2 , of motifels less similar to the focus motifel. The goal is to assign to the focus motifel an inhomogeneity value based on Ω n,1 alone disregarding motifels in Ω n,2 . In order to identify motifels in Ω n,1 we use the method described in Deng et al. (1999) but applied to motifels instead of color pixels. The neighborhood motifels are first sorted in ascending dissimilarity from the focus motifel. Second, the criterion proposed by Deng et al. (1999) is used to find the peer group Ω n,1 . Note, however, that whereas Deng et al. (1999) used an absolute maximum of their criterion to determine the peer group, we use the first maximum (the smallest possible peer group). Fig.5 shows four examples of neighborhoods around focus motifels divided by heavy orangecolored polygons into Ω n,1 (inside the orange-colored polygon) and Ω n,2 .
Once the peer group is determined, its inhomogeneity µ(Ω n,1 ) is calculated as µ(Ω n,1 ) = D(M n , Ω n,1 ), where M n is the focus motifel. Note, that for the purpose of In each example, the focus motifel is indicated by the heavy black-colored outline and its peer group is indicated by the heavy orange-colored outline. The numbers beneath neighborhoods give values of peer group inhomogeneity and standard variation of dissimilarities between the focus motifel and members of its peer group, respectively. Legend to NLCD is given in Fig.1.
forming the priority queue we measure inhomogeneity as an average dissimilarity of the peer group to the focus motifel and not as an average mutual dissimilarity between peer group motifels. This is because we want the nucleus of the forming segment to be as similar to the seed as possible. Standard variation, σ(Ω n,1 ), of all dissimilarity values between motifels in Ω n,1 is also calculated; the sum T n = µ(Ω n,1 ) + σ(Ω n,1 ) is a threshold for stopping the growth of the segment started using the focus motifel as a seed. The two numbers given in Fig.5 beneath examples of neighborhoods are µ(Ω n,1 ) and σ(Ω n,1 ), respectively. We observe, that, as designed, the value of T n increases with increasingly inhomogeneous neighborhoods.

Segment growing
The process of growing segments proceeds according to the priority queue. A new segment starts from a motifel currently at the top of the queue. It grows by incorporating available motifels (those that have not been previously incorporated into segments) from among its current perimeter until growth stopping criterion is met.
For all available motifels belonging to the current segment's perimeter, the algorithm calculates the value of D(motifel, current segment) and selects the one with the smallest value, D min , as a candidate for possible absorption into the segment. The absorption occurs if D min < T , where T is a segment growth threshold. Recall from the previous subsection that each seed is assigned its own growth threshold T n . However, the algorithm also features two static thresholds, T min and T max whose values need to be set before segmentation starts. The value of actual threshold T is as follows The purpose of T min is to control (to some degree) sizes of the segments, if T min > T n it overwrites the seedspecific threshold relaxing the growth stopping criterion and allowing the formation of larger segment. The purpose of T max is to prevent the growth of the segment from the seed surrounded by an exceptionally inhomogeneous neighborhood. Note that a segment grows one motifel at the time with its perimeter recalculated after every absorption (see Fig.6). Once the segment is fully grown, its constituent motifels are deleted from the priority queue and the next segment starts to be grown starting from the new top of the queue. The segment growing procedure ends when all motifels are incorporated into segments.
At that point, the segmentation process could be considered as finished but our algorithm also includes three optional steps designed to further improve the quality of the segmentation, they are: segment merging, removal of small segments, and adjustment to boundaries between segments.

Segment merging
The segment growing procedure described in the previous subsection may result in the over-segmentation of the raster. This means that some adjacent segments are sufficiently similar to justify merging them into a single segment. This situation occurs because growing procedure considers for accretion only motifels from the perimeter of the growing segment possibly missing a large set of very similar motifels located just outside the perimeter. Segment merging procedure is designed to fix this problem thus reducing the over-segmentation of the raster (Shih and Cheng, 2005). Segment merging is based on an undirected graphbased process. The segments are graph's nodes while graph's edges indicate segments adjacencies. The weight of each edge is set by the value of dissimilarity between the two nodes (adjacent segments) calculated using the linkage D(S 1 , S 2 ) (eqn.3), where, in this context, S 1 and S 2 indicate adjacent segments. The merging of segments is tantamount with graph pruning. Only edges fulfilling the criterion are considered for pruning. In this criterion T S 1 and T S 2 are the growth stopping thresholds (possibly modified by T min or T max ) used to growth segments S 1 and S 2 , respectively. The pruning is performed iteratively, at each iteration the edge with the smallest weight is identified and eliminated (the segments are merged). Next, the graph is recalculated resulting in the new set of edges with the new values of weights and the procedure is repeated until no more pruning (segment merging) is possible. Fig.7 illustrates the hypothetical case of merging segments. On the left of the figure the part of the segmented raster is shown with segments labeled by letters from A to H highlighted for consideration. The adjacencies between segments that fulfills the criterion given in eq.6 are indicated by pairs of arrows and the values of the linkage. The smallest linkage is between the segments A and B. On the right of the figure the several iterations of merging are illustrated in the form of the graphs. It can be seen that segments A and B are merged in the first iteration, the segments AB and D are merged in the second iteration, and the segments ABD and F are merged in the third iteration. The segments ABDF and G are merged next (not shown). No more merges are possible because the edges ABDFG-C, ABDFG-E, and ABDFG-H have the weight too large to fulfill the criterion given by eq.6.

Removal of small segments
After completion of growing and merging processes further (optional) improvement of the segmentation's utility can be achieved by removal of small segments. Some small segments are relevant and need to be kept while others are more of a nuisance and could be incorporated into a neighboring, larger segment without much loss of information but resulting in a gain in the segmentation's utility. Our algorithm removes segments which are smaller than T A motifels, where T A is a free parameter. The removal of small segments is performed using an undirected graph-based process (just like in the merging process). Only edges fulfilling the criterion A S 1 or A S 2 ≤ T A and D(S 1 , S 2 ) ≤ √ min[T S 1 , T S 2 ] (7) are considered for pruning, where A S 1 and A S 2 are areas (in the units of motifels) of segments S 1 and S 2 , respectively. If the criterion above is fulfilled the smaller segment is incorporated into the larger segment. The removal is performed iteratively, at each iteration the small segment edge with the smallest weight is identified and eliminated. Next, the graph is recalculated and the procedure is repeated until there are no more eligible small segments to remove.

Adjustment to boundaries between segments
The final segmentation quality enhancement procedure is an adjustment to boundaries between segments. Because the segment growing is a greedy process there is a chance that some motifels located at segments' boundaries show more affinity to motifels across the boundary than to motifels in its own segment.
To alleviate this problem we adjust the boundaries between the segments. For this procedure, only border If ∆D = D par − D neigh is positive (the linkage to the segment across the boundary is smaller than the linkage to parent segment) the motifel is marked for possible reassignment (resulting in an adjustment to the boundary). All marked motifels are sorted in an ascending order of ∆D and the top motifel in this ordering is reassigned. Linkages affected by the change to the two segments are recalculated resulting in a possible change to the sorted list of marked motifels. The procedure is repeated until no more marked motifels are left.

Key parameters
Table I summarizes the key parameters of our segmentation algorithm. The last three columns give the possible range of each parameter, the default value, and values which have special meaning. Default values are established on the basis of experimentation with four datasets discussed in section 4. For each dataset, we run our segmentation algorithm several times with different set of parameters. Each segmentation was visually and quantitatively assessed (see section 4 for details). Default values represent optimal choices derived from these experiments. Note that no single set of parameters is optimal for all conceivable datasets. Most of parameters in Table I (k, T min , T max , and T A ) have been already discussed in previous sections. The raster, and thus some of its motifels, may contain cells with nodata values. If the motifel contains too many such cells the entire motifel is considered as null, but if it contains a relatively small number of nodata cells, it is considered a normal motifel and its histogram is calculated from available cells. The parameter n min controls this transition, the default value is 0.5 (50%). Calculating linkage between two segments (eq.3) could be computationally expensive if the segments contain large numbers of motifels. When the number of motifels in the segment exceeds N lev we resort to so-called leveraging to speed up the calculation. This means that only N lev motifels, randomly selected from the segment, are used to calculate the linkage. Finally, we adjust the boundary between segments only when ∆D > ∆D min , or when the motifel has significantly more affinity to the neighboring segment than it has to its own segment.

Testing segmentation algorithm
In this section, we demonstrate the performance of our algorithm by applying it to four pertinent datasets representing different thematic meanings. Note that our algorithm is not meant to be a general purpose, instead, it suppose to be applied only to large categorical datasets in which categories form spatial patterns having physiographic, environmental, or social meaning. For examples of such datasets, we demonstrate that our algorithm can delineate patterns in a relatively short computational time and that the resulting segments are of high quality using both, visual and quantitative assessment.
The first dataset is the NLCD2011 (land cover theme) over the entire conterminous U.S. NLCD cells have 30 m resolution and there are 16 different land cover categories. The second is the NED (topography theme)a 30 m resolution DEM covering the entire conterminous U.S. We converted NED to a 30 m resolution map of ten most common landform elements using the geomorphons method (Jasiewicz and Stepinski, 2013). The third is the EOSD (land cover theme), the 25 m resolution, 23 categories land cover map that covers forest over the entire Canada. The fourth is the high resolution (1m) RGB image (image theme) quantized to 26 color-classes and covering the Los Angeles, CA area. Table 2 gives the sizes of these datasets as well as the number of resultant segments depending on the selected size of motifel. The last column in Table 2 gives the total time of calculation. All calculations were performed on a computer with Intel 3.4GHz, 4-core processor and 16 GB of memory running the Linux operating system. Default parameters were used. For NLCD and topography, the decomposition signature of patterns was used, whereas for the EOSD and for high-resolution image the co-occurrence signature was used.
As can be seen from Table 2 the total computation time is short considering size of datasets. For a given dataset the major determinant of computation time is the motifel's size; the larger the motifel the smaller the grid of motifels and the shorter the computation time.
The short computation times are due to relatively small sizes of grids of motifels. For example, for the 128× 128 motifel the size of the grid of motifels for the NLCD is 812×1256. Thus, although pattern-based segmentation is applied to very large datasets, grids that are actually segmented are not that large because they consist of motifels -objects constructed from thousands to hundreds of thousands of cells. Fig.8 shows small portions of the NED, EOSD, and high-resolution RGB image datasets with the boundaries of segments superimposed to allow for a visual assessment of a degree to which segments delineate homogeneous patterns. Fig.1(B) can be used for visual assessment of NLCD segmentation. Readers can also visually assess segmentations of the entire NLCD and NED by downloading the full results from http://sil.uc.edu/. Note that visual assessment, although subjective and tedious, is the most common method of evaluating the effectiveness of a segmentation (Zhang et al., 2008), and, if practical, should always be taken into consideration. Fig.8(A) shows a segmentation of a portion of the map of geomorphons (calculated from a NED). Ten different colors indicate ten common landforms (Jasiewicz and Stepinski, 2013). We don't include here the geomorphons legend as the precise meaning of different landforms is not relevant for assessing the segmentation results. Segments ( indicated by black lines) delineate homogeneous topographic landscapes or physiographic units. We observe that segments enclose homogeneous patterns of colors, which translates into enclosure of homogeneous topographic landscapes because patterns of colors visualize the patterns of landscape elements. Thus, segments in the panel A delineate various physiographic units. We stress that, in this example, and throughout the rest of the paper, we don't use color for any calculations, we only use it to visualize the raster's categories. Fig.8(B) shows a segmentation of a portion of the land cover map (EOSD) of Canada. EOSD concentrates on the forest, thus most of its 23 categories pertain to various types of forest; for a legend see Wulder et al. (2008). Segments in Fig.8(B) delineate regions exhibiting homogeneous patterns of forest types (forest units).
Finally, Fig.8(C) shows a segmentation of (quantized) high-resolution RGB image depicting a small portion of the Los Angeles area. This is an actual image but with only 26 colors. The quantized colors are used as abstract categories (no similarity between different colors are taken into account). Segments delineate urban structure types (USTs) -distinct spatial patterns of the urban structure at the neighborhood scale, which can be interpreted in terms of the type of activity or of a residential pattern.

Dependence on the growth threshold
The two parameters of our segmentation algorithm which are most likely to be adjusted are the size of the motifel (k) and the lower growth threshold (T min ). In this subsection we discuss dependence of segmentation results on the choice of T min . Recall that T min can relax the seed-specific growth stopping criterion, in effect making it easier for a segment to grow to a larger size. Larger segments are desirable because a map with too many small segments starts to loose its utility. On the other hand, larger segments are inevitably less homogeneous. We calculated segmentations of a 16896×16192 cells subset of NLCD2011 (denoted in Table 2 as NLCD (exmp.)) using 128×128 motifel size and two values of T min , 0.15 and 0.25. This subset covers the entire state of Wisconsin and portions of Iowa and Minnesota. With our choice of motifel size the region is tiled into 132×126=16632 local landscapes, each having a scale of ∼4 km.
The isolation of a segment is calculated as an average linkage (eq.3) between the focus segment and all of its neighbors. It measures how much the segments stand out from its surroundings. The ratio of inhomogeneity to isolation (δ/γ) captures both metrics of segment's "goodness" in a single number -segment quality; we actually use 1 − (δ/γ) so increased values indicate increased quality. The metrics δ, γ and s are calculated for each segment in the region so their spatial dependence is mappable. Average values of these metrics over all segments are used to compare different segmentations. Fig.9 shows the segmentation of the NLCD example region using T min = 0.15 together with the maps show-ing spatial variation of the metrics' values. Fig.9(A) shows segments boundaries (in black) overlaying the land cover map; this is intended for visual assessment of segmentation. There are 605 segments. Fig.9(B) shows the segments color-coded according to their values of δ, smaller values (greener colors) indicate more homogeneous segments. The range of the values of δ is between 0 and 0.32 with the mean over all segments ⟨δ⟩ = 0.13. Fg.9(C) shows the segments color-coded according to their values of γ, larger values (greener colors) indicate larger isolation. The range of the values of γ is between 0.1 and 0.84 with the mean over all segments ⟨γ⟩ = 0.3. Fig.9(D)   Comparing the metrics for the two segmentations we find that the segmentation with T min = 0.15 is characterized by slightly better values of homogeneity while the segmentation with T min = 0.25 is characterized by slightly better values of isolation, so the overall segmentation quality, as measured by the metrics, is about the same. Spatial distributions of metrics' values are also very similar in both segmentations. This indicates that the segmentation algorithm is stable with respect to this important parameter. It also means that a choice of what value of T min to use is likely to come down to subjective, visual examination.

Dependence on motifel's size
Size of motifels sets the spatial scale over which pattern is perceived. Thus, segmenting with different values of motifel's size leads to segments defined by homogeneity of different patterns. Because of this, the segmentations with different values of k do not form a hierarchy. To demonstrate how segmentations change with the value of k we calculated segmentations for a subset of the region featured in Fig.9. The only reason for using a relatively small region for this demonstration is so the illustration fits on the journal page. We calculated segmentations using motifel's size of 64×64, 128×128, 256×256 cells, and the same set of remaining parameters. Fig.11 has six panels. The first column (Fig.11(A) and Fig.11(D)) shows segmentation boundaries and assessment of segments' inhomogeneities, respectively, for segmentation with k = 64. Local landscapes (patterns of land cover) are assessed on the scale of ∼2 km. There is a lot of variability in the landscape on such small scale resulting in the large number of segments. Many segments are small because many patterns are spatially restricted as they are location-specific. However, there are also larger segments when the pattern on the ∼2 km scale persists over larger regions. The homogeneity of segments ( Fig.11(D)) is well-controlled, but it does not appear to be proportional to segments' sizes. This is again due to high variability of landscapes on the scale of ∼2 km.
The third column (Fig.11(C) and Fig.11(F)) show segmentation boundaries and assessment of segments' inhomogeneity, respectively, for segmentation with k = 256. In this case local landscapes are assessed on the scale of ∼8 km (16 times bigger area than in the case of ∼2 km landscapes). In this segmentation the segments are larger not only because they are built from larger motifels but also because, at the larger scale, there is less variability between local landscapes within the region. The homogeneities of segments are high and in apparent inverse proportion to their sizes, which is what we expect.
The second column (Fig.11(B) and Fig.11(E)) show segmentation boundaries and assessment of segments' inhomogeneity, respectively, for segmentation with k = 128. The appearance of this segmentation is somewhere in-between the two segmentations calculated for k = 64 and k = 256. The inhomogeneity of segments ( Fig.11(E)) is well-controlled with small segments having predominantly high homogeneity, but the least homogeneous segments do not coincide with the largest segments. What size of motifel to chose for segmentation depends on its purpose and needs to be selected in consultations with domain experts.

Summary and Conclusions
Although segmentation is one of the most widely researched and used tools in image analysis, the available algorithms are all tuned to data structures associated with images as they use pixels as basic units of analysis. Our purpose in this paper was to describe a unique algorithm especially designed for a niche application of segmenting a grid where basic units of analysis are larger tracts of land containing patterns of pixels. Such algorithm enables an automatic delineation of regions characterized by stationary patterns of the variable of interest. To the best of our knowledge, no other presently available algorithm can perform such task. Although a number of datasets most-suited for being segmented by our algorithm is limited, they all are of high importance for analyzing large-scale environments. Depending on an input data, delineated segments may be interpreted as physiographic units, ecological units, urban structure units, etc.
We choose to base our algorithm on the SRG concept because it has been already thoroughly investigated in applications to images and because it is computationally efficient. Modifications necessary to make the SRG concept applicable to the grid of patterns were described in details in section 3. Some of these modifications involved the development of original concepts. One example of such concept is the brick-wall grid topology (section 2.1). This simple concept significantly improves segmentation quality at small computational cost. Another example is the method for seeds ranking and for setting an individual growing threshold for each segment based on the degree of pattern cohesion in a seed's neighborhood.
Because of the niche character of the algorithm, its validation (section 4) is carried out on subset of an actual datasets for which it was designed. Four datasets covering three pertinent themes (land cover, topography, high-resolution urban image) were segmented multiple times using different sets of parameters. The results of these experiments were threefold. First, we demonstrated that our algorithm is computationally efficient (Table II). For example, the entire NLCD is segmented in 12 minutes to 2 hours depending on the size of the motifel. Second, we demonstrated that our algorithm produces a good quality segmentation (see below). Third, we use the results of our experiments to set values of default parameters (Table I).
Segmentation is an unsupervised procedure, thus objective evaluation of its results is an issue (Zhang et al., 2008). Most segmentations are evaluated subjectively by an analyst who visually reviews the result and passes judgment based on domain knowledge. Another version of subjective evaluation is to compare algorithmic segmentation to a manual segmentation by an expert. This last method is not practical for the types of datasets we aim at because of their size and their character that makes manual segmentation impossible (see Figs. 1 and  8). Thus, we visually inspected segmentations (calculated with a default set of parameters) superimposed on the map depicting patterns of categorical variable and inspect a degree to which segments enclose discernible patterns. Based on such inspection we came to a conclusion that the algorithm successfully regionalizes different patterns. A reader can observe this in Figs. 1 and 8 or make a more thorough inspection by downloading (http://sil.uc.edu/) our segmentations of the entire NLCD and superimposing it on the NLCD map. Segmentations of the U.S. topography (as represented by geomorphons derived from the NED) and an associated map of geomorphons are also available for inspection.
We also performed a quantitative assessment of calculated segmentations using two metrics, inhomogeneity and isolation. These metrics were calculated for each segment. We found an average value of segment inhomogeneity to be 0.13. To get an indication of the meaning of this value a reader may visually inspect a segment enclosed by an orange-colored outline in Fig. 4(B) whose inhomogeneity is 0.1. We found an average value of segment isolation to be 0.3, and we found the average ratio of inhomogeneity to isolation to be about 2, indicating that motifels (local patterns) within a segment are twice as similar to each other than the segment as a whole is similar to its neighboring segments. Overall, both, visual and quantitative assessments validate that the algorithm successfully performs the task to which it was designed.
Future work will extend the functionality of the algorithm by allowing segmentation of geographical space based on multiple layers of information. Presented version of the algorithm segments geographical space on the basis of patterns of a single variable, for example, land cover or topography. Future work on our algorithm will aim at segmentation on the basis of patterns of multiple variables, for example, land cover and topography. Such an algorithm would be able to delineate segments characterized by homogeneous patterns of land cover as well as homogeneous patterns of topography. This would make it well-suited for delineation of ecoregions. evaluation: A survey of unsupervised methods. Computer Vision and Image Understanding 110 (2), 260-280.