Umbralización de Otsu (Otsu’s thresholding)

INFORMACIÓN DE CONTACTO:

• Correo electrónico: renzoangeldelacruz@gmail.com
  
• GitHub: https://github.com/delacruz-renzo
  
• LinkedIn: https://www.linkedin.com/in/renzo-delacruz/

I. Introducción

La segmentación de imágenes constituye una etapa fundamental en la visión por computadora, ya que permite transformar información continua en clases discretas con significado temático. En el ámbito de la teledetección (remote sensing), esta tarea adquiere especial relevancia, pues facilita la delimitación de coberturas terrestres como vegetación, cuerpos de agua o áreas antrópicas a partir de imágenes satelitales. En este sentido, la elección del umbral de separación suele ser determinante para la calidad del resultado. No obstante, la selección del umbral puede introducir subjetividad cuando se basa únicamente en valores reportados en la bibliografía o en criterios visuales, ya que dichos valores no son necesariamente transferibles entre sensores, resoluciones espaciales, condiciones atmosféricas o características propias de cada escena.

En este contexto, el método de umbralización de Otsu (Otsu’s thresholding) ofrece una alternativa estadística, reproducible y basada en datos para estimar automáticamente un umbral óptimo a partir de la distribución de intensidades de la imagen. Propuesto por Nobuyuki Otsu (1979), este método se fundamenta en maximizar la separación entre dos clases mediante la optimización de la varianza interclase, lo que permite obtener una segmentación binaria robusta bajo un marco probabilístico.

El objetivo de este proyecto es desarrollar el método de Otsu aplicado a un raster NDVI derivado de imágenes satelitales CBERS-4A sobre el Humedal de Santa Rosa. Para ello, se elaboro un flujo de trabajo que incluye la preparación de los datos, la discretización del índice en niveles de intensidad, el cálculo del umbral óptimo mediante una función implementada manualmente y la comparación de resultados con un método estándar (paquete EBImage). Finalmente, se presenta la segmentación binaria obtenida y se discuten aspectos prácticos como la influencia del rango de reescalamiento, las diferencias numéricas entre implementaciones y las limitaciones del método en escenarios donde el histograma no es aproximadamente bimodal.

II. Materiales

2.1. Área de Estudio

El Humedal de Santa Rosa, también conocido como Humedal del Cascajo, se localiza en el distrito de Chancay, provincia de Huaral, departamento de Lima, en la costa central del Perú. Este ambiente palustre cumple un rol fundamental como refugio y área de alimentación para numerosas especies de aves residentes y migratorias del Pacífico sudamericano, además de albergar peces, invertebrados y flora halófila adaptada a condiciones salobres. A pesar de su alta importancia ecológica y paisajística, el humedal se encuentra sometido a fuertes presiones antrópicas derivadas de la expansión urbana, las actividades agrícolas y la disposición inadecuada de residuos, lo que resalta la necesidad de su conservación y manejo sostenible como parte del patrimonio natural costero.

library(sf)
library(sp)
library(leaflet)
library(mapview)
# Definir el directorio de trabajo del proyecto.
setwd("d:/Lenovo/Desktop/GEO/PROJECT_16")
# Importar el archivo shapefile del área de estudio.
limHumedal <- read_sf("shp/aoi_humedal.shp")
# Visualizar el área de estudio.
mapview::mapview(limHumedal, layer.name = "Área de Estudio", alpha.regions = 0, color = "blue", 
                 lwd = 2, map.types = c("Esri.WorldImagery"))

2.2. Datos Satelitales

Las imágenes satelitales utilizadas corresponden al satélite CBERS‑4A (China–Brazil Earth Resources Satellite), una misión conjunta entre Brasil y China orientada a la observación de la Tierra. En este estudio se emplearon productos del sensor WPM (Wide-Field Panchromatic and Multispectral Camera), el cual proporciona cuatro bandas multiespectrales en el rango del visible e infrarrojo cercano, con una resolución espacial de 8 m, así como una banda pancromática de mayor resolución espacial (2 m). Las imágenes analizadas corresponden a productos de Nivel 4, los cuales se encuentran geométricamente corregidos y ortorrectificados, lo que las hace adecuadas para análisis espaciales y diversas aplicaciones de teledetección.

library(raster)
library(terra)
# Obtener la lista de archivos raster multiespectrales en formato .tif.
CBERS_4A_files <- sort(list.files(path = "raster/MS", pattern = "\\.tif$", full.names = TRUE))
# Crear un raster multibanda a partir de los archivos individuales.
CBERS_4A_MS <- rast(CBERS_4A_files)
CBERS_4A_MS

## class       : SpatRaster 
## size        : 14584, 14360, 4  (nrow, ncol, nlyr)
## resolution  : 8, 8  (x, y)
## extent      : 194388, 309268, 8628390, 8745062  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18S (EPSG:32718) 
## sources     : CBERS_4A_WPM_20230319_245_128_L4_BAND1.tif  
##               CBERS_4A_WPM_20230319_245_128_L4_BAND2.tif  
##               CBERS_4A_WPM_20230319_245_128_L4_BAND3.tif  
##               CBERS_4A_WPM_20230319_245_128_L4_BAND4.tif  
## names       : CBERS_4~4_BAND1, CBERS_4~4_BAND2, CBERS_4~4_BAND3, CBERS_4~4_BAND4

# Asignar nombres a las bandas del raster multiespectral.
names(CBERS_4A_MS) <- c("CBERS4A_BAND1", "CBERS4A_BAND2", "CBERS4A_BAND3", "CBERS4A_BAND4")
# Recortar el raster multibanda al límite del área de estudio.
CBERS_4A_MS <- terra::crop(CBERS_4A_MS, limHumedal)
CBERS_4A_MS <- terra::mask(CBERS_4A_MS, limHumedal)
# Visualizar las bandas multiespectrales resultantes.
plot(CBERS_4A_MS)

# Importar la banda pancromática.
CBERS_4A_PAN <- rast("raster/PAN/CBERS_4A_WPM_20230319_245_128_L4_BAND0.tif")
CBERS_4A_PAN

## class       : SpatRaster 
## size        : 58338, 57440, 1  (nrow, ncol, nlyr)
## resolution  : 2, 2  (x, y)
## extent      : 194388, 309268, 8628386, 8745062  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18S (EPSG:32718) 
## source      : CBERS_4A_WPM_20230319_245_128_L4_BAND0.tif 
## name        : CBERS_4A_WPM_20230319_245_128_L4_BAND0

# Asignar nombre a la banda pancromática.
names(CBERS_4A_PAN) <- c("CBERS4A_BAND0")
# Recortar la banda pancromática al límite del área de estudio.
CBERS_4A_PAN <- terra::crop(CBERS_4A_PAN, limHumedal)
CBERS_4A_PAN <- terra::mask(CBERS_4A_PAN, limHumedal)
# Visualizar la banda pancromática en escala de grises.
plot(CBERS_4A_PAN, col = gray.colors(256, start = 0, end = 1))

III. Preprocesamiento y Cálculo de NDVI

3.1. Corrección Radiométrica

La información registrada por los sensores satelitales puede verse alterada por diversas interferencias durante el proceso de adquisición, lo que hace necesario aplicar procedimientos de corrección antes de realizar cualquier análisis. Entre estas interferencias se encuentran las distorsiones radiométricas, asociadas principalmente a las características y limitaciones del sensor. La calibración radiométrica es el proceso mediante el cual se corrigen estas imperfecciones, transformando los niveles digitales (Digital Numbers, DN) de los píxeles de la imagen en una magnitud física cuantificable, como la radiancia, la reflectancia o la temperatura, permitiendo así el análisis cuantitativo de la información espectral.

En este estudio, la calibración radiométrica se realizó mediante la conversión de los niveles digitales a radiancia a nivel del sensor, utilizando los coeficientes de calibración provistos en los metadatos de las imágenes CBERS-4A. La relación empleada se expresa de la siguiente manera:

\[ L_\lambda = DN \cdot C_{\text{abs}} \] Donde:

• \(L_\lambda\): Radiancia espectral registrada por el sensor para una banda espectral específica (W·m⁻²·sr⁻¹·µm⁻¹).

• \(DN\): Nivel digital (Digital Number) del píxel correspondiente.

• \(C_{\text{abs}}\): Coeficiente de calibración radiométrica absoluta, específico para cada banda espectral y obtenido de los metadatos del producto CBERS-4A.

Este procedimiento se realiza a nivel de Top of Atmosphere (TOA), es decir, corrige únicamente los efectos instrumentales asociados al sensor, sin considerar aún las influencias de la atmósfera, tales como la dispersión y absorción atmosférica, las cuales requieren métodos adicionales de corrección atmosférica. No obstante, dado que en el presente trabajo se analiza una única escena satelital y no se realiza una comparación multitemporal, se optó por aplicar exclusivamente la calibración radiométrica, la cual resulta suficiente para los objetivos del análisis desarrollado.

# Definir coeficientes de conversión a radiancia por banda.
coef_rad <- c(CBERS4A_BAND0 = 0.184471, CBERS4A_BAND1 = 0.29107, CBERS4A_BAND2 = 0.297832,
              CBERS4A_BAND3 = 0.232504, CBERS4A_BAND4 = 0.178993)
# Aplicar corrección radiométrica a bandas multiespectrales.
bandas_MS <- CBERS_4A_MS * coef_rad[names(CBERS_4A_MS)]
# Visualizar las bandas multiespectrales con valores de radiancia calibrados.
plot(bandas_MS)

# Aplicar la corrección radiométrica a la banda pancromática.
banda_PAN <- CBERS_4A_PAN * coef_rad[names(CBERS_4A_PAN)]
# Visualizar la banda pancromática con valores de radiancia calibrados.
plot(banda_PAN, col = gray.colors(256, start = 0, end = 1))

3.2. Pan-Sharpening

El pan-sharpening, también conocido como fusión pancromática, es el proceso mediante el cual se integra la información espectral de las bandas multiespectrales de menor resolución espacial con la información espacial de una banda pancromática de mayor resolución, con el objetivo de generar una imagen multiespectral de alta resolución espacial, mejorando el detalle espacial sin perder completamente su contenido espectral.

En el ejemplo presentado, se incrementa la resolución espacial de las bandas multiespectrales del satélite CBERS-4A, pasando de 8 m a 2 m, mediante la técnica Brovey, un método clásico de pan-sharpening basado en una relación aritmética entre la intensidad de la banda pancromática y la suma de las intensidades de las bandas multiespectrales. Este enfoque prioriza la mejora del contraste y la nitidez espacial, y ha sido ampliamente utilizado y evaluado en la literatura, tanto en estudios comparativos de algoritmos de fusión como en aplicaciones prácticas de imágenes satelitales de alta resolución.

# Visualizar la composición en color natural (RGB: Red, Green y Blue).
plotRGB(bandas_MS, r = 3, g = 2, b = 1, axes = TRUE, box = TRUE, stretch = "lin")

# Visualizar la composición en falso color infrarrojo (NIR, Red y Green).
plotRGB(bandas_MS, r = 4, g = 3, b = 2, axes = TRUE, box = TRUE, stretch = "lin")

# Aplicar Pan-Sharpening a la composición en falso color infrarrojo.
library(RStoolbox)
imagen_fusion <- RStoolbox::panSharpen(bandas_MS, banda_PAN, r = 4, g = 3, b = 2, method = "brovey")
imagen_fusion

## class       : SpatRaster 
## size        : 700, 700, 3  (nrow, ncol, nlyr)
## resolution  : 2, 2  (x, y)
## extent      : 252142, 253542, 8716498, 8717898  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18S (EPSG:32718) 
## source(s)   : memory
## names       : CBERS4A_BAND2_pan, CBERS4A_BAND3_pan, CBERS4A_BAND4_pan 
## min values  :          8.873853,          6.017082,          2.692099 
## max values  :         77.268342,         61.753420,         51.159028

# Visualizar la imagen fusionada con mayor resolución espacial.
plotRGB(imagen_fusion, r = 3, g = 2, b = 1, axes = TRUE, box = TRUE, stretch = "lin")

3.3. NDVI

El índice NDVI (Normalized Difference Vegetation Index, por sus siglas en inglés) es uno de los índices espectrales más utilizados en estudios de ecología y monitoreo ambiental, debido a su capacidad para estimar el vigor y la densidad de la cobertura vegetal. Su principal aplicación es la identificación y delimitación de áreas con vegetación, ya que aprovecha el marcado contraste espectral entre la alta reflectancia de la vegetación sana en el infrarrojo cercano (NIR) y la fuerte absorción de la radiación en la región del rojo del espectro electromagnético, asociada a la clorofila.

La expresión matemática del NDVI es: \[ NDVI = \frac{(NIR - Red)}{(NIR + Red)} \] Donde:

• \(NIR\): Reflectancia correspondiente a la banda del infrarrojo cercano del espectro electromagnético (CBERS-4A WPM: Banda 04).

• \(Red\): Reflectancia correspondiente a la banda del rojo del espectro visible (CBERS-4A WPM: Banda 03).

El NDVI toma valores en el intervalo \([-1, 1]\). Valores cercanos a 1 indican vegetación densa y saludable, mientras que valores próximos a 0 corresponden a superficies con escasa o nula cobertura vegetal, como suelos desnudos o áreas urbanas. Por el contrario, valores negativos suelen asociarse con cuerpos de agua, sombras o superficies artificiales, donde la reflectancia en el infrarrojo cercano es baja en comparación con la banda roja. Debido a estas características, el NDVI constituye una herramienta fundamental para el análisis espacial de la vegetación y la delimitación de áreas verdes en distintos ecosistemas.

# Extraer las bandas del infrarrojo cercano (NIR) y del rojo (Red) desde la imagen fusionada.
NIR <- imagen_fusion["CBERS4A_BAND4_pan"]
Red <- imagen_fusion["CBERS4A_BAND3_pan"]
# Calcular el Índice de Vegetación de Diferencia Normalizada (NDVI).
NDVI <- (NIR - Red)/(NIR + Red)
# Visualizar el NDVI.
plot(NDVI, main = "NDVI", col = colorRampPalette(c("red", "yellow", "green"))(255))

IV. Método de Otsu

El método de umbralización de Otsu (Otsu’s thresholding), propuesto originalmente por Nobuyuki Otsu en 1979, es una técnica estadística ampliamente utilizada en visión artificial y procesamiento de imágenes para la selección automática de un umbral óptimo que permite segmentar una imagen en dos clases. Inicialmente fue desarrollado con el objetivo de separar el primer plano del fondo; este método maximiza la separación entre las clases al optimizar la varianza interclase, a partir del análisis de la distribución de intensidades de la imagen, lo que da lugar a una segmentación estadísticamente robusta.

Este método aborda de manera explícita la problemática asociada a la selección del umbral de segmentación, al evitar la dependencia de enfoques subjetivos, como la adopción de valores propuestos en la bibliografía sin una justificación específica para el conjunto de datos analizado o para el contexto del área de estudio, así como la elección manual del umbral tras la inspección visual de múltiples resultados obtenidos con distintos valores.

El método de Otsu inicia disponiendo de una imagen monobanda (en nuestro ejemplo el raster NDVI), donde cada píxel está representado por un único valor numérico.

El objetivo es dividir la imagen en dos clases utilizando un único umbral “𝑡”:

• Clase \(C_{\text{0}}\): Píxeles con valores bajos.

• Clase \(C_{\text{1}}\): Píxeles con valores altos.

Formalmente, la asignación de clases se define como:

\[ C = \begin{cases} C_0, & \text{si } x \le t \\ C_1, & \text{si } x > t \end{cases} \]

La pregunta central del método es: ¿Cuál es el mejor valor del umbral “𝑡”?

El método no trabaja píxel por píxel, sino sobre la distribución estadística de los valores de la imagen. Por ello, el primer paso conceptual consiste en construir un histograma.

library(ggplot2)
# Convertir el raster NDVI a un data frame para su análisis estadístico.
NDVI_df <- as.data.frame(NDVI)
colnames(NDVI_df) <- "NDVI"
# Calcular métricas descriptivas del NDVI (media y mediana).
media_NDVI   <- mean(NDVI_df$NDVI, na.rm = TRUE)
mediana_NDVI <- median(NDVI_df$NDVI, na.rm = TRUE)
# Graficar el histograma de la distribución de valores NDVI.
ggplot(NDVI_df, aes(x = NDVI)) +
  geom_histogram(aes(fill = after_stat(x)), color = "white", linewidth = 0.1, bins = 30) +
  scale_fill_gradient(low = "#e5f5e0", high = "#006d2c") +
  geom_vline(xintercept = media_NDVI, color = "red", linewidth = 0.7, linetype = "dashed") +
  geom_vline(xintercept = mediana_NDVI, color = "blue", linewidth = 0.7, linetype = "dashed") +
  theme_grey(base_size = 9) +
  annotate("text", x = media_NDVI, y = 3.5e4, label = "Media", color = "red", angle = 90, vjust = -0.5, 
           size = 3.4) +
  annotate("text", x = mediana_NDVI, y = 3.5e4, label = "Mediana", color = "blue", angle = 90, vjust = -0.5, 
           size = 3.4) +
  labs(title = "Histograma de valores NDVI", x = "NDVI", y = "Frecuencia") +
  scale_x_continuous(breaks = seq(-1, 1, by = 0.2)) +
  scale_y_continuous(labels = scales::scientific) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 9.5), legend.position = "none")

El histograma se obtiene dividiendo el rango de valores de la imagen en un conjunto de intervalos o bins.

• \(L\): Número total de bins (niveles) del histograma.

Resulta importante aclarar que el parámetro \(L\) no constituye una propiedad física inherente al fenómeno analizado, sino una decisión numérica adoptada para aproximar la distribución de los valores mediante su discretización. En la presente actividad se asigna un valor de \(L = 256\) (valores entre 0 y 255), correspondiente a una imagen de 8 bits, lo cual representa un estándar ampliamente utilizado en teledetección debido a su adecuado equilibrio entre precisión radiométrica y costo computacional.

De manera alternativa, podría emplearse una correspondencia de 10 bits, equivalente a \(L = 1024\) (valores entre 0 y 1023); sin embargo, esta elección solo resulta recomendable cuando se requiere una mayor precisión radiométrica, ya que implica un incremento en el costo computacional y en el almacenamiento, por lo que debe ser debidamente justificada.

Previo a la aplicación del método, los valores de NDVI deben ser reescalados a una profundidad radiométrica de 8 bits (0–255), con el objetivo de discretizar el índice continuo y permitir la construcción del histograma correspondiente, utilizando así \(L = 256\).

Definiciones fundamentales:

• \(h(i)\): Histograma de la imágen. \[ h(i) = \text{numero de pixeles que caen en el bin } i \]
• \(N\): Número total de píxeles válidos de la imagen. \[ N = \sum_{i=0}^{L-1} h(i) \]
• \(P(i)\): Distribución de probabilidad empírica. En este contexto, \(P(i)\) representa la probabilidad de que un píxel elegido al azar de la imagen tenga intensidad \(i\).

\[ P(i) = \frac{h(i)}{N} \]

A partir de esta normalización, se cumple la propiedad fundamental de toda distribución de probabilidad discreta, según la cual la suma de las probabilidades de todos los eventos posibles es igual a 1:

\[ \sum_{i=0}^{L-1} P(i) = 1 \]

En consecuencia, el método deja de operar sobre conteos absolutos y pasa a trabajar en un marco probabilístico, interpretando los valores de intensidad de los píxeles como realizaciones de una variable aleatoria discreta. Esto permite formular el problema de segmentación en términos de probabilidades, lo que aporta ventajas teóricas y prácticas, como la invariancia del método al tamaño de la muestra y la posibilidad de definir las clases como eventos probabilísticos, es decir, calcular la probabilidad de que un píxel elegido al azar pertenezca a cada clase.

Despues se selecciona un umbral candidato “𝑡” (un bin del histograma), lo que divide la distribución en dos clases: \[ C_0 = \{0,1,\ldots,t\}, \quad C_1 = \{t+1,\ldots,L-1\} \]

• Probabilidad de la clase \(C_{\text{0}}\): \[ \omega_0(t) = \sum_{i=0}^{t} P(i) \]

• Probabilidad de la clase \(C_{\text{1}}\): \[ \omega_1(t) = \sum_{i=t+1}^{L-1} P(i) = 1 - \omega_0(t) \]

Estas probabilidades indican qué fracción de la imagen pertenece a cada clase, y actúan como factores de ponderación.

No es suficiente conocer el tamaño de cada clase; también es necesario conocer dónde se concentran sus valores. Por eso se calculan las medias (centros) de las clases.

• \(\mu_T\): Media global de la imagen. \[ \mu_T = \sum_{i=0}^{L-1} i\,P(i) \] Representa el centro de masa del histograma completo.

• \(\mu_0(t)\): Media de la clase \(C_0\). \[ \mu_0(t) = \frac{1}{\omega_0(t)} \sum_{i=0}^{t} i\,P(i) \] Es el valor promedio de los píxeles menores o iguales al umbral.

• \(\mu_1(t)\): Media de la clase \(C_1\). \[ \mu_1(t) = \frac{1}{\omega_1(t)} \sum_{i=t+1}^{L-1} i\,P(i) \] Es el valor promedio de los píxeles mayores al umbral.

El método de Otsu se fundamenta en el principio de que un umbral óptimo es aquel que permite separar la distribución de los datos en dos clases bien diferenciadas y con tamaños representativos.

Este criterio se cuantifica mediante la varianza interclase, la cual mide el grado de separación entre las medias de las dos clases generadas por un determinado umbral. Matemáticamente, la varianza interclase se expresa mediante la siguiente formulación:

\[ \sigma_b^2(t) = \omega_0(t)\,\omega_1(t)\,[\mu_0(t) - \mu_1(t)]^2 \] En donde:

• \((\mu_0 - \mu_1)^2\): Mide la separación entre las clases.

• \(\omega_0\,\omega_1\): Penaliza clases demasiado pequeñas.

El método evalúa todos los valores posibles de “𝑡” y selecciona el que maximiza la varianza inter-clase: \[ t^* = \arg \max_{0 \le t \le L} \, \sigma_b^2(t) \]

El valor \(t^*\) es el umbral óptimo según el método de Otsu.

# Definir función para cálcular el umbral óptimo según el método de Otsu.
otsu_threshold <- function(x_rast, L = 256, range = NULL, eps = 1e-10) {
  # Verificar que la entrada sea un objeto "SpatRaster" (paquete "terra").
  if (!inherits(x_rast, "SpatRaster")) {x_rast <- terra::rast(x_rast)}
  # Asegurar que el número de niveles de intensidades ("L") sea entero válido.
  if (length(L) != 1 || !is.finite(L)) stop("L debe ser un único valor numérico finito.")
  L <- as.integer(L)
  if (L < 2) stop("L debe ser un entero >= 2.")
  # Definir el rango de los valores (ignorando NA).
  if (is.null(range)) {
    r_min <- terra::global(x_rast, "min", na.rm = TRUE)[1, 1]
    r_max <- terra::global(x_rast, "max", na.rm = TRUE)[1, 1]
  } else {
    stopifnot(length(range) == 2)
    r_min <- as.numeric(range[1])
    r_max <- as.numeric(range[2])
  }
  # Validaciones del rango.
  if (!is.finite(r_min) || !is.finite(r_max)) {
    stop("El rango contiene valores no finitos (NA/Inf).")
  }
  if (r_max <= r_min) {
    stop("Rango degenerado: r_max <= r_min. Otsu no es informativo.")
  }
  # Reescalamiento lineal a niveles discretos [0, L-1] (por defecto 8 bits: [0, 255]).
  rast_nbits <- round((x_rast - r_min) / (r_max - r_min) * (L - 1))
  # Acotar los valores al rango válido ante ruido o errores numéricos.
  rast_nbits <- terra::clamp(rast_nbits, lower = 0, upper = L - 1, values = TRUE)
  # Extraer los valores del raster como vector numérico.
  vals_nbits <- terra::values(rast_nbits, mat = FALSE)
  # Eliminar valores faltantes (NA) y no validos.
  vals_nbits <- vals_nbits[is.finite(vals_nbits)]
  if (length(vals_nbits) == 0) {
    stop("No hay píxeles válidos para calcular el histograma.")
  }
  # Construcción del histograma.
  h <- tabulate(vals_nbits + 1, nbins = L) 
  # Número total de píxeles válidos.
  N <- sum(h)
  # Distribución de probabilidad empírica.
  P <- h / N
  # Niveles de intensidad discretos.
  i <- 0:(L - 1)
  # Probabilidad de la clase C0.
  omega0 <- cumsum(P)
  # Probabilidad de la clase C1.
  omega1 <- 1 - omega0
  # Media global.
  muT <- sum(i * P)
  # Media de la clase C0.
  mu0 <- cumsum(i * P) / pmax(omega0, eps)
  # Media de la clase C1.
  mu1 <- (muT - cumsum(i * P)) / pmax(omega1, eps)
  # Varianza interclase, según el criterio de maximización del método de Otsu.
  sigma_b2 <- omega0 * omega1 * (mu0 - mu1)^2
  sigma_b2[omega0 < eps | omega1 < eps] <- -Inf
  # Umbral óptimo.
  threshold_bits <- which.max(sigma_b2) - 1
  # Reconversión del umbral a la escala original.
  threshold_scaled <- r_min + (threshold_bits / (L - 1)) * (r_max - r_min)
  # Salida.
  return(list(
    threshold_scaled = threshold_scaled,
    threshold_bits   = threshold_bits,
    sigma_b2         = sigma_b2,
    levels           = i,
    range_used       = c(r_min, r_max)
  ))
}
# Calcular el umbral de Otsu para el raster NDVI.
otsu_ndvi <- otsu_threshold(NDVI)
cat("Umbral óptimo de Otsu (t*) en niveles discretos:", otsu_ndvi$threshold_bits, "\n")

## Umbral óptimo de Otsu (t*) en niveles discretos: 137

cat("Umbral óptimo de Otsu (t*) en escala original:", otsu_ndvi$threshold_scaled, "\n")

## Umbral óptimo de Otsu (t*) en escala original: 0.0632096

cat("Rango utilizado para cálcular el umbral óptimo de Otsu (t*):", otsu_ndvi$range_used, "\n")

## Rango utilizado para cálcular el umbral óptimo de Otsu (t*): -0.6167652 0.6488813

# Graficar el histograma de NDVI y el umbral estimado por Otsu.
threshold_otsu_NDVI <- otsu_ndvi$threshold_scaled
ggplot(NDVI_df, aes(x = NDVI)) +
  geom_histogram(aes(fill = after_stat(x)), color = "white", linewidth = 0.1, bins = 30) +
  scale_fill_gradient(low = "#e5f5e0", high = "#006d2c") +
  geom_vline(xintercept = media_NDVI, color = "red", linewidth = 0.7, linetype = "dashed") +
  geom_vline(xintercept = mediana_NDVI, color = "blue", linewidth = 0.7, linetype = "dashed") +
  geom_vline(xintercept = threshold_otsu_NDVI, color = "black", linewidth = 0.7, linetype = "dashed") +
  theme_grey(base_size = 9) +
  annotate("text", x = media_NDVI, y = 3.5e4, label = "Media", color = "red", angle = 90, vjust = -0.5, 
           size = 3.4) +
  annotate("text", x = mediana_NDVI, y = 3.5e4, label = "Mediana", color = "blue", angle = 90, vjust = -0.5, 
           size = 3.4) +
  annotate("text", x = threshold_otsu_NDVI, y = 3.5e4, label = "Otsu (t*)", color = "black", angle = 90, vjust = -0.5, 
           size = 3.4) +
  labs(title = "Histograma de valores NDVI", x = "NDVI", y = "Frecuencia") +
  scale_x_continuous(breaks = seq(-1, 1, by = 0.2)) +
  scale_y_continuous(labels = scales::scientific) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 9.5), legend.position = "none")

# Binarizar el NDVI con el umbral calculado mediante el método de Otsu.
mask_otsu <- NDVI > threshold_otsu_NDVI
# Aplicar un filtro espacial mayoritario (modal) para suavizar el ruido.
w <- matrix(1, 5, 5)
mask_otsu <- focal(mask_otsu, w, fun = modal, na.rm = TRUE)
# Visualizar la comparación entre el NDVI y la segmentación binaria obtenida con Otsu.
mapview::mapview(NDVI,
                 layer.name = "NDVI",
                 col.regions = c("red", "yellow", "green"),
                 map.types = "Esri.WorldImagery",
                 legend = FALSE) +
  mapview::mapview(mask_otsu,
                   layer.name = "Otsu",
                   col.regions = c("gray90", "darkgreen"),
                   legend = FALSE)

library(dplyr)
# Guardar el umbral óptimo de Otsu en niveles discretos (0–255).
threshold_bits_NDVI <- otsu_ndvi$threshold_bits
# Reescalar linealmente el NDVI a niveles discretos para análisis por intensidades.
NDVI_df <- NDVI_df %>%
  mutate(NDVI_bits = round((NDVI - min(NDVI)) / (max(NDVI) - min(NDVI)) * (256 - 1))) %>%
  # Etiquetar clases binarias según el umbral de Otsu.
  mutate(class = ifelse(NDVI_bits <= threshold_bits_NDVI, "0", "1")) %>%
  # Contar píxeles por nivel de intensidad y clase.
  count(NDVI_bits, class, name = "count")
# Reescalar la varianza interclase para superponerla sobre el histograma de conteos.
scale_factor  <- as.numeric(max(NDVI_df$count, na.rm = TRUE) / max(otsu_ndvi$sigma_b2, na.rm = TRUE))
varInterclases_plot <- data.frame(
  level = otsu_ndvi$levels[is.finite(otsu_ndvi$sigma_b2)],
  sigma = otsu_ndvi$sigma_b2[is.finite(otsu_ndvi$sigma_b2)] * scale_factor)
# Identificar el punto de corte (t*) y su σ_b^2 correspondiente en la curva de varianza interclase.
point_otsu <- varInterclases_plot %>%
  slice_min(abs(level - threshold_bits_NDVI), n = 1, with_ties = FALSE) %>%
  select(level, sigma)
# Graficar el histograma de intensidades y la curva de varianza interclase.
ggplot() +
  geom_col(data = NDVI_df, aes(x = NDVI_bits, y = count, fill = class), color = "NA", width = 1) +
  geom_line(data = varInterclases_plot, aes(x = level, y = sigma), color = "red", linewidth = 0.9, alpha = 0.7) +
  geom_vline(xintercept = threshold_bits_NDVI, color = "black", linewidth = 0.7, linetype = "dashed") +
  geom_point(data = point_otsu, aes(x = level, y = sigma), color = "red", size = 3) +
    theme_minimal(base_size = 9) +
  scale_fill_manual(values = c("0" = "gray65", "1" = "#006d2c")) +
  scale_x_continuous(breaks = seq(0, 255, by = 50)) +
  scale_y_continuous(name = "Conteo de Píxeles", labels = scales::scientific, breaks = seq(0, 1e4, by = 2.38e3),
                     sec.axis = sec_axis(~ . / scale_factor, name = "Varianza Interclase")) +
  annotate("text", x = threshold_bits_NDVI, y = 4.6e3, label = "Otsu (t*)", color = "grey10", angle = 90, 
           vjust = -0.5, size = 3.4) +
  labs(x = "Intensidad NDVI (bins 0–255)", 
       title = "Histograma de intensidades NDVI y varianza interclase") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 9.5),
        panel.grid.minor = element_blank(), legend.position = "none")

Como se observó en el resultado obtenido, el umbral óptimo según el método de Otsu fue de 0.0632, calculado mediante la función desarrollada, la cual utilizó como entradas únicamente el raster del índice NDVI y el número de bins representado por la variable “L”, con un valor de 256. Sin embargo, otra variable importante a definir es el rango de valores. En este caso, dicho input fue calculado automáticamente por la misma función a partir de los valores mínimo y máximo del raster (–0.6168, 0.6489), lo cual resulta adecuado cuando se analiza una escena individual. No obstante, cuando se trabaja con una serie temporal de escenas, se recomienda utilizar un rango teórico fijo para que los resultados sean comparables entre escenas.

Para nuestro ejemplo desarrollado, el rango teórico de los índices espectrales es [–1, 1]; a continuación, se procede a calcular el umbral utilizando este rango teórico.

# Calcular el umbral de Otsu para el raster NDVI utilizando el rango teórico [-1, 1].
otsu_ndvi2 <- otsu_threshold(NDVI, range = c(-1, 1))
cat("Umbral óptimo de Otsu (t*) utilizando rango teórico, en niveles discretos:", 
    otsu_ndvi2$threshold_bits, "\n")

## Umbral óptimo de Otsu (t*) utilizando rango teórico, en niveles discretos: 136

cat("Umbral óptimo de Otsu (t*) utilizando rango teórico, en escala original:", 
    otsu_ndvi2$threshold_scaled, "\n")

## Umbral óptimo de Otsu (t*) utilizando rango teórico, en escala original: 0.06666667

Al calcular el umbral utilizando el rango teórico, se obtiene un valor de 0.0666, lo que representa una ligera diferencia respecto al umbral obtenido al emplear los valores mínimo y máximo del raster (0.0632). Esta diferencia se debe principalmente a la forma en que se distribuyen los niveles de intensidad al reescalar los valores del NDVI.

Cuando se utiliza el rango mínimo–máximo de la imagen, los valores del raster se reescalan para ocupar todo el espectro de intensidades discretas (0–255), maximizando el uso del histograma y adaptándolo específicamente a la distribución de la escena analizada. En cambio, al emplear un rango teórico fijo, los valores del NDVI de la imagen no necesariamente cubren todo ese rango, por lo que los niveles de intensidad se concentran en un subconjunto del histograma. Esta diferencia en la distribución del histograma afecta el cálculo de las probabilidades de clase y, por ende, el valor del umbral óptimo estimado por el método de Otsu, produciendo pequeñas variaciones entre ambos enfoques.

El umbral fue calculado de forma manual mediante el desarrollo de una función propia, basada en la formulación matemática previamente descrita. No obstante, existen paquetes especializados que permiten estimar este umbral de manera directa. Entre ellos destaca el paquete EBImage, orientado al procesamiento y análisis de imágenes. Este paquete forma parte del proyecto Bioconductor, cuya misión es promover el análisis estadístico y la comprensión de ensayos biológicos de alto rendimiento actuales y emergentes (https://www.bioconductor.org/).

En este proyecto se desarrollo la función “otsu”, la cual devuelve el valor del umbral y requiere como entrada los valores de la imagen en un objeto de tipo matriz (array). Además, es necesario definir el rango de valores de la imagen, ya que por defecto la función asume un rango normalizado entre 0 y 1. Finalmente, el número de niveles de intensidad se establece por defecto en 256.

library(BiocManager)
library(EBImage)
# Convertir el raster NDVI a una matriz para su procesamiento con EBImage.
NDVI_array <- as.array(NDVI)
# Calcular el umbral óptimo de Otsu utilizando la implementación del paquete EBImage.
otsu_ndvi3 <- otsu(x = NDVI_array, range = c(-1, 1), levels = 256)
cat("Umbral óptimo de Otsu (t*) cálculado mediante el paquete EBImage:", otsu_ndvi3, "\n")

## Umbral óptimo de Otsu (t*) cálculado mediante el paquete EBImage: 0.06640625

Se observa que, al calcular el umbral mediante el paquete EBImage, se obtuvo un valor de 0.0664, mientras que el calculado con la función desarrollada manualmente fue de 0.0666, presentando una diferencia mínima en el cuarto decimal. Esto ocurre aun cuando ambos cálculos se realizaron bajo las mismas condiciones: mismos valores de NDVI, mismo número de niveles de intensidad (bins) y el mismo rango teórico.

No obstante, ambas estimaciones son conceptualmente consistentes, y la pequeña diferencia observada responde principalmente a efectos numéricos de cuantización y a detalles de implementación, como el manejo de los decimales, el criterio de redondeo o truncamiento de los valores al discretizarlos en niveles de intensidad, entre otros factores. Estas diferencias no reflejan discrepancias metodológicas del algoritmo de Otsu, sino variaciones de la implementación computacional.

4.1. Relación entre la Varianza Interclase y la Varianza Intraclase

Resulta pertinente establecer la relación entre la varianza intraclase y la varianza interclase desarrollada previamente. La varianza intraclase cuantifica el grado de dispersión de los valores dentro de cada una de las clases generadas por un umbral dado; en otras palabras, mide cuán homogéneas son las observaciones al interior de cada clase. Un valor bajo de varianza intraclase indica que los píxeles asignados a una misma clase presentan intensidades similares, lo que sugiere una partición coherente desde el punto de vista estadístico. Formalmente, la varianza intraclase se define como:

\[ \sigma_w^2(t) = \omega_0(t)\,\sigma_0^2 + \omega_1(t)\,\sigma_1^2 \]

Donde \(\omega_0(t)\) y \(\omega_1(t)\) representan las probabilidades de ocurrencia de cada clase, mientras que \(\sigma_0^2(t)\) y \(\sigma_1^2(t)\) corresponden a las varianzas condicionales asociadas a cada clase para un umbral 𝑡. Esta expresión muestra que la dispersión total dentro de las clases es una combinación ponderada de la variabilidad interna de cada grupo

En el contexto de un histograma discretizado, las varianzas condicionales por clase pueden expresarse en términos de los momentos de primer y segundo orden de la distribución de intensidades condicionada a cada clase:

\[ \sigma_k^2(t) = \mathbb{E}\!\left[i^2 \mid C_k\right] - \mu_k^2(t), \quad k \in \{0,1\} \]

Donde \(\mu_k(t)\) es la media de la clase \(C_k\), correspondiente al primer momento de la distribución condicional, y \(\mathbb{E}\!\left[i^2 \mid C_k\right]\) representa el segundo momento condicional de los niveles de intensidad \(i\) dentro de dicha clase. En la práctica computacional, estos términos se calculan de forma eficiente a partir de sumas acumuladas del primer y segundo momento del histograma discretizado.

Por otra parte, en el caso de una imagen discretizada en niveles de intensidad, la varianza total del conjunto de datos puede descomponerse como la suma de la varianza intraclase y la varianza interclase:

\[ \sigma_T^2 = \sigma_w^2(t) + \sigma_b^2(t) \] Donde la varianza total se define como:

\[ \sigma_T^2 = \mathbb{E}\!\left[i^2\right] - \mu_T^2. \] Debido a que la varianza total \(\sigma_T^2\) depende únicamente del histograma global de la imagen y no del umbral 𝑡, maximizar la varianza interclase \(\sigma_b^2(t)\) es matematicamnte matemáticamente equivalente a minimizar la varianza intraclase \(\sigma_w^2(t)\). En consecuencia, ambos criterios conducen al mismo umbral óptimo \(t^*\), lo que explica por qué el método de Otsu puede formularse indistintamente como un problema de maximización de la separación entre clases o como uno de minimización de la dispersión interna dentro de ellas.

# Definir función para cálcular el umbral óptimo según el método de Otsu.
# (Ahora usando el criterio equivalente: MINIMIZAR la varianza intraclase).
otsu_threshold_min <- function(x_rast, L = 256, range = NULL, eps = 1e-10) {
  # Verificar que la entrada sea un objeto "SpatRaster" (paquete "terra").
  if (!inherits(x_rast, "SpatRaster")) {x_rast <- terra::rast(x_rast)}
  # Asegurar que el número de niveles de intensidades ("L") sea entero válido.
  if (length(L) != 1 || !is.finite(L)) stop("L debe ser un único valor numérico finito.")
  L <- as.integer(L)
  if (L < 2) stop("L debe ser un entero >= 2.")
  # Definir el rango de los valores (ignorando NA).
  if (is.null(range)) {
    r_min <- terra::global(x_rast, "min", na.rm = TRUE)[1, 1]
    r_max <- terra::global(x_rast, "max", na.rm = TRUE)[1, 1]
  } else {
    stopifnot(length(range) == 2)
    r_min <- as.numeric(range[1])
    r_max <- as.numeric(range[2])
  }
  # Validaciones del rango.
  if (!is.finite(r_min) || !is.finite(r_max)) {
    stop("El rango contiene valores no finitos (NA/Inf).")
  }
  if (r_max <= r_min) {
    stop("Rango degenerado: r_max <= r_min. Otsu no es informativo.")
  }
  # Reescalamiento lineal a niveles discretos [0, L-1] (por defecto 8 bits: [0, 255]).
  rast_nbits <- round((x_rast - r_min) / (r_max - r_min) * (L - 1))
  # Acotar los valores al rango válido ante ruido o errores numéricos.
  rast_nbits <- terra::clamp(rast_nbits, lower = 0, upper = L - 1, values = TRUE)
  # Extraer los valores del raster como vector numérico.
  vals_nbits <- terra::values(rast_nbits, mat = FALSE)
  # Eliminar valores faltantes (NA) y no validos.
  vals_nbits <- vals_nbits[is.finite(vals_nbits)]
  if (length(vals_nbits) == 0) {
    stop("No hay píxeles válidos para calcular el histograma.")
  }
  # Construcción del histograma.
  h <- tabulate(vals_nbits + 1, nbins = L) 
  # Número total de píxeles válidos.
  N <- sum(h)
  # Distribución de probabilidad empírica.
  P <- h / N
  # Niveles de intensidad discretos.
  i <- 0:(L - 1)
  # Probabilidad de la clase C0.
  omega0 <- cumsum(P)
  # Probabilidad de la clase C1.
  omega1 <- 1 - omega0
  # Sumas acumuladas de primer y segundo momento.
  cum_ip  <- cumsum(i * P)
  cum_i2p <- cumsum((i^2) * P)
  # Medias por clase.
  mu0 <- cum_ip / pmax(omega0, eps)
  muT <- sum(i * P)
  mu1 <- (muT - cum_ip) / pmax(omega1, eps)
  # Segundo momento condicional por clase.
  Ei2_total <- sum((i^2) * P)
  E2_0 <- cum_i2p / pmax(omega0, eps)
  E2_1 <- (Ei2_total - cum_i2p) / pmax(omega1, eps)
  # Varianzas por clase.
  sigma0_2 <- E2_0 - mu0^2
  sigma1_2 <- E2_1 - mu1^2
  # Varianza intraclase.
  sigma_w2 <- omega0 * sigma0_2 + omega1 * sigma1_2
  sigma_w2[omega0 < eps | omega1 < eps] <- Inf
  # Umbral óptimo (minimiza varianza intraclase).
  threshold_bits <- which.min(sigma_w2) - 1
  # Reconversión del umbral a la escala original.
  threshold_scaled <- r_min + (threshold_bits / (L - 1)) * (r_max - r_min)
  # Salida.
  return(list(
    threshold_scaled = threshold_scaled,
    threshold_bits   = threshold_bits,
    sigma_w2         = sigma_w2,
    levels           = i,
    range_used       = c(r_min, r_max)
  ))
}

# Calcular el umbral de Otsu minimizando la varianza intraclase.
otsu_min_ndvi <- otsu_threshold_min(NDVI)
cat("Umbral óptimo por varianza intraclase:", otsu_min_ndvi$threshold_scaled, "\n")

## Umbral óptimo por varianza intraclase: 0.0632096

cat("Umbral óptimo por varianza interclase:", otsu_ndvi$threshold_scaled, "\n")

## Umbral óptimo por varianza interclase: 0.0632096

4.2. Criterio Discriminante de Fisher

El método de Otsu no introduce un criterio de optimización nuevo; su fundamento teórico se basa en la aplicación del criterio discriminante de Fisher (Fisher’s discriminant) al caso particular de una variable unidimensional, correspondiente a la discretización del histograma de intensidades de una imagen.

El criterio de Fisher cuantifica la separabilidad entre clases en relación con su dispersión interna mediante el cociente entre la varianza interclase y la varianza intraclase. Matemáticamente, el criterio se define como:

\[ J(t) = \frac{\sigma_b^2(t)}{\sigma_w^2(t)} \]

Donde:

• \(\sigma_b^2(t)\): Corresponde a la varianza interclase, la cual mide el grado de separación entre las clases generadas por el umbral 𝑡.

• \(\sigma_w^2(t)\): Corresponde a la varianza intraclase, asociada a la dispersión interna de los valores dentro de cada clase.

En un problema de clasificación o segmentación, la separación óptima se obtiene al seleccionar el umbral 𝑡 que maximiza este cociente, favoreciendo simultáneamente una alta separación entre clases y una baja dispersión interna dentro de cada una. En el método de Otsu, este criterio aparece de forma implícita y se simplifica mediante la maximización directa de la varianza interclase. En particular, el umbral óptimo se define como:

\[ t^* = \arg \max_t \; \sigma_b^2(t) \]

Esta formulación es matemáticamente equivalente al criterio discriminante de Fisher para el caso de dos clases. En el contexto del histograma de intensidades, la varianza total del conjunto de datos es constante y puede expresarse como la suma de la varianza interclase y la varianza intraclase. Dado que la varianza total no depende del umbral 𝑡, maximizar la varianza interclase es equivalente a minimizar la varianza intraclase y, en consecuencia, a maximizar el criterio de Fisher. Por esta razón, el método de Otsu no requiere utilizar explícitamente el cociente completo del criterio de Fisher, ya que todos estos criterios conducen al mismo umbral óptimo.

Si bien el criterio discriminante de Fisher está formalmente definido para problemas multiclase y espacios de características de mayor dimensión, el método de Otsu puede interpretarse como una especialización de este principio al caso unidimensional del histograma de intensidades, aplicada al problema de segmentación binaria de imágenes.

# Definir función para cálcular el umbral óptimo según el método de Otsu.
# (Ahora usando el criterio equivalente: MAXIMIZAR el criterio discriminante de Fisher).
otsu_threshold_fisher <- function(x_rast, L = 256, range = NULL, eps = 1e-10) {
  # Verificar que la entrada sea un objeto "SpatRaster" (paquete "terra").
  if (!inherits(x_rast, "SpatRaster")) {x_rast <- terra::rast(x_rast)}
  # Asegurar que el número de niveles de intensidades ("L") sea entero válido.
  if (length(L) != 1 || !is.finite(L)) stop("L debe ser un único valor numérico finito.")
  L <- as.integer(L)
  if (L < 2) stop("L debe ser un entero >= 2.")
  # Definir el rango de los valores (ignorando NA).
  if (is.null(range)) {
    r_min <- terra::global(x_rast, "min", na.rm = TRUE)[1, 1]
    r_max <- terra::global(x_rast, "max", na.rm = TRUE)[1, 1]
  } else {
    stopifnot(length(range) == 2)
    r_min <- as.numeric(range[1])
    r_max <- as.numeric(range[2])
  }
  # Validaciones del rango.
  if (!is.finite(r_min) || !is.finite(r_max)) stop("El rango contiene valores no finitos (NA/Inf).")
  if (r_max <= r_min) stop("Rango degenerado: r_max <= r_min. Otsu no es informativo.")
  # Reescalamiento lineal a niveles discretos [0, L-1].
  rast_nbits <- round((x_rast - r_min) / (r_max - r_min) * (L - 1))
  rast_nbits <- terra::clamp(rast_nbits, lower = 0, upper = L - 1, values = TRUE)
  # Extraer valores.
  vals_nbits <- terra::values(rast_nbits, mat = FALSE)
  vals_nbits <- vals_nbits[is.finite(vals_nbits)]
  if (length(vals_nbits) == 0) stop("No hay píxeles válidos para calcular el histograma.")
  # Histograma.
  h <- tabulate(vals_nbits + 1, nbins = L)
  N <- sum(h)
  P <- h / N
  # Niveles de intensidad.
  i <- 0:(L - 1)
  # Probabilidades por clase.
  omega0 <- cumsum(P)
  omega1 <- 1 - omega0
  # Sumas acumuladas de primer y segundo momento.
  cum_ip  <- cumsum(i * P)
  cum_i2p <- cumsum((i^2) * P)
  # Medias por clase.
  mu0 <- cum_ip / pmax(omega0, eps)
  muT <- sum(i * P)
  mu1 <- (muT - cum_ip) / pmax(omega1, eps)
  # Segundo momento condicional por clase.
  Ei2_total <- sum((i^2) * P)
  E2_0 <- cum_i2p / pmax(omega0, eps)
  E2_1 <- (Ei2_total - cum_i2p) / pmax(omega1, eps)
  # Varianzas por clase.
  sigma0_2 <- E2_0 - mu0^2
  sigma1_2 <- E2_1 - mu1^2
  # Varianza intraclase.
  sigma_w2 <- omega0 * sigma0_2 + omega1 * sigma1_2
  sigma_w2[omega0 < eps | omega1 < eps] <- NA_real_
  # Varianza interclase.
  sigma_b2 <- omega0 * omega1 * (mu0 - mu1)^2
  sigma_b2[omega0 < eps | omega1 < eps] <- NA_real_
  # Criterio de Fisher.
  J <- sigma_b2 / pmax(sigma_w2, eps)
  # Umbral óptimo (maximiza J).
  threshold_bits <- which.max(J) - 1
  # Reconversión del umbral a la escala original.
  threshold_scaled <- r_min + (threshold_bits / (L - 1)) * (r_max - r_min)
  # Salida.
  return(list(
    threshold_scaled = threshold_scaled,
    threshold_bits   = threshold_bits,
    J_fisher         = J,
    sigma_w2         = sigma_w2,
    sigma_b2         = sigma_b2,
    levels           = i,
    range_used       = c(r_min, r_max)
  ))
}
# Calcular el umbral de Otsu maximizando el criterio discriminante de Fisher (J).
otsu_fisher_ndvi <- otsu_threshold_fisher(NDVI)
cat("Umbral óptimo por varianza interclase:", otsu_ndvi$threshold_scaled, "\n")

## Umbral óptimo por varianza interclase: 0.0632096

cat("Umbral óptimo por varianza intraclase:", otsu_min_ndvi$threshold_scaled, "\n")

## Umbral óptimo por varianza intraclase: 0.0632096

cat("Umbral óptimo por criterio de Fisher:", otsu_fisher_ndvi$threshold_scaled, "\n")

## Umbral óptimo por criterio de Fisher: 0.0632096

# Convertir el raster NDVI a un data frame para análisis por intensidades.
NDVI_df2 <- as.data.frame(NDVI)
# Renombrar la columna que contiene los valores de NDVI.
colnames(NDVI_df2) <- "NDVI"
# Reescalar linealmente los valores de NDVI a niveles discretos [0, 255].
NDVI_df2 <- NDVI_df2 %>%
  mutate(NDVI_bits = round((NDVI - min(NDVI)) / (max(NDVI) - min(NDVI)) * (256 - 1))) %>%
  # Construir el histograma: conteo de píxeles por nivel de intensidad (bin).
  count(NDVI_bits, name = "count") %>%
  # Normalizar el histograma para compararlo con funciones objetivo normalizadas.
  mutate(count_norm = count / max(count, na.rm = TRUE))
# Normalizar las funciones objetivo (σ²_B, σ²_W y J) para graficarlas en la misma escala.
varCriterios_plot <- data.frame(
  level = otsu_fisher_ndvi$levels,
  sigma_b2 = otsu_fisher_ndvi$sigma_b2 / max(otsu_fisher_ndvi$sigma_b2, na.rm = TRUE),
  sigma_w2 = otsu_fisher_ndvi$sigma_w2 / max(otsu_fisher_ndvi$sigma_w2, na.rm = TRUE),
  J = otsu_fisher_ndvi$J / max(otsu_fisher_ndvi$J, na.rm = TRUE)) %>%
  # Conservar únicamente valores válidos.
  dplyr::filter(is.finite(sigma_b2), is.finite(sigma_w2), is.finite(J))
# Graficar el histograma de intensidades NDVI y las funciones objetivo del método de Otsu.
ggplot() +
  geom_col(data = NDVI_df2, aes(x = NDVI_bits, y = count_norm), fill = "gray75", color = NA, width = 1) +
  geom_line(data = varCriterios_plot, aes(x = level, y = sigma_b2, color = "Interclase"), linewidth = 1, 
            alpha = 0.65) +
  geom_line(data = varCriterios_plot, aes(x = level, y = sigma_w2, color = "Intraclase"), linewidth = 1, 
            alpha = 0.65) +
  geom_line(data = varCriterios_plot, aes(x = level, y = J, color = "Fisher"), linewidth = 1, 
            alpha = 0.65) +
  scale_color_manual(values = c("Interclase" = "red", "Intraclase" = "blue", "Fisher" = "darkgreen"),
                     labels = c("Interclase" = "Varianza\ninterclase (σ²_B)",
                                "Intraclase" = "Varianza\nintraclase (σ²_W)",
                                "Fisher"     = "Discriminante\nde Fisher (J)")) +
  scale_x_continuous(breaks = seq(0, 255, by = 50)) +
  geom_vline(xintercept = threshold_bits_NDVI, color = "black", linewidth = 0.7, linetype = "dashed") +
  geom_point(data = varCriterios_plot %>% 
               dplyr::filter(level == threshold_bits_NDVI), aes(x = level, y = sigma_b2 + 0.004),
             color = "red", size = 3) +
  geom_point(data = varCriterios_plot %>% 
               dplyr::filter(level == threshold_bits_NDVI), aes(x = level, y = sigma_w2),
             color = "blue", size = 3) +
  geom_point(data = varCriterios_plot %>% 
               dplyr::filter(level == threshold_bits_NDVI), aes(x = level, y = J - 0.004),
             color = "darkgreen", size = 3) +
  theme_minimal(base_size = 9) +
  annotate("text", x = threshold_bits_NDVI, y = 0.5, label = "Otsu (t*)", color = "grey10", angle = 90, 
           vjust = -0.5, size = 3.4) +
  labs(x = "Intensidad NDVI (bins 0–255)", y = "Normalización",
       title = "Histograma de intensidades NDVI y funciones objetivo del método de Otsu",
       color = "Función objetivo") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 9.5),
        legend.title = element_text(face = "bold"),
        panel.grid.minor = element_blank())

4.3. Limitaciones y Variaciones del Método Otsu

En este proyecto se implementó el método de Otsu clásico, también conocido como Otsu binivel, el cual presenta un desempeño óptimo cuando la distribución de intensidades de la imagen es aproximadamente bimodal, es decir, cuando existen dos grupos bien diferenciados de valores que corresponden a dos clases espectralmente separables (por ejemplo, vegetación y no vegetación en el caso del NDVI). En este contexto, el umbral óptimo tiende a ubicarse en el valle entre ambos modos del histograma, maximizando la separación entre clases.

Sin embargo, en el caso de imágenes satelitales y productos derivados de sensores remotos, donde es común encontrar alta heterogeneidad espectral y espacial, el histograma suele presentar distribuciones multimodales o un solapamiento significativo entre clases; en estos escenarios, el método puede producir umbrales subóptimos desde el punto de vista semántico, aun cuando el criterio estadístico sea formalmente correcto. Por esta razón, la inspección del histograma y del contexto espectral de la escena constituye una etapa fundamental para evaluar la idoneidad del método en cada aplicación.

Debido a estas limitaciones, diversos autores han propuesto extensiones al método original, ya sea mediante ajustes en la formulación matemática del criterio de Otsu o a través de la incorporación de etapas de procesamiento adicionales.

Una extensión directa del método clásico es el Otsu multinivel (multi-threshold Otsu), el cual generaliza el criterio de maximización de la varianza interclase para estimar múltiples umbrales de forma simultánea. Esto permite segmentar la imagen en tres o más clases, abordando la limitación binaria del método clásico cuando el histograma presenta múltiples modos asociados a distintos rangos de intensidad. No obstante, su principal limitación es el costo computacional, ya que el tiempo de procesamiento y el consumo de memoria crecen de forma exponencial con el número de umbrales a determinar, especialmente en imágenes de gran tamaño o con alta profundidad radiométrica.

Otra variante es el Otsu local o adaptativo (también denominado por ventanas, bloques o tiles), en la cual la imagen se divide en regiones y el umbral se calcula de manera independiente para cada una de ellas. Este enfoque permite adaptarse a variaciones locales de iluminación, contraste o distribución de intensidades que no pueden ser capturadas mediante un único umbral global, constituyendo una de las principales alternativas para escenas con alta heterogeneidad de coberturas. Sin embargo, su principal desventaja es la aparición de discontinuidades entre regiones adyacentes, ya que cada bloque se segmenta de forma independiente. Este efecto puede mitigarse mediante el uso de ventanas solapadas, interpolación entre umbrales vecinos o procesos de suavizado espacial posteriores.

Adicionalmente, otra variante relevante es el Otsu ponderado (Weighted Otsu) modifica la formulación original de la varianza interclase introduciendo factores de ponderación en las probabilidades de las clases. El objetivo de esta variante es controlar el sesgo del método clásico hacia las clases dominantes del histograma, favoreciendo la separación de clases minoritarias que podrían quedar absorbidas por una clase mayoritaria. Este enfoque resulta útil cuando se desea priorizar la detección de estructuras u objetos que ocupan una proporción reducida del total de píxeles, ya que el criterio de optimización deja de depender únicamente de la distribución global de frecuencias.

Es común complementar el método de Otsu con etapas de preprocesamiento y posprocesamiento. Por ejemplo, el uso de filtros de suavizado (como el filtro gaussiano) permiten estabilizar el histograma y reducir el ruido de la imagen antes del cálculo del umbral; este enfoque suele denominarse Robust Otsu thresholding o Preprocessed Otsu thresholding. Asimismo, después de la segmentación es habitual aplicar operaciones espaciales o morfológicas para mejorar la coherencia espacial del output, eliminando regiones pequeñas aisladas y regularizando bordes. Estas operaciones no modifican el criterio de Otsu, pero sí mejoran la calidad geométrica de las regiones segmentadas y la interpretabilidad del resultado final.

Finalmente, es importante mencionar que estas variantes del método de Otsu requieren un desarrollo metodológico más detallado y evaluaciones comparativas específicas para analizar su comportamiento y desempeño, por lo que se abordarán como una línea de trabajo independiente en proyectos posteriores dentro del portafolio.

V. Conclusiones

En este proyecto se implementó el método de umbralización de Otsu aplicado a un raster NDVI derivado de imágenes CBERS-4A, consolidando un flujo de trabajo reproducible desde el preprocesamiento hasta la obtención de una segmentación binaria del área de estudio. El umbral óptimo estimado mediante la función desarrollada manualmente fue consistente con el obtenido usando una implementación estándar (EBImage), observándose diferencias mínimas (en el cuarto decimal), lo que confirma que ambas soluciones responden al mismo principio estadístico y que las discrepancias se deben principalmente a efectos numéricos propios de la implementación computacional.

Asimismo, se evidenció que la selección del rango utilizado para discretizar el NDVI (mínimo–máximo de la escena vs rango teórico [−1,1]) influye directamente en la forma del histograma y, por ende, en el valor del umbral resultante. Para el análisis de una sola escena, el rango mínimo–máximo permite aprovechar la totalidad del histograma; no obstante, en análisis multitemporales se recomienda emplear un rango teórico fijo para garantizar la comparabilidad entre escenas. De manera complementaria, la validación del método desde tres perspectivas equivalentes —maximización de la varianza interclase, minimización de la varianza intraclase y maximización del criterio discriminante de Fisher— refuerza el sustento teórico del enfoque y justifica la formulación adoptada por Otsu.

Si bien el método de Otsu proporciona un umbral objetivo basado en la distribución estadística de los valores, su aplicación directa sobre índices espectrales continuos como el NDVI debe interpretarse con cautela, ya que la separación estadística entre clases no garantiza una correspondencia directa con clases semánticas bien definidas en el terreno. Diferentes tipos de cobertura vegetal pueden presentar rangos de NDVI similares y quedar agrupados en una misma clase, mientras que superficies no vegetadas con firmas espectrales intermedias pueden solaparse parcialmente con la clase de vegetación, generando ambigüedades en la interpretación temática del resultado.

Finalmente, la aplicación de un filtro espacial mayoritario permitió mejorar la coherencia geométrica de la máscara binaria, reduciendo píxeles aislados y facilitando su interpretación cartográfica. No obstante, el desempeño del método es óptimo cuando el histograma es aproximadamente bimodal; en escenas con alta heterogeneidad espectral o solapamiento entre clases, el umbral obtenido puede ser estadísticamente correcto pero conceptualmente limitado desde el punto de vista temático. En contextos operativos, este enfoque se complementa con variantes adaptativas o multinivel, así como con métodos supervisados o modelos de aprendizaje automático para una discriminación temática más precisa.

Referencias

• Alkhimova, S., & Sliusar, S. (2019). Analysis of effectiveness of thresholding in perfusion ROI detection on T2-weighted MR images. arXiv. Ver recurso 🔗

• Navarrete, P. H. (2016). Umbralización múltiple utilizando el método de Otsu para reconocer la luz roja en semáforos. Scientia, 17(17), 247–262. Ver recurso 🔗

• Niño-Rondón, C. V., Castro-Casadiego, S. A., Medina-Delgado, B., Guevara-Ibarra, D., & Camargo-Ariza, L. L. (2021). Comparativa entre la técnica de umbralización binaria y el método de Otsu para la detección de personas. Revista UIS Ingenierías, 20(2), 65–74. Ver recurso 🔗

• Otsu, N. (1979). A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man, and Cybernetics, 9(1), 62–66. Ver recurso 🔗

• Shaikh, A., Botcha, S., & Krishna, M. (2022). Otsu based differential evolution method for image segmentation. arXiv. Ver recurso 🔗

• Shapiro, L. G., & Stockman, G. C. (2001). Computer vision. Prentice Hall. Ver recurso 🔗

• Triana, N., Jaramillo, A. E., Gutiérrez, R. M., & Rodríguez, C. A. (2016). Técnicas de umbralización para el procesamiento digital de imágenes de GEM-Foils. Scientific and Technological Advances, 21(4), 352–361. Ver recurso 🔗