Revista Mexicana de Ingeniería Biomédica

Registro Elástico de Imágenes Médicas Multimodales

Basado en Medidas de Variabilidad Local

INTRODUCCIÓN

El registro de imágenes se ha convertido en una tarea muy importante en el área de procesamiento digital de imágenes médicas, ya que puede ser empleado en varios procesos como: caracterización de los cambios anatómicos del corazón en un ciclo cardíaco o la atrofia gradual del cerebro en el envejecimiento; comparación o modelado de estructuras anatómicas (morfometría); segmentación de estructuras o tejidos por medio de átlases médicos; corrección de los artefactos causados por el movimiento en imágenes fetales; detección de las discrepancias entre imágenes de un mismo paciente que fueron adquiridas durante distintas etapas de algún tratamiento; entre muchas otras [1, 2, 3, 4]. De acuerdo a los tipos de deformaciones que se deseen aproximar con el registro, éste puede dividirse en dos: en el registro paramétrico y en el elástico o no paramétrico [5]. En la literatura, las técnicas basadas en el descenso (o ascenso) de gradiente [6], [7] son las más comúnmente utilizadas para optimizar una medida de similitud (ej. Información Mutua [8, 9, 10]) y obtener la transformación espacial (ej. afín o perspectiva [11]) que alinee las imágenes candidata y patrón [5]. Otras opciones basadas en métodos de optimización global como los algoritmos genéticos [12] y filtrado de partículas [13], han comenzado a ganar terreno dentro de este campo [14, 15, 16]. Por otro lado, el registro elástico es un problema más complejo que el registro paramétrico [3], especialmente cuando se intentan registrar imágenes provenientes de distintas fuentes (registro multimodal). El registro elástico es un área de investigación relativamente reciente donde la mayoría de los algoritmos desarrollados se encuentran en sus primeras etapas de evaluación y validación [17], y pocas veces buscan resolver el problema para imágenes multimodales, y de hacerlo, generalmente se centran en parametrizar el espacio de deformaciones para optimizar una medida de similitud obligando en ocasiones al usuario a ubicar puntos de interés en las imágenes, lo que resulta en procesos con muy alto costo computacional y/o semi-autónomos [18, 19, 20, 21, 22, 23]; razones por las que se ha dificultado su uso en tareas de imagenología médica.

Una propuesta reciente para resolver el registro elástico se enfoca en el uso del flujo óptico (FO) [24] para encontrar los desplazamientos faltantes, después de haber realizado una aproximación inicial del campo de deformación por medio de un registro paramétrico basado en el uso del filtro de partículas (FP), donde se estiman los parámetros de una matriz de transformación [25]; este método muestra resultados prometedores en [26] y [27]. Por otro lado, una limitante de dicho algoritmo es la restricción para uso solo en imágenes unimodales o la necesidad de contar con una función de transferencia inyectiva de intensidades entre las imágenes a registrar, lo que no se tiene en los casos médicos de registro multimodal.

Dentro del contexto descrito en los párrafos anteriores, en este artículo se propone una metodología para resolver dos problemas importantes en la imagenología médica, pero poco abordados simultáneamente en la literatura: el registro elástico de imágenes y a la vez multimodal. En este estudio lo que se propone es mapear las imágenes a un espacio donde sus intensidades puedan ser comparadas, y entonces sea posible abordar el problema sin parametrizar el espacio de deformaciones al aplicar iterativamente una técnica de FO como en [25], bajo un esquema de espacio de escalas [28].

El artículo está organizado de la siguiente manera: en la sección 2 se describe la teoría necesaria para la comprensión de la metodología propuesta; en las subsecciones 2.1 y 2.2 se revisa el registro paramétrico basado en el FP y el registro elástico basado en FP más FO, respectivamente, y los detalles del algoritmo propuesto se describen en la subsección 2.3. Después, en la sección 3 se explican los experimentos realizados para llevar a cabo la evaluación del algoritmo de registro elástico multimodal propuesto, empleando datos generados sintéticamente. A continuación, en la sección 4 se presentan los resultados obtenidos de los experimentos realizados para la evaluación (cuantitativos y cualitativos), y en la sección 5 se realiza el análisis y discusión de dichos resultados. Finalmente, en la sección 6 se muestran las conclusiones derivadas de este trabajo, así como el posible trabajo a futuro sobre la metodología presentada.

Metodología

El registro elástico puede ser visto como encontrar el campo vectorial de desplazamientos V (r) que alinee la imagen candidata I_C(r) con la patrón I_P (r),

(1)

donde r = (x,y)^⊤ denota una posición dentro del dominio rectangular Ω ⊂ ℝ² de las imágenes, y F[⋅] representa una relación entre las intensidades de ambas imágenes. De acuerdo con la ecuación (1), si el registro es unimodal, es decir si F es la identidad, el problema puede ser visto como obtener el FO entre las imágenes. Sin embargo, para simplificar la estimación del FO es mejor restringir la búsqueda a desplazamientos pequeños entre los pixeles correspondientes. Por esta razón, antes de calcular el FO, es conveniente acercar los pixeles entre ambas imágenes por medio de una transformación inicial usando algún método de registro paramétrico, por ejemplo el basado en el FP [25].

Registro paramétrico basado en el filtro de partículas

La idea principal del registro paramétrico basado en el FP consiste en estimar el vector de parámetros desconocidos θ de una transformación, afín o perspectiva, por medio de una búsqueda estocástica sobre una superficie de optimización (función de costo) [29]. Este objetivo se logra usando un conjunto de N puntos de prueba llamados partículas (θ₁,…,θ_N), y sus pesos asociados W_j, calculados por medio de una función de verosimilitud entre las imágenes,

(2)

donde σ_η² es la varianza del ruido en la medición z, T(r|θ) es una transformación paramétrica dependiente del vector θ, y NMI(⋅,⋅) representa la información mutua normalizada [30] entre dos imágenes; ésta medida de similitud ofrece el mejor desempeño para este método como se reporta en [31].

Los pesos W_j son utilizados para aproximar de forma iterativa la función de distribución de probabilidad a posteriori P(θ|z) de los parámetros desconocidos θ con respecto a una medición z [13, 32]. Ésta distribución se aproxima considerando una ventana de k mediciones (z¹,…,z^k), la cual es usada para estimar el vector de parámetros de la transformación ^k, tomando por ejemplo el valor esperado de la distribución a posteriori:

(3)

Para ver más detalles teóricos y de implementación de este algoritmo, el lector puede dirigirse a [13, 14, 16, 26, 31].

Registro elástico basado en flujo óptico

Una vez logrado el registro paramétrico empleando el FP, los desplazamientos existentes entre los pixeles correspondientes de las imágenes, en teoría deberían ser pequeños. Además, sí las imágenes son unimodales, entonces es posible encontrar dichos desplazamientos sin parametrizar el espacio de deformaciones empleando una técnica de estimación de FO. Esto se puede lograr bajo un enfoque solamente espacial dado que no se cuenta con una secuencia de imágenes variante en el tiempo como en los problemas de FO tradicionales [24].

Sea w = ^⊤ el vector de desplazamientos a calcular, donde u y v son los desplazamientos en x y en y respectivamente, I(x + u,y + u,t + 1) = I(x,y,t), con t = 0 para nuestro problema en particular; la mayoría de los métodos de FO de la literatura pueden ser representados por medio de la siguiente función de energía [33]:

(4)

donde ψ₁(⋅) es una función que penaliza el error en el término de datos, ψ₂(⋅) penaliza el término de regularización pesado por una constante positiva α, ∇I = ^⊤ es el gradiente de la imagen donde los subíndices representan las derivadas parciales, y dado un kernel Gaussiano K_ρ(r) con desviación estándar ρ; la función J_ρ(⋅) está dada por

(5)

donde * representa el operador convolución. Las técnicas de FO pueden clasificarse en métodos locales que optimizan una expresión local de energía (más robustos al ruido), y las estrategias globales que minimizan una función global de energía (más aptos para estimar flujos densos) [33]. De la ecuación (4), variando los valores de ρ, α y las funciones ψ₁(⋅) y ψ₂(⋅), es posible derivar por ejemplo, el método local propuesto en [34] sí ψ₁ es la identidad, ρ > 0 y α = 0; o la técnica tanto global como local en [33] si ρ > 0 y α > 0.

Para este trabajo, nuestro interés es observar el desempeño del algoritmo propuesto en la estimación del campo denso de desplazamientos, dado que no se está enfocando a una aplicación médica en específico, por lo que un método global resulta lo más viable para probar el desempeño de la propuesta. Entonces, tomando ψ₁ como la identidad, ψ₂(w) = |∇u|² + |∇v|², ρ = 0 y α > 0, se obtiene la técnica global de FO propuesta en [35], y de la cual se basa la mayoría de las propuestas globales existentes [24]. La solución al problema de optimización en (4) de acuerdo al planteamiento en [35], se encuentra resolviendo el sistema de ecuaciones en (6) para cada pixel, empleando para ello un método iterativo (por ejemplo Gauss-Seidel [36]):

un+1 = un - I x[Ixun + I yvn + I t]∕[α2 + I x2 + I y2],
vn+1 = vn - I y[Ixun + I yvn + I t]∕[α2 + I x2 + I y2],		(6)

donde u y v son el promedio de los desplazamientos en el vecindario en torno al pixel de interés. Para más detalles teóricos y de implementación del algoritmo de FO empleado, el lector interesado puede consultar [35, 37].

Por otro lado, una evaluación del método de registro propuesto empleando diversas técnicas de FO debe ser contemplada en un futuro de acuerdo a la aplicación médica que se pretenda abordar; ya que dependerá del problema el tipo de evaluación requerida, así como las características de las imágenes. A diferencia de las técnicas que parametrizan el espacio de deformaciones, abordar el problema de registro elástico desde esta perspectiva, convierte a la estimación del campo vectorial denso en un problema de optimización sobre una superficie convexa, cuya solución está dada por un sistema lineal de ecuaciones, lo que garantiza un cómputo rápido en la estimación de dicho campo. Además, el hecho de no parametrizar el espacio de deformaciones incrementa la gama de deformaciones que pueden ser logradas (cumpliendo ciertas restricciones de homogeneidad), debido a que éstas no se ven limitadas por los dominios de parámetros.

Así pues, de acuerdo con [25] es posible aplicar de forma iterativa un registro paramétrico y el FO, y acumular los campos vectoriales obtenidos en cada iteración hasta satisfacer cierto criterio de paro. Para el caso del registro multimodal, proponen aproximar la función de transferencia de intensidades entre las imágenes usando el histograma conjunto, pero este enfoque sólo es factible si el mapeo de intensidades F[⋅] en (1) es inyectivo. Por lo tanto, en este trabajo se propone una extensión de esta metodología de registro elástico con la finalidad de solventar las limitaciones mencionadas anteriormente, para lo cual después de un registro paramétrico se aplica un mapeo de intensidades para las imágenes candidata y patrón usando medidas de variabilidad local, y enseguida se realiza un proceso iterativo por medio del FO en un espacio de escalas.

Algoritmo propuesto para el registro elástico multimodal

La idea principal del algoritmo propuesto es el mapear las imágenes a registrar, después de haber sido alineadas por una transformación paramétrica, a un espacio donde la intensidad de cada pixel en una de las imágenes pueda ser comparada con la de su correspondiente en la otra imagen, a pesar de sus características multimodales, es decir, establecer un mapeo G tal que

(7)

Para definir dicho mapeo, sugerimos el uso de medidas que no dependan del nivel de gris del pixel, pero si de la variabilidad de las intensidades en los pixeles vecinos. Dos medidas que cumplen con lo anterior son la entropía [38] y la varianza, calculadas sobre una ventana centrada en el pixel de interés, y a las cuales nosotros llamamos Medidas de Variabilidad Local (MVL). Así, la metodología propuesta para el registro elástico multimodal está compuesta por los siguientes pasos:

1. Registro paramétrico. Encontrar el vector de parámetros (usando el FP) de una transformación de perspectiva T(r|) que proporcione la mejor alineación entre las dos imágenes, I_P e I_C; obtener la imagen registrada I_R(r) = I_C(T(r|)), y el correspondiente campo vectorial de desplazamientos V _R(r) = T(r|) - r.
2. Mapeo basado en MVL. Aplicar el mapeo G basado en una MVL (entropía, varianza o una combinación de ambas) sobre todos los pixeles en ambas imágenes (I_P e I_R). Esto es, obtener Ĩ_P (r) = G[I_P (r)] e Ĩ_R(r) = G[I_R(r)] de acuerdo a las 4 siguientes propuestas:

   (8)

   (9)

   (10)

   (11)

   donde MVL₁ representa la medida de variabilidad local usando la entropía, MVL₂ empleando la varianza, MVL₃ una ponderación Euclideana entre la varianza y la entropía en cada pixel, y MVL₄ considerando el valor máximo entre la varianza y la entropía nuevamente por pixel. Además, N_r representa el conjunto de pixeles de una ventana centrada en r y tamaño n × n, p_r(I(s)) es la distribución de probabilidad de las intensidades de la imagen I(s) asociada a la ventana N_r, y μ_r es el valor esperado de las intensidades en la misma ventana, dado por

   (12)

3. Ecualización. Después de aplicar el mapeo G, las intensidades de las imágenes resultantes (Ĩ_P e Ĩ_R) tienen valores pequeños concentrados en un rango dinámico muy estrecho. Por esta razón, es necesario escalar dichas intensidades en un rango de 0 a 255 (escala de grises), y aplicar una ecualización de histograma [11].
4. Espacio de escalas. Con la intención de poder aproximar deformaciones más complejas encontrando desplazamientos más grandes, el FO se aplica empleando una técnica de espacio de escalas [28] (ver Figura 1):

   Dadas las dos imágenes ecualizadas Ĩ_P e Ĩ_R de tamaño 2^K × 2^K, y una subescala m (escala más pequeña sobre la que se desea calcular el FO):

   1. Tomar como cero el valor inicial de los desplazamientos para el FO en la escala m, V ₀^m = 0.
   2. Escalar las imágenes a un tamaño 2^K-m × 2^K-m; aplicando un suavizado con un kernel Gaussiano, de tamaño n × n y desviación estándar σ_G, a la imagen que se desea escalar y después realizando un proceso de submuestreo.
   3. Estimar el FO entre Ĩ_P ^m e Ĩ_R^m, y establecer t = 0.
      1. Obtener una imagen auxiliar, dada por I_A = Ĩ_R^m(r + V _t^m(r)).
      2. Incrementar t, t = t + 1.
      3. Calcular el FO V _t^m entre Ĩ_P ^m e I_A.
      4. Acumular los desplazamientos obtenidos en esta escala, V _t^m = V _t^m + V _t-1^m.
      5. Si ∑ _r|V _t^m(r) - V _t-1^m(r)| < ϵ, ir a d), si no regresar a (i).
   4. Ir a la siguiente escala m = m - 1.
   5. Si m ≥ 0, propagar los desplazamientos en la siguiente escala V ₀^m(r) = 2B(V ^m-1(r)) e ir a b), si no V _E(r) = V ⁰(r) y continuar con el paso 5. Aquí B(⋅) representa una función de interpolación (ej. bilineal [11]).
5. Obtener el registro elástico. Finalmente, se suman los campos vectoriales obtenidos en los pasos 1 y 4, V (r) = V _R(r) + V _E(r), y después es posible obtener la imagen candidata registrada elásticamente I_E(r) = I_C(r + V (r)), tal que I_P (r) ≈ F[I_E(r)].

---

---

En el algoritmo descrito anteriormente, se utilizó una escala inicial de m = 5 y para suavizar las imágenes antes de realizar el submuestro, se empleó un kernel Gaussiano de tamaño n = 5 con desviación estándar σ_G = 0,3(n∕2 - 1) + 0,8. También se debe tener en cuenta el umbral de convergencia ϵ usado en el paso 4.c).(v), y para el cual se puede establecer un valor con base a la resolución deseada en el FO y el tamaño de las imágenes, o puede ser remplazado por algún otro criterio de parada, por ejemplo calculando una medida de similitud como la entropía conjunta entre las imágenes, que sirva como indicador del estado del registro. Por otro lado, el tamaño de ventana en pixeles establecido para el cálculo de la entropía es de 7 × 7 y de 13 × 13 para la varianza, donde éstos valores fueron obtenidos por medio de una evaluación numérica del método con respecto al tamaño de ventana. Sin embargo, el tamaño óptimo de ventana puede variar dependiendo de las características de las texturas y del contraste con respecto al fondo de las imágenes a registrar. Adicionalmente, debido a la factibilidad del FP para ser paralelizado, éste fue implementado en una arquitectura de procesamiento en paralelo, lo que incrementa en gran medida la velocidad de cómputo del algoritmo.

Experimentos

Con el propósito de analizar el comportamiento del algoritmo propuesto con datos reales y para mostrar su eventual desempeño en el registro de átlases médicos (donde se registran los datos de un paciente con una muestra típica o promedio de varias imágenes de diferentes pacientes) [39], comenzamos resolviendo el problema de registro elástico de dos imágenes de resonancia magnética T1 en corte transversal de distintos pacientes tomadas de [40] (ver Figura 2).

Además, con la finalidad de obtener datos numéricos que nos proporcionen más información sobre el desempeño del algoritmo, se simularon 20 pares de imágenes alineadas de resonancia magnética T1 y tomografía computarizada (MR/CT), empleando el software de simulación reportado en [41] y disponible en línea en http://www.bic.mni.mcgill.ca/brainweb/. Siete de los pares fueron en corte sagital, siete en axial y seis en coronal, con un tamaño de pixel de 1×1 mm y con un nivel de inhomogeneidad de intensidades de entre 0% y 40%. También se adquirieron una tomografía computarizada y una resonancia magnética T1 en corte transversal (tamaño de pixel de 1×1 mm), Figuras 3.(d) y 3.(g), tomadas de diferentes pacientes con un tumor cerebral en el Hospital San Raffaele y proporcionadas por el Institute of Molecular Bioimaging and Physiology, ambos en Milan, Italia. Con éstas muestras se simularon dos pares de imágenes multimodales alineadas usando dos polinomios de grado 6 y 9 como mapeo de intensidades (ver Figura 4). Después, a una imagen de cada par simulado se le aplicó un campo vectorial de deformación elástico generado a partir de Thin Plate Splines [42], los cuales se tomaron como campo de referencia (ground truth GT, ver Figura 3) y se utilizaron para calcular el error del algoritmo de registro en la estimación del campo vectorial de la deformación, usando como métrica la norma Euclideana.

Figura 2: Registro elástico entre dos pacientes diferentes, el cual presenta un escenario similar al problema de registrar una imagen con la de un atlas médico. (a) Imagen candidata tomada de [40], (b) Imagen patrón tomada de [40], (c) Imágenes superpuestas antes del registro en RGB. En la tercera columna la imagen patrón se encuentra en el canal verde y la candidata en el rojo.

Dado que conocemos el mapeo de intensidades en dos pares de imágenes simulados, también se midió el error del algoritmo obteniendo las diferencias de intensidades entre estos pares.

Otro factor importante a tener en cuenta dentro de los algoritmos computacionales que emplean técnicas de optimización, es garantizar la convergencia del método. Por lo tanto, para establecer un criterio de convergencia para el algoritmo propuesto, se realizó un análisis usando la entropía conjunta entre las imágenes con respecto al número de iteraciones.

Finalmente, cabe mencionar que existen algoritmos de registro elástico que no parametrizan el espacio de deformaciones [43], pero no encontramos en la literatura, técnicas multimodales que sigan este mismo enfoque aplicado en imágenes médicas, motivo por el cual no se compara esta propuesta con otros algoritmos de registro elástico multimodal.

Resultados

En la Figura 5 se presentan los resultados visuales de registrar las imágenes de las Figuras 2.(a) y 2.(b), utilizando las cuatro MVL propuestas para el algoritmo. De igual forma en la Figura 6 se muestra el resultado del registro entre las imágenes en las Figuras 3.(g) y 3.(h), empleando las cuatro MVL. Calculando el error usando la norma Euclideana entre el GT y los vectores de desplazamiento obtenidos para cada pixel al registrar la imagen de la Figura 3.(a) con 3.(b), se obtuvieron los errores que se muestran en la Figura 7, usando la entropía como MVL en la Figura 7.(a), la varianza en 7.(b), la ponderación Euclideana entre entropía y varianza en 7.(c), y el máximo entre ellas mismas en 7.(d). También se calculó el error en las intensidades para el registro de las imágenes de las Figuras 3.(d) y 3.(e) (dado que conocemos su transformación sintética de intensidades tonales), el cual se puede observar en la Figura 8.(a) para la MVL₁, 8.(b) con la MVL₂, 8.(c) para la MVL₃ y 8.(d) con la MVL₄.

En la Tabla 1 se presentan la media del error y su desviación estándar, calculado con la norma Euclideana, del campo vectorial estimado por el algoritmo de registro elástico con las cuatro MVL, donde también se incluye la media del error obtenida con la versión unimodal del algoritmo que nos servirá como punto de referencia, ya que idealmente los resultados de las MVL tendrían que ser muy cercanos a estos valores y que es factible calcular dado que se utilizaron transformaciones de intensidad sintéticas. De la misma forma, en la Tabla 1 se muestra la media y desviación estándar del error en las intensidades obtenidas por las cuatro MVL con los dos pares de imágenes de los cuales se conoce su mapeo de intensidades que se muestran en la Figura 4.

Por último, como parte de los resultados en la Figura 9 se muestran las gráficas de la entropía conjunta y la media del error en el campo vectorial de las imágenes registradas, ambas con respecto al número de iteraciones del algoritmo, con la intensión de mostrar gráfica y numéricamente la convergencia del método presentado.

Figura 3: Tres ejemplos de problemas sintéticos de registro generados para evaluar el algoritmo. (a), (d) y (g) imagen candidata. (b), (e) y (h) imagen patrón (creada artificialmente a partir de la candidata). (c), (f) e (i) deformación sintética para generar la imagen patrón.

Figura 4: Mapeos de intensidades empleados para generar las imágenes sintéticas. Para la imagen mostrada en la Figura 3.(h) se empleó el mapeo 1, y para la mostrada en la Figura 3.(e) el mapeo 2.

Discusión

Observando la Figura 5, se puede apreciar que las estructuras del cráneo, la tráquea, el tracto nasal, el bulbo y la médula espinal, se alinean congruentemente (en tonalidades amarillo verdosas) en tres de las cuatro MVL (entropía, varianza y obteniendo el máximo); con respecto a otras estructuras o tejidos, es difícil para alguien sin formación médica y experiencia en análisis de imágenes de este tipo, tener la certeza de que se encuentran alineadas correctamente todas las estructuras y tejidos. Para realizar un análisis cuantitativo sobre estos resultados, se requeriría del uso de técnicas de segmentación así como de índices de evaluación sobre ciertos estándares médicos, además de que la evaluación dependería de los tejidos y/o estructuras que un médico especialista quisiera analizar en las imágenes [17, 27], asunto que no era de interés para este trabajo por el momento, ya que el estudio está orientado a probar la eficiencia del algoritmo en la estimación del campo vectorial denso de la deformación. En este punto cabe señalar que las imágenes registradas en este experimento a pesar de ser ambas resonancias magnéticas T1 (ver Figura 2), las intensidades en escala de grises no representan un mapeo biyectivo entre ambas, debido a que pertenecen a diferentes pacientes y fueron adquiridas bajo diferentes condiciones; además de que la deformación elástica requerida para lograr alinearlas no es fácil de estimar y es muy poco probable que se presente en casos reales de registro dentro de la imagenología médica, sin embargo se optó por probar con dichas imágenes para tener una idea de la capacidad del algoritmo presentado para estimar este tipo de deformaciones, además de ser un problema similar al que se tiene en el registro de átlases médicos [39].

Figura 5: Resultados obtenidos de aplicar el algoritmo de registro elástico multimodal sobre las imágenes de la Figura 2, empleando las cuatro MVL propuestas. (a) Entropía o MVL₁, (b) varianza o MVL₂, (c) ponderación Euclideana entre entropía y varianza o MVL₃, y (d) el máximo entre ellas mismas o MVL₄. Las imágenes se muestran superpuestas en RGB, donde la imagen patrón está en el canal verde y la registrada en el rojo.

Figura 6: Resultados obtenidos de aplicar el algoritmo de registro elástico multimodal sobre las imágenes de la Figura 3.(g) y 3.(h), empleando las cuatro MVL propuestas. (a) Entropía o MVL₁, (b) varianza o MVL₂, (c) ponderación Euclideana entre entropía y varianza o MVL₃, y (d) el máximo entre ellas mismas o MVL₄. Las imágenes se muestran superpuestas en RGB, donde la imagen patrón está en el canal verde y la registrada en el rojo.

Figura 7: Error en la estimación del campo vectorial para el registro elástico multimodal entre las imágenes de las Figuras 3.(a) y 3.(b), usando las MVL. (a) Entropía o MVL₁, (b) varianza o MVL₂, (c) ponderación Euclideana entre entropía y varianza o MVL₃, y (d) el máximo entre ellas mismas o MVL₄.

Figura 8: Error en las intensidades para el registro elástico multimodal entre las imágenes de las Figuras 3.(d) y 3.(e), usando las MVL. (a) Entropía o MVL₁, (b) varianza o MVL₂, (c) ponderación Euclideana entre entropía y varianza o MVL₃, y (d) el máximo entre ellas mismas o MVL₄.

Figura 9: Ejemplo de gráfica de la entropía conjunta (izquierda) y del error en la estimación del campo vectorial (derecha) con respecto al número de iteraciones del algoritmo. Gráficas obtenidas de registrar el par de imágenes de las Figuras 3.(d) y 3.(e).

Por otra parte, observando más a detalle la Figura 5.(a) donde se utilizó la entropía como MVL, se puede apreciar que las zonas donde existen alineamientos poco congruentes (pixeles totalmente rojos o verdes), son principalmente en los bordes del cráneo, y en la Figura 5.(b) donde la MVL utilizada fue la varianza, dichas zonas están principalmente en la materia gris.

Esto nos hace suponer que los errores de alineamiento en los bordes del cráneo son debidos a que el cálculo de la entropía se ve afectado por el fondo negro que cubre una gran parte de la ventana (N

r)cuandoseencuentracentradaenestaszonas;encontraparte,esteefectonosepresentaconlavarianzaquedefinemejorlosbordesdebidoalcontrasteconelfondo.Demanerasimilarsetieneunefectocontrarioconlamateriagris,dondenoexistemuchocontrasteentretejidos,perositexturasmuymarcadas.DebidoaestefenómenoobservadoenMV L_1yenMV L_2,esquesedecidióprobarcombinandoambasmedidas(MV L_3yMV L_4)conlaintensióndecomplementarlainformaciónentreellaseintentarreducirloserroresenlaalineación.EnlosresultadosobtenidoscombinandoambasMV LmostradosenlaFigura5.(c)usandounaponderaciónEuclideana,yenla5.(d)empleandoelmáximovalorentreellas,visualmentenoparecenmuyfavorablesparalaMV L_3,perolaMV L_4sipresentaunregistromuysimilaralasotrasdos.

Analizando un caso más cercano a los problemas reales de registro en imagenología, como lo es el correspondiente a los resultados mostrados en de la Figura 6, se puede observar que las regiones que presentan un mayor error, es decir zonas no traslapadas correctamente (regiones totalmente rojas o verdes) son mayores para la MVL₂ y la MVL₃ (Figuras 6.(b) y 6.(c)); contrario a lo que presentan la MVL₁ y la MVL₄ (Figuras 6.(a) y 6.(d)), donde se aprecia que en su mayoría todas las regiones están correctamente traslapadas (tonos amarillos verdosos). Por lo que en este escenario, la información proveniente de la entropía, permitió tener un correcto alineamiento, es decir con MVL₁ y MVL₄.

Con respecto a los resultados numéricos, se comprueba lo descrito en el párrafo anterior, ya que analizando las Figuras 7 y 8, así como la Tabla 1, se observa que las MVL que presentan el menor error para los casos sintéticos analizados, son la entropía (MVL₁ con un error medio menor a 1 mm) y el máximo entre las medidas (MVL₄ con un error medio de 1 mm), ambas con desviaciones estándar menores a 0.5 mm; a diferencia de la varianza y la ponderación Euclideana que presentan errores mayores a 1 mm con una desviación estándar mayor a 0.5 mm (ver Tabla 1). Además, los errores en las intensidades son para MVL₄ los más cercanos a los obtenidos por la versión unimodal del algoritmo (ver Tabla 1). Una diferencia que se puede observar entre MVL₁ y MVL₄ es que los errores son más homogéneos para el máximo entre la entropía y varianza, que para la entropía por si sola (ver Figuras 7.(a), 7.(d), 8.(a) y 8.(d)), probablemente debido a la información adicional que proporciona la varianza y entropía combinadas en esta MVL.

Por último, debido a que el algoritmo propuesto es un método iterativo, es necesario establecer las condiciones que deben cumplirse para detener el proceso, además de asegurar su convergencia a una solución. Así pues, observando las gráficas de la Figura 9 se puede apreciar que la entropía conjunta entre las dos imágenes tiene su mínimo aproximadamente en la misma iteración en que se alcanza el mínimo error por el algoritmo (más visible para la MVL₁ y MVL₄), y después tiende a incrementarse al igual que el error. Por lo tanto, se puede establecer como criterio de parada para el método, el fijar un umbral para la magnitud de la derivada de la entropía conjunta y de esta forma se asegura la convergencia del algoritmo.

Conclusiones

Después de analizar los resultados obtenidos, cualitativos y cuantitativos, de los experimentos realizados para evaluar el algoritmo propuesto, se puede concluir que este método ofrece un buen desempeño para el registro elástico multimodal empleando como medida de variabilidad local la entropía, obteniendo en la evaluación cuantitativa un error medio menor a 1 mm. Además, ofrece como ventajas un cálculo numérico más rápido y con una mayor libertad en las deformaciones que se puedan lograr, debido a que el campo vectorial de las deformaciones no es parametrizado. Sin embargo, es deseable buscar mejoras para el método, intentando igualar los errores obtenidos por la versión unimodal, con un error medio menor a 0.5 mm. Dado que los experimentos fueron realizados empleando imágenes médicas simuladas y reales, la técnica propuesta puede ser considerada como una alternativa viable para su aplicación en la imagenología médica, aunque por otro lado es necesaria una evaluación de acuerdo a los estándares médicos dependiendo de la aplicación hacia la que se decida enfocar el método.

Como trabajo futuro, se propone realizar modificaciones al algoritmo para reducir el error en la estimación del campo vectorial e intentar que sea similar al reportado por técnicas de registro unimodal. También se contempla una implementación del algoritmo para el registro de volúmenes en la herramienta de análisis de imágenes médicas ITK (Insight Segmentation and Registration Toolkit [44]), además de realizar una comparación exhaustiva, de acuerdo a estándares médicos y a una aplicación específica, con otros métodos de registro elástico multimodal reportados en la literatura. Finalmente, buscamos realizar una evaluación del método propuesto para una aplicación médica específica.

Agradecimientos

Los autores agradecen al Institute of Molecular Bioimaging (IBFM), Milan, Italia, por proporcionar dos de las imágenes utilizadas en este trabajo y en especial a Aldo Mejía Rodríguez. Los autores agradecen a CONACyT por el apoyo financiero otorgado a través de los proyectos #168140 y #154146. I. Reducindo agradece a CONACyT por el apoyo financiero brindado a través de una beca de posgrado (CVU/Becario 266553/218513).

Referencias

1. Hill D. L. G., Batchelor P. G., Holden M., y Hawkes D. J., “Medical image registration,” Physics in Medicine and Biology, vol. 46, no. 3, 2001.
2. Pluim J. P. y Fitzpatrick J. M., “Image registration,” IEEE Transactions on Medical Imaging, vol. 22, pp. 1341–1343, Noviembre 2003.
3. Rueckert D. y Aljabar P., “Nonrigid registration of medical images: Theory, methods, and applications,” IEEE Signal Processing Magazine, vol. 27, pp. 113–119, Julio 2010.
4. Fischer B. y Modersitzki J., “Ill–posed medicine–an introduction to image registration,” Inverse Problems, vol. 24, Mayo 2008.
5. Zitová B. y Flusser J., “Image registration methods: a survey,” Image and vision computing, vol. 21, pp. 977–1000, 2003.
6. Nocedal J. y Wright S. J., Numerical Optimization. Springer, 2 ed., 2006.
7. Pluim J., Maintz J. y Viergever M., “Image registration by maximization of combined mutual information and gradient information,” Medical Imaging, IEEE Transactions on, vol. 19, pp. 809–814, Agosto 2000.
8. Viola P. A., Wells III W. M., Atsumi H., Nakajima S., y Kikinis R., “Multi–modal volume registration by maximization of mutual information,” Medical Image Analysis, vol. 1, no. 1, pp. 35–51, 1996.
9. Maes F., Collignon A., Vandermeulen D., Marchal G. y Suetens P., “Multimodality image registration by maximization of mutual information,” IEEE Transactions on Medical Imaging, vol. 16, pp. 187–198, Abril 1997.
10. Pluim J. P. W., Maintz J. B. A. y Viergever M. A., “Mutual-information-based registration of medical images: a survey,” IEEE Transactions on Medical Imaging, vol. 22, pp. 986–1004, Agosto 2003.
11. Gonzalez R. C., Woods R. E., y Eddins S. L., Digital Image Processing Using MATLAB. Gatesmark Publishing, 2 ed., 2009.
12. Man K., Tang K., y Kwong S., Genetic Algorithms.Springer, 3 ed., 2001.
13. Arulampalam M. S., Maskell S., Gordon N. y Clapp T., “A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking,” IEEE Transactions on Signal Processing, vol. 50, pp. 174–188, Febrero 2002.
14. Das A. y Bhattacharya M., “Affine-based registratiom of ct and mr modality images of human brain using multiresolution approaches: comparative study on genetic algorithm and particle swarm optimization,” Neural Computing and Applications, Mayo 2010.
15. Arce-Santana E. R., Campos-Delgado D. U. y Alba A., “Affine image registration guided by particle filter,” IET Image Processing, vol. 6, no. 5, pp. 455–462, 2002.
16. Reducindo I., “Registro rígido de imágenes guiado por filtro de partículas,” Tesis de Maestría, Universidad Autónoma de San Luis Potosí, Diciembre 2010.
17. Klein A., Andersson J., Ardekani B. A., Ashburner J., et al., “Evaluation of 14 nonlinear deformation algorithms applied to human brain mri registration,” Neuroimage, vol. 46, pp. 786–802, Julio 2009.
18. Xuan J., Wang Y., Freedman M. T., Adali T., y Shields P., “Nonrigid medical image registration by finite-element deformable sheet-curve models,” International Journal of Biomedical Imaging, pp. 1–9, 2006.
19. Fei W. y Vemuri B. C., “Non-rigid multi-modal image registration using cross-cumulative residual entropy,” International Journal of Computer Vision, vol. 74, no. 2, pp. 201–215, 2007.
20. Wang X. y Feng D. D., “Automatic elastic medical image registration based on image intensity,” International Journal of Image and Graphics, vol. 5, no. 2, pp. 351–369, 2005.
21. Serifović-Trbalić A., Demirović D., Prljaca N., Szekely G., y Cattin P. C., “Intensity-based elastic registration incorporating and isotropic landmark erros and rotational information,” International Journal of Computer Assisted Radiology and Surgery, vol. 4, pp. 463–468, Junio 2009.
22. Lange T., Papenberg N., Heldmann S., Modersitzki J., et al., “3d ultrasound-ct registration of the liver using combined landmark-intensity information,” International Journal of Computer Assisted Radiology and Surgery, vol. 4, pp. 79–88, 2008.
23. Joshi A., Shattuck D., Thompson P., y Leahy R., “Brain image registration using cortically constrained harmonic mappings,” in Information Processing in Medical Imaging, Lecture Notes in Computer Science, pp. 359–371, Springer–Verlag, 2007.
24. Baker S., Scharstein D., Lewis J., Roth S., Black M., y Szeliski R., “A database and evaluation methodology for optical flow,” International Journal of Computer Vision, vol. 92, pp. 1–31, 2011.
25. Arce-Santana E. R., Campos-Delgado D. U., y Alba A., “A non-rigid multimodal image registration method based on particle filter and optical flow,” in Advances in Visual Computing, vol. 6453 of Lecture Notes in Computer Science, pp. 35–44, Springer–Verlag, 2010.
26. Reducindo I., Arce-Santana E. R., Campos-Delgado D. U., y Alba A., “Evaluation of multimodal medical image registration based on particle filter,” Int. Conf. on Electrical Eng., Computing Science and Automatic Control, Septiembre 2010.
27. Mejia-Rodriguez A., Arce-Santana E., Scalco E., Tresoldi D., Mendez M., Bianchi A., Cattaneo G. y Rizzo G., “Elastic registration based on particle filter in radiotherapy images with brain deformations,” in Engineering in Medicine and Biology Society,EMBC, 2011 Annual International Conference of the IEEE, pp. 8049–8052, Septiembre 2011.
28. ter Haar Romeny B., Florack L., Koenderink J. y Viergever M., “Scale space: Its natural operators and differential invariants,” in Information Processing in Medical Imaging, vol. 511 of Lecture Notes in Computer Science, pp. 239–255, Springer Berlin / Heidelberg, 1991.
29. Simon D., Optimal State Estimation.Wiley–Interscience, 1 ed., 2006.
30. Studholme C., Hill D. L. G., y Hawkes D. J., “An overlap invariant entropy measure of 3d medical image alignment,” Pattern Recognition, vol. 32, pp. 71–86, Enero 1999.
31. Reducindo I., Arce-Santana E. R., Campos-Delgado D. U., Alba A. y Vigueraz-Gómez J. F., “An exploration of multimodal similarity metrics for parametric image registration based on particle filtering,” Int. Conf. on Electrical Eng., Computing Science and Automatic Control, Octubre 2011.
32. Trees H. L. V., Detection, Estimation and Modulation Theory: Part 1.John Wiley & Sons, 2001.
33. Bruhn A., Weickert J. y Schnörr C., “Lucas/kanade meets horn/schunck: Combining local and global optic flow methods,” International Journal of Computer Vision, vol. 61, pp. 211–231, 2005.
34. Lucas B. D. y Kanade T., “An iterative image registration technique with an application to stereo vision,” Proceedings of Imaging Understanding Workshop, pp. 121–130, 1981.
35. Horn B. K. y Schunck B. G., “Determining optical flow,” Artificial Intelligence, vol. 17, 1981.
36. Biswal P. C., Numerical Analysis.Prentice-Hall of India Pvt. Ltd., 2008.
37. Zikic D., Kamen A. y Navab N., “Revisiting horn and schunck: Interpretation as Gauss-Newton optimisation,” in In Proceedings of BMVC’2010, 2010.
38. Shannon C. E., “A mathematical theory of communication,” The Bell System Technical Journal, vol. 27, pp. 379–423, 623–656, Julio 1948.
39. Gefen S., Kiryati N. y Nissanov J., “Atlas-based indexing of brain sections via 2-d to 3-d image registration,” Biomedical Engineering, IEEE Transactions on, vol. 55, pp. 147 –156, Enero 2008.
40. Periaswamy S., General–Purpose Medical Image Registration. Tesis de Doctorado, Dartmouth College, Abril 2003.
41. Kwan R.-S., Evans A. y Pike G., “Mri simulation-based evaluation of image-processing and classification methods,” Medical Imaging, IEEE Transactions on, vol. 18, pp. 1085 –1097, Noviembre 1999.
42. Bookstein F., “Principal warps: thin-plate splines and the decomposition of deformations,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 11, no. 6, 1989.
43. Vercauteren T., Pennec X., Perchant A. y Ayache N., “Diffeomorphic demons: efficient non-parametric image registration” Neuroimage, vol. 45, pp. 61–72, Marzo 2009.
44. Ibanez L., Schroeder W., Ng L. y Cates J., The ITK Software Guide. Kitware, Inc. ISBN 1-930934-15-7, http://www.itk.org/ItkSoftwareGuide.pdf, second ed., 2005.