1 INTRODUCCION

Este apartado muestra los fundamentos de una técnica de clasificación de imágenes, denominada segmentación, a partir de una fotografía en color natural RGB, es decir, una imagen con tres bandas que representan el espectro visible. En una imagen en color natural, la radiación en las longitudes de onda verdes se muestra en verde, la radiación en las longitudes de onda rojas se muestra en rojo y la radiación en las longitudes de onda azules se muestra en azul.

En el procesamiento de imágenes digitales, la segmentación es el proceso de dividir una imagen digital en múltiples segmentos (conjuntos de píxeles, también conocidos como objetos de imagen o superpíxeles). El objetivo de la segmentación es simplificar y/o cambiar la representación de una imagen en algo que sea más significativo y más fácil de analizar. Más precisamente, la segmentación de imágenes es el proceso de asignar una etiqueta a cada píxel de una imagen de modo que los píxeles con la misma etiqueta compartan ciertas. características. La segmentación de imágenes se utiliza normalmente para localizar objetos y límites (líneas, curvas, etc.) en imágenes.

Los contenidos provienen en gran medida de las páginas web https://rpubs.com/ials2un/segment_rgb y de https://www.r-bloggers.com/2020/03/analyzing-remote-sensing-data-using-image-segmentation/

Los paquetes necesarios para su ejecución son los siguientes:

#install.packages("rasterVis")
#install.packages('OpenImageR')
#install.packages('RImagePalette')
#library("rasterVis")
#library('OpenImageR')
#library('RImagePalette')

La imagen utilizada proviene del sitio web NASA Earth Observatory (https://earthobservatory.nasa.gov/images/49352/agricultural-fields-near-perdizes-minas-gerais-brazil). Muestra el paisaje agrícola en la parte occidental del estado de Minas Gerais en Brasil, al SW de la ciudad de Perdizes, consistente en una combinación de campos poligonales cuadriculados regularmente y campos circulares, acompañados de pequeños riachuelos del río Araguarí. Los cultivos son muy variados: girasoles, trigo, patatas, café, arroz, soja y maíz. Los campos en barbecho, que no tienen un uso agrícola activo, muestran los suelos violetas, rojizos y tostados claros, comunes en esta parte de Brasil. Los suelos más oscuros suelen ser ricos en óxidos de hierro y aluminio, y son típicos de suelos muy erosionados que se forman en climas cálidos y húmedos.

cultivos <- stack("D:/G174_Teledeteccion_2022/OCW/Lab_7_Clasificacion_imagenes/figuras/segmentation/cultivos.jpg")
plotRGB(cultivos, 1, 2, 3)

2 ÍNDICES DE VEGETACIÓN OBTENIDOS SOBRE IMÁGENES EN COLOR VERDADERO

Un Índice de Vegetación (VI) es una transformación espectral de dos o más bandas diseñada para mejorar la contribución de las propiedades de la vegetación y permitir la comparación entre áreas y durante diferentes momentos del tiempo, mostrando variaciones en la actividad fotosintética y variaciones estructurales del dosel.

Hay muchos índices de vegetación (IV), y muchos de ellos son funcionalmente equivalentes. Muchos de los índices utilizan la relación inversa entre la reflectancia del rojo y del infrarrojo cercano asociada con una vegetación verde saludable.

También existen algunos índices que combinan las bandas espectrales del visible. Jian et al. (2018) propusieron TCVI, un nuevo índice de vegetación de color natural RGB, y compararon su rendimiento con otros índices basados en RGB.

library(knitr)
include_graphics("D:/G174_Teledeteccion_2022/OCW/Lab_7_Clasificacion_imagenes/figuras/segmentation/VIs.png")

A partir de la tabla anterior, probaremos aquí dos índices: el ExG y el NGRDI.

Para calcular el índice ExG se puede utilizar la siguiente función

ExG <- function(img, i, j, k) {
  bi <- img[[i]]
  bj <- img[[j]]
  bk <- img[[k]]
  vi <- (2*bj-bi-bk) / (bj + bk +bi)
  return(vi)
}

Para calcular el índice NGRDI:

NGRDI <- function(img, i, j, k) {
  bi <- img[[i]]
  bj <- img[[j]]
  bk <- img[[k]]
  vi <- (bj-bi) / (bj + bi)
  return(vi)
}

Gracias a Matthias Forkel podemos utilizar una bonita paleta para índices de vegetación:

colores <- function(n) {
    .fun <- colorRampPalette(c("chocolate4", "orange", "yellow", "grey", "green", "green3", "darkgreen"))
    col <- .fun(n)
    return(col)
}

Cálculo del índice ExG: en el caso de la imagen de los cultivos la banda 1 = red, la banda 2 = green y la banda 3 = azul.

egvi <- ExG(cultivos, 1, 2,3)
cortes <- seq(-0.35, 1, by=0.1)
cols <- colores(length(cortes)-1)
##plot(ndvitrend, 2, col=cols, breaks=classbreaks)
plot(egvi, 2, col = cols, breaks=cortes, 
     main = "Índice de vegetación ExG")

Índice EXG

En el caso del índice NGR

ngvi <- NGRDI(cultivos, 1, 2,3)

Visualización

cortes <- seq(-0.35, 1.0, by=0.1)
cols <- colores(length(cortes)-1)
plot(ngvi, 2, col = cols, breaks=cortes, main = "NGRDI Vegetation Index")

Índice NGRDi

3 SEGMENTACIÓN DE IMÁGENES

Hay varios paquetes de análisis de imágenes en R para realizar la segmentación. Los objetos de imagen creados por el paquete OpenImageR, desarrollado por Mouselimis (2018), tienen una estructura de datos muy simple y transparente que los hace ideales para el análisis de datos y, en particular, para la integración con el paquete raster.

A su vez, el paquete OpenImageR utiliza agrupación iterativa lineal simple. El algoritmo SLIC mide la distancia total entre píxeles como la suma de dos componentes:

Específicamente, la distancia Ds se calcula como

\[D_s = d_lab + \frac{m}{S} * d_{xy}\]

donde \(S\) es el intervalo de cuadrícula de los píxeles y \(m\) es un parámetro de compacidad que permite enfatizar o restar énfasis a la compacidad espacial de los superpíxeles. El algoritmo comienza con un conjunto de K centros de conglomerados dispuestos regularmente en la cuadrícula (x,y) y realiza una agrupación de K-medias de estos centros hasta que se alcanza una medida de convergencia predefinida.

Para cargar la imagen se utiliza la función readImage() del paquete OpenImageR:

RGB.Color <- readImage("D:/G174_Teledeteccion_2022/OCW/Lab_7_Clasificacion_imagenes/figuras/segmentation/cultivos.jpg")

Usemos ahora SLIC para la segmentación de imágenes. Para obtener una explicación detallada del código, consulte la documentación correspondiente del citado paquete.

Region.slic <- superpixels(input_image = RGB.Color, 
                           method = "slic", 
                           superpixel = 80,
                           compactness = 30, 
                           return_slic_data = TRUE,
                           return_labels = TRUE, 
                           write_slic = "",
                           verbose = FALSE)

imageShow(Region.slic$slic_data)

Índice NGVI

La función imageShow() es parte del paquete OpenImageR. La función superpixels() crea superpíxeles, pero en realidad no crea una segmentación de imagen en el sentido que definimos al principio de esta sección. Esto se hace mediante funciones en el paquete SuperPixelImageSegmentation cuya discusión está más allá del alcance de este apartado. Para obtener más información, revise la documentación del paquete.

4 ANÁLISIS DE DATOS

La estructura simple de los objetos generados por la función superpixels() hace que sea fácil utilizar estos objetos para nuestros propósitos. En la última sección leemos los datos de la imagen en el objeto RGB.Color para su análisis. Echemos un vistazo a la estructura de este objeto.

str(RGB.Color)

El resultado es una matriz tridimensional. El “alto” y el “ancho” del cuadro son el número de filas y columnas respectivamente en la imagen y en cada valor de alto y ancho la “profundidad” es un vector cuyos tres componentes son el color rojo, verde y azul. valores de esa celda en la imagen, escalados a [0,1] dividiendo los valores RGB de 8 bits, que van de 0 a 255, por 255. Esta es la fuente del mensaje de advertencia que acompaña a la salida de la función superpíxels(), que decía que los valores se multiplicarán por 255. Sin embargo, para nuestros propósitos, significa que podemos analizar fácilmente imágenes que consisten en transformaciones de estos valores, como el NGRDI.

Consigamos una matriz con valores NGRDI:

ngvi.mat <- matrix(ngvi@data@values,
            nrow = ngvi@nrows,
            ncol = ngvi@ncols, byrow = TRUE)

La función imageShow() trabaja con datos que están en el rango de ocho bits 0 – 255 o en el rango [0,1] (es decir, el rango de x entre 0 y 1 inclusive). Sin embargo, no funciona con valores NGRVI si estos valores son negativos. Por lo tanto, escalaremos los valores de NGRDI a [0,1].

normalize <- function(x) {
  min <- min(x)
  max <- max(x)
  return(1* (x - min) / (max - min))
}
ngvi.mat1 <- normalize(ngvi.mat)
imageShow(ngvi.mat1)

ngvi.data <- RGB.Color
#ngvi.data[,,1] <- ngvi.mat1
#ngvi.data[,,2] <- ngvi.mat1
#ngvi.data[,,3] <- ngvi.mat1

En la siguiente sección describiremos cómo crear una segmentación de imágenes de los datos del NGRDI y cómo utilizar el análisis de conglomerados para crear una clasificación de usos de suelo.

5 Clasificación de la imagen

Aquí se muestra una aplicación de la función superpixels() para la clasificación de usos de suelo con los valores de NGRDI.

ngvi.80 <- superpixels(input_image = ngvi.data,
                      method = "slic", 
                      superpixel = 80,
                      compactness = 30, 
                      return_slic_data = TRUE,
                      return_labels = TRUE, 
                      write_slic = "",
                      verbose = FALSE)

imageShow(ngvi.80$slic_data)

Obsérvese la estructura del objeto ngvi.80.

str(ngvi.80)

Es una lista con dos elementos:

Veamos si esto es cierto.

sort(unique(as.vector(ngvi.80$labels)))

Hay 75 etiquetas únicas. Aunque al llamar a superpixels() se especificaron 80 superpíxeles, la función sólo generó 76. Podemos ver qué píxeles tienen el valor de etiqueta 0 estableciendo los valores de todos los demás píxeles en [255, 255, 255], que los trazará como blanco. .

R0 <- ngvi.80
for (i in 1:nrow(R0$label))
    for (j in 1:ncol(R0$label))
       if (R0$label[i,j] != 0)
          R0$slic_data[i,j,] <- c(255,255,255)
imageShow(R0$slic_data)

La operación aísla el superpíxel en la esquina superior derecha de la imagen, junto con la porción correspondiente del límite general. Podemos utilizar fácilmente este enfoque para determinar qué valor de ngvi.80$label corresponde a qué superpíxel.

Ahora ocupémonos del límite. Una pequeña exploración del objeto ngvi.80 sugiere que los píxeles en el límite tienen los tres componentes iguales a cero. Aislamos y trazamos todos esos píxeles coloreando todos los demás píxeles de blanco.

Bdry <- ngvi.80
  for (i in 1:nrow(Bdry$label))
   for (j in 1:ncol(Bdry$label))
     if (!(Bdry$slic_data[i,j,1] == 0 &
     Bdry$slic_data[i,j,2] == 0 &
       Bdry$slic_data[i,j,3] == 0))
       Bdry$slic_data[i,j,] <- c(255,255,255)
Bdry.norm <- NormalizeObject(Bdry$slic_data)
imageShow(Bdry$slic_data)

La figura anterior muestra que efectivamente hemos identificado los píxeles límite. Tenga en cuenta que la función imageShow() muestra estos píxeles en blanco con un borde negro, en lugar de negro puro.

Habiendo realizado un análisis preliminar, se puede organizar el proceso de segmentación en dos pasos.

Se utilizará la función make.segments() para realizar la segmentación. El primer argumento es el objeto de superpíxeles y el segundo es la medida funcional de tendencia central. Aunque en este caso cada uno de los tres colores del objeto ngvi.80 tiene los mismos valores, esto puede no ser cierto para todas las aplicaciones, por lo que la función analiza cada color por separado.

La función, escrita por Richard Plant, es la siguiente:

5.0.1 PASO 1. Identificar una medida de tendencia central de cada superpixel.

make.segments <- function(x, ftn){
# El argumento ftn es cualquier medida funcional de tendencia central
   z <- x
# Para cada superpixel identificado, calcula una medida de tendencia central
   for (k in unique(as.vector(x$labels))){
# Identificar los miembror del superpixel que tienen una etiqueta 
      in.super <- matrix(0, nrow(x$label), ncol(x$label))
      for (i in 1:nrow(x$label))
         for (j in 1:ncol(x$label))
            if (x$label[i,j] == k)
               in.super[i,j] <- 1
# Identificar las celdas situadas en el los límites como aquellas que tienen todos los valores como 0
      on.bound <- matrix(0, nrow(x$label), ncol(x$label))
      for (i in 1:nrow(x$label))
         for (j in 1:ncol(x$label))
            if (in.super[i,j] == 1){
               if (x$slic_data[i,j,1] == 0 & x$slic_data[i,j,2] == 0 
                  & x$slic_data[i,j,3] == 0)
                     on.bound[i,j] <- 1
         }
# Identificar los píxeles del superpixel situados fuera del límite
      sup.data <- matrix(0, nrow(x$label), ncol(x$label))
         for (i in 1:nrow(x$label))
         for (j in 1:ncol(x$label))
            if (in.super[i,j] == 1 & on.bound[i,j] == 0)
               sup.data[i,j] <- 1
# Calcular la medida de tendencia central de los píxeles en R, G, B
      for (n in 1:3){
# Crea una matriz M del mismo tamaño que la matriz de valores superpixel
         M <- matrix(0, dim(x$slic_data)[1], dim(x$slic_data)[2]) 
         for (i in 1:nrow(x$label))
            for (j in 1:ncol(x$label))
# Asigna los valores M en el superpixel.
               if (sup.data[i,j] == 1) M[i,j] <- x$slic_data[i,j,n]
          if (length(M[which(M > 0 & M < 255)]) > 0)
# Calcula la medida de tendencia central
            ftn.n <- round(ftn(M[which(M > 0 & M < 255)]), 0)
         else
            ftn.n <- 0
         for (i in 1:nrow(x$label))
            for (j in 1:ncol(x$label))
               if (in.super[i,j] == 1) z$slic_data[i,j,n] <- ftn.n
           }
      }
   return(z)
   }

Aquí está la aplicación de esa función al objeto ngvi.80 con el segundo argumento configurado como el promedio. Podemos ver el resultado.

ngvi.means <- make.segments(ngvi.80, mean)
imageShow(ngvi.means$slic_data)

5.0.2 PASO 2. El siguiente paso es desarrollar grupos que representen tipos de cobertura terrestre identificables.

En un proyecto real, el procedimiento sería recopilar un conjunto de datos reales del sitio, pero esa opción no está disponible para nosotros. En su lugar, trabajaremos con la reproducción en color real de la foto del astronauta.

La cobertura del suelo se puede subdividir utilizando K-medias en cualquier número de tipos. Tenga en cuenta que no he interpretado aquí la relación entre los grupos de K-medias y las clases de cobertura del suelo existentes en Perdizes. Sin embargo, en una aplicación real, esta tarea es imprescindible.

set.seed(123)
ngvi.clus <-
 kmeans(as.vector(ngvi.means$slic_data[,,1]), 10)
vege.class <- matrix(ngvi.clus$cluster,
    nrow = ngvi@nrows,
   ncol = ngvi@ncols, byrow = FALSE)
class.ras <- raster(vege.class, xmn = 0,
    xmx = 1440, ymn = 0, ymx = 960, crs = crs(cultivos))

A continuación podemos usar la función ráster ratify() para asignar niveles de factores descriptivos a los grupos.

class.ras <- ratify(class.ras)
rat.class <- levels(class.ras)[[1]]
rat.class$landcover <- c("LC 10", "LC 09", "LC 08", "LC 07", "LC 06", "LC 05", "LC 04", "LC 03", "LC 02", "LC 01")

Crea una paleta de 10 colores.

ncultivos <-  jpeg::readJPEG("./cultivos.jpg")
cultivos.pal <- image_palette(ncultivos, n=10)
levels(class.ras) <- rat.class
levelplot(class.ras, margin=FALSE,
   col.regions= cultivos.pal[c(5,7,9,4,8,10,2,4,1,6)],
    main = "Tipos de usos de suelo en Perdizes")

Eche un vistazo a los resultados y evalúe si son una representación fiel de las clases de cobertura terrestre existentes o no.