Máster Universitario en Fı́sica y Matemáticas

Seminario BIOMAT 2017

Trabajo presentado por

Raúl Enrı́quez Delgado

Curso 2016/17

Seminario de Invitados BIOMAT

1
 ANÁLISIS MULTIESCALA DE SIMILITUDES ENTRE IMÁGENES EN
VARIEDADES DE RIEMANN

La tercera intervención por parte de la doctora Coloma Ballester en el seminario BIOMAT

2017 sobre modelos matemáticos en imágenes biomédicas, tuvo como objetivo principal la
comparación de dos parches (patches) de una imagen definida en una variedad de Riemann,
que a su vez puede ser definida por el dominio de la imagen con una métrica adecuada que
dependerá de esta imagen.

Este enfoque, se basa principalmente en saber que la mayorı́a de los dispositivos prácticos
para el análisis de imágenes asumen un paso preliminar en el procesamiento de imágenes
que consiste en pasar de la imagen original a versiones suavizadas, que contiene información
totalmente significativa. En esta etapa del procesamiento, la imagen se transforma en nuevas
imágenes utilizando por ejemplo un parámetro que mida el grado de suavizado, o dicho de
otra forma, el tamaño de zonas vecinas que se utilicen para dar una estimación del brillo de
la imagen en un punto dado. Además, la comparación de imágenes es importante desde un
punto de vista computacional, ya que deriva muchas aplicaciones como el reconocimiento de
objetos, visión estéreo, interpolación de imágenes, eliminación de imágenes, imágenes basa-
das en ejemplares, etc.

El tamaño del parche no se determina a priori, y se lo identifica con una escala variable .
Se realiza la comparación de dos puntos, es decir desde un enfoque considerado como una
extensión no local, mediante un análisis multiescala, como se denomina generalmente en la
teorı́a del procesamiento de imágenes, por lo tanto tiende a dar información más fiable sobre
el nivel de gris en el pixel original fluctuante.

El análisis multiescala define una familia de transformaciones (Tt )t≥0 , las que se aplican
a la imagen original denotada como f (x), y éstas producen una secuencia de imágenes
u(t, x) = (Tt f )(x). Se podrı́a considerar Tt f como una versión semi-local de f donde las
zonas vecinas de tamaño t alrededor de x serán escanneadas para determinar el valor de
Tt f (x). Como se comento anteriormente, el ”tamaño” t no se define aún.

El estudio sistemático detallado del análisis multiescala (aproximación axiomática) para

imágenes fue desarrollado por (Alvarez et al., 1993); donde a partir de varios axiomas que
definen la estructura del espacio multiescala y un conjunto de invariantes geométricos y
fotométricos, se definieron soluciones de una ecuación en derivadas parciales parabólica de-
generada.

Matemáticamente, la idea fundamental de la aproximación axiomática es que dadas dos

imágenes u, v definidas en sus respectivos dominios (por ejemplo en R2 ), se quiere comparar
sus zonas vecinas en los puntos x, y ∈ R2 , respectivamente. La forma más sencillas de
compararlas serı́a usando la distancia Euclı́dea, por lo tanto se tiene
Z
D(t, x, y) = gt (h)(u(x + h) − v(y + h))2 dh, (1)
R2
donde gt es una ventana data que se asume Gaussiana de varianza t, y además se asume que
el dominio de la imagen es el plano Euclideano. Sin embargo, se deben definir medidas de
similaridad en el caso de imágenes sobre variedades de Riemann, entonces ya no serı́a posible

2
una fórmula explı́cita como (1). Se muestra que estas medidas están dadas por la solución
de una EDP elı́ptica degenerada en las variables (x, y),
2
Dt = ∆x D + 2Tr(Dxy D) + ∆y D, (2)
que representa un caso de una EDP lineal que expresa de forma simple la comparación multi-
escalar de dos parches de una imagen. En el caso de comparar parches de imagen definidos en
variedades Riemannianas, aparecerá una gran familia de posibilidades, derivada del enfoque
axiomático, donde se puede conseguir una especificación adicional incluyendo suposiciones
lineales o morfológicas. Además, se podrı́an considerar invarianza geométrica (traslación o
rotación) para imágenes definidas en Rn .

Uno de los ejemplos del análisis lineal multiescala de una medida de similitud es el modelo
D
= Tr(G1 (x)−1 Dx2 D) + 2Tr(G1 (x)−1/2 G2 (x)−1/2 Dxy
2
D) + Tr(G2 (y)−1 Dy2 D), (3)
t
donde Mi = (RN , Gi (x)), i = 1, 2, son dos variedades de Riemann, y x ∈ M1 , y ∈ M2 . En
el caso particular, si asumimos que las dos métricas son constantes en ambas imágenes, el
modelo es

Dt = Tr(At ADx2 D) + 2Tr(AB t Dxy D) + Tr(B t BDy2 D),

Rdonde A, B son dos matrices N × N . La medida de similitud multiescala es C(t, x, y) =

g (z)C(0, x + Ah, y + Bh)dh, donde gt es la Gaussiana de escala t, y C(0, x, y) = (I(x) −
R2 t
J(y)2 ) satisfacen la ecuación

Ct = Tr(At ADx2 C) + 2Tr(AB t Dxy C) + Tr(B t BDy2 C).

Por simplicidad se escribirá las variedades en el análisis multiescala de similaridades como
N = M1 × M2 . La métrica se denota por g = g 1 × g 2 , y G(ξ), G1 (ξ), G2 (ξ) las matrices co-
rrespondientes, donde ξ = (ξ1 , ξ2 ) ∈ N . Y (κ) := κn es una secuencia creciente de constantes
nonegativas.

Q((κ)) := {C ∈ Cb∞ (N ) : k Dα C k∞ ≤ κn ∀n ≥ 0 ∀|α| ≤ n}.

Se asume que Tt : Cb (N ) → Cb (N ) es un operador nolineal para cualquier t ≥ 0. Se
denotará C(t, x, y) = Tt C(x, y), C ∈ Cb (N ), asumiendo una conexión a priori de P sobre N .

La aproximación axiomática en el procesamiento de imágenes está inspirada en la derivación

sistemática de modelos en la mecánica del continuo desde los primeros principio, y se basa
en tres tipos de axiomas.

El primer tipo de axiomas se ocupa de la arquitectura del análisis de imágenes multiescala

(causalidad, recursividad, regularidad, localidad), llamados axiomas arquitectónicos.

[Recursividad] T0 (C) = C, Ts (Tt C) = Ts+t C ∀s, t ≥ 0, ∀C ∈ Cb (N ).

[Causalidad] Tt+s,s Ts para todo 0 ≤ s, t < ∞.

3
El axioma de recursividad es una versión fuerte de la causalidad que implica que la medida
de similitud a una escala más robusta se puede deducir de una más fina, que es una propiedad
natural en el análisis de la imagen.

Cuando M = RN con la métrica Euclideana el axioma es

Th C − C
[Generador infinitesimal] → A(C)ash → 0 + .
h
para C ∈ Cb∞ (RN × RN ).
En cierto sentido el sistema de coordenadas alrededor de cada punto es siempre el mismo
sistema canónico, por lo tanto están relacionados por la identidad

[Regularidad] kTt (C + hC̃) − (Tt (C) + hC̃)k∞ ≤ M ht ∀h, t ∈ [0, 1], ∀C, C̃ ∈ Q((κ)),

donde la constante M depende de Q((κ)).

[Localidad] Tt (C)(x) − Tt (C̃)(x) = o(t), t → 0+, x ∈ RN ∀C, C̃ ∈ Cb (N )

tal que Dα C(x) = Dα C̃(x) para todos los multi-ı́ndices α.

[Principio de comparación] Tt C ≤ Tt C̃ ∀t ≥ 0, ∀C, C̃ ∈ Cb∞ (N ), C ≤ C̃.

El principio de comparación es una propiedad que preserva el orden, lo que significa que si
una medida de similitud es siempre menor que otra, se puede aplicar el análisis multiescala
sin invertir esta relación. El análisis multiescala produce versiones de baja resolución de las
medidas de similitud, que deben ser consistentes con las iniciales.

Los axiomas morfológicos indican que el análisis de la imagen debe ser invariante bajo
las fluctuaciones de la luz y bajo los cambios de posición, orientación y escala de la forma
planar. La comprensión de las imágenes y formas están ligadas con la arbitrariedad de la
escala de grises en las imágenes perceptivas o digitales. Si algún punto de una región es más
brillante que otro en la imagen original, este orden deberı́a preservarse en la imagen final. Si
consideramos el análisis de la imagen independiente del nivel de grises, tenemos que

[Grey-level-shift invariance] Tt (0) = 0, Tt (C + κ) = Tt (C) + κ

∀t ≥ 0, ∀C ∈ Cb∞ (N ), ∀κ ∈ R.

Este axioma significa que no se hace una suposición a priori sobre el alcance de la claridad
de una imagen a observar, además que la interpretación de una fotografı́a es ampliamente
independiente del tiempo de exposición, es decir que puede ser oscura o clara y aún ası́
identificarse como esencialmente la misma imagen.

La forma fuerte de este axioma es

[Gray scale invariance] Tt (f (C)) = f (Tt (C)) ∀t ≥ 0, ∀C ∈ Cb∞ (N ),

para cualquier función estrictamente creciente f : R → R.

4
El caso lineal.

Es un caso particular para el análisis multiescala sobre similitud de funciones que satisfacen
el axioma de linealidad Tt (aC1 + bC2 ) = aTt (C1 ) + bTt (C2 ), ∀a, b ∈ R, ∀C1 , C2 ∈ Cb (N ).
Siendo Tt un análisis multiescala que satisface todos los axiomas arquitectónicos, el principio
de comparación, los nivel de grises de desplazamiento invariante y además se asume lineal,
entonces
∂C 2
= F (DN C, ξ, G),
∂t
donde

F (X, ξ, G) = c11 (ξ)Tr((G1 )−1 (ξ1 )X11 ) + 2c1,2 (ξ, G)Tr(D̃12 I 1 (ξ1 )−1 X12 )
+c22 (ξ)Tr((G2 )−1 (ξ2 )X22 ),

donde D̃12 es de la forma isométrica (Tξ1 M1 , G1 (ξ1 )) → (Tξ2 M2 , G2 (ξ2 )). La elipticidad de
F implica que c11 , c22 ≥ 0.

El caso de (Mr , g r ) = (RN , g r ).

Se considera M1 = M2 = M = RN y gijr (x) métricas generales en RN , r = 1, 2. Además, se

sabe que ei = GT (x)−1/2 fi es una base ortonormal de (Tx Mr , g r (x)) si fi es una base ortonor-
mal Euclideana. Dado I r (x) : (RN , g r (x)) → (RN , (g r )−1 (x)) viene dada por I r (x)ei = e∗i .
Entonces

I r (x) = Gr (x).
Se podrı́a definir P (x, y)(v) = G2 (y)−1/2 G1 (x)1/2 v, v ∈ RN , como una conexión a priori de
x y y. Dado que D̃1,2 : (RN , g 1 (x)) → (Rn , g 2 (y)) es una isometrı́a, en este caso se tiene que
D̃1,2 = G2 (y)−1/2 G1 (x)1/2 . Entonces

2c12 (x, y)D̃1,2 I 1 (x)−1 = 2c12 (x, y)G2 (y)−1/2 G1 (x)−1/2 .

La EDP obtenida es
∂C
= a(x, y)∆Mx C + 2c12 (x, y)Tr(G2 (y)−1/2 G1 (x)−1/2 Dxy C) + c(x, y)∆My C,
∂t
donde
∆Mx C = Tr(G1 (x)−1 (Dxx C(x) − Γ(1) (Dx C)(x))),
y una expresión similar para el operador ∆My .

Aproximaciones WKB (Wentzel-Kramers-Brillouin).

Este método de aproximación desarrolla dos fórmulas prácticas para aproximar la solución
correspondiente a la EDP lineal (3).
Z
C(t, p) = gt (h)C(0, x + G1 (x)−1/2 h, y + G2 (y)−1/2 h)dh, (4)

5
considerando el caso C(0, x, y) = (u(x) − v(y))2 .

Por otro lado, si se considera una distancia geodésica d(x, x0 ) en la variedad M1 , y x0 =

x + A(x)h, entonces
Z
H0 2
C(t, p) = √ e−d(x,x+A(x)h) /t C(0, x + A(x)h, y + B(y)h)dh. (5)
t

Experimento sobre la medida de similitud multiescalar

Con el fin de mostrar de cambiar t en toda la imagen, se calculan los mapas de disparidad
para un el par estéreo de la Figura 1. Las disparidades se obtuvieron mediante la búsqueda
exhaustiva de la mejor coincidencia en el rango de posibilidades de x ∈ [−55, −5], usando
pesos Gaussianos o pesos geodésicos mediante las fórmulas (4) y (5) respectivamente. Para
este experimento el parámetro de escala fue seleccionado para t ∈ {16.6, 25, 50, 150, 15000}.
El parámetro r se fijó en 150. Para los pesos geodésicos aproximados, los coeficientes se
establecieron como κspatial = 1.0 y κcolor = 3.0.

Figura 1: Par estéreo usado para el cálculo del mapa de disparidad

(a) Base (b) t = 16,6 (c) t = 15000

Figura 2: Mapas de disparidad

6
 (a) Máscara de oclusión (b) t = 16,6 (c) t = 15000

Figura 3: Errores en las disparidades

Como se puede ver, el cálculo de valores de similitud en una escala gruesa da lugar a ma-
pas de disparidad lisos que se ajustan estrechamente a la verdad del terreno dentro de los
objetos, pero que se superponen en los lı́mites. En contraste, las similitudes de escala fina
hacen hincapié en los lı́mites agudos de los objetos, pero producen una gran cantidad de des-
ajustes en regiones planas. En la Figura 2 los colores más brillantes corresponde a mayores
compensaciones, y en la Figura 3 los colores más oscuros corresponden a errores más grandes.

La estructura afı́n de un campo vector se utilizará como una métrica anisotrópica en

el dominio de la imagen y nos permitirá construir un mapa de conexión a priori para com-
para adecuadamente dos imágenes. El cálculo de los tensores covariantes está estrechamente
relacionado con el problema de estimar regiones covariantes afines. Dado dos puntos de coin-
cidencia, la transformación afı́n se estima entonces hasta una rotación.

Para la construcción de tensores estructurales y sus vecinos correspondientes los cuales son
afines covariantes, se considera u una imagen dada, u RN → R. Sea GL(N ) un conjunto de
matrices invertibles en RN . Sea A ∈ GL(N ) y denotamos uA (x) := u(Ax).

Por definición, Hu sea un tensor (1, 1) definido en RN tal que para cada x ∈ RN , está
representado por una matriz N × N , Hu (x) mapea un vector en RN a otro vector en RN . Se
dice que Hu es un tensor afı́n covariante si satisface

HuA (x) = At Hu (Ax)A,

donde uA (x) := u(Ax) para A ∈ GL(N ).

Por ejemplo, dada una imagen u : R2 → R, y considerando los tensores covariantes afines
N T (k) (u)(x) y vecindades covariantes afines BN T (k) (u) (x, r) para x ∈ R2 , k ∈ N, r > 0. Se
observa que para x en el dominio de la imagen, cualquiera de estos tensores está garantizado
de ser un tensor covariante afı́n. El propósito del algoritmo iterativo no es asegurar la co-
varianza afı́n, sino disminuir la dependencia de la primera iteración, que depende solamente
de un solo gradiente y por lo tanto es muy sensible al ruido. Entonces, dada una imagen u,
se calculan los tensores en cada punto x con un número diferente de iteraciones k. Luego se
(k)(u)(x) (k+1)(u)(x)
calculan las normas de Frobenius de la diferencia kNT − NT kF para dos valores
consecutivos de k. Finalmente, se promedia las normas sobre todo x.

7
La Figura 4 muestra estos cambios promedio sobre k para tres imágenes seleccionadas. El
hecho de que los valores vayan hacia un valor distinto de cero se explica por la alternancia
ocasional entre dos tensores covariantes afines. Este comportamiento es tı́pico para todas las
imágenes probadas.

Figura 4: Cambio promedio de los tensores desde una iteración a otra, dependiendo del número de
iteraciones.

En el caso de la Figura 5, los vecindarios covariantes afines BN T (k) (u) (x, r) son elipses para
N = 2, y en el experimento se tomó k = 30.

Figura 5: Imagenes que muestran algunas vecindades afines.

El parámetro r es libre y controla el tamaño del vecindario covariante afı́n en un punto dado.
Por otro lado, el tamaño del vecindario también se ve afectado por la textura en la vecindad
de ese punto.

El método iterativo propuesto para el cálculo de tensores estructurales covariantes afines nos
permite estimar la conexión a priori entre dos puntos dados y por lo tanto una afinidad local.

8
La Figura 6 ilustra los valores C(t, x, y) de esta medida de similitud calculada entre un punto
fijo x en la imagen de referencia y todos los puntos y de la imagen secundaria, que se toma
desde un punto de vista diferente.

Figura 6: Ilustración de una medida de similitud.

Los valores de la medida de similitud se calculan entre un punto x de la imagen de refe-

rencia (primera fila, izquierda) y todos los puntos y de la imagen secundaria (primera fila,
derecha), que se toma desde un punto de vista diferente. En la segunda fila se muestra un
primer plano que contiene el punto x y el valor de la medida de similitud. Los pı́xeles rojos
indican valores más bajos de C(t, x, y). Se observa además que los mı́nimos de la medida de
similitud ocurren en los puntos donde la estructura de la imagen secundaria es similar al
parche de referencia. Las ventanas de comparación son elipses que corresponden a la bola
unitaria asignada de acuerdo con A y B. La tercera fila ilustra el paisaje de similitud para
una escala mayor. La variación de la escala del análisis corresponde en este caso al aumento
del tamaño de la ventana.

9
Referencias
[1] Alvarez, L., Guichard, F., Lions, P., Mores, J., 1993. Axioms and fundamental equations
of image processing, Arch. Ration. Mech. Anal., 123.
[2] Ballester, C., Calderero, F., Caselles, V., Facciolo, G., 2014. Multiscale analysis of si-
milarities between images on Riemannian manifolds, Multiscale Model. Simul.
[3] Fedorov, V., Arias, P., Sadek, R., Facciolo, G., Ballester, C., 2015. Linear Multiscale
Analysis of Similarities between Images on Riemannian Manifolds: Practical Formula
and Affine Covariant Metrics, Society for Industrial and Applied Mathematics, SIAM
J. Imaging Sciences, Vol. 8, No 3.

10