Spatial segmentation in R using the {supercells} package

class: inverse, left, nonum, clear
background-image: url("figs/title-cover.jpg")
background-size: cover

.titlestyle2[Spatial segmentation in R] <br>
.titlestyle2[using the supercells package]

.captionstyle[Jakub Nowosad, nowosad.github.io]
.pull-right2[.captionstyle[OpenGeoHub Summer School 2021, 2021-09-02]]

---
# OBIA

Considering cells' values without its spatial context is often not enough to understand underlying objects or processes.

<div class="figure" style="text-align: center">
<img src="figs/blaschke.png" alt="Blaschke (2010), doi:10.1016/j.isprsjprs.2009.06.004" width="60%" />
<p class="caption">Blaschke (2010), doi:10.1016/j.isprsjprs.2009.06.004</p>
</div>

----

Generalization of spatial information is often associated with (geographic) **object-based image analysis (OBIA)**.
OBIA consists of two major stages:

1. **Segmentation** - partition of space to identify homogeneous objects
2. Classification - classification of the obtained objects

---
# Spatial segmentation

.pull-left[
OBIA applies many techniques, including the multiresolution segmentation (MRS) that uses cells as the underlying representation.

Recently, **an approach of superpixel become considered as a promising alternative**.
]

.pull-right[
<div class="figure" style="text-align: center">
<img src="figs/esch.png" alt="Esch (2008), doi:10.1109/LGRS.2008.919622" width="75%" />
<p class="caption">Esch (2008), doi:10.1109/LGRS.2008.919622</p>
</div>
]

---
# Superpixels

.pull-left[
<div class="figure" style="text-align: center">
<img src="figs/csillik.png" alt="Csillik (2017), doi:10.3390/rs9030243" width="100%" />
<p class="caption">Csillik (2017), doi:10.3390/rs9030243</p>
</div>
]

.pull-right[
**Main idea:** create groupings of cells with similar values, which result in an over-segmentation.

Each superpixel represents a desired level of homogeneity while at the same time maintains boundaries and structures.

Superpixels also carry more information than each cell alone, and thus they can speed up the subsequent processing efforts.

A large number of methods for creating superpixels were developed (Stutz et al (2018), doi:10.1016/j.cviu.2017.03.007).

The SLIC algorithm (Achanta et al. (2012), doi:10.1109/TPAMI.2012.120) is the most prominent. 
]

---
# SLIC

**SLIC starts with regularly located cluster centers** spaced by the interval of `$S$`.

.pull-left[
<img src="figs/slic1.png" width="90%" style="display: block; margin: auto;" />
]

.pull-right[
<img src="figs/slic2.png" width="90%" style="display: block; margin: auto;" />
]

---
# SLIC

.pull-left[

Next, **each cell is assigned to the nearest cluster center**, and **the distance `$D$` is calculated between the cluster centers and cells** in the `$2S \times 2S$` region.

$$
D = \sqrt{\left(\frac{d_c}{m}\right)^2 + \left(\frac{d_s}{S}\right)^2}
$$

where `$d_c$` is the color (spectral) distance, `$m$` is the compactness parameter, `$d_s$` is the spatial (Euclidean) distance, and `$S$` is the interval between the initial cluster centers.

]

.pull-right[
<img src="figs/slic2.png" width="90%" style="display: block; margin: auto;" />
]

---
# SLIC

The color (spectral) distance is calculated between values `$I(x_i,y_i,s_p)$` and `$I(x_j,y_j,s_p)$` for a spectral band `$s_p$` in the set of spectral bands `$B$`:

$$
dc = \sqrt{\sum{p \in B}{(I(x_i,y_i,s_p)-I(x_j,y_j,s_p))^2}}
$$

The spatial (Euclidean) distance between cells represents spatial proximity:

$$
d_s = \sqrt{(x_j - x_i)^2 + (y_j - y_i)^2}
$$

**The color distance controls the homogeneity of superpixels.**

**The spatial distance is related to spatial contiguity.**

---
# SLIC

.pull-left[

Afterward, **new cluster centers (centroids) are updated for the new superpixels**, and their color values are the average of all the cells belonging to the given superpixel.

**The SLIC algorithm works iteratively**, repeating the above process until it reaches the expected number of iterations.

Experiments of Achanta et al. (2012) showed that **between 4 and 10 iterations suffices in the case of RGB images**.
]
.pull-right[
<img src="figs/ga.gif" width="90%" style="display: block; margin: auto;" />
]

---
# SLIC

The last, optional, step **enforces the 4-connectivity of the cell belonging to the same superpixels** by reassigning disjoint cells.

.pull-left[
<img src="figs/slic3.png" width="90%" style="display: block; margin: auto;" />
]

.pull-right[
<img src="figs/slic4.png" width="90%" style="display: block; margin: auto;" />
]

---
# SLIC - two main parameters

.pull-left[
<img src="figs/slic5.png" width="65%" style="display: block; margin: auto;" />
`$S$` **controls the size of each superpixel**, which is directly related to the number of output superpixels.
]
.pull-right[
<img src="figs/slic6.png" width="65%" style="display: block; margin: auto;" />
`$m$`**, compactness, controls the influence of the spectral distance on the results.**

Its large values result in more regularly shaped superpixels, while lower values create more spatially adapted, irregularly shaped, superpixels.

]

---
# SLIC - software

.pull-left[
SLIC exist in [many](https://scikit-image.org/docs/dev/api/skimage.segmentation.html?highlight=slic#skimage.segmentation.slic) [image](https://www.mathworks.com/help/images/ref/superpixels.html) [processing](https://docs.opencv.org/3.4/df/d6c/group__ximgproc__superpixel.html) [libraries](https://cran.r-project.org/package=OpenImageR).

**SLIC has also been implemented in a number of GIS software**, including:

<a href="https://developers.google.com/earth-engine/apidocs/ee-algorithms-image-segmentation-snic" target="_blank"><img src="figs/gee.png" width="70%" style="display: block; margin: auto;" /></a>
]

.pull-right[
## SLIC's limitation

**The SLIC algorithm**, in principle, **enables creation of superpixels based on less and more than three variables**.

However, **both color (spectral) and spatial distances are based on the Euclidean distance**.

Using the Euclidean distance to calculate color distances is adequate in many cases, however, **it limits the possible usability of the SLIC algorithm for spatial raster data**.
]

---
# supercells <img src="figs/logo.png" align="right" width="400" />

**The supercells package aims to utilize the concept of superpixels to a variety of spatial data.**

This package works on spatial data with **one variable** (e.g., continuous raster), **many variables** (e.g., RGB rasters), and **spatial patterns** (e.g., areas in categorical rasters).

In this extension, **the distance between values can be replaced with any applicable dissimilarity measure**.

<br><hr><br>
**supercells** is not, yet, on CRAN

```r
install.packages("supercells", repos = "https://nowosad.r-universe.dev")
```

(Any ideas, comments, etc. are welcomed!)

---
# Basic example (I)

.pull-left[

```r
library(supercells)    # superpixels for spatial data
library(terra)         # spatial raster data
library(sf)            # spatial vector data
```

```r
vol = rast(system.file("raster/volcano.tif",
                       package = "supercells"))
plot(vol)
```
]
.pull-right[
<img src="index_files/figure-html/index-4-1.png" style="display: block; margin: auto;" />
]

---
# Basic example (II)

```r
vol_slic1 = supercells(vol, k = 50, compactness = 1) # creating superpixels
```

```r
vol_slic1
```

```
## Simple feature collection with 49 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 2667400 ymin: 6478705 xmax: 2668010 ymax: 6479575
## Projected CRS: NZGD49 / New Zealand Map Grid
## First 10 features:
##    supercells       x       y elevation                       geometry
## 1           1 2667491 6479499  95.59932 POLYGON ((2667400 6479575, ...
## 2           2 2667472 6479376 100.08730 POLYGON ((2667560 6479445, ...
## 3           3 2667464 6479282 105.08633 POLYGON ((2667550 6479395, ...
## 4           4 2667462 6479182 109.57724 POLYGON ((2667510 6479305, ...
## 5           5 2667487 6479142 117.17857 POLYGON ((2667530 6479255, ...
## 6           6 2667460 6479004 125.84211 POLYGON ((2667490 6479115, ...
## 7           7 2667447 6478873 111.01835 POLYGON ((2667400 6478995, ...
## 8           8 2667462 6478759 107.58741 POLYGON ((2667400 6478875, ...
## 9           9 2667575 6478816 137.95283 POLYGON ((2667460 6478935, ...
## 10         10 2667546 6478784 119.21839 POLYGON ((2667460 6478895, ...
```

---
# Basic example (III)

.pull-left[

```r
plot(vol)
plot(st_geometry(vol_slic1),
     add = TRUE, lwd = 0.5)
```
]
.pull-right[
<img src="index_files/figure-html/bep2-1.png" width="672" style="display: block; margin: auto;" />
]

---
# Basic example (IV)

.pull-left[

```r
vol_slic2 = supercells(vol, k = 50, compactness = 20)
```
<img src="index_files/figure-html/index-9-1.png" style="display: block; margin: auto;" />
]
.pull-right[

```r
vol_slic3 = supercells(vol, k = 500, compactness = 1)
```
<img src="index_files/figure-html/index-11-1.png" style="display: block; margin: auto;" />
]

---
# RGB example (I)

```r
ortho = rast(system.file("raster/ortho.tif", package = "supercells"))
plot(ortho)
```

---
# RGB example (II)

```r
ortho_slic1 = supercells(ortho, k = 2000, compactness = 10, transform = "to_LAB")
```

```r
ortho_slic1
```

```
## Simple feature collection with 1784 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 358900.8 ymin: 505399.5 xmax: 359774.8 ymax: 505799.5
## Projected CRS: ETRS89 / Poland CS92
## First 10 features:
##    supercells        x        y   X735_1   X735_2   X735_3
## 1           1 358907.6 505792.5 171.2376 172.1786 164.9187
## 2           2 358907.8 505768.7 199.4737 198.9204 188.7943
## 3           3 358910.4 505757.9 137.2053 146.8596 135.2775
## 4           4 358916.0 505745.3 167.4763 170.6358 168.8508
## 5           5 358903.8 505733.7 170.8143 172.7115 168.2873
## 6           6 358908.7 505711.2 176.1361 176.4549 170.7990
## 7           7 358904.7 505691.0 167.7467 167.2519 157.6308
## 8           8 358908.7 505669.7 191.9513 193.7822 183.4240
## 9           9 358909.8 505657.2 157.3911 158.7681 152.9694
## 10         10 358906.0 505649.5 158.7732 134.6349 137.2844
##                          geometry
## 1  POLYGON ((358900.8 505799.5...
## 2  POLYGON ((358900.8 505777.5...
## 3  POLYGON ((358912.8 505767.5...
## 4  POLYGON ((358914.8 505753.5...
## 5  POLYGON ((358900.8 505747.5...
## 6  POLYGON ((358906.8 505723.5...
## 7  POLYGON ((358900.8 505699.5...
## 8  POLYGON ((358900.8 505681.5...
## 9  POLYGON ((358906.8 505663.5...
## 10 POLYGON ((358900.8 505655.5...
```

---
# RGB example (III)

```r
plot(ortho); plot(st_geometry(ortho_slic1), add = TRUE, lwd = 0.5)
```

---
# RGB example (IV)

.pull-left[
<img src="index_files/figure-html/index-17-1.png" style="display: block; margin: auto;" />

<img src="index_files/figure-html/index-18-1.png" style="display: block; margin: auto;" />
]

.pull-right[
*Converting average colors of each superpixel from RGB to a hexadecimal representation:*

```r
rgb_to_hex = function(x){
  apply(t(x), 2, function(x) rgb(x[1], x[2], x[3],
                                 maxColorValue = 255))
}
ortho_slic1_df = st_drop_geometry(ortho_slic1[4:6])
avg_colors = rgb_to_hex(ortho_slic1_df)
```
]

---
# Categorical example (I)

```r
library(motif)
lc = rast(system.file("raster/landcover2015.tif", package = "motif"))
plot(lc)
```

---
# Categorical example (II)

*Calculating proportions of land cover categories in spatial windows of 20 by 20 cells:*

```r
lc_comp = lsp_signature(lc, type = "composition", window = 20)
lc_comp = lsp_restructure(lc_comp)
lc_comp = lsp_add_terra(lc_comp)
lc_comp
```

```
## class       : SpatRaster 
## dimensions  : 191, 368, 9  (nrow, ncol, nlyr)
## resolution  : 6000, 6000  (x, y)
## extent      : -1091676, 1116324, -1184556, -38556.49  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=cea +lat_ts=5.5 +lon_0=140.8 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs 
## source      : memory 
## names       :           id,      na_prop,           X1,           X2,           X3,           X4, ... 
## min values  :            1,            0,            0,            0,            0,            0, ... 
## max values  : 7.028800e+04, 9.000000e-01, 1.000000e+00, 1.000000e+00, 1.000000e+00, 4.690909e-01, ...
```

---
# Categorical example (III)

```r
lc_comp_t = subset(lc_comp, 3:9); plot(lc_comp_t)
```

---
# Categorical example (IV)

```r
lc_slic1 = supercells(lc_comp_t, k = 1000, compactness = 0.03,
                      dist_fun = "jsd", clean = TRUE)
```

```r
lc_slic1
```

```
## Simple feature collection with 197 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -1091676 ymin: -1184556 xmax: 1116324 ymax: -38556.49
## Projected CRS: unnamed
## First 10 features:
##    supercells          x          y           X1        X2           X3
## 1           1 -1069695.3 -397344.95 0.2032499379 0.6577740 0.000000e+00
## 2           2 -1039509.4 -840389.82 0.1156097630 0.8139993 0.000000e+00
## 3           3 -1005004.2 -169267.42 0.0388159189 0.9086585 0.000000e+00
## 4           4  -931269.2  -69487.52 0.0006973719 0.9650979 0.000000e+00
## 5           5  -935945.8 -202073.34 0.0434398541 0.9414739 0.000000e+00
## 6           6  -946390.4 -311801.38 0.0036186524 0.8534272 0.000000e+00
## 7           7  -915787.2 -127156.49 0.0310925926 0.9638333 0.000000e+00
## 8           8  -868075.3 -360931.98 0.0315480475 0.9143126 5.928854e-05
## 9           9  -866937.0 -251643.44 0.0708165724 0.8008102 8.152174e-04
## 10         10  -875506.3 -628764.03 0.2544807610 0.5618048 0.000000e+00
##              X4 X5           X6          X7                       geometry
## 1  1.066861e-02  0 0.0030427372 0.125264729 MULTIPOLYGON (((-1061676 -8...
## 2  3.429355e-04  0 0.0005851882 0.069462786 MULTIPOLYGON (((-1025676 -7...
## 3  2.089414e-04  0 0.0009771445 0.051339467 MULTIPOLYGON (((-983676.1 -...
## 4  1.354058e-04  0 0.0000000000 0.034069367 MULTIPOLYGON (((-947676.1 -...
## 5  1.825843e-04  0 0.0002808989 0.014622766 MULTIPOLYGON (((-965676.1 -...
## 6  1.268359e-03  0 0.0000000000 0.141685780 MULTIPOLYGON (((-965676.1 -...
## 7  1.203704e-03  0 0.0024074074 0.001462963 MULTIPOLYGON (((-899676.1 -...
## 8  7.905138e-05  0 0.0075660597 0.046434931 MULTIPOLYGON (((-869676.1 -...
## 9  3.811406e-04  0 0.0489953686 0.078181515 MULTIPOLYGON (((-869676.1 -...
## 10 1.685905e-03  0 0.0023474186 0.179681117 MULTIPOLYGON (((-851676.1 -...
```

---
# Categorical example (V)

```r
plot(lc, reset = FALSE)
plot(st_geometry(lc_slic1), add = TRUE, border = "red")
```

---
# Segmentation quality - external data

.pull-left[
`vol_slic1`
<img src="index_files/figure-html/index-25-1.png" style="display: block; margin: auto;" />
]
.pull-right[
`vol_regions`
<img src="index_files/figure-html/index-26-1.png" style="display: block; margin: auto;" />
]

---
# Segmentation quality - external data

.pull-left[
**sabre** - **S**patial **A**ssociation **B**etween **RE**egionalizations (Nowosad and Stepinski (2018), doi:doi:10.1080/13658816.2018.1511794)

```r
library(sabre)
vm = vmeasure_calc(vol_slic1, vol_regions,
                   "supercells", "group")
vm
```

```
## The SABRE results:
## 
##  V-measure: 0.6 
##  Homogeneity: 0.49 
##  Completeness: 0.77 
## 
##  The spatial objects can be retrieved with:
##  $map1 - the first map
##  $map2 - the second map
```
]
.pull-right[
.font130[Global measures:]

1. **Homogeneity** is a measure of how well regions from the first map fit inside of regions from the second map

2. **Completeness** measures how well regions from the second map fit inside of regions from the first map

3. The final value of **V-measure** is calculated as the weighted harmonic mean of homogeneity and completeness

All of these metrics have **values between 0 and 1, where larger values indicate better spatial agreement**.
]

---
# Segmentation quality - external data

.font130[Local measures:] **regions' inhomogeneities** (rih)

.pull-left[

```r
plot(vm$map1["rih"])
```

<img src="index_files/figure-html/index-28-1.png" width="80%" style="display: block; margin: auto;" />
]

.pull-right[

```r
plot(vm$map2["rih"])
```

<img src="index_files/figure-html/index-29-1.png" width="80%" style="display: block; margin: auto;" />
]

---
# Segmentation quality

The **regional** packages calculates intra-regional and inter-regional similarities based on user-provided spatial vector objects (regions) and spatial raster objects (cells with values)

**regional** is also not, yet, on CRAN

```r
install.packages("regional", repos = "https://nowosad.r-universe.dev")
```

(Any ideas, comments, etc. are also very welcomed!)

```r
library(regional)
```

<hr>

Implemented metrics include:

1. **Inhomogeneity** - a degree of mutual dissimilarity between values of all cells in a region

2. **Isolation** - an average distance between the focus region and all of its neighbors (Haralick and Shapiro (1985), doi:10.1016/S0734-189X(85)90153-7, Jasiewicz et al. (2018), doi:10.1016/j.cageo.2018.06.003)

3. **Distinction** -  an average distance between the focus region and all of the other regions (Nowosad (2021), doi:10.1080/13658816.2021.1893324)

---
# Segmentation quality - inhomogeneity

```r
vol_slic1$inh = reg_inhomogeneity(vol_slic1, vol)
```

.pull-left[
<img src="index_files/figure-html/index-31-1.png" style="display: block; margin: auto;" />
]

.pull-right[
<img src="index_files/figure-html/index-32-1.png" style="display: block; margin: auto;" />
]

---
# Segmentation quality - isolation

```r
vol_slic1$iso = reg_isolation(vol_slic1, vol)
```

.pull-left[
<img src="index_files/figure-html/index-34-1.png" style="display: block; margin: auto;" />
]

.pull-right[
<img src="index_files/figure-html/index-35-1.png" style="display: block; margin: auto;" />
]