Redes Markovianas en Procesamiento de Imágenes

Una aplicación importante de redes markovianas es en procesamiento de imágenes, en esta diciplina las redes markovianas típicamente se conocen como Campos Aleatorios de Markov (Markov Random Fields) y se utilizan con distintos objetivos, por ejemplo, en la eliminación de ruido (eliminar lo borroso de una imagen), segmentación de imagenes o reconocimiento de objetos.

En la mayor parte de estas aplicaciones el modelo sigue una estructura de red Markoviana a pares donde los nodos corresponden a los pixeles y las aristas a interacciones entre pixeles adyacentes, por lo tanto cada pixel interior tiene exactamente cuatro vecinos. La definición de los potenciales depende de la aplicación, sin embargo es usual que los modelos se especifiquen en términos de energías (log-potenciales negativos) de tal manera que los valores representan penalizaciones y un valor menor corresponde a una configuración de mayor probabilidad.

En el caso de elimincación de ruido el objetivo es recuperar el verdadero valor de los pixeles dado que algunos pixeles presentan ruido. Denotemos por \(X_i\) el verdadero valor del pixel \(i\) (este es el valor que buscamos inferir) y por \(Y_i\) el valor del pixel que observamos (que presenta ruido). En este problema suponemos una correlación fuerte entre \(X_i\) y \(Y_i\), además de haber correlación positiva entre pixeles vecinos \(X_i\) y \(X_j\). La gráfica correspondiente a este modelo se muestra a continución.

La gráfica muestra dos tipos de cliques, cada uno con dos variables. Los cliques de la forma \(\{X_i, Y_i\}\) a la que asociaremos una función de energía que cumple con el efecto de asignar menor energía (y por tanto alentando mayor probabilidad) cuando el valor de las variables es cercano. Por otra parte, los cliques de la forma \(\{X_i,X_j\}\) donde \(i\) y \(j\) son índices de pixeles vecinos, tienen asociada una función de energía que penaliza diferencias entre pixeles vecinos.

Veamos un ejemplo sencillo donde la imagen ruidosa esta descrita por una matriz de pixeles que toman uno de dos valores \(y_i \in \{-1,1\}\), el índice \(i\) indica el pixel. Supongamos que la imagen ruidosa se obtiene de una imagen libre de ruido que se descibe por una matriz donde los pixeles toman los valores \(x_i \in \{-1,1\}\), el ruido se añade cambiando el signo de los pixeles de la imagen original con una probabilidad baja. La figura de abajo es un ejemplo donde se cambio el signo de los pixeles con probabilidad 0.10.

Capturamos la relación entre \(X_i\) y \(Y_i\) asociando la función de energía \(-\eta x_i y_i\) donde \(\eta\) es una constante positiva, observemos que la energía es menor cuando \(x_i\) y \(y_i\) tienen el mismo signo y mayor cuando tienen el signo opuesto, ¿Cómo se traduce esto en los potenciales? En el caso de cliques que relacionan pixeles vecinos \(X_i\), \(X_j\) elegimos la función de energía \(-\beta x_i x_j\) donde \(\beta\) es una constante positiva, el comportamiento de esta función es análogo a la función que relaciona \(X_i\) y \(Y_i\).

Finalmente, añadimos el término \(hx_i\) para cada pixel \(i\) cuya función es sesgar el modelo hacia valores de pixeles que compartan el mismo signo. La función de energía completa para el modelo toma la forma: \[E(x,y) = h\sum_i x_i - \beta \sum_i x_i x_j - \eta \sum_i x_i y_i\] que determina la distribución conjunta: \[p(x,y)=\frac{1}{Z}exp(-E(x,y)).\] Una vez definidas las funciones de energía se asignan a \(y\) los valores de los pixeles observados en la imagen ruidosa que define implícitamente una distribución condicional \(p(x|y)\) sobre los pixeles libres de ruido. Buscamos entonces encontrar una imagen \(x\) de probabilidad alta (idealmente máxima) condicional a la imagen observada.
Para recuperar la imagen utilizamos una técnica iterativa llamada Modas Condicionales Iteradas (o ICM por sus siglas en inglés) que es una aplicaciónn de ascenso de gradiente:

  1. Se inicializan las variables \(\{x_i\}\) igualando \(x_i=y_i\) para toda \(i\).
  2. Para cada nodo \(X_j\) evaluamos la energía total para cada estado \(x_j = 1\) y \(x_j = -1\).
  3. Actualizamos el valor que toma \(X_j\) al estado que minimiza la energía. Repetimos el procedimiento en un nuevo nodo hasta que se satisfaga algún criterio de paro. Este paso puede mantener la probabilidad de constante o incrementarla si cambiamos el valor que toma \(X_j\).
  4. Repetimos el paso 2 y 3 hasta satisfacer algún criterio de paro.

La función de energía es en principio arbitraria, sin embargo, existen varias alternativas estándar para procesamiento de imagenes, por ejemplo, si la imagen se encuentra en escala de grises, los pixeles pueden tomar cualquiera de 256 valores (8 bits) y una función de energía apropiada es la siguiemte:

\[E(x,y) = \sum_i \frac{(x_i - y_i)^2}{2\sigma^2} + \gamma \sum_{\{i,j\}} min((x_i - x_j)^2, \beta)\]

donde \(\sigma^2\) representa la creencia de que la imagen está corrompida por ruido con una varianza \(\sigma^2\).

library(bmp)
library(pixmap)
library(jpeg)

perro <- readJPEG("imagenes/perrito.jpg", native = FALSE)[,,1]
pr <- pixmapGrey(perro)
plot(pr)

perro_mat <- 255 * perro

# Añadir ruido
noise <- function(image, sdev = 20){
  dimension <- dim(image)
  a <- round(c(image) + rnorm(n = length(image), 0, sd = sdev))
  a[a < 0] <- 0
  a[a > 255] <- 255
  image_blur <- array(data = a, dim = dimension)
}
perro_n <- noise(perro_mat) 
pr_n <- pixmapGrey(perro_n / 255)
plot(pr_n)


icm <- function(image, sdev = 25, max.diff = 200, weight.diff = 0.15, 
                iterations = 20){
  
  # sdev: error estándar del ruido gaussiano.
  # max_diff: contribución máxima al potencial de la diferencia entre los valores
  #   de dos pixeles vecinos 
  # weight_diff: es la ponderación asociada al componente del potencial 
  #   debida a la diferencia entre los valores de pixeles vecinos.
  # iter: es el número de iteraciones.

  # Siempre tengo dos imágenes, en cada iteración se alternan entre imagen 
  # fuente e imagen destino.
  
  dimension <- dim(image)
  buffer <- array(0, dim = c(dimension, 2)) 
  buffer[, , 1] <- image
  s <- 2
  d <- 1 
  # Este valor siempre será mayor que el potencial de cualquier configuración
  # de valores de pixeles
  V.max = (dimension[1] * dimension[2]) * ((256 ^ 2) / (2 * sdev) + 
    4 * weight.diff * max.diff)
  for(i in 1:iterations){
    # Switch source and destination buffers.
    if(s == 1){
      s = 2;
      d = 1;
    }
    else{
      s = 1;
      d = 2;
    }
    # Variamos cada pixel individualmente para encontrar los valores que
    # minimizan los potenciales locales.
    for(r in 1:dimension[1]){
      for(c in 1:dimension[2]){
        V.local = V.max
        min.val = -1
        for(val in 0:255){
          # Componente del potencial correspondiente a los datos observados.
          V.data = (val - image[r,c])^2 / (2 * sdev) 
          # Componente del potencial correspondiente a la diferencia entre los
          # pixeles vecinos
          V.diff = 0
          if(r > 1){
            V.diff = V.diff + min((val - buffer[r-1,c,s])^2, max.diff)
          }
          if(r < dimension[1]){
            V.diff = V.diff + min((val - buffer[r+1,c,s])^2, max.diff)
          }
          if(c > 1){
            V.diff = V.diff + min((val - buffer[r,c-1,s])^2, max.diff)
          }
          if(c < dimension[2]){
            V.diff = V.diff + min((val - buffer[r,c+1,s])^2, max.diff)
          }
          V.current = V.data + weight.diff * V.diff
          if(V.current < V.local){
            min.val = val
            V.local = V.current
          }
        }
        buffer[r, c, d] = min.val
      }    
    }
  }
  buffer[,,d]
}

perro_c <- icm(perro_n)
#writeJPEG(perro_n / 255, "imagenes/p_ruido.jpg")
#writeJPEG(perro_c / 255, "imagenes/p_recuperado_2.jpg")

pr_c <- pixmapGrey(perro_c/255)
plot(pr_c)

El siguiente es un ejemplo de clasificación de pixeles usando redes markovianas. En este ejemplo se tienen dos imágenes satelitales (1000x1000 pixeles cada una) del mismo lugar, tomadas en diferentes tiempos y se busca encontrar cambios. Los resultados se muestran abajo, la tercera imagen es la imagen de diferencias (usando IMAD_MAF) y las imágenes de abajo corresponden a dos algoritmos de detección de cambios, el primero es mezclas gaussianas y da resultados muy ruidosos. El segundo utiliza campos aleatorios de Markov y resulta en cambios menos ruidosos.

Referencias

Comparación entre modelos gráficos dirigidos y no dirigidos

Hemos estudiado dos tipos de modelos gráficos, cada uno tiene distintos puntos fuertes y débiles:

Ahora veremos como se relacionan las redes bayesianas y las redes markovianas, para ello veremos como pasar de un tipo de modelo al otro.

Redes Bayesianas a Redes Markovianas

Podemos ver la relación entre redes Bayesianas y Markovianas desde dos perspectivas: 1) Dada una red Bayesiana \({\mathcal B}\) como representar la distribución \(p_{\mathcal B}\) como una parametrización correspondiente a una red markoviana, o 2) Dada una gráfica dirigida \({\mathcal G}\) como represento las indepencias en \({\mathcal G}\) usando una gráfica no dirigida \({\mathcal H}\).

Respecto al primer punto, es fácil notar que las densidades marginales y condicionales que definen una red bayesiana son potenciales, por tanto una factorización de una gráfica dirigida en densidades condicionales corresponde a una factorización de una distribución de Gibbs donde la constante de normalización es \(Z=1\).

En cuanto a la representación gráfica, para transformar una red bayesiana en una red markoviana definimos un clique sobre cada familia (cada nodo y sus padres) en la red bayesiana. Definir un clique sobre cada familia corresponde en la gráfica a conectar los padres de cada nodo (en caso de que no exista esta arista) y eliminar la dirección de las aristas. El proceso de transformar una red bayesiana en una red markoviana se conoce como moralización debido a que se casan (o conectan) los padres de una variable, la definición formal es:

La gráfica moral \({\mathcal M[G]}\) de una red Bayesiana con estructura \({\mathcal G}\) sobre \(V\) (conjunto de nodos/variables aleatorias) es una gráfica no dirigida sobre \(V\) que contiene una arista entre \(X\) y \(Y\) si: a) hay una arista dirigida entre ellas (sin importar la dirección), ó b) \(X\) y \(Y\) son padres del mismo nodo.

Un corolario de la definición anterior es que si \({\mathcal G}\) es la estructura de una red bayesiana, entonces para cualquier distribución \(p_{\mathcal B}\) tal que \({\mathcal B}\) es una parametrización de \({\mathcal G}\), tenemos que \({\mathcal M[G]}\) es un mapeo de las independencias de \(p_{\mathcal B}\).

Ejemplo

Notemos que todas las independencias condicionales que se representan en la gráfica no dirigida también se leen en la gráfica dirigida, sin embargo, hay independencias condicionales en el modelo dirigido que no se representan en el modelo no dirigido. ¿Qué independencias implicadas en la red Bayesiana hemos perdido al moralizar la gráfica?

Del ejemplo anterior concluímos que el proceso de moralizar una gráfica dirigida puede conllevar que perdamos información de independencias; sin embargo, la siguiente proposición implica que moralizar es un mecanismo adecuado para transformar una gráfica dirigida en una no dirigida.

Sea \({\mathcal G}\) una gráfica asociada a una red Bayesiana. La gráfica moralizada \({\mathcal M[G]}\) es un mapeo mínimo de las independencias de \({\mathcal G}\).

La proposición anterior nos dice que si eliminamos una arista de \({\mathcal M[G]}\) estaríamos implicando relaciones de independencia que no se leen de la gráfica \({\mathcal G}\) y si añadimos aristas estaríamos perdiendo información de independencias de \({\mathcal G}\) que si se representan en \({\mathcal M[G]}\).

Vale la pena destacar que no siempre hay pérdida de información al moralizar una gráfica dirigida \({\mathcal G}\). Intuitivamente, la pérdida de información ocurre cuando se añaden aristas para conectar nodos, decimos que una red Bayesiana es moral si para cada par de variables \(X\), \(Y\) que comparten un hijo, existe una arista que une a \(X\) y \(Y\).

Si la gráfica dirigida \({\mathcal G}\) es moral, entonces su gráfica moralizada \({\mathcal M[G]}\) es un mapeo perfecto de \({\mathcal G}\). Esto es, todas las independencias que leemos en \({\mathcal G}\) se leen también en \({\mathcal M[G]}\).

Otra manera de leer este resultado es que las independencias en \({\mathcal G}\) que no estan presentes en la gráfica no dirigida que contiene las mismas aristas son las correspondientes a estructuras \(v\), a menos que la estructura \(v\) este protegida, en este último caso la gráfica dirigida no induce indepencias que no se lean también de la gráfica no dirigida.

Redes Markovianas a Redes Bayesianas

Consideremos ahora el problema de encontrar una red Bayesiana que se un mapeo minimal de las independencias de la red Markoviana. Veremos que en general la transformación en esta dirección es considerablemente más difícil.

Ejemplo

Consideremos la estructura de la red Markoviana de la figura de arriba y supongamos que buscamos una mapeo de esta red en una red bayesiana. Una manera es enumerar los nodos de la red y definir los padres de cada nodo en términos de las relacones de independencia que se leen de la red markoviana. Ordenemos los nodos de la siguiente manera: \(A,B,C,D,E,F\), la relación entre \(A\) y \(B\) es fácil, pero veamos que ocurre cuando añadimos a \(C\), introducimos \(A\) como padre de \(C\) pues no son independientes; sin embargo \(C\) tampoco es independiente de \(B\) condicional a \(A\) por lo que debemos añadir a \(B\) como padre de \(C\). Similarmemte consideramos \(D,E,F\) hasta obtener la gráfica dirigida de la derecha.

Es claro que hemos introducido aristas hasta formar una gráfica cordal (todos los ciclos están particionados en triángulos). Nos podemos preguntar si otro ordenamiento de los nodos podría conllevar introducir menos aristas, pero la respuesta es no: cualquier I-mapeo de una red Bayesiana que represente a la red markoviana debe introducir aristas que triangulan la gráfica produciendo una gráfica cordal.

Sea \({\mathcal H}\) una estructura de red Markoviana, y sea \({\mathcal G}\) una red bayesiana tal que es un mapeo minimal de las independencias de \({\mathcal H}\). Entonces \({\mathcal G}\) no puede tener inmoralidades.

Un corolario del teorema anterior es que \({\mathcal G}\) es necesariamente cordal. Esto se debe al proceso de triangulación que se requiere para convertir la red markoviana en red bayesiana.

En la transformación de una red no dirigda en una dirigida también perdemos información de independencias al introduci aristas. Volviendo al ejemplo anterior ¿Que relaciones de independencia se leen en la gráfica markoviana que no podemos leer en la red bayesiana?

Gráficas cordales

Hemos visto que la conversión entre redes bayesianas y markovianas puede resultar en la introducción de aristas que conlleva pérdida de informacióon de independencias implicada por la estructura de la gráfica original. Es interesante preguntarnos, ¿Cuándo un conjunto de supuestos de independencia se puede representar de manera perfecta usando ya sea una red markoviana o una red bayesiana? Resulta que esta es la clase de las gráficas cordales no dirigidas.

Sea \({\mathcal H}\) una red markoviana no cordal. Entonces no existe una red Bayesiana \({\mathcal G}\) que sea un mapeo perfecto de \({\mathcal H}\), esto es, no existe una red bayesiana de cuya gráfica se puedan leer todas las independencias que se leen de la red markoviana.