Repositorio Institucional de Documentos

Trabajo Fin de Máster DESARROLLO Y EVALUACIÓN DE UNA HERRAMIENTA DE REGISTRO DIFEOMÓRFICO POR LANDMARKS 3D DE ESTRUCTURAS CEREBRALES EN IMÁGENES MRI Autor David Bretón Domínguez Director Salvador Olmos Gassó Escuela de Ingeniería y Arquitectura (EINA) 2011 - 2012 Desarrollo y evaluación de una herramienta de registro difeomórco por landmarks 3D de estructuras cerebrales en imágenes MRI RESUMEN En la disciplina de Anatomía Computacional se requiere una herramienta de registro o normalización espacial de imágenes que sea capaz de modelar grandes deformaciones (ya que existe gran variabilidad anatómica entre los cerebros humanos, especialmente en la corteza) garantizando suavidad en las transformaciones espaciales. Con este n, recientemente se ha propuesto el paradigma de los difeomorsmos, tanto en imágenes volumétricas como en registro de puntos anatómicos o landmarks. Este proyecto versa sobre el registro difeomórco de estructuras cerebrales mediante landmarks 3D a partir de imágenes de resonancia magnética. El objetivo nal consiste en desarrollar y evaluar las prestaciones de una herramienta que proporcione una transformación espacial que consiga una correspondencia entre los surcos y giros más importantes de la corteza cerebral. En concreto, se trata de un problema de registro no rígido basado en una transformación difeomórca entre los puntos de un atlas cerebral (imagen y sus correspondientes etiquetas de segmentación) y los puntos de unas imágenes objetivo. Dicha transformación se ha utilizado además como inicialización de una herramienta ya existente de registro volumétrico de todo el cerebro, con el objetivo de dotarla de mayor robustez. La evaluación de prestaciones se ha realizado sobre un conjunto de imágenes MRI provenientes del estudio Alzheimer's Disease Neuroimaging Initiative (ADNI). Los puntos anatómicos se han obtenido mediante la herramienta BrainVisa. Aunque estos landmarks son realmente tridimensionales, a lo largo de esta memoria y por simplicidad en la visualización de los resultados, se presentan en dos dimensiones para un corte axial representativo. Los créditos prácticos se han realizado en el Servicio de Radiodiagnóstico del Grupo Hospitalario Quirón, Zaragoza, bajo la tutela del Jefe de Servicio, el Dr. Nicolás Fayed Miguel. El objetivo fue conocer el equipamiento de radiología (resonancia magnética nuclear, tomografía computerizada, ecografías, rayos X, etc.) así como el sistema PACS de transferencia, gestión y almacenamiento de imágenes médicas. Índice general 1. Introducción 1 1.1. Interésdelproyecto................................... 1 1.2. Estadodelarte..................................... 2 1.3. Objetivodelproyecto ................................. 3 1.4. Contenidosdelamemoria............................... 3 2. Registro difeomórco de landmarks 5 2.1. Planteamiento del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2. Registrodifeomórco.................................. 6 2.3. Función de energía y gradiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3. Resultados 9 3.1. Implementación numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1.1. Campo de velocidades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1.2. Condiciones de contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1.3. Integración del campo de velocidades . . . . . . . . . . . . . . . . . . . . . 11 3.1.4. Operador L y parámetro de suavizado (α) .................. 11 3.1.5. Parámetro de equilibrio (λ) .......................... 12 3.1.6. Elección del campo inicial v0(x) ....................... 13 3.1.7. Estrategia de optimización . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2. Resultados en imagen médica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2.1. Material y procedimientos . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2.2. Resultados ................................... 18 3.2.3. Discusión .................................... 20 4. Conclusiones y líneas futuras 21 4.1. Resumen del proyecto y conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.2. Líneasfuturas...................................... 21 Bibliografía 22 A. Cálculo de la derivada de un campo vectorial respecto a una perturbación 27 B. Fórmula BakerCampbellHausdor (BCH) 31 C. Parte práctica del TFM 33 i ii Índice general Índice de guras 2.1. Modelado de un difeomorsmo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.1. Efecto de las condiciones de contorno. . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2. Efecto del dominio de denición. . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.3. Efecto de la discretización del campo de velocidades. (a) Valor mínimo normalizado de la energía en función del paso de discretización. (b) Tiempo de ejecución en función del paso de discretización. . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.4. Efecto del parámetro de regularización en el registro por landmarks. Para cada valor de α se muestra ϕv y el entorno de deformación asociado. . . . . . . . . . . 12 3.5. Efecto del parámetro de equilibrio λ en el registro por landmarks. Los números junto a la gráca indican el valor de λ . Los valores de los términos de ajuste y de regularización han sido normalizados al intervalo [0,1]. . . . . . . . . . . . . . . . 13 3.6. Campo vectorial asociado a una transformación afín. La orientación y el tamaño de las echas representan la dirección y el módulo del campo. Puntos p (aspas), T · p (cuadrados) y q (círculos). ............................ 14 3.7. (a) Velocidades que afectan a la trayectoria. (b) Campo completado minimizando la energía de regularización. Nótese la continuidad en los bordes a diferencia de la gura3.6......................................... 15 3.8. BCHdecampos. .................................... 16 3.9. (a) Campo nal optimizado. (b) Proceso de optimización. En rojo, energía. En azul,tiempo. ...................................... 16 3.10. (a),(b): Imagen fuente y objetivo, respectivamente. Los contornos en blanco delimitan las estructuras. Cada círculo numerado representa un landmark. (c) Difeomorsmo de landmarks. Puntos origen (aspas) y destino (círculos). (d) Relación de las estructuras cerebrales asociadas a los landmarks. . . . . . . . 17 3.11. Distancia media entre landmarks calculada sobre los 120 sujetos. . . . . . . . . . 18 3.12. (a) Distribución del promedio espacial de la distancia entre landmarks. (b) Distribución del solapamiento de regiones corticales. . . . . . . . . . . . . . . 19 3.13. Mapa de intersección en las regiones corticales (valores desde 1, no hay intersección, hasta 120, intersección de todos los sujetos). . . . . . . . . . . . . . 19 3.14. Distribución del logaritmo del jacobiano. . . . . . . . . . . . . . . . . . . . . . . . 20 C.1.Workstation. ...................................... 33 C.2. Adquisición de Imagen de Resonancia Magnética (MRI). . . . . . . . . . . . . . . 34 C.3. Equipo conectado al PACS desde el que el Dr. Fayed realiza las labores de diagnósticoporimagen. ................................ 35 iii iv Índice de guras Capítulo 1 Introducción 1.1. Interés del proyecto El registro de imágenes es el proceso de poner en correspondencia espacial dos o más imágenes de un mismo objeto, tomadas en diferentes momentos, desde diferentes puntos de vista o mediante distintos sensores, de manera que queden correctamente alineadas geométricamente. Consiste principalmente en establecer una correspondencia biunívoca entre los puntos de las diferentes imágenes, donde esta correspondencia puede considerarse como una transformación entre una imagen fuente y otra objetivo. Para calcular dicha transformación, normalmente se plantea un problema de optimización, en el que se busca, dentro de un espacio de transformaciones posibles, aquella que minimice una métrica establecida. En imagen médica, el registro se utiliza normalmente para alinear imágenes del mismo sujeto adquiridas con diferentes técnicas (MRI, PET, etc.), alinear imágenes del mismo tipo pero tomadas en diferentes instantes temporales (detección de cambios o el control de un tumor), o bien imágenes de diferentes sujetos que se desean comparar. A menudo, implican un registro no rígido para hacer frente a la deformación (debido a la respiración, a los cambios anatómicos, etc.) o a la variabilidad anatómica entre distintos sujetos. En el caso de neuroimagen, el interés del registro de imagen viene dado por su aplicación en la comparación de una misma estructura cerebral entre individuos diferentes, y en el desarrollo de una teoría estadística que permita estudiar la forma de dichas estructuras en diferentes poblaciones. Esta rama de la investigación médica, que se conoce como Anatomía Computacional, tiene su base en las diferencias que existen entre las estructuras cerebrales de diferentes grupos naturales, o bien entre una población control y otra afectada por enfermedad, fármaco, etc. En este sentido, se puede aprender mucho de una enfermedad estudiando las estructuras a las que afecta, y en último término, diagnosticar o caracterizar el estado de una enfermedad por la forma de una estructura anatómica en particular. Por ejemplo, en las fases tempranas de la enfermedad de Alzheimer ya se pueden apreciar atroas signicativas en la corteza entorrinal y en el giro hipocámpico [Braak and Braak, 1995], mientras que la disminución progresiva del estado cognitivo está relacionada con atroa de regiones del hipocampo y la amígdala [Bossa et al., 2010]. El registro de imagen también se utiliza en el análisis funcional de la corteza cerebral. La hipótesis de que el patrón de pliegues de la corteza está relacionado con la arquitectura neuronal subyacente y con la organización funcional ha sido rearmada por una serie de estudios recientes [Regis et al., 2005, Fischl et al., 2008]. En este caso, resulta crucial alinear los pliegues anatómicamente homólogos, tarea nada sencilla debido a la alta variabilidad de forma y topología entre los valles y surcos de la corteza de diferentes individuos. Así pues, el registro de imagen se considera una herramienta vital en el análisis de imágenes médicas cerebrales, ya que permite registrar las estructuras internas y alinear correctamente los pliegues corticales homólogos, asegurando así la consistencia y la sensibilidad de todas las medidas anatómicas y funcionales posteriores. 1 8 2.3. Función de energía y gradiente Para reducir el valor de esta energía, se sigue un proceso de optimización basado en un método de descenso. Este tipo de métodos están basados en algoritmos iterativos que suelen utilizar información sobre el gradiente de la función, que en este caso se calcula como: ∂vE(ϕv, p, q) = λ∂v(E1(ϕv, p, q)) + ∂v(E2(ϕv)) (2.6) Desarrollando el gradiente del término de ajuste: ∂v(E1(ϕv, p, q)) = −2 n X i=1 (qi−ϕv(pi)) ∂ϕv(pi) ∂v (2.7) Para poder resolver este gradiente primero es necesario calcular como varía ϕv al perturbar el campo de velocidades v en la dirección h . De acuerdo con [Beg et al., 2005] esta variación viene dada por: ∂hϕv= l´ım →0 ϕv+h −ϕv =Dϕvˆ1 0 (Dϕuv)−1h(ϕuv)du (2.8) donde Dϕ es la matriz jacobiana del difeomorsmo. Tanto la demostración de la expresión (2.8) como el cálculo de Dϕ pueden verse en el Apéndice A. De esta manera, conocida la variación de ϕv ante una perturbación cualquiera, se puede calcular ∂ϕv(pi) ∂v aplicando (2.8) respecto a las perturbaciones h que afectan a ϕv(pi) (es decir, las perturbaciones h que afectan a la trayectoria impuesta por el campo de velocidades v en el punto pi ) ∂ϕv(pi) ∂v =∂ϕv(pi) ∂h h(ϕv(pi)) =Dϕv(pi)ˆ1 0 (Dϕuv (pi))−1hu(ϕuv (pi)) du (2.9) y sustituyendo en la expresión (2.7): ∂v(E1(ϕv, p, q)) = −2 n X i=1 (qi−ϕv(pi)) Dϕv(pi)ˆ1 0 (Dϕuv (pi))−1hu(ϕuv (pi)) du (2.10) Por otro lado, el gradiente del término de suavidad es : ∂v(E2(ϕv)) = ∂v(hLv, Lvi) = ∂vL†Lv, v= 2L†Lv (2.11) donde L† el operador adjunto de L . Juntando las dos partes del gradiente se obtiene: ∂vE(ϕv, p, q) =−2λ n X i=1 (qi−ϕv(pi)) Dϕv(pi)ˆ1 0 (Dϕuv (pi))−1h(ϕuv (pi)) du + 2L†Lv (2.12) Este gradiente tiene componentes de alta frecuencia, por lo que resulta recomendable ltrarlo por ejemplo con K=L†L−1, resultando un gradiente nal: ∂vE(ϕv, p, q) = −2Kλ n X i=1 (qi−ϕv(pi)) Dϕv(pi)ˆ1 0 (Dϕuv (pi))−1h(ϕuv (pi)) du + 2v (2.13) Capítulo 3 Resultados 3.1. Implementación numérica 3.1.1. Campo de velocidades Aunque la herramienta de registro se ha desarrollado y evaluado para un conjunto de landmarks 3D, por simplicidad en las expresiones y en la visualización de las guras, se presentan aquí los resultados obtenidos para un problema 2D. Así pues, nuestro problema de registro parte de dos conjuntos de landmarks {pi∈Ω|i= 1,2, . . . , n} y {qi∈Ω|i= 1,2, . . . , n} , con Ω⊆R2 y un campo de velocidades v denido sobre una rejilla estructurada {xi}a i=1 × {yj}b j=1 , donde para cada uno de los puntos (xi, yj)∈Ω se conoce el valor de las componentes vx y vy . Es necesario que el dominio de denición Ω sea lo sucientemente amplio para que, además de englobar a los puntos origen y objetivo, el campo de velocidades sea capaz de decaer hasta valores sucientemente pequeños en los extremos del dominio. Sin embargo, hay que tener en cuenta que un dominio demasiado grande conlleva medir la energía de suavizado E2(ϕv) en un número elevado de puntos, lo que genera un gasto computacional adicional. En cuanto a la discretización del campo de velocidades, cuanto más pequeños sean los pasos de discretización (xi+1 −xi) y (yi+1 −yi) , mayor número de velocidades guiarán al difeomorsmo. Al aumentar los grados de libertad, el registro es capaz de alcanzar soluciones con mejor conguración energética, pero a cambio, se eleva considerablemente el tiempo de optimización. 3.1.2. Condiciones de contorno Para calcular el valor de v en cualquier otro punto del dominio distinto de la rejilla sobre la que está denido, se dene una regla de interpolación, que por simplicidad se considera lineal. Si se desea calcular el valor del campo en un punto que no pertenece a Ω, entonces es necesario establecer unas condiciones de frontera, que pueden ser: Valor jo: El valor del campo en el exterior de Ω es igual a un valor arbitrario vout . Simétrico: Fuera del dominio, el valor del campo se considera simétrico al del interior . Más cercano: El valor del campo en un punto pf en el exterior de Ω toma el valor del campo en el punto pi del interior más cercano a pf . Circular: El campo de velocidades es periódico en la extensión de su dominio. Decaimiento: El valor del campo decae linealmente desde el valor en los extremos de Ω hasta cero, en función de un parámetro que regula la pendiente del decaimiento. En la gura 3.1 puede verse el efecto de las condiciones de contorno en el perl de un campo unidimensional con Ω = {x∈R|0≤x≤30} . 9 10 3.1. Implementación numérica −30 −15 0 15 30 45 60 0 0.5 1 1.5 (a) Sin condición de contorno −30 −15 0 15 30 45 60 0 0.5 1 1.5 (b) Valor jo vout =0.1 −30 −15 0 15 30 45 60 0 0.5 1 1.5 (c) Simétrico −30 −15 0 15 30 45 60 0 0.5 1 1.5 (d) Más cercano −30 −15 0 15 30 45 60 0 0.5 1 1.5 (e) Circular −30 −15 0 15 30 45 60 0 0.5 1 1.5 (f) Decaimiento Figura 3.1: Efecto de las condiciones de contorno. Para garantizar la invertibilidad del difeomorsmo, el campo de deformaciones ϕv(x) debe ser suave, lo que implica también la suavidad del campo de velocidades v(x) . Si se utiliza la condición de frontera circular se asegura además un contorno periódico, lo que unido a la suavidad de v(x) permite poder calcular correctamente la derivada del campo vectorial en cualquier punto del dominio. En la gura 3.2 se puede ver el efecto que tiene el dominio del campo de velocidades en un ejemplo sintético de registro de landmarks con condiciones de contorno periódicas. Dados cuatro puntos origen {pi∈Ω|i= 1,2,3,4} y destino {qi∈Ω|i= 1,2,3,4} , ambos con Ω⊆R2 , la gura muestra la deformación que sufre el espacio por la acción del campo de velocidades. Como se puede observar, para dominios demasiado ajustados a los datos, las velocidades en un extremo del dominio afectan a las velocidades del extremo contrario. Figura 3.2: Efecto del dominio de denición. Capítulo 3. Resultados 11 Del mismo modo, se puede observar como afecta la discretización del campo en el valor mínimo de la energía (gura 3.3a) y en el tiempo del proceso de optimización (gura 3.3b). 0.01 0.1 1 10 0.1 0.2 0.5 1 Paso (ud) Energia normalizada (a) 10−2 10−1 100101102 100 101 102 103 104 Paso (ud) Tiempo (s) (b) Figura 3.3: Efecto de la discretización del campo de velocidades. (a) Valor mínimo normalizado de la energía en función del paso de discretización. (b) Tiempo de ejecución en función del paso de discretización. 3.1.3. Integración del campo de velocidades Dado un campo de velocidades y unos puntos iniciales pi , es necesario calcular las posiciones ϕtv (pi) que denen la trayectoria del punto pi a lo largo del tiempo, las cuales están determinadas por el difeomorsmo solución a la ecuación (2.2). Siguiendo el esquema de integración explícito de Euler, las nuevas posiciones se calculan como: ϕ(t+k)v(pi) = ϕtv (pi) + k · vϕtv (pi) (3.1) Este método de integración tiene la ventaja de ser muy simple de implementar, aunque es necesario ajustar empíricamente el paso k elegido para la integración. 3.1.4. Operador L y parámetro de suavizado (α) El operador L fue elegido como L=Id −α∆ , donde ∆ es el operador Laplaciano 4=∂2 ∂x ² +∂2 ∂y ² +∂2 ∂z ² y el parámetro α penaliza las derivadas de segundo orden del campo de velocidades ([Beg et al., 2005]). Concretamente se utiliza un operador Laplaciano con esquema de diferencias centradas que, asumiendo condiciones de contorno periódicas, es autoadjunto ( L†=L ) y tanto Lv como L†Lv pueden ser calculados en el dominio de Fourier. Para pequeños valores de α , L se asemeja a la identidad, y la energía de regularización penaliza simplemente valores altos del campo de velocidades E2(ϕv) = ´ΩkLv(x)k2dx ≃´Ωkv(x)k2dx . Esto conduce a que el campo óptimo tiene la mayoría de velocidades nulas, pero permite variar abruptamente sus valores para poder ajustar los puntos. En cambio, si α es grande, en el operador L predomina el laplaciano, la energía de regularización penaliza las derivadas de segundo orden del campo de velocidades, y el campo de velocidades resultante es suave. En la gura 3.4 se puede ver este efecto para diferentes valores del parámetro de suavizado α , además de la trayectoria que siguen los puntos p en su deformación hasta registrarse con los puntos q . 12 3.1. Implementación numérica (a) α= 0,1 (b) α= 1 (c) α= 5 (d) α= 20 Figura 3.4: Efecto del parámetro de regularización en el registro por landmarks. Para cada valor de α se muestra ϕv y el entorno de deformación asociado. Estos resultados pueden interpretarse también desde otro punto de vista. De acuerdo con [Marsland and Twining, 2004], para minimizar la energía de regularización la mejor forma de mover un punto diferencialmente, es deformar un entorno de dicho punto. El tamaño de este entorno depende del operador L escogido. Como se puede observar en la gura 3.4, para valores pequeños de α , el tamaño del entorno que se deforma es pequeño; para valores mayores de α , se deforma un entorno mayor, que puede llegar a englobar varios landmarks. 3.1.5. Parámetro de equilibrio (λ) En la función de energía, el parámetro de equilibrio λ pondera la importancia relativa entre el ajuste de los landmarks y la suavidad del campo de velocidades (expresión (2.4)). Un valor alto de λ da mayor importancia al término de ajuste sobre el de suavizado y conduce a conguraciones en las que los landmarks origen son deformados hasta coincidir perfectamente con los landmarks objetivo, aunque el campo de velocidades resultante puede no ser suave. Por el contrario, valores pequeños de λ preponderan el término de regularización y llevan a campos de velocidades suaves, pero que no logran ajustar completamente los landmarks. Hay que destacar que los valores de λ que logran un ajuste de landmarks aceptable dependen de cada problema concreto (número de landmarks, dominio y discretización del campo de velocidades). Como criterio, se elige el menor λ que devuelva un ajuste de landmarks prácticamente exacto. Por ejemplo, para el caso de la gura 3.5, el valor escogido para el parámetro de equilibrio sería λ= 10000 . Capítulo 3. Resultados 13 Figura 3.5: Efecto del parámetro de equilibrio λ en el registro por landmarks. Los números junto a la gráca indican el valor de λ . Los valores de los términos de ajuste y de regularización han sido normalizados al intervalo [0,1]. 3.1.6. Elección del campo inicial v0(x) El proceso de optimización necesita un valor inicial a partir del cual iniciar el descenso. A continuación se muestran diferentes estrategias para establecer el campo de velocidades inicial v0(x) . 1. Campo nulo Consiste sencillamente en establecer v0(x) = 0 . La conguración energética inicial es fácil de determinar en esta situación. Si v= 0 en todo el dominio, entonces el término de energía de regularización o suavizado E2(ϕv) = ´ΩkLv(x)k2dx es también nulo. Por otro lado, ϕv(pi)|v=0=pi , y el término de ajuste E1(ϕv, p, q) = Pn i=1 kqi−ϕv(pi)k2=Pn i=1 kqi−pik2 es simplemente la suma de la norma al cuadrado de las distancias entre los puntos iniciales y los puntos nales. Esta conguración se corresponde con un punto singular de la función de energía, y dependiendo del valor del parámetro de equilibrio λ , el proceso de optimización puede detenerse desde el inicio en un mínimo local. Por otro lado, y aun cuando esto no ocurriese, iniciar la optimización desde v0(x) = 0 conlleva una convergencia lenta, ya que la optimización parte de una conguración inicial carente de información. Así pues, aunque el campo inicial nulo es el más sencillo de todos, su conguración energética asociada no es la más apropiada, lo que hace necesario el diseño de otras conguraciones de campos iniciales. 14 3.1. Implementación numérica 2. Campo afín Sean p=  px 1px 2. . . px n py 1py 2. . . px n 1 1 . . . 1  y q=  qx 1qx 2. . . qx n qy 1qy 2. . . qx n 1 1 . . . 1  las coordenadas homogéneas en 2D de los puntos origen y destino respectivamente. Se puede calcular la matriz de la transformación afín T que minimiza kq−T · pk2 . kq−T · pk2= tr (q−T · p)T(q−T · p) = tr qTq−2 tr TpqT+ tr TppTTT (3.2) Derivando esta expresión respecto a T e igualando a cero se obtiene: T=qpTppT+ (3.3) donde ( · )+ denota pseudoinversa. Esta transformación existe para cualquier p y q , y siempre que exista su logaritmo, se podrá utilizar el campo vectorial asociado a dicha transformación como campo inicial para nuestra optimización. Figura 3.6: Campo vectorial asociado a una transformación afín. La orientación y el tamaño de las echas representan la dirección y el módulo del campo. Puntos p (aspas), T · p (cuadrados) y q (círculos). Sea vT un campo vectorial estacionario que se pueda expresar de forma lineal como vaf (x) = A · x∀x∈Ω, (3.4) donde v y x están expresadas en coordenadas homogéneas. Entonces la solución a la ecuación (2.2) se demuestra fácilmente que es ϕtv (x) = eAt · ϕ0(x) = eAt · x donde eAt es la exponencial de matrices y ϕtv (x) está expresado también en coordenadas homogéneas. Particularmente, para t= 1 se obtiene ϕv(x) = eA · x (3.5) A partir de esta expresión, y teniendo en cuenta que los puntos T·p son la transformación de los puntos p , identicando términos se obtiene que eA=T , y por tanto A= logm T . Sustituyendo el valor de A en la ecuación (3.4), el campo vectorial asociado a la transformación afín queda de la forma: vaf (x)=( logm T) · x (3.6) Capítulo 3. Resultados 15 En la gura 3.6 se puede observar que el campo afín obtenido mediante la expresión 3.6 no cumple la condición de continuidad y periodicidad en los bordes establecida en la subsección 3.1.2, por lo que el campo no se podría derivar correctamente en todo su dominio. Para solucionarlo, lo que se hace a continuación es ltrar dicho campo. Para ello, se aprovecha la medida de regularización denida en la expresión (2.5). Preservando el valor de las velocidades que afectan directamente a la trayectoria de los puntos (gura 3.7a), se completa el resto del dominio con velocidades que minimizan dicho término de regularización (gura 3.6). (a) (b) Figura 3.7: (a) Velocidades que afectan a la trayectoria. (b) Campo completado minimizando la energía de regularización. Nótese la continuidad en los bordes a diferencia de la gura 3.6. 3. BCH de campos Como se puede comprobar en la gura 3.7b, la transformación afín no asegura que los puntos T·p coincidan con los puntos q . Una transformación afín implica conservar las líneas paralelas, y por tanto es insuciente para poder proyectar cualquier conjunto de n puntos en otro. Sería interesante por tanto, poder deformar a continuación los puntos ya transformados (que estarán muy próximos a los puntos objetivo) de manera que se aproximen todavía más a los puntos objetivo. Una posible solución es establecer una trayectoria lineal desde T·p hasta q y calcular el campo vlin en los puntos asociados a dicha trayectoria: vlin = (q−T·p) Es necesario entonces componer el campo vaf asociado a la transformación afín y el campo vlin asociado a la trayectoria lineal. Es decir, hay que determinar un campo de velocidades v=Z(vlin, vaf ) tal que exp (v)≈exp (vlin)◦exp (vaf ) (3.7) ϕv(x)≈(ϕvlin ◦ϕvaf ) (x) = ϕvlin (ϕvaf (x)) Esta composición de campos no es trivial. Para resolverla se utiliza la fórmula de BakerCampbellHausdor (BCH). Para más detalles, ver el Apéndice B. Como puede apreciarse en la gura 3.8, al combinar el campo afín inicial, con uno lineal, se produce una mejora apreciable en el ajuste de los landmarks. Sin embargo, este ajuste sigue sin ser del todo exacto, ya que la solución de una BCH de campos implica una sucesión de innitos términos, que se debe truncar en algún momento. En cualquier caso, estos pequeños desajustes se verán reducidos posteriormente durante el proceso de optimización. 16 3.1. Implementación numérica (a) Campo vectorial afín (ltrado) (b) Campo vectorial lineal (detalle) (c) Combinación de (a) y (b) Figura 3.8: BCH de campos. 3.1.7. Estrategia de optimización El proceso de optimización se lleva a cabo mediante un método de gradiente conjugado para funciones no lineales ([Hager and Zhang, 2005, Nocedal and Wright, 1999]). Este método genera una secuencia vk , k≥1 a partir del campo inicial v0 , utilizando la ley de recurrencia: vk=vk−1+γkdk donde el paso positivo γk se obtiene mediante line search . La dirección de búsqueda es una actualización de la dirección de descenso anterior con la dirección del gradiente negativo (la dirección opuesta al gradiente calculado según la expresión (2.13)) que depende de un parámetro de actualización βk . dk=βkdk−1−∂vE(p, q, ϕv) En la gura 3.9 se muestra la optimización del campo de la imagen 3.8c. (a) (b) Figura 3.9: (a) Campo nal optimizado. (b) Proceso de optimización. En rojo, energía. En azul, tiempo. Capítulo 3. Resultados 17 3.2. Resultados en imagen médica 3.2.1. Material y procedimientos La herramienta de registro basado en landmarks se ha evaluado en un conjunto de imágenes MRI obtenidas de sujetos reales. Para ello, se seleccionaron aleatoriamente 120 sujetos (40 sanos, 40 enfermos de Alzheimer y 40 con deterioro cognitivo leve) de la base de datos de ADNI ( Alzheimer's Disease Neuroimaging Initiative ). La imagen referencia fue una imagen MRI de un sujeto sano diferente a los anteriores. El primer paso fue el establecimiento de landmarks en todas las imágenes. Con la herramienta BrainVisa (http://brainvisa.info) se segmentaron y etiquetaron un total de 22 estructuras de la corteza cerebral. Para establecer los landmarks se calculó el punto anatómico centroide de cada estructura. Posteriormente estas posiciones fueron vericadas y corregidas manualmente (guras 3.10a y 3.10b). (a) (b) (c) 1. Ínsula Izquierda 2. Occipital Izquierdo 3. Surco Cerebral Lateral (posterior) Izquierdo 4. Surco Cerebral Lateral (anterior) Izquierdo 5. Surco Cerebral Lateral (ascendente) Izquierdo 6. Surco Calcarino (anterior) Izquierdo 7. Giro Frontal Inferior (anterior) Izquierdo 8. Giro Frontal Interno Izquierdo 9. Giro Temporal Inferior (posterior) Izquierdo 10. Giro Temporal Superior Izquierdo 11. Ínsula Derecha 12. Occipital Derecho 13. Surco Cerebral Lateral (posterior) Derecho 14. Surco Cerebral Lateral (ascendente) Derecho 15. Surco Calloso−Marginal (anterior) Derecho 16. Surco Calcarino (anterior) Derecho 17. Giro Frontal Inferior (anterior) Derecho 18. Giro Frontal Intermedio Derecho 19. Giro Frontal Marginal Derecho 20. Giro Intralingual (posterior) Derecho 21. Giro Temporal Inferior (posterior) Derecho 22. Giro Temporal Superior Derecho (d) Figure 3.10: (a),(b): Imagen fuente y objetivo, respectivamente. Los contornos en blanco delimitan las estructuras. Cada círculo numerado representa un landmark. (c) Difeomorsmo de landmarks. Puntos origen (aspas) y destino (círculos). (d) Relación de las estructuras cerebrales asociadas a los landmarks. 24 Bibliografía [Fischl et al., 2008] Fischl, B., Rajendran, N., Busa, E., Augustinack, J., Hinds, O., Yeo, B. T. T., Mohlberg, H., Amunts, K., and Zilles, K. (2008). Cortical folding patterns and predicting cytoarchitecture. Cerebral Cortex , 18(8):19731980. [Hager and Zhang, 2005] Hager, W. W. and Zhang, H. (2005). A survey of nonlinear conjugate gradient methods. Science , 2(0203270):3558. [Hernandez et al., 2009] Hernandez, M., Bossa, M. N., and Olmos, S. (2009). Registration of anatomical images using paths of dieomorphisms parameterized with stationary vector eld ows. International Journal of Computer Vision , 85(3):291306. [I. and Havel, 1994] I., N. and Havel, T. F. (1994). Derivatives of the matrix exponential and their computation. Technical Report TR-33-94, Biological Chemistry and Molecular Pharmacology, Harvard Medical School, Boston, MA 02115-5718. [Jaccard, 1901] Jaccard, P. (1901). Distribution de la ore alpine dans le bassin des dranses et dans quelques régions voisines. Bulletin de la Société Vaudoise des Sciences Naturelles , 37. [Joshi et al., 2007] Joshi, A. A., Shattuck, D. W., Thompson, P. M., and Leahy, R. M. (2007). Surface-constrained volumetric brain registration using harmonic mappings. IEEE TRANS. MED. IMAG , 26(12):16571669. [Liu et al., 2004] Liu, T., Shen, D., and Davatzikos, C. (2004). Deformable registration of cortical structures via hybrid volumetric and surface warping. NeuroImage , 22(4):1790801. [Lyttelton et al., 2007] Lyttelton, O., Boucher, M., Robbins, S., and Evans, A. (2007). An unbiased iterative group registration template for cortical surface analysis. NeuroImage , 34(4):15351544. [Marsland and Twining, 2004] Marsland, S. and Twining, C. J. (2004). Constructing dieomorphic representations for the groupwise analysis of nonrigid registrations of medical images. IEEE Trans. Med. Imaging , 23(8):10061020. [Miller et al., 2006] Miller, M. I., Trouvé, A., and Younes, L. (2006). Geodesic shooting for computational anatomy. J. Math. Imaging Vis. , 24(2):209228. [Nocedal and Wright, 1999] Nocedal, J. and Wright, S. J. (1999). Numerical Optimization , volume 43. Springer. [Pennanen et al., 2004] Pennanen, C., Kivipelto, M., Tuomainen, S., Hartikainen, P., Hanninen, T., Laakso, M. P., Hallikainen, M., and Vanhanen, M. (2004). Hippocampus and entorhinal cortex in mild cognitive impairment and early ad. Neurobiology of Aging , 25(3):303  310. [Postelnicu et al., 2009] Postelnicu, G., Zollei, L., and Fischl, B. (2009). Combined volumetric and surface registration. IEEE TRANS. MED. IMAG , 28(4):508522. [Qiu et al., 2009] Qiu, A., Fennema-Notestine, C., Dale, A. M., and Miller, M. I. (2009). Regional shape abnormalities in mild cognitive impairment and alzheimer's disease. NeuroImage , 45(3):656  661. [Regis et al., 2005] Regis, J., Mangin, J.-F., Ochiai, T., Frouin, V., Riviere, D., Cachia, A., Tamura, M., and Samson, Y. (2005). Sulcal root generic model: a hypothesis to overcome the variability of the human cortex folding patterns. Neurologia medicochirurgica , 45(1):117. [Rueckert et al., 2003] Rueckert, D., Frangi, A. F., and Schnabel, J. A. (2003). Automatic construction of 3-d statistical deformation models of the brain using nonrigid registration. IEEE TRANS. MED. IMAG , 22(8):10141025. [Shen and Davatzikos, 2002] Shen, D. and Davatzikos, C. (2002). Hammer: hierarchical attribute matching mechanism for elastic registration. Bibliografía 25 [Styner et al., 2003] Styner, M., Gerig, G., Lieberman, J., Jones, D., and Weinberger, D. (2003). Statistical shape analysis of neuroanatomical structures based on medial models. Medical Image Analysis , 7(3):207  220. Functional Imaging and Modeling of the Heart. [Subsol et al., 1997] Subsol, G., Roberts, N., Doran, M., Thirion, J. P., and Whitehouse, G. H. (1997). Automatic analysis of cerebral atrophy. Magnetic Resonance Imaging , 15(8):917927. [Tosun and Prince, 2008] Tosun, D. and Prince, J. L. (2008). A geometry-driven optical ow warping for spatial normalization of cortical surfaces. IEEE TRANS. MED. IMAG , 27(12):17391753. [Van Essen and Dierker, 2007] Van Essen, D. C. and Dierker, D. L. (2007). Surface-based and probabilistic atlases of primate cerebral cortex. Neuron , 56(2):209225. [Vercauteren et al., 2009] Vercauteren, T., Pennec, X., Perchant, A., and Ayache, N. (2009). Dieomorphic demons: Ecient non-parametric image registration. NeuroImage , 45(1, Supplement 1):S61  S72. [Wang et al., 2007] Wang, L., Beg, F., Ratnanather, T., Ceritoglu, C., Younes, L., Morris, J. C., and Csernansky, J. G. (2007). Large deformation dieomorphism and momentum based hippocampal shape discrimination in dementia of the alzheimer type. IEEE TRANS. MED. IMAG , 26(4):462470. 26 Bibliografía Apéndice A Cálculo de la derivada de un campo vectorial respecto a una perturbación La variación de ϕtv cuando v es perturbada con h , es de la forma: ∂hϕtv = l´ım →0 ϕt(v+h)−ϕtv =Dϕtv ˆt 0 (Dϕuv)−1h(ϕuv)du (A.1) Demostración: Asumiendo que en (A.1) la derivada con respecto a  existe (la prueba de existencia puede realizarse mediante ecuaciones diferenciales ordinarias), se procede a su identicación. Partiendo de dϕt(v+h) dt = (v+h)ϕt(v+h)=v◦ϕt(v+h)+h ◦ϕt(v+h) (A.2) se calcula su derivada respecto a  : ∂d dt ϕt(v+h)=d dt ∂ϕt(v+h)=∂vϕt(v+h)+∂h ϕt(v+h) =Dv ϕt(v+h)·∂ϕt(v+h)+hϕt(v+h)+·Dh ϕt(v+h)·∂ϕt(v+h) (A.3) siendo ∂ϕt(v+h)= lim δ→0 ϕt(v+(+δ)h)−ϕt(v+h) δ (A.4) Si se sustituye = 0 en (A.3) y (A.4) d dt ∂ϕt(v+h)=Dv ϕtv·∂ϕtv+hϕtv (A.5) ∂ϕtv=0 = l´ım δ→0 ϕt(v+hδ)−ϕtv δ=∂hϕtv (A.6) y por lo tanto d dt ∂hϕtv=Dv ϕtv·∂hϕtv +hϕtv (A.7) 27 28 Esta ecuación tiene la forma de una ecuación diferencial no homogénea: ˙w(t) = Dv ϕtv·w(t) + g (A.8) Una forma de obtener la solución de este tipo de ecuaciones diferenciales es encontrando primero una solución w(t) de la ecuación homogénea asociada ˙w(t) = Dv ϕtv·w(t) (A.9) .En este caso, partiendo de la denición de la velocidad ˙ ϕtv =vϕtv (A.10) se calcula el jacobiano de ambos miembros ˙ Dϕtv =Dv ϕtv·Dϕtv (A.11) y al compararlo con la expresión (A.9) se obtiene w(t) = Dϕtv como solución de la ecuación diferencial homogénea. Para calcular una solución particular por variación de parámetros, se necesita una base del espacio de soluciones de la ecuación homogénea. En este caso, dado que al ser ϕtv una deformación guiada por velocidades ( |D(ϕtv)|>0 ), las las de D(ϕtv) son linealmente independientes y forman una base del espacio de soluciones. A continuación, se ensaya como solución de la ecuación inhomogénea una combinación lineal de los elementos de la base w(t)c(t) , determinando c(t) a a partir de las condiciones iniciales. De esta forma: (w(t)c(t))0=Dv ϕtv·w(t)c(t) + hϕtv (A.12) que comparado con la derivada de (w(t)c(t)) (w(t)c(t))0=˙ w(t)c(t) + w(t)˙ c(t) (A.13) resulta w(t)˙ c(t) = hϕtv (A.14) ˙ c(t) = w−1(t)hϕtv (A.15) con solución c(t) = c(0) + ˆt 0 w−1(u)h(ϕuv)du =c(0) + ˆt 0 (Dϕuv)−1h(ϕuv)du (A.16) donde aplicando la condición inicial para calcular el valor de c(0) : ∂hϕtv |t=0=w(t)c(t)|t=0= 0 (A.17) Id · c(0) + ˆ0 0 (Dϕuv)−1h(ϕuv)du= 0 (A.18) c(0) = 0 (A.19) Apéndice A. Cálculo de la derivada de un campo vectorial respecto a una perturbación 29 Finalmente c(t) = ˆt 0 (Dϕuv)−1h(ϕuv)du (A.20) ∂hϕtv =Dϕtv ˆt 0 Dv (ϕuv)−1h(ϕuv)du (A.21) quedando así demostrado (A.1). Sustituyendo t= 1 se obtiene el caso particular de la expresión (2.8) ∂hϕv=Dϕvˆ1 0 (Dϕuv)−1h(ϕuv)du (A.22)  La expresión anterior depende de la matriz jacobiana del campo de deformaciones Dϕv . Esta matriz tiene tamaño d x d y determinante positivo. Expresa, dado un punto determinado del dominio, cómo se deforma su entorno al aplicarle el difeomorsmo. Partiendo de las ecuaciones A.10 y A.11, renombrando Jt=Dϕtv e introduciendo la condición inicial se obtiene (˙ Jt=Dv (ϕtv)·Jt J0=Id (A.23) ecuación diferencial que puede resolverse mediante el esquema de integración de Euler. Jt+k=Jt+k · Dv ϕtvJt (A.24) Pero siguiendo este método, puede darse el caso de que al pasar del instante t al instante t+k y avanzar un paso k · Dv (ϕtv)Jt , el jacobiano resultante deje de tener determinante positivo. Es recomendable por tanto encontrar un esquema diferente de resolución. Partiendo de la suposición de que en el intervalo (t , t +k) el valor de Dv (ϕtv) es constante, entonces ˙ Jt=C·Jt (A.25) que, de acuerdo a [I. and Havel, 1994] tiene como solución Jt+k=ekC Jt (A.26) Sustituyendo por el valor de C Jt+k=ekDv(ϕtv )Jt=1 + kDv ϕtv+. . .Jt (A.27) se recupera el esquema de Euler y se obtiene nalmente una expresión que permite calcular la matriz jacobiana del difeomorsmo al mismo tiempo que se va calculando la propia deformación. 30 Apéndice B Fórmula BakerCampbellHausdor (BCH) En matemáticas, la fórmula BakerCampbellHausdor es la solución a Z= log eXeY (B.1) para X e Y elementos no conmutativos de un espacio de dimensión nita. Esta ecuación debe su nombre a Henry Frederick Baker, John Edward Campbell, y Felix Hausdor. Fue descrita inicialmente por Campbell; elaborada por Henri Poincaré y Baker; y sistematizada geométricamente por Hausdor. Concretamente, si G es un grupo de Lie simple, con g el Álgebra de Lie asociada X, Y, Z ∈g , exp : g→G la función exponencial, y exp : G→g la función logaritmo (inversa de la exponencial), entonces: Z= log(eXeY) = X n>0 (−1)n−1 nX ri+si>0 1≤i≤n (Pn i=1(ri+si))−1 r1!s1!· · · rn!sn![Xr1Ys1Xr2Ys2. . . XrnYsn] (B.2) donde sn y rn son enteros no negativos según la siguiente notación: [Xr1Ys1. . . XrnYsn] = [X, [X, . . . [X | {z } r1 ,[Y, [Y, . . . [Y | {z } s1 , . . . [X, [X, . . . [X | {z } rn ,[Y, [Y, . . . Y | {z } sn ]] . . .]] (B.3) De acuerdo con [Casas and Murua, 2008], los primeros términos de esta fórmula general combinatoria son: Z(X, Y ) = log(eXeY) =X+Y+1 2[X, Y ] +1 12[[X, Y ], Y ]−1 12[X, [X, Y ]] + 1 24[X, [[X, Y ], Y ]] +1 720[[[[X, Y ], Y ], Y ], Y ] + 1 360[[X, [X, Y ],[X, Y ], Y ]] +1 120[[X, Y ],[[X, Y ], Y ]] + 1 180[X, [[[X, Y ], Y ], Y ]] +1 180[X, [X, [[X, Y ], Y ]]] −1 720[X, [X, [X, [X, Y ]]]] +. . . (B.4) 31 32 Esta fórmula es válida para espacios dimensionalmente nitos, pero se ha demostrado ([Bossa et al., 2007]) que puede ser utilizada con éxito en la composición de difeomorsmos. Así pues, dada una transformación gobernada por un campo de velocidades w y a continuación otra deformación de campo v la fórmula BCH permite calcular el campo z=BCH (v, w) tal que exp (z)≈exp (v)◦exp (w) (B.5) ϕz(x)≈(ϕv◦ϕw) (x) = ϕv(ϕw(x)) El operador [v, w] se dene como [v, w] = wDv −vDw (B.6) donde D · signica matriz jacobiana del campo de velocidades y wDv se calcula como: (wDv)j=X i wi∂ivj (B.7) siendo wi la componente i-ésima de w y ∂i la derivada con respecto a dicha coordenada. El término vDw se calcula análogamente. Concretando para un espacio bidimensional: ([v, w]1= (w1∂1v1+w2∂2v1)−(v1∂1w1+v2∂2w1) [v, w]2= (w1∂1v2+w2∂2v2)−(v1∂1w2+v2∂2w2) (B.8) Apéndice C Parte práctica del TFM Los créditos prácticos de este proyecto se han llevado a cabo en el Servicio de Radiodiagnóstico del Grupo Hospitalario Quirón, Zaragoza, bajo la tutela del Jefe de Servicio, el Dr. Nicolás Fayed Miguel. Uno de los objetivo de esta parte práctica ha sido conocer el equipamiento de radiología (resonancia magnética nuclear, tomografía computerizada, ecografías, rayos X, etc.) así como el sistema PACS ( Picture Archiving and Communication System ) de transferencia, gestión y almacenamiento de imágenes médicas. Cada equipo de radiología está conectado con una consola desde la que se controlan el proceso de adquisición de imagen (posición del paciente, secuencias en el caso de MRI,...) y las imágenes resultantes se almacenan en un servidor, que se encuentra situado en la clínica de La Floresta, al que tienen acceso los workstation de ambas clínicas mediante el software desarrollado por General Electrics. Figura C.1: Workstation. Además de este programa, se utiliza Cosmosalud para las labores administrativas (establecer sesiones para el uso del equipamiento del servicio de radiología) y Klinic para manejar el historial y los datos de hospitalización del (consultas externas, urgencias, quirófano, etc.), así como otros programas más especícos como por ejemplo LCModel6.2, que se utiliza para realizar espectrografías detalladas. 33