Desarrollo de algoritmo rápido de esqueletonización de estructuras óseas en 3D a partir de imágenes médicas

Escuela Técnica Superior de Ingeniería Informática Universitat Politècnica de València Desarrollo de algoritmo rápido de esqueletonización de estructuras óseas en 3D a partir de imágenes médicas Proyecto Final de Carrera Ingeniería Informática Junio 2012 Autor: Mª Isabel Mompeán Escartí Directores: José Vicente Manjón Herrera Ángel Alberich Bayarri Página I Resumen La esqueletonización pretende obtener de un objeto, un patrón continúo que lo represente con la menor cantidad de datos posibles. Con el fin de llevar a cabo dicho proceso existen diferentes métodos. Todos ellos tienen por objetivo la eliminación de puntos del objeto hasta obtener la curva o superficie media de menor espesor que define su forma. Sin embargo, el proceso utilizado para conseguirlo es diferente y por lo tanto, los principales métodos para la esqueletonización se pueden clasificar en tres clases: basados en Adelgazamiento, basados en Transformada de la Distancia y basados en Diagramas de Voronoi. En el presente proyecto se comparan estas tres clases y se analizan sus ventajas y desventajas. Para ello, se estudian los resultados obtenidos al aplicar una serie de algoritmos pertenecientes a cada clase con la finalidad de identificar la metodología más adecuada y desarrollar su versión 3D. Así pues, se presenta un nuevo algoritmo de esqueletonización 3D basado en la transformada de la distancia que extrae el esqueleto de estructuras óseas a partir de imágenes médicas adquiridas mediante resonancia magnética. La intención es conseguir modelos más simples de la estructura del tejido trabecular que nos permitan evaluar la respuesta mecánica y la influencia de la erosión en el hueso mediante simulaciones de estas estructuras. Página II Agradecimientos Ante todo, me gustaría agradecer al Hospital Quirón Valencia la oportunidad que me ha brindado al posibilitarme participar en un proyecto realmente interesante, pudiendo a su vez, utilizarlo como caso de estudio para el desarrollo de mi proyecto final de carrera. Especialmente, me gustaría expresar mi gratitud al Dr. Ángel Alberich-Bayarri, coordinador científico técnico del departamento Cuantificación Quirón, por todo su tiempo invertido y su contribución decisiva al resultado final del proyecto. También quisiera darle las gracias a mi supervisor del proyecto, José Vicente Manjón, profesor del departamento de Física Aplicada de la Universidad Politécnica de Valencia, por su soporte académico, consejo y contribución a esta investigación. Finalmente, quiero agradecer el ánimo recibido por mi familia, novio y amigos, sin ellos este proyecto no habría sido posible. Página III Motivación Desde hace décadas se ha estudiado la complejidad de la forma tratando de reducir su estructura con el fin de obtener un modelo simplificado que contenga a su vez, toda la información suficiente para detallar la geometría y estructura de los objetos que representan. Este modelo de representación se obtiene mediante el proceso de esqueletonización. La extracción del esqueleto es de gran interés en cualquier área científico-técnica por la capacidad del mismo para representar, de forma simplificada, las estructuras complejas de los objetos. Sus aplicaciones prácticas son muchas y muy variadas, desde el análisis gráfico, reconocimiento de formas o compresión de datos, hasta el análisis de imágenes médicas. De particular interés para esta investigación es su aplicación en la imagen médica donde se utiliza para el análisis de imágenes extraídas de pacientes. El esqueleto permite la extracción de la información médica relevante, como las características morfológicas y parámetros para el análisis clínico. Esto posibilita el diagnóstico y tratamiento de enfermedades con el uso de técnicas no invasivas. Debido a las razones expuestas y la necesidad en el área de la medicina, la investigación en este campo ha sido muy amplia y, por lo tanto, es extensa la literatura sobre esqueletonización en 2D, las metodologías existentes y los procesos para su aplicación en imágenes 2D de estructuras muy variadas. Sin embargo, el proceso de esqueletonización en 3D ha recibido menos atención. Es por ello necesario estudiar nuevas metodologías de esqueletonización 3D y desarrollar nuevos algoritmos capaces de obtener esqueletos a partir de estructuras 3D. Los esqueletos 3D de la estructura trabecular ósea (hueso esponjoso) representan mejor las características topológicas y la disposición microestructural de las trabéculas, y resultan de gran ayuda para el diagnóstico y caracterización de pacientes con osteoporosis. En la actualidad, los algoritmos de esqueletonización 3D aplicados a hueso trabecular son escasos, por lo que existe la necesidad de desarrollar nuevos algoritmos para su aplicación en el tejido óseo. De entre los insuficientes ejemplos de algoritmos en 3D aplicados a hueso trabecular, son pocos los eficientes. Por norma general, los algoritmos actuales son muy lentos debido a la complejidad de la estructura trabecular y el nivel de detalle que se obtiene con altas resoluciones de IRM o TC. Sin embargo, el algoritmo presentado en este proyecto está específicamente diseñado para resolver estructuras complejas de forma eficiente, con un tiempo de computación muy bajo en comparación con el resto de propuestas existentes hasta el momento aplicadas sobre la misma materia. Página IV Finalmente, el desarrollo de nuevos algoritmos de esqueletonización 3D más rápidos y eficientes, se pueden integrar perfectamente en una plataforma de postprocesamiento clínico. Como resultado, se facilita la investigación de la caracterización in vivo de la microarquitectura del hueso trabecular y la evaluación morfométrica y mecánica del mismo, pudiendo mejorar el diagnostico y tratamiento en enfermedades que afectan directamente al tejido óseo. Página V Objetivos El propósito de este proyecto es desarrollar un algoritmo de procesado de imagen basado en esqueletonización 3D para su aplicación en el estudio de la microarquitectura del hueso trabecular. Este objetivo puede subdividirse en cuatro objetivos más concretos que deben alcanzarse para completar este proyecto: 1. Analizar el estado del arte de las metodologías de esqueletonización e identificar aquella más adecuada para la aplicación a estructuras 3D de hueso trabecular obtenidas a partir de imágenes médicas. 2. Desarrollar un algoritmo de esqueletonización 3D basado en la metodología escogida. 3. Aplicar el algoritmo desarrollado a estructuras óseas extraídas a partir del procesado de imágenes de pacientes. 4. Evaluar las ventajas frente a otros algoritmos de similar fundamento utilizados en la actualidad en el área de la imagen médica. Página VI Índice de contenidos Resumen ...................................................................................................................................I Agradecimientos ......................................................................................................................II Motivación ............................................................................................................................. III Objetivos ................................................................................................................................. V 1. Introducción ..................................................................................................................... 1 1.1. Antecedentes de la informática médica...................................................................... 2 1.2. Adquisición de imágenes médicas .............................................................................. 5 1.3. Reconstrucción de imágenes de hueso trabecular .................................................... 16 1.3.1. Reconstrucción 2D ........................................................................................... 16 1.3.2. Reconstrucción 3D ........................................................................................... 17 2. Morfología y Esqueletonización...................................................................................... 19 2.1 Formas básicas en 2D Y 3D ....................................................................................... 20 2.2 Características Morfológicas..................................................................................... 22 2.2.1 Homotopía (Preservar la topología).................................................................. 22 2.2.2 Invarianza bajo transformaciones isométricas .................................................. 22 2.2.3 Reconstrucción ................................................................................................ 23 2.2.4 Delgadez .......................................................................................................... 23 2.2.5 Conectividad .................................................................................................... 23 2.2.6 Centrado .......................................................................................................... 24 2.2.7 Robustez .......................................................................................................... 24 2.2.8 Detección de componentes .............................................................................. 24 2.2.9 Suavidad .......................................................................................................... 25 2.2.10 Complejidad Estructural ................................................................................... 25 2.2.11 Eficiencia.......................................................................................................... 25 2.2.12 Exactitud .......................................................................................................... 26 2.2.13 Aplicabilidad de las características en los métodos analizados.......................... 26 2.3 Metodologías de Esqueletonización ......................................................................... 27 2.3.1 Adelgazamiento ............................................................................................... 27 2.3.2 Transformada de la distancia (TD) .................................................................... 31 Página VII 2.3.3 Métodos Geométricos ..................................................................................... 34 2.3.4 Funciones de Campo-General .......................................................................... 37 3. Esqueletonización en 2D ................................................................................................. 39 3.1 Algoritmos de esqueletonización 2D ........................................................................ 40 3.2 Discusión.................................................................................................................. 51 4. Desarrollo e implementación de la esqueletonización en 3D ......................................... 53 4.1 Estructura del algoritmo ........................................................................................... 54 4.2 Aplicación a hueso trabecular................................................................................... 63 5. Resultados y Discusión ................................................................................................... 67 5.1. Tiempos de ejecución ............................................................................................... 70 6. Conclusiones ................................................................................................................... 72 7. Futuras líneas de investigación y recomendaciones ....................................................... 73 Referencias............................................................................................................................. 75 Abreviaturas ........................................................................................................................... 81 Apéndices ............................................................................................................................... 82 A) Algoritmo de Esqueletonización 3D .......................................................................... 83 B) Función Procesado del esqueleto final ..................................................................... 84 C) Función Detección puntos extremos y uniones en el esqueleto ................................ 85 Mª Isabel Mompeán Escartí Página 7 IMAGEN POR RESONANCIA MAGNÉTICA La Imagen por Resonancia Magnética (IRM) es una de las modalidades de imagen médica de más reciente aparición y que ha experimentado un avance más rápido. El nacimiento de los fundamentos físicos de esta técnica de obtención de imágenes tiene sus orígenes en las primeras décadas del siglo XX, más concretamente en 1929, cuando el físico Isidor Isaac Rabi y su equipo de científicos investigaban con haces de moléculas haciéndolas atravesar una cámara de vacío bajo la acción de distintos campos magnéticos y así, poder estudiar las propiedades magnéticas y la estructura interna de las moléculas, los átomos y los núcleos de estas. Por estas investigaciones fue galardonado con el premio Nobel de física en 1944. El desarrollo más importante de la RM fue descubierto en 1946, de forma simultánea pero independiente por Felix Bloch de la universidad de Stanford y Edward Purcell de la universidad de Harvard. Ambos descubrieron “el fenómeno de la resonancia magnética nuclear en materia condensada” y compartieron el premio Nobel de física en 1952 (Ferreirós Domínguez J, 2004). A principios de la década de los 70, un experimento realizado por Raymonde Damadian, en el que se habían implantado en ratas una serie de tumores de rápido crecimiento, mostró que los tejidos malignos extirpados presentaban una señal diferente a la de los sanos al ser observados mediante RM. Estos experimentos no demostraron la suficiente fiabilidad clínica para la detección del cáncer, con lo que Paul Lauterbur ideó una técnica mejorada para obtener la información suficiente en el diagnóstico de tumores, observando los tiempos de relajación en tejidos sanos y tumorales. La clave consistía en superponer un campo magnético más débil que variaba de posición, creando lo que se conoce hoy en día como gradiente de campo magnético (Vlaardingerbroek MT, Den Boer JA, 2003). El mismo año en que Godfrey Hounsfield introdujo el TAC (1973), Paul Lauterbur publicó la utilización del fenómeno de la RM para producir imagen. Combinó la idea del gradiente con la idea del escáner tomográfico para realizar varias proyecciones y así, obtuvo la primera imagen por RM en 1975. A su vez, en este mismo año, se instalaba el primer TAC en España. También, en este mismo año, Richard Ernst propuso el procedimiento actual de la técnica IRM, utilizando gradientes de campo magnético en los tres ejes del espacio, codificando la frecuencia y fase, y utilizando la transformada de Fourier para obtener las imágenes de RM, lo que le valió el Premio Nobel de Química de 1991 (Ferreirós Domínguez J, 2004). Mª Isabel Mompeán Escartí Página 8 Principios Físicos Al someter un haz de moléculas a ondas electromagnéticas de una frecuencia fija, mientras se varía la intensidad del campo magnético aplicado, se observa que los momentos magnéticos de los núcleos se pueden invertir con respecto a la dirección del campo magnético a una frecuencia muy precisa. Explicado formalmente, los núcleos de los átomos de la materia se componen en protones y neutrones. Estas cargas presentan la propiedad de espín, el giro de las partículas alrededor del propio eje. Dadas las propiedades magnéticas de las cargas en movimiento, los espines tienen asociado un momento magnético µ, representado por un vector con dirección al eje de giro. Se muestra un ejemplo gráfico en la figura 1.3. µ 𝑩𝒁 s Fig. 1. 3. Diagrama del movimiento de precesión del spin del electrón S, dirigido en dirección opuesta a la dirección del vector momento magnético µ, bajo la aplicación de un campo magnético Bz. Cuando el número de protones o de electrones es par, los espines tienen sentidos opuestos y sus momentos magnéticos se anulan. Por eso, la RM es efectiva sólo en núcleos con un número impar de protones y/o neutrones, como el Hidrógeno (1H), Carbono (3C), Flúor (9F), Sodio (23Na) y Fósforo (31P). Debido a que el 1H es el elemento más abundante en el cuerpo humano, se estudia el comportamiento de sus espines en RM. Mª Isabel Mompeán Escartí Página 9 En ausencia de un campo magnético, los espines de los protones están orientados al azar, pero tras aplicar un campo magnético B0, estos tienden a alinearse en sentido paralelo (estado energético bajo) o en sentido antiparalelo (estado energético más alto), siempre en dirección a B0. Ver figura 1.4. Este fenómeno se explica por la mecánica cuántica. Debido al campo magnético y al movimiento de rotación de los espines, se produce un movimiento giroscópico o de precesión. Este movimiento lleva asociado una frecuencia de Larmor o de precesión directamente proporcional al campo B0, y se define como: (Ec. 1.1) Donde  es la constante giromagnética y para el protón de Hidrógeno, bajo un campo magnético B0 de 1T, la frecuencia de Larmor f0 es  H = 42.58 MHz/T. B µ µ µ µ µ µ µ µ µ µ µ Fig. 1. 4. Diagrama de los momentos magnéticos. A la izquierda, en ausencia de campo magnético, los espines se orientan aleatoriamente. A la derecha, tras aplicar un campo magnético B, los espines se alinean con el campo magnético, en sentido paralelo o antiparalelo. El vector resultante de sumar todos los vectores de las partículas, es un vector de magnetización M0, paralelo a B0 y proporcional al número de átomos de 1H. Se produce entones el estado de equilibrio, en el que la componente longitudinal MZ=M0 y la componente trasversal se anula MXY=0. El proceso de RM mide la relación que emite una señal denominada FID (Free Induction Decay) con respecto al tiempo. Mediante un proceso de excitación, provocado por un pulso de radiofrecuencia (RF), se rompe el estado de equilibrio. El tiempo que tardan los núcleos en volver a su estado, se llama proceso de relajación, e induce la señal FID que da lugar a la imagen de RM. Mª Isabel Mompeán Escartí Página 10 Secuencias de Pulsos Una secuencia de pulsos es la distribución en el tiempo de las distintas señales de RF que se envía al paciente para generar ecos de la FID necesarios para la correcta medición de esta señal. La obtención de IRM en sus comienzos era un proceso lento. Por ello, el diseño de nuevas secuencias de pulsos RF es crucial, y su finalidad aparte de conseguir mejores contrastes, es la de disminuir el tiempo de exploración en los pacientes, que se consigue disminuyendo el tiempo en que se genera un eco (TE). Existen muchas configuraciones de pulsos, pero se consideran dos secuencias de pulsos como la base de todas las demás. Estas son: un sistema RF en secuencias Espín-Eco (SE) y un sistema de gradientes en secuencias Gradiente-Eco (GE). Secuencias SE: Las secuencias SE empiezan con un pulso de 90 para excitar los espines de los núcleos. Seguidamente se emite un pulso de 180 para compensar el desfase producido en la precesión de los espines por inhomogeneidades del campo magnético. El tiempo que tardan los espines en relajarse tras el pulso inicial de 90, es de TE/2 al perder la coherencia de fase. El pulso de 180 es importante para que los espines recuperen su coherencia de fase y se produzca el eco de la señal FID. El proceso se repite cada TR segundos. La figura 1.5 muestra el diagrama de las secuencias SE. Señales Pulsos RF 90 90 180 S-Eco FID TR TE/2 FID TE Fig. 1. 5. Secuencia pulsos SE. El tiempo que tarda en aparecer el eco de la señal FID tras el pulso de 180 es el mismo tiempo entre el pulso de 90 y el de 180 Secuencias GE: Este tipo de secuencias se caracterizan por reducir el tiempo de respuesta debido a que se envía un pulso menor de 90, lo que permite una recuperación más rápida. Para poder generar un eco que sea medible al igual que en las secuencias SE, se sustituye el pulso de 180 por un gradiente de campo magnético bipolar (desfase, refase), que provoca un primer desfase en los espines y luego una coherencia de fase que da lugar a la señal eco. Se muestra a continuación en la figura 1.6 un diagrama de las secuencias GE. < 90 Pulsos RF < 90 Señales FID G-Eco FID TE TR Fig. 1. 6. Secuencia pulsos GE. Mª Isabel Mompeán Escartí Página 11 Fig. 1. 7. a) Imagen obtenida con una secuencia SE, potenciada en T1. Con este tipo de secuencia de pulsos no hay mucho contraste entre la sustancia blanca y la sustancia gris del cerebro. b) Imagen de SE, potenciada en T2. En ella brillan los líquidos y también la grasa subcutánea. El parénquima encefálico es más oscuro que en la potenciación T2*. c) Imagen obtenida con secuencias GE, potenciada en T2*. Aparece con una tonalidad más clara que la de b). También brillan intensamente los líquidos. La grasa subcutánea se oscurece de manera suave, en comparación con T2. El cráneo es completamente negro en este tipo de secuencias y no se aprecia la banda central del díploe craneal que sí se puede observar en las imágenes potenciadas en T2 obtenidas con secuencias SE. Las ventajas de la RM son muchas y muy variadas. La principal es que no utiliza radiación ionizante, por lo que es inocua para el organismo, a diferencia de los Rayos X que se utilizan en el diagnóstico convencional o el TAC, que son radiaciones ionizantes con efectos biológicos. Otra de las ventajas de la RM es que proporciona una buena resolución espacial y temporal, así como un excelente contraste. La RM puede proporcionar un gran abanico de imágenes morfológicas y funcionales. Se ha convertido en la principal herramienta de diagnostico, desplazando a otros procedimientos como TC y técnicas radiológicas, cuya base es la radiación ionizante y puede provocar alteraciones cromosómicas en el organismo. Los avances más importantes en IRM se derivan de la introducción de nuevos mecanismos de contraste que producen nuevas imágenes morfológicas o estructurales. Con estas técnicas también se puede estudiar el flujo sanguíneo (angiografía por RM) u obtener imágenes funcionales, por ejemplo, de la activación cerebral. Tras describir los métodos de adquisición de imagen, y dado que las aplicaciones de cada uno de ellos son muy amplias, nos centraremos en la aplicación a hueso trabecular con el fin de realizar un estudio, en líneas generales, del conjunto de información que se obtiene con cada una de estas modalidades. Sin embargo, lo primero es describir la estructura y morfología del tejido óseo para facilitar la mejor comprensión de sus aplicaciones. a) T1 b) T2 c) T2* Mª Isabel Mompeán Escartí Página 12 CARACTERIZACIÓN DEL HUESO El hueso es un tejido conjuntivo mineralizado muy vascularizado e inervado, que está estructurado en laminillas de matriz osteoide calcificada (Fernández-Tresguerres HernándezGil I et al., 2006). Se distinguen dos clases de tejido que desempeñan una función específica en cada sección en la cual se encuentran; estos son: el hueso trabecular y el cortical. Fig. 1. 8. El gráfico muestra la estructura del hueso compacto o cortical y una visión del interior en el que se encuentran las trabéculas, hueso esponjoso. (Visto en: Natural Medicine Seminars) El hueso cortical es una masa sólida y continua con canales únicamente microscópicos. Aproximadamente, el 80% de la masa esquelética de un ser humano adulto es hueso cortical, el cual se forma como una pared exterior de los huesos y sus funciones principales son de soporte y protección. Su matriz ósea mineralizada está dispuesta en laminillas, formada por osteoblastos. Las laminillas se disponen de 3 formas, aunque la gran mayoría se dispone concéntricamente alrededor de un canal longitudinal vascular llamado conducto de Havers, que contiene capilares y vénulas postcapilares formando estructuras cilíndricas llamadas osteonas o sistemas haversianos visibles únicamente al microscopio óptico. Los canales haversianos comunican entre sí con la superficie o la cavidad medular por canales transversales u oblicuos llamados canales perforantes o de Volkman, que poseen vasos que vienen del periostio y del endostio. El periostio es la superficie externa de los huesos y está en contacto con el tejido blando adyacente. El endostio es la superficie interna frente a la cavidad medular. El hueso trabecular también conocido como poroso o esponjoso, está compuesto de una red interconectada de placas y barras que delimitan las cavidades areolares, en cuyo interior se encuentra la médula ósea, recibiendo estas el nombre de trabéculas. A densidades relativamente bajas, las celdas forman una red abierta de barras. Conforme la densidad relativa aumenta, más material se acumula en las paredes de las celdas y la estructura se transforma en una red cerrada de placas. Hueso compacto Cavidad medular Cartílago articular Trabécula s Hueso Esponjoso Trabéculas Hueso Cortical Cartílago Mª Isabel Mompeán Escartí Página 13 El hueso esponjoso se encuentra principalmente cerca de los terminales de los huesos largos y en los cuerpos vertebrales. En los huesos largos se distribuyen las cargas mecánicas desde las superficies articulares hasta el hueso cortical, mientras que en las vértebras, es el constituyente principal que absorbe la carga. Fig. 1. 9. Esquema de la estructura del hueso. Se observa en la parte exterior del hueso compacto las osteonas, formadas por un conducto central y longitudinal llamado el Canal de Havers, que está rodeado por anillos concéntricos (laminillas). Estos canales se interconectan mediante canales transversales llamados los canales de Volkman. En el interior del hueso se encuentran las trabéculas formadas por placas y barras y por la medula ósea. (Visto en: SEER Training) La densidad, arquitectura y función del hueso en una posición localizada, dependen de la magnitud de las cargas a las que está sometido. En las imágenes, las estructuras con bajas densidades y celdas abiertas se desarrollan en regiones de bajos esfuerzos, mientras que las de mayores densidades y celdas cerradas se desarrollan en regiones sometidas a altos esfuerzos. Existe una creciente evidencia de que, además de la fracción volumétrica de los huesos, a menudo cuantificados en términos de la densidad ósea (Gomberg BR et al., 2000), las propiedades mecánicas de este tipo de hueso se pueden determinar no sólo por su porosidad, sino también por la disposición estructural de las trabéculas. Por lo tanto, se ha sugerido que la predicción de las propiedades mecánicas del hueso, que hasta hace poco se basaba sólo en la densidad ósea, se puede aumentar si se estudia también la estructura trabecular (Ulrich D et al., 1999). Es decir, la relación entre la disposición de placas y barras, y las propiedades elásticas que determinan el comportamiento mecánico del hueso esponjoso (Pothuaud L et al., 2004). TRABECULAS DEL HUESO ESPONJOSO OSTEONA DEL HUESO COMPACTO CONDUCTO DE HARVERS CANAL DE VOLKMANN LAMELAS INTERSTICIALES OSTEONA PERIOSTIO CANALÍCULO OSTEOCITOS Mª Isabel Mompeán Escartí Página 14 Fig. 1. 10. Microfotografía de las trabéculas del tejido óseo. En a) de un hueso sano, en estado normal, y en b) se observa el deterioro de las trabéculas por ostoporosis. (Visto en: Science photo library) Tras comentar la caracterización del hueso trabecular, que es en definitiva el objeto de estudio de este proyecto, nos centraremos en explicar cada metodología de obtención de imagen aplicada al tejido óseo. La técnica utilizada debe contar con la especificidad y precisión suficiente para estudiar individualmente los distintos elementos articulares y hacer estimaciones cuantitativas de las alteraciones. Además, la técnica debe ser tan inocua como para poder utilizarse repetidamente durante la prolongación de una enfermedad en estudio: La Radiografía convencional proporciona excelentes detalles de las características óseas, pero en general, es pobre en la definición de las estructuras de los tejidos blandos. El objeto más pequeño registrado en la película va de 0,1 mm a 0,2 mm, excepto con la ayuda de un procedimiento invasivo que implica la introducción de un medio de contraste (BucklandWright JC, 1994). TC es más sensible que la radiografía convencional. Obtiene una imagen del tejido óseo de gran calidad y puede proporcionar reconstrucciones 3D. Se ha demostrado su eficacia para evaluar cambios precoces en el hueso trabecular en modelos experimentales de artrosis (Dequeker J et al., 1997). Como técnica basada en rayos X, su principal inconveniente reside en que se trata de una prueba que conlleva una cantidad no despreciable de radiación ionizante, lo que limita su empleo en estudios secuenciales. Con la ecografía sólo pueden valorarse las alteraciones en la superficie y el grosor del cartílago. Debido a la densidad del hueso, las ondas ultrasónicas no pueden penetrar en el hueso y obtener así la estructura de las trabéculas. La obtención de imágenes mediante RM constituye un método de diagnóstico ideal para la evaluación del tejido óseo en general. Debido a la alta resolución espacial proporciona una definición anatómica y clara de los tejidos blandos, y obtiene imágenes en múltiples planos (Wehrli FW, 2007). También posibilita la cuantificación de sustancias y variables funcionales implicadas en la fisiopatología del aparato locomotor. Además, esta técnica no tiene efectos adversos y se tolera lo suficientemente bien como para utilizarla en exámenes secuenciales, incluso en pacientes asintomáticos. Desde el punto de vista experimental no es destructiva, por lo que puede utilizarse para estudiar simultáneamente varias características en la misma zona. (a) (b) Mª Isabel Mompeán Escartí Página 15 Una vez analizados los métodos de adquisición de imagen, se puede concluir lo siguiente. Debido al excelente contraste entre el hueso y los tejidos blandos (Hildebrand T et al., 1999), parece que la técnica de TC sea adecuada para visualizar el tejido óseo de manera macroscópica. Sin embargo, dadas las limitaciones por la dosis de radiación ionizante y el fenómeno de emborronamiento por la función de punto disperso (PSF) producida por los rayos X, dificulta su aplicación al estudio del tejido óseo de manera microscópica con buena resolución, incluso en regiones anatómicas periféricas como brazos y piernas. Como resultado, la RM parece ser el método más eficaz para estudiar la caracterización in vivo de la microarquitectura del hueso trabecular. Los actuales sistemas de RM permiten una gran resolución espacial en cualquier orientación del plano, garantizando una visión clara de las características más relevantes de la estructura de interés. Asimismo, proporcionan un alto contraste entre el hueso y la medula ósea, además de ser una técnica inocua por su radiación no ionizante. Las ondas electromagnéticas utilizadas en RM tienen baja energía y parecen ser seguras en condiciones normales (Vlaardingerbroek MT, Den Boer JA, 2003). En (Alberich-Bayarri A, 2010) se demuestra, mediante la configuración de un protocolo de adquisición RM optimizado con un sistema de 3 Teslas, que la utilización de RM consigue obtener una alta resolución espacial, una satisfactoria relación señal-ruido y un contraste bien definido, para realizar un análisis mecánico de las imágenes in vivo del hueso trabecular. Mª Isabel Mompeán Escartí Página 16 1.3. Reconstrucción de imágenes de hueso trabecular Las imágenes 2D de TC o RM son de gran utilidad para realizar diagnósticos. Sin embargo, tienen el inconveniente de que las estructuras u órganos internos no pueden verse en su forma original. Para solucionar este problema se utiliza la reconstrucción 3D. La mejora que aportan las imágenes en 3D es que sirven tanto para el diagnostico, como para la planificación y simulación quirúrgica. Para ello, las técnicas de reconstrucción 3D utilizan los cortes transversales de las imágenes obtenidas con TC o RM. Dichos cortes proporcionan toda la información necesaria para crear vistas 3D de la estructura de estudio y sus partes internas siendo, por lo tanto, la base para poder visualizar un volumen en 3D. En este apartado se comentarán las técnicas de reconstrucción y visualización tanto para 2D como 3D. 1.3.1. Reconstrucción 2D La utilidad de la imagen 2D depende de la orientación física del plano de la imagen con respecto a la estructura de interés. La mayoría de sistemas de obtención de imagen tienen una capacidad limitada para obtener una imagen óptima directamente, debido a que el posicionamiento y la orientación del escáner están generalmente restringidos. Por lo tanto, es importante un proceso para generar y mostrar imágenes 2D óptimas a partir de volúmenes reales. Los pasos comunes que se siguen en el proceso de reconstrucción de imágenes médicas de cualquier región de la anatomía humana, con pequeñas variaciones dentro de cada fase dependiendo de la calidad y resolución espacial de los datos obtenidos, son los siguientes: Fig. 1. 11. Proceso de reconstrucción de imágenes médicas. Este proceso es equivalente en la reconstrucción de imágenes médicas 3D. Así pues, se explicará en el siguiente apartado, a grandes rasgos, la reconstrucción de imágenes 3D. Concretamente, la reconstrucción 3D de IRM de hueso trabecular, necesaria a la hora de aplicar algoritmos de post-procesamiento clínico, como el algoritmo de esqueletonización 3D desarrollado en este proyecto. Detección Región de Interes Segmentación Umbralización Interpolación Binarización Mª Isabel Mompeán Escartí Página 23 2.2.3 Reconstrucción La reconstrucción se refiere a la capacidad de recuperar el objeto original a partir de la curva-esqueleto. Teniendo en cuenta la relación de la curva-esqueleto y la definición del esqueleto como el conjunto de centros de las esferas máximas inscritas, un buen método de reconstrucción consiste en calcular la unión de todas las esferas máximas centradas en cada punto de la curva-esqueleto (Blum H, 1967; Gagvani N, Silver D, 2001). El radio de cada esfera viene dado por el valor de la transformada de la distancia, la cual especifica la distancia al punto más cercano del contorno del objeto. Esta propiedad tiene una aplicación inmediata en la compresión de la forma y la animación del volumen. La reconstrucción puede ser mejorada almacenando más información de cada punto de la curva-esqueleto o incrementando el número de ramificaciones. La habilidad de reconstruir el objeto desde su esqueleto, podría ser utilizado como un indicador de la calidad de la forma en tareas de análisis. Si el grado de reconstrucción es muy bajo significa que el esqueleto no representa bien el objeto original. Sin embargo, recientes investigaciones niegan esta afirmación. En (Shilane P et al., 2004) se muestra que usando algunos de los mejores procedimientos para la representación de formas, no se consigue reconstruir el objeto original en ningún caso. 2.2.4 Delgadez La curva-esqueleto debe ser unidimensional, lo que significa un vóxel de espesor en cada dirección con excepción en las uniones, donde el grosor puede ser mayor para asegurar la conectividad entre las diferentes ramas. Para los objetos 3D, si calculamos la superficieesqueleto, pueden existir regiones que también tengan más de un vóxel de espesor debido a la presencia de áreas que se representan como superficies. La delgadez y la reconstrucción son dos propiedades en conflicto. Incluso para objetos cuya superficie media realmente contiene sólo curvas, por ejemplo, objetos tubulares, una curva-esqueleto de un vóxel de espesor podría no contener a todas las esferas máximas necesarias para reconstruir el objeto con precisión, ya que la superficie media discreta es normalmente más espesa que un vóxel, debido a la naturaleza discreta del objeto (Cornea ND et al., 2007). 2.2.5 Conectividad Esta característica es una consecuencia directa de la homotopía. Si la curva-esqueleto corresponde a un objeto con una sola componente conexa, entonces, para mantener la topología de este objeto, la curva-esqueleto tendría que estar formada por una sola componente conexa también. Mª Isabel Mompeán Escartí Página 24 2.2.6 Centrado Una característica importante de la curva-esqueleto es el hecho de que se encuentre centrada respecto al objeto. Para conseguir un centrado perfecto, es necesario que la curvaesqueleto recaiga sobre la superficie media ya que la superficie media está centrada respecto del objeto. En el reconocimiento de formas u otras aplicaciones científicas, el centrado exacto de la curva-esqueleto es esencial. Sin embargo, en la mayoría de los casos el centrado del esqueleto no es necesario. Debido a la sensibilidad de la superficie media a pequeñas perturbaciones en el límite del objeto, limitar el esqueleto a que coincida con la superficie media podría provocar que esta fuese igualmente sensible a este tipo de cambios. En lugar de un centrado exacto, un centrado aproximado es probablemente suficiente para muchas aplicaciones como la navegación virtual o la animación. Por ejemplo, para colonoscopías virtuales, son características más importantes la precisión y la suavidad en el camino de navegación. 2.2.7 Robustez Una propiedad deseable de la curva-esqueleto es que la sensibilidad al ruido sea mínima para que resulten similares un esqueleto procedente de un objeto libre de ruido y el de un objeto con ruido. Como se muestra en la figura 2.4, el eje medio es muy sensitivo a pequeños cambios en los contornos del objeto. Una curva-esqueleto robusta no puede ser exactamente centrada, ya que el centrado exacto limita la curva-esqueleto a la superficie media la cual es extremadamente sensible a las perturbaciones de los contornos del objeto. 2.2.8 Detección de componentes El esqueleto debe ser capaz de distinguir las diferentes componentes del objeto original reflejando las partes de su estructura. Esto significa que las componentes lógicas del objeto deben tener una correspondencia unívoca con las componentes lógicas de la curva-esqueleto. Las componentes lógicas del objeto se pueden hallar con los puntos de unión del esqueleto. Una partición del objeto original puede ser realizada para producir una correspondencia uno a uno entre las diferentes componentes del objeto y de la curvaesqueleto. Comprobar si una curva-esqueleto satisface esta propiedad es una tarea complicada porque la definición de las componentes del objeto no es suficientemente precisa, ya que envuelve la percepción humana que es subjetiva. Sin embargo, definiciones específicas para cada aplicación pueden ser usadas para este propósito. Por ejemplo, el esqueleto de la forma de una persona, para la animación digital, podría ser modelado tal y como un niño dibuja a una persona, es decir, con barras y círculos para cada extremidad del cuerpo. Mª Isabel Mompeán Escartí Página 25 2.2.9 Suavidad La suavidad no es solamente una característica estética, sino que es realmente útil en algunas aplicaciones. Por ejemplo, en la navegación virtual, que usa la curva-esqueleto como la trayectoria del camino para una cámara virtual, la curvatura debe ser tan suave como se pueda para evitar cambios abruptos en la imagen mostrada. Se puede definir la suavidad de un segmento del esqueleto como la variación de la tangente sobre la dirección de la curva a lo largo de todo el segmento. Para ello, se pueden medir los ángulos entre direcciones tangentes en localizaciones sucesivas a lo largo de la curva y tomar la desviación estándar como una medida de variación. Con el fin de asegurar una navegación suave, la variación en las direcciones tangentes a medida que nos desplazamos de un punto al siguiente debería ser tan pequeña como fuese posible. 2.2.10 Complejidad Estructural Debido al hecho de que el esqueleto es una aproximación de las componentes complejas de un objeto, el proceso de esqueletonización y la propia curva-esqueleto deberían reflejar la complejidad estructural del objeto de estudio. Un esqueleto que se aproxime exactamente a la complejidad del objeto original se requiere en aplicaciones de compresión para poder reconstruir su forma en el análisis de imágenes médicas o en el reconocimiento de formas. Algunas metodologías no permiten el control total de la complejidad estructural de los esqueletos. Por ejemplo, los métodos basados en adelgazamiento topológico funcionan muy bien para objetos lisos y regulares, ya que el resultado de su esqueleto es una línea central muy simplificada. Sin embargo, la mayoría de los objetos del mundo real presentan estructuras complejas y con contornos muy abruptos, lo que requiere estrategias capaces de mantener la complejidad de la forma. 2.2.11 Eficiencia Esta propiedad es necesaria para todas las aplicaciones por ordenador, no sólo para las comentadas a lo largo del proyecto. Es cierto que muchas aplicaciones necesitan una computación en tiempo real para lo que la eficiencia es crucial como, por ejemplo, en el ámbito de la medicina o en la simulación por ordenador. En nuestro caso determinadas metodologías de extracción de esqueletos, que se verán a continuación, presentan un mayor tiempo de ejecución en la extracción del esqueleto de estructuras complejas. Mª Isabel Mompeán Escartí Página 26 2.2.12 Exactitud La exactitud, se refiere a la propiedad de la curva-esqueleto por la que todos los puntos del contorno de la imagen son visibles desde al menos un punto de la curva-esqueleto. En otras palabras, para cualquier punto de la superficie existe una línea recta que lo conecta con un punto del esqueleto y, además, no cruza ningún otro contorno. Esta característica es de obligado cumplimiento en endoscopias virtuales, donde se debe asegurar que la superficie interior del órgano sea plenamente examinada por el facultativo cuando realiza el procedimiento virtual. Un algoritmo de fuerza bruta que comprueba la exactitud de la curva-esqueleto, revisa la visibilidad de cada punto de la frontera con una línea recta a cada punto de la curva-esqueleto. Aquellos puntos de la frontera que no puedan ser conectados sin intersecar otra superficie, son no visibles. 2.2.13 Aplicabilidad de las características en los métodos analizados A continuación se muestra una tabla con las propiedades descritas, indicando si se cumplen o no en los 4 métodos que se describirán en la siguiente sección. Cuando una propiedad se puede cumplir para algún algoritmo de la clase pero no para todos, la celda correspondiente se presenta en blanco. Adelgazamiento Transformada de la distancia Métodos Geométricos Campo General Homotopía S S N Invariancia bajo transformaciones S S Reconstrucción N N N Delgadez S Conectividad S Centrado Robustez N N N S Suavizado S Jerarquía N S Eficiencia S S S N Exactitud Tabla 2. 1. Resumen de las características morfológicas. Mª Isabel Mompeán Escartí Página 27 2.3 Metodologías de Esqueletonización Una revisión del estado del arte de las metodologías de esqueletonización podría clasificar los tipos de métodos en base a la implementación del paso inicial para la extracción del esqueleto en: Adelgazamiento Transformada de la distancia (TD) Métodos Geométricos: Diagramas de Voronoi (DV) Funciones de campo general 2.3.1 Adelgazamiento El proceso de adelgazamiento consiste en eliminar puntos de la imagen desde los límites del objeto y continuar hacia el interior, hasta que no se pueden eliminar más puntos, mediante una serie de restricciones. El objetivo es obtener la línea o superficie media que preserva la topología y mantiene la forma del objeto tanto como sea posible. Esta clase de algoritmos dependen de la clasificación de los puntos del objeto introducido por Morgenthaler (Morgenthaler DG, 1981). Antes de seguir con la explicación del método de adelgazamiento, introduciremos unos conceptos básicos sobre la teoría de la topología que utiliza las relaciones de adyacencia propuestas por Kong y Rosenfeld (Kong TY, Rosenfeld A, 1989).  Dos vóxeles son 6-adyacentes si sólo comparten cualquiera de sus 6 caras.  Dos vóxeles son 18-adyacentes si sólo comparten cualquiera de sus caras y aristas.  Dos vóxeles son 26-adyacentes si sólo comparten cualquiera de sus caras, aristas y esquinas (vértices). Si dos puntos son -adyacentes donde ∈ {6, 18, 26} entonces estos dos puntos también se pueden llamar -vecinos. Ver figura 2.5. Fig. 2. 5. Representación de 6-adyacencia (a), 18-adyacencia (b) y 26-adyacencia (c). (a) (b) (c) Mª Isabel Mompeán Escartí Página 28 Así pues, para extraer el esqueleto de la imagen debemos detectar todos aquellos puntos simples. Una característica importante de los mismos es que se pueden identificar solamente estudiando sus vecinos locales. Un punto simple, dentro del concepto de volumen, es un vóxel que puede ser eliminado sin cambiar la topología del objeto (Kong TY, Rosenfeld A, 1989). Esto significa que el número de componentes conexas, de agujeros y de cavidades del objeto no cambian (Saha PK, Chaudhuri BB, 1994-1996). Teniendo en cuenta que se utiliza la nomenclatura de punto o vóxel negro, para referirnos a vóxeles del objeto, y punto o vóxel blanco, para vóxeles que no pertenecen al objeto, definiremos un punto P como punto-simple, si satisface las siguientes condiciones (Saha PK et al., 1997):  P tiene al menos un punto negro 26-vecino.  P tiene al menos un punto blanco 6-vecino.  El conjunto de puntos negros 26-vecinos de P está 26-conectado.  El conjunto de puntos blancos 6-vecinos de P está 6-conectado al conjunto de puntos blancos 18-vecinos de P. En cualquier otro caso la eliminación de un punto puede formar nuevas cavidades o túneles, la cual alteraría la topología de la forma original. En la literatura encontramos diversos diseños para el algoritmo de adelgazamiento. El esquema propuesto por She consiste en estudiar cada vóxel en cada iteración e ir eliminándolos capa a capa (She FH et al., 2009): Repetir 1. Marcar todos los puntos negros que son 26-adyacentes a un punto blanco. 2. Repetir  Borrar del conjunto de puntos marcados, aquellos que cumplen las condiciones anteriores,  Hasta que se han borrados todos los puntos marcados, o no se puede borrar ninguno más. 3. Desmarcar todos los puntos sin eliminarlos y volver al paso 1. Hasta que no se pueden marcar más puntos. Del mismo modo, Saha propone un algoritmo de adelgazamiento que consta de dos pasos principales: “adelgazamiento-primario” y “adelgazamiento-final”. Como resultado de los mismos, se obtiene un “esqueleto-primario” y un “esqueleto-final” utilizando los conceptos de s-puntos, e-puntos y v-puntos para analizar la topología del objeto (Saha PK et al., 1997). A continuación, se define cada uno de estos conceptos y en la figura 2.6 se incluye una descripción gráfica de ellos. Mª Isabel Mompeán Escartí Página 29 s-puntos: conjunto de puntos que son 6-adyacentes entre ellos. Ver figura 2.6 (a). e-puntos: conjunto de puntos que son 18-adyacentes entre ellos, pero no 6-adyacentes. Ver figura 2.6 (b). v-puntos: conjunto de puntos que son 26-adyacentes entre ellos, pero no 18-adyacentes. Ver figura 2.6 (c). Fig. 2. 6. Clasificación de los puntos s, e y v En el primer proceso, el “adelgazamiento-primario”, se realizan tres escáneres de los puntos del objeto para clasificarlos como s-puntos, e-puntos o v-puntos. Para cada uno de estos grupos, se borrarán aquellos que sean puntos simples. Esto se repite en cada iteración hasta que no se pueden clasificar más puntos. En el último proceso, el “adelgazamiento-final”, que contiene una única iteración, se detecta aquellos puntos que forman una superficie errónea en el esqueleto-primario para ser eliminados cumpliendo las condiciones mencionadas anteriormente. En este proceso de adelgazamiento, Saha et al. emplean dos versiones de la imagen; una contiene los puntos binarios de la imagen anterior a la iteración actual, mientras que la segunda contiene la información de la imagen actual. Los puntos simples siempre se detectan en la versión de la imagen actual, mientras que las condiciones para preservar la topología se comprueban en la imagen de la iteración anterior. Esto mejora considerablemente la obtención del esqueleto, respetando la forma inicial del objeto sin añadir ramificaciones no deseadas al esqueleto, provenientes de pequeñas irregularidades en los contornos de la imagen. Ver figura 2.7 y figura 2.8. Fig. 2. 7. Imagen original Esqueleto ideal Esqueleto con excesivas ramificaciones Fig. 2. 8. Imagen original Esqueleto ideal Esqueleto con excesivas ramificaciones a) s-puntos b) e-puntos c) v-puntos Mª Isabel Mompeán Escartí Página 30 Otro posible enfoque para el algoritmo de adelgazamiento se describe en (Ma CM, Sonka M, 1996) que está basado en máscaras predefinidas de tamaño 3x3x3 que a priori contenían todas las combinaciones posibles para que un punto pudiera ser borrado. El método consistía en hacer coincidir el punto central de las máscaras con el punto que se estaba estudiando. Si el conjunto de sus 26-vecinos coincidía con cualquiera de ellas, el punto ya se podía borrar de la imagen sin alterar la topología del objeto. Sin embargo, como se explica detalladamente en (Wang T, Basu A, 2007), una de las configuraciones de estas máscaras, para casos muy concretos, no preservaba la topología del objeto ya que lo desconectaba. Así pues, una solución que da Wang y Basu fue modificar esta configuración en un conjunto de nuevas máscaras, que ahora sí, conservaban las condiciones de topología. A su vez, estos algoritmos de adelgazamiento se pueden dividir en subclases según la detección de los puntos simples. Estas son: métodos de adelgazamiento basados en direcciones, métodos de adelgazamiento secuenciales basados en subconjuntos y métodos de adelgazamiento completamente paralelos. Basados en direcciones: El método de adelgazamiento basado en direcciones consiste en eliminar los vóxeles desde una dirección particular, por ejemplo, norte, sur, este u oeste (Gong W, Bertrand G, 1990; Bertrand G, 1995; Lohou C, Bertrand G, 2004; Palagyi K, Kuba A, 1999). Estos métodos son sensibles en cuanto al orden de procesamiento de las direcciones, y puede causar que el esqueleto no esté centrado en el objeto debido a la orientación inicial escogida. Secuenciales basados en subconjuntos: Este método divide el espacio discreto del objeto en subconjuntos, llamados subespacios, y en cada subiteración sólo considera los vóxeles que pertenecen a un subespacio. El número de subespacios que se pueden utilizar en 3D es dos (Ma CM, Wan SY, 2001), cuatro (Ma CM et al., 2002) u ocho (Saha PK et al., 1997). Métodos completamente paralelos: Esta clase de algoritmos considera todos los puntos del contorno a borrar en una sola iteración (Ma CM, Sonka M, 1996; Borgefors G et al., 1999; Manzanera A et al., 1999). Aquí podemos observar que estas estrategias son más rápidas que las técnicas basadas en nsubiteraciones, ya que el número de iteraciones para obtener el resultado es menor (Perroton L, Aktouf Z, 1997). Sin embargo, el proceso de eliminación de puntos de la imagen se complica, y con el fin de mantener la topología del objeto, según Cornea el conjunto de vecinos que se deben inspeccionar para borrar un vóxel debe ser superior a los 26 vecinos locales, es decir, con máscaras de 5x5x5 como expone Palágyi (Palágyi K, 2008). Mª Isabel Mompeán Escartí Página 31 2.3.2 Transformada de la distancia (TD) La metodología basada en TD normalmente es sólo aplicable a imágenes binarias (Breu H et al, 1995), ya que estas imágenes presentan una clara distinción entre el objeto de estudio y el fondo de la imagen. Esto posibilita la detección de los contornos que es esencial en el proceso de cálculo de TD. El método TD se define para cada punto interior de un objeto y especifica la mínima distancia al punto más cercano del contorno del objeto. Existen distintas formas de calcular la TD dependiendo de las métricas de distancia utilizadas para determinar la distancia entre dos puntos. Las métricas de distancia son muy utilizadas en el procesamiento de imágenes. De entre ellas se consideran destacables las siguientes: Distancia Euclidea: La distancia Euclídea es la distancia más utilizada para el cálculo de la línea central, debido a la exactitud que se requiere en aplicaciones de procesamiento de imágenes médicas. Sean dos puntos p1 y p2 con coordenadas (x1, y1) y (x2, y2) entonces, la distancia Euclidea viene dada por la ecuación 2.1: (Ec. 2.1) Distancia "City Block": La distancia “City Block”, también conocida como “Manhattan”, consiste en que ir desde un punto a otro es sólo posible de forma horizontal y/o vertical pero nunca diagonal. Los movimientos diagonales, por lo tanto, no están permitidos. La distancia “City Block” viene dada por la ecuación 2.2: 𝑩 (Ec. 2.2) Distancia "Chebyshev": La distancia de “Chebyshev”, también conocida como “Chessboard” o tablero de ajedrez, permite los mismos movimientos que la distancia “City Block" pero, además, añade la posibilidad de realizar movimientos diagonales. Dichos movimientos diagonales contarán de la misma forma que uno horizontal o vertical. La distancia “Chebyshev” viene dada por la ecuación 2.3: 𝑩 (Ec. 2.3) La métrica de “City Block” y “Chebyshev” son una aproximación a la distancia Euclidea y consiguen computaciones más rápidas al calcular la TD. Sin embargo, usando la distancia Euclidea se consiguen esqueletos más robustos bajo rotaciones del objeto. Mª Isabel Mompeán Escartí Página 32 En un mapa de distancias calculado con el método TD, los puntos centrales del objeto presentan valores máximos de TD, por consiguiente, estos máximos locales pertenecerán al esqueleto y no podrán ser eliminados. Ver figura 2.9. Fig. 2. 9. a) Imagen original. b) Resultado tras aplicar TD sobre (a). (Rumpf M, Telea A, 2002) Al aplicar la TD sobre la figura 2.9 (a), obtenemos un mapa de distancias que da como resultado la figura 2.9 (b), donde podemos observar que los puntos más brillantes representan el esqueleto de la imagen original. Si representásemos todas las distancias de cada punto como una función, los picos que se observan representarían los puntos que se encuentran en el centro del objeto. La metodología basada en TD revisada en la literatura, se centra en encontrar estos puntos, candidatos para formar el esqueleto. Los métodos que se usan para encontrar estos puntos se basarán en los conceptos descritos a continuación.  Búsqueda de máximos locales controlados por parámetros de adelgazamiento. Gagvani y Silver (Gagvani N, Silver D, 1997) utilizan la media de la TD de los 26vecinos de un punto. Si la media de las distancias de los 26-vecinos es mayor o igual al valor de su propia distancia, dicho punto no pertenecerá al esqueleto. Sin embargo, el conjunto de puntos que se obtiene no es mínimo para poder extraer el esqueleto, así que acota dicho conjunto mediante un parámetro de adelgazamiento para formar el esqueleto final.  Ordenación previa de distancias. El método se basa en la eliminación secuencial, en orden ascendente, de los puntos del mapa de distancia hasta que no es posible seguir eliminando puntos. Durante este proceso debe tenerse en consideración que un punto sólo puede eliminarse si no afecta a las condiciones topológicas (Pudney C, 1998).  Cálculo del gradiente y la divergencia (Dimitrov P et al., 2003). Propone calcular el gradiente del mapa de distancias y aplicar una alternativa al teorema de la divergencia. Muy similar a este método es el descrito en (Bouix S et al., 2005), que se basa en el método propuesto por Siddiqi y Bouix (Siddiqi K, Bouix S, 2002), donde se utiliza la media del flujo del vector gradiente para distinguir los puntos del esqueleto combinando este método, con un procedimiento de adelgazamiento para conservar las propiedades topológicas. Página 39 3. Esqueletonización en 2D En el siguiente capítulo se comparan los distintos algoritmos de esqueletonización 2D enunciados en el apartado anterior. Estos algoritmos son: adelgazamiento, Diagramas de Voronoi (DV) y transformada de la distancia (TD). Tras ejecutarlos con dos imágenes distintas, se va a analizar los pros y contras de los resultados obtenidos con cada uno de ellos. Mª Isabel Mompeán Escartí Página 40 3.1 Algoritmos de esqueletonización 2D Partiendo de las imágenes originales 3.1 y 3.2, se va a realizar un estudio entre algoritmos basados en adelgazamiento, DV y TD. Estas imágenes presentan estructuras diferentes que mostrarán, para cada uno de los algoritmos, ventajas y desventajas al obtener sus respectivos esqueletos. Fig. 3. 1. Imagen binarizada de un corte de RM vascular de pulmón humano Fig. 3. 2. Imagen binarizada de un corte del cráneo humano adquirido mediante TAC. Mª Isabel Mompeán Escartí Página 41 ADELGAZAMIENTO El primero de los algoritmos probados está basado en adelgazamiento. Los algoritmos clasificados en el adelgazamiento de la forma se basan en ir eliminando sucesivamente puntos de la imagen desde los contornos de la misma, hasta llegar a un conjunto de puntos que se encuentran en el centro de la forma y contienen al esqueleto. Estos métodos como se ha comentado en el capitulo anterior, se pueden clasificar según su programación, en secuenciales o en paralelos. Los métodos de adelgazamiento secuencial, consideran capa a capa el conjunto de puntos a eliminar. Para ello, empiezan desde los límites de la imagen considerando sólo aquellos que cumplen la condición de puntos simples. Se considerará como “simples” aquellos puntos cuya eliminación no afecta la topología de la imagen. De manera inmediata se prevé que la conectividad también se preservará. Los métodos de adelgazamiento en paralelo, hacen una primera división de la imagen con el fin de realizar un borrado simultáneo del conjunto de puntos simples mediante iteraciones. En esta clase de algoritmos, el estudio de los puntos se basa únicamente en el resultado de la iteración anterior. Todos los puntos pueden ser examinados independientemente y de forma paralela en cada iteración. Para alcanzar el esqueleto de una imagen se deben eliminar tantos puntos como sea necesario hasta dejar sólo aquellos que se consideran puntos finales. Esta consideración es ambigua desde el punto de vista de su aplicación. Por ejemplo, la figura 3.3 muestra dos esqueletos para la imagen de la figura 3.2, obtenidos con versiones distintas de adelgazamiento. Como se observa en la imagen de la izquierda, está difiere de la imagen de la derecha en los puntos finales. Cada rama del esqueleto induce un punto extremo que se considera final. En ambas imágenes, los puntos marcados en amarillo son considerados puntos finales y, por ello, no se pueden eliminar. Mª Isabel Mompeán Escartí Página 42 Fig. 3. 3. Esqueletos de la imagen 3.2, obtenidos mediante versiones distintas de adelgazamiento. Los esqueletos obtenidos con estos algoritmos pueden estar influenciados erróneamente, al final de una rama, por la presencia de una protuberancia en los contornos de la imagen. Una posible mejora, sería aplicar heurísticas para clasificar los puntos del contorno como ruido y eliminarlos para no producir ramificaciones no deseadas. En la figura 3.4, se observan dos esqueletos para la imagen 3.1, calculados con versiones distintas de adelgazamiento. En la imagen de la derecha se observa una menor ramificación con respecto a la imagen de la izquierda. También se podrían eliminar las ramas no deseadas del esqueleto mediante algoritmos de podas. Fig. 3. 4. Esqueletos de la imagen 3.1 mediante distintas versiones del algoritmo basado en adelgazamiento, donde se observan ramificaciones no deseadas debido a irregularidades en el contorno de la imagen original. En amarillo están marcados los puntos finales de cada rama. Mª Isabel Mompeán Escartí Página 43 En general, cuando los píxeles se procesan de forma secuencial, no existe problema en preservar la conectividad. Por lo tanto, el requisito de preservación topológica se cumple para estos algoritmos. Sin embargo, conservar las propiedades geométricas es un problema más difícil. La principal dificultad es que para minimizar la complejidad del algoritmo es deseable considerar sólo vecinos locales, por ejemplo con tamaños de ventanas de 3x3. Estudiar dichas agrupaciones de vecinos no nos proporciona información global de la forma para poder distinguir entre ruido y los auténticos puntos finales. Además, la aplicación de ventanas de 3x3 de forma secuencial, podría producir un exceso de erosión o grave acortamiento de las ramas. Para evitar la excesiva erosión y la creación de falsos puntos finales al mismo tiempo, se necesitaría considerar información más global de la imagen. Esta idea puede permitir una mejor correlación entre la forma del esqueleto y el contorno exterior, que no se conseguiría mediante operaciones locales solamente. Por supuesto, esta consideración más global, llevaría a un aumento en el tiempo de cálculo y procedimientos más complejos. Para estudiar el comportamiento del algoritmo de adelgazamiento explicado en (Lam L. et al., 1992), se ha considerado la figura 3.5 que contiene 6 formas simples, todas ellas diferentes entre sí. Fig. 3. 5. Ejemplo de 6 formas distintas En cada forma de la figura 3.6 está dibujado su esqueleto correspondiente calculado mediante el algoritmo anterior. Además se muestra con un punto verde, los puntos extremos del esqueleto. Como se puede observar, las figuras son todas distintas, sin embargo, sus esqueletos son para todas ellas el mismo. El hecho de que distintas formas tengan un esqueleto idéntico no es apropiado. Debido a esta igualdad entre los esqueletos, se hace imposible la reconstrucción de la forma, ya que podría ser cualquiera de ellas la forma original. Así pues un algoritmo basado en los métodos de adelgazamiento, no garantiza conservar las propiedades geométricas ya que no es posible la reconstrucción exacta de la imagen original. 10 20 30 40 50 60 5 10 15 20 25 30 35 40 45 Mª Isabel Mompeán Escartí Página 44 Fig. 3. 6. Esqueleto de la imagen 3.5, obtenido mediante adelgazamiento, donde se muestran los puntos extremos mediante puntos verdes Este proyecto desea crear un algoritmo capaz de obtener el esqueleto de volúmenes de datos formados por imágenes médicas 3D, que se caracterizan por su gran tamaño y la elevada cantidad de puntos a procesar. Los métodos de adelgazamiento en una imagen 2D, presentan un tiempo de ejecución aceptable pero al convertir estos métodos para que sean capaces de resolver el esqueleto de una imagen 3D, su coste computacional aumenta considerablemente. Esto se debe, principalmente, al incremento de las numerosas condiciones que se deben implementar. Para solventar este problema se puede decidir paralelizar este tipo de algoritmos, de lo contrario el coste computacional es excesivamente elevado. Un algoritmo de adelgazamiento en paralelo, elimina simultáneamente un conjunto de puntos simples. Sin embargo, la eliminación simultánea de estos puntos puede alterar la topología de la imagen. Por este hecho, se deben agregar nuevas estrategias a las condiciones iniciales para procesar imágenes en 3D, lo que supone una complejidad añadida para mantener la topología de la imagen original. Mª Isabel Mompeán Escartí Página 45 DIAGRAMAS DE VORONOI (DV) Como se ha comentado en el capitulo anterior, el DV se forma a partir de la intersección de rectas perpendiculares que equidistan de dos puntos del contorno. El algoritmo aplicado crea un subgrafo, donde cada arista se traza como una recta perpendicular y equidistante a dos puntos del contorno de la imagen. Donde se intercepten estas rectas se formará un punto del esqueleto. Este método, como cualquier otro, es fácil de aplicar en formas relativamente simples, sin demasiadas protuberancias en el contorno de la forma o en imágenes que se pueden clasificar como tubulares. Esto se debe a que las paredes de una forma tubular son apreciablemente paralelas y se pueden trazar fácilmente líneas rectas que se interceptarán en un punto determinado dando una línea central limpia, sin ramificaciones y claramente visible como se observa en la figura 3.7. Fig. 3. 7. Esqueleto de la imagen 3.2 mediamente un algoritmo basado en DV Por el contrario, aun siendo una imagen tubular, es delicadamente curva y se necesitan muchos segmentos de rectas para aproximar el contorno de la imagen a su forma real (ver figura 3.9) y así, generar el DV para obtener el esqueleto que representa a la línea central de la imagen. Esto repercute en el tiempo de ejecución del algoritmo. Mª Isabel Mompeán Escartí Página 46 Fig. 3. 9. Representación del DV de la imagen 3.2 Otra de las limitaciones de este algoritmo reside en la intersección de las rectas que forman el DV. Especialmente relevante es cuando se trata con imágenes que presentan una gran cantidad de líneas convexas en su forma, tales como los vértices de un polígono, induciéndose en estos casos una rama adicional en el esqueleto final. Por ello, es muy difícil cumplir con los requisitos de invariabilidad y robustez mediante un algoritmo basado en DV. Al igual que en las metodologías de adelgazamiento, algunos algoritmos basados en DV resuelven dicho problema mediante técnicas de suavizado en los contornos de la imagen o técnicas de poda de las ramas adicionales en el esqueleto final. Fig. 3. 11. Representación del DV de la imagen original 3.1 Fig. 3. 8. Detalle del DV de la figura 3.9, donde se observan los vértices (puntos rojos), centrados en cada celda del diagrama. Fig. 3. 10. Detalle de la figura 3.11, donde se observan los vértices (puntos rojos), centrados en cada celda del diagrama. Mª Isabel Mompeán Escartí Página 47 Como se puede observar en la imagen 3.12, al trazar el esqueleto partiendo del DV, se crean numerosas ramas adicionales. Se puede ver que las partes del esqueleto que se encuentran centradas en el interior de la forma son menos sensibles a los cambios en los puntos del contorno, mientras que si nos alejamos del centro y llegamos a los contornos, las ramas son influenciadas por la variabilidad de estos. Por lo tanto, se obtiene una representación que no se corresponde con la forma real de la imagen. Fig. 3. 12. Esqueleto resultado de la imagen 3.1 mediante el procesado del DV. En el caso del espacio 3D, la aproximación adecuada de los contornos de un volumen a una forma poligonal concreta, se convierte en crucial para la complejidad y la topología del esqueleto. Por ejemplo, una aproximación mediante superficies rectangulares, en casos donde los contornos no varían demasiado, daría una buena aproximación de la forma. Sin embargo, para volúmenes muy irregulares con superficies curvas sería mejor una aproximación mediante superficies triangulares gracias al ajuste que se puede establecer con los ángulos de sus vértices. Mª Isabel Mompeán Escartí Página 48 TRANSFORMADA DE LA DISTANCIA (TD) Por último, los algoritmos basados en el método TD computan la mínima distancia Euclídea de cada punto de la imagen a un punto del contorno más cercano. Como se ha visto hasta ahora, una de las características más importantes que tiene que tener un algoritmo de esqueletonización es la propiedad de conservación de la topología. Los algoritmos basados en TD cumplen dicha propiedad. Estos algoritmos crean un mapa de distancias para todos los puntos de la imagen. Las figuras 3.13 y 3.14 representan los mapas de distancia de sus imágenes originales 3.1 y 3.2 respectivamente. En dichas imágenes, se puede observar mediante una escala de grises, con rango normalizado entre [0-1], como los puntos centrales de la imagen, los más lejanos al contorno, presentan valores más altos, es decir, más cercanos al color blanco, y como los puntos más oscuros son los puntos que están más cerca del contorno y, por lo tanto, con menor valor de distancia. Estos valores de distintas intensidades, hacen referencia a la distancia Euclídea del punto a su contorno más cercano. Fig. 3. 13. Representación en escala de grises del mapa de distancias de la imagen 3.1 Fig. 3. 14. Representación en escala de grises del mapa de distancias de la imagen 3.2 50 100 150 200 250 300 50 100 150 200 250 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 50 100 150 200 250 300 350 400 450 50 100 150 200 250 300 350 400 450 500 550 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Mª Isabel Mompeán Escartí Página 55 DETECCIÓN DE LOS PASOS POR 0 EN LOS TRES PLANOS X, Y y Z MAPA DE DISTANCIAS GRADIENTE ∑ DERIVADAS PARCIALES x PROCESAR ESQUELETO PRIMARIO A B E D C HG F I J Fig. 4. 1. Diagrama del algoritmo de esqueletonización 3D. Mª Isabel Mompeán Escartí Página 56 Cálculo del mapa de distancia: TD Partiendo del volumen original, mediante TD se obtiene una matriz 3D de distancias. Ésta contiene, para cada vóxel del volumen, la distancia Euclidea al vóxel más cercano del contorno. Expresado formalmente, la TD (ecuación 4.1) de un vóxel p = {x1, y1, z1} se define como (Ec. 4.1) Donde dt es la distancia Euclídea (ecuación 4.2) desde el vóxel (x1, y1, z1) a cualquier vóxel (x2, y2, z2) que pertenece al conjunto de vóxeles del contorno OC. (Ec. 4.2) Además del conjunto de vóxeles del contorno OC, clasificamos la imagen como el conjunto de vóxeles que definen en el objeto O que se quiere esqueletonizar (OC ⊂O) y como el conjunto de vóxeles que forman el fondo de la imagen. Entonces el resultado de TDP se clasifica como: (Ec. 4.3) La figura 4.2, muestra a modo de ejemplo, el mapa de distancias de los puntos del contorno de la imagen, calculado con la mínima distancia Euclidea. Como podemos comprobar la distancia en los contornos es 1, excepto para los puntos del contorno que son 26-vecinos a su punto más cercano en . Fig. 4. 2. Ejemplo de TD en los contornos. 000000000000000 0110000000001 1 1 01 1,4 1 01 1 1 1 1 1 1 1 0 0 01 1,4 1 1,4 1 0 0 0 01 1 1,4 1,4 1 0 0 0001 1 1,4 1 0 0 000001 1,4 1 0 00001 1,4 1,4 1 1 1,4 1,4 1 0001 1,4 1 001 1 1 1 1 0 001 1,4 1 0 0 0 0 0 0 0 0 001 1 0 0 0 0 0 0 0 0 001 1,4 1 0 0 0 0 0 0 0 001 1 1 1,4 1,4 1 0 0 0 0 0 0 000001110 0 0 0 0 0 0 000000000000000 Mª Isabel Mompeán Escartí Página 57 Cálculo de los gradientes: Mediante mapa de distancia Tras obtener el mapa de distancias TD, se obtienen las derivadas parciales en cada dirección del espacio, mediante el cálculo del gradiente. El gradiente de una función se define como el campo vectorial cuyas funciones direccionales son las derivadas parciales del campo escalar, esto es: (Ec. 4.4) El vector gradiente de evaluado en un punto genérico del dominio de , , indica la dirección en la cual el campo varía más rápidamente y su módulo representa el ritmo de variación de en la dirección de dicho vector gradiente. Dado un campo escalar definido en el espacio, análogamente a la ecuación 4.4, se tendrá que: (Ec. 4.5) Donde corresponde a la diferencia en la dirección x (horizontal), corresponde a la diferencia en la dirección y (vertical) y corresponde a la diferencia en la dirección Z (profundidad). De forma geométrica, el gradiente es un vector normal (perpendicular) a la curva de nivel en el punto sometido a estudio. Aplicando la función gradiente a la matriz 3D de distancias, obtenemos 3 matrices (ver figura 4.3), una para cada dirección del espacio, figura A (dirección eje X), figura B (dirección eje Y), y figura C (dirección eje Z). C BA Fig. 4. 3. Representación del flujo del vector gradiente aplicado en los tres planos X (A), Y (B) y Z (C) en una sección de trabécula. Mª Isabel Mompeán Escartí Página 58 Si interpretamos los resultados del gradiente, al aplicarlo sobre el mapa de distancias, lo que obtenemos es como varían las distancias de la imagen en cada plano con respecto a la forma de sus contornos y, es más, hacia donde varían. Dicho de otro modo, en cada plano del espacio, la componente direccional que se está estudiando se anulará en un punto, que se sitúa a la máxima distancia de los contornos, es decir en el centro de la imagen. Así pues, combinando los tres planos en una única matriz, como la multiplicación lógica de los 0’s de cada componente, se obtiene una imagen única, ver figura 4.1 (I), que contiene los máximos de las derivadas parciales de cada dirección del espacio. Como último paso de este proceso, la imagen obtenida se multiplica por la imagen inicial originando el esqueleto inicial Si, que se muestra en la imagen de la figura 4.4 (A) como resultado del proceso inicial de esqueletonización. Procesado del esqueleto inicial: En este momento se obtiene un primer esqueleto Si, centrado en la imagen, pero que no es definitivo. Antes de obtener el esqueleto final Sf, se realiza un paso intermedio para hallar el esqueleto adelgazado de la imagen 3D de la que partimos. Como se observa en la figura 4.4, dado que tenemos el mapa de distancias (B) de la imagen original y el esqueleto Si (A), al multiplicar ambas imágenes, obtenemos en cada punto del esqueleto Si, la distancia a los contornos (C). x B A C Fig. 4. 4. Esquema del cálculo de las distancias del esqueleto inicial. Donde (A) es el esqueleto Si, (B) es el mapa de distancias de la imagen original y (C) el mapa de distancias del esqueleto Si. Mª Isabel Mompeán Escartí Página 59 Combinación TD con adelgazamiento: Finalmente, tras hallar el mapa del esqueleto inicial, (figura 4.6, A), se procesa el mismo para obtener un esqueleto adelgazado. Este paso es una combinación de la técnica de mapas de distancias TD con la técnica de adelgazamiento. La idea principal consiste en evaluar los tramos gruesos del esqueleto Si, y desde los límites del mismo, ir eliminando todos aquellos vóxeles que no varían la topología de la imagen original. En este sentido, es importante la utilización del mapa de distancias, y no una imagen binaria de ésta. Esto es debido a la necesidad de conocer de ante mano dónde se centraría el esqueleto definitivo Sf en la imagen que se está tratando. El algoritmo aplicado para esta sección se muestra a continuación mediante pseudocódigo: Esqueleto procesado desde (el vóxel x ∈ Si con la menor TD) … … hasta (el vóxel y ∈ Si con la máxima TD) hacer si (x es simple) entonces Borrar x en Sf; fin si Fig. 4. 5. Esqueleto final Sf, como resultado de adelgazar la imagen 4.4 (C) mediante el método descrito. Mª Isabel Mompeán Escartí Página 60 Detección de extremos, uniones y vóxeles aislados en el esqueleto Un proceso complementario al de esqueletonización, consiste en detectar cuatro tipos de puntos en el esqueleto. Estos tipos son: puntos extremos, son los puntos finales de cada rama; puntos de curva, que son los que forman cada rama; puntos de unión entre ramas; y por último, los puntos aislados, que son un caso particular de los puntos extremos. En la figura 4.6, se muestra esta clasificación. Fig. 4. 6. Representación de la clasificación de los vóxeles en modelo real de trabecular El método consiste en estudiar únicamente el conjunto de los 26-vecinos para cada vóxel del esqueleto. De esta forma, la función desarrollada analiza la conectividad de los 3x3x3 vecinos para cada vóxel PSf del esqueleto de la siguiente forma:  Si PSf no tiene ningún 26-vecino, el vóxel se clasifica como Aislado.  Si PSf tiene un sólo 26-vecino, el vóxel se clasifica como Extremo.  Si PSf tiene exactamente dos 26-vecinos, el vóxel se clasifica como Curva.  Si PSf tiene tres o más 26-vecinos, el vóxel se clasifica como Unión. Fig. 4. 7. Aplicación de la función descrita sobre una muestra de la estructura trabecular tras el proceso de esqueletonización 3D. Se puede observar, la detección de los extremos de cada rama y las uniones entre ellas. Aislado Unión Extremo Curva Mª Isabel Mompeán Escartí Página 61 A continuación se muestran las representaciones gráficas en 3D de hueso trabecular y su esqueleto como resultado del algoritmo descrito en el presente proyecto. Mediante técnicas de visualización del volumen, se pueden conseguir representaciones 3D diversas que aportan información visual muy variada. Fig. 4. 8. Imágenes obtenidas por el mismo método de reconstrucción 3D de la microarquitectura del hueso trabecular. El conjunto de datos está formado por una colección de cortes de IRM que se apilan en profundidad en una matriz 3D. Mediante la extracción de la capa superficial y el suavizado de la misma, combinada con la extracción de las capas interiores de las trabéculas y aplicando un mapa de color y focos de luz sobre la geometría dando la sensación de volumen, se puede visualizar el volumen completo del tejido trabecular (a). En (b), se muestra el esqueleto de (a) al modelar con la misma técnica, la matriz resultante de aplicar el algoritmo de esqueletonización 3D desarrollado en este proyecto. Fig. 4. 9. Se muestra una única representación por la fusión de ambas imágenes (figura 4.8), aplicando sobre el volumen original del hueso trabecular (a), una capa de transparencia para apreciar en el interior de las trabéculas, el esqueleto representado por la línea central (azul). (a) (b) Mª Isabel Mompeán Escartí Página 62 Otra forma de visualización del volumen, consiste en representar la matriz de datos 3D mediante la técnica de elementos finitos (EF). En el caso particular del hueso trabecular, esta técnica consiste en definir la estructura del hueso, como una red de pequeños elementos del orden de micras, por ejemplo, tetraedros o ‘cubos’ y las conexiones entre sus nodos. El desarrollo del método se explica en (Alberich-Bayarri A, 2010). Fig. 4. 10. Reconstrucción del hueso trabecular como una red de EF basada en elementos de cubo. Los cubos se muestran apilados formando la estructura trabecular. Fig. 4. 11. Reconstrucción de la estructura trabecular tras el proceso de esqueletonización 3D sobre la imagen de la figura 4.10. Modelando el esqueleto con EF basados en cubos se puede observar el grosor de la línea central de un vóxel como máximo. El numeró de cubos en la imagen corresponde con el número de vóxeles que componen el esqueleto del hueso trabecular. Mª Isabel Mompeán Escartí Página 63 4.2 Aplicación a hueso trabecular Como en muchas otras ciencias, al transformar el objeto de estudio en su forma más simple, se consigue observar la esencia del mismo, mejorando con ello su comprensión, al centrarse únicamente en la naturaleza de su forma. Esto mismo sucede en las estructuras trabeculares del tejido óseo. Debido a la complejidad de su organización trabecular, el ideal reside en obtener un modelo simplificado que represente su morfología y la disposición de las trabéculas en el interior de su estructura. Al reducir la forma hasta obtener su esqueleto 3D, se consigue este propósito. El hueso trabecular se compone de un entramado de placas y barras. Ver figura 4.12. A través del esqueleto se puede observar esta estructura de la red. Muchos estudios han demostrado que la pérdida de masa ósea está acompañada por los cambios en la arquitectura que se refieren tanto a la magnitud como a la topología de la red. Por ello, la cuantificación de la composición estructural del hueso trabecular es de principal interés (Saha PK et al., 2000). Fig. 4. 12. Representación 3D del tejido trabecular, en base a una red de placas y barras. Las propiedades del tejido óseo trabecular pueden ser analizadas a partir de una gran variedad de parámetros, extraídos mediante adquisiciones IRM in vivo con una resolución espacial elevada. Se ha reconocido que la competencia mecánica del hueso trabecular, además de medirse en función de los parámetros morfológicos, grosor y densidad de masa ósea, (Majumdar S et al., 1997; Laib A et al., 2002), también se relaciona con la estructura de su arquitectura (Feldkamp LA et al., 1989), y la orientación de las trabéculas e índices relacionados con la elasticidad (Gomberg BR et al., 2003). El estudio combinado de estas tres clases de parámetros (Wehrli FW, 2007) da una nueva visión de la calidad del hueso trabecular en comparación con el cálculo de las características morfológicas por sí solas (Pothuaud L et al., 2002). Placas Barras Mª Isabel Mompeán Escartí Página 64 Gracias al modelado del hueso trabecular y su representación esqueletonizada 3D, se puede obtener cuantitativamente una caracterización estructural del hueso trabecular en base a sus parámetros morfológicos. Los estudios de Laib et al. (Laib A et al., 2002) y Alberich-Bayarri (Alberich-Bayarri A, 2010), muestran que los parámetros morfológicos más relevantes en la caracterización del hueso trabecular son: la relación del volumen óseo frente al volumen total (BV/TV), el número trabecular (Tb.N), la separación trabecular media (Tb.Sp), y el grosor trabecular medio (Tb.Th). Este último parámetro Tb.Th es muy fácil de calcular gracias a la metodología de esqueletonización elegida en el presente proyecto. Con el algoritmo basado en TD, al proporcionar para todos los vóxeles la mínima distancia a los contornos, posibilita obtener a su vez la distancia de los vóxeles que forman el esqueleto, y si duplicamos esta distancia, obtenemos el grosor trabecular en ese punto. Esto se muestra en el diagrama de la figura 4.11. Fig. 4. 13. Representación de una sección de trabécula (gris) y su esqueleto (negro), para el cálculo del grosor trabecular Tb.th. Mediante TD obtenemos la distancia a los contornos. En cada vóxel del esqueleto tenemos la distancia mínima, de forma que duplicando esta distancia, obtenemos el grosor trabecular local. El segundo grupo de parámetros a estudiar está formado por aquellos basados en el análisis topológico. Dicho análisis versa sobre la disposición de la red ósea, con el fin de caracterizar las propiedades topológicas del hueso trabecular. En esta línea, Feldkamp et al (Feldkamp LA et al., 1989) mostró que la conectividad de la red trabecular puede expresarse en términos de sus características topológicas. El adelgazamiento uniforme de la estructura trabecular, afecta a los parámetros de escala, pero esto no implica un cambio en la topología de la estructura, siempre y cuando no se rompan las conexiones o no se perfore una placa trabecular (Wehrli FW, 2007). Enfermedades tales como la osteoporosis, se caracteriza por ser un trastorno que resulta de la erosión ósea, e implican tanto la pérdida neta de masa ósea, como el deterioro de su arquitectura (OMS, 1994). d 2d Mª Isabel Mompeán Escartí Página 71 El nuevo algoritmo de esqueletonización 3D se ha comparado en tiempos de ejecución con un método basado únicamente en adelgazamiento descrito en (Alberich-Bayarri A, 2010). Como se ha argumentado en los capítulos anteriores, los métodos basados en adelgazamiento, procesan cada vóxel del volumen secuencialmente; aplicando una amplia gama de condiciones a fin de preservar la topología de la imagen original. Debido a esta característica, para imágenes con un gran volumen de datos, el tiempo aumenta considerablemente. En la siguiente gráfica (ver figura 5.3.) se muestran los tiempos de ejecución anteriormente presentados para el algoritmo que se ha desarrollado, comparado con los tiempos de ejecución obtenidos tras aplicar el algoritmo de adelgazamiento de AlberichBayarri sobre el mismo conjunto de imágenes de la tabla 5.1. Fig. 5. 2. Comparativa de los tiempos de ejecución (en segundos) entre el algoritmo desarrollado basado en TD y el basado en adelgazamiento. Como muestra la gráfica, para tamaños de imágenes de un volumen pequeño y mediano, los tiempos de ejecución son similares. Ahora bien, conforme aumenta el tamaño de las imágenes, el tiempo de ejecución entre ambos algoritmos es muy distinto. En el citado escenario, es claramente apreciable la reducción del tiempo de ejecución obtenida mediante la aplicación del nuevo algoritmo basado en TD. En resumen, se puede concluir que el algoritmo de esqueletonización 3D presentado en este proyecto, ha obtenido unos resultados muy buenos en la extracción del esqueleto en imágenes 3D de estructuras trabeculares obtenidas mediante RM. Haciendo esto, además, con una mejora notable en los tiempos de ejecución gracias a la metodología escogida. 4,36 19,41 150,78 1.208,68 3.201,38 8.953,44 0,02 0,06 1,48 48,37 346,57 2.073,55 -2.000,00 0,00 2.000,00 4.000,00 6.000,00 8.000,00 10.000,00 Tiempo (Seg.) Comparativa Tiempos de Ejecución Thinning Algoritm Nuevo algoritmo de Esqueletonización 3D Algoritmo de adelgazamiento de Alberich-Bayarri Mª Isabel Mompeán Escartí Página 72 6. Conclusiones Teniendo en cuenta las expectativas y objetivos establecidos en este trabajo, podemos extraer las siguientes conclusiones:  Las diferentes metodologías estudiadas en este proyecto han demostrado su viabilidad y aplicabilidad para calcular el esqueleto del hueso trabecular. La metodología basada en TD permite la representación detallada de estructuras complejas. Para el caso concreto del tejido trabecular, en el cual encontramos estructuras muy complejas que necesitan de una representación lo más exacta posible, el método más apropiado será TD frente a los métodos basados en adelgazamiento y DV que obtienen un esqueleto simplificado de la estructura a la que representan.  Se ha conseguido utilizar un método de alta eficiencia computacional como es la TD para realizar la esqueletonización 3D de volúmenes de datos formados por imágenes médicas 3D, que se caracterizan por su gran tamaño y la elevada cantidad de puntos a procesar. Es importante remarcar que la aplicación de metodologías de adelgazamiento o DV sobre estas estructuras tan características supondría un elevado coste computacional.  La metodología elegida puede ser aplicada de manera eficiente a estructuras óseas extraídas a partir del procesado de imágenes de pacientes, con el fin de caracterizar su forma y estructura volumétrica. Además, garantiza que el esqueleto resultante esté centrado en la estructura, por la definición de la línea central como los máximos locales del mapa de distancias. Como consecuencia, se obtuvo un esqueleto que representaba fehacientemente las características topológicas de las trabéculas del hueso.  El algoritmo implementado es completamente automático y no necesita la interacción con el usuario, frente a otros algoritmos que requieren la selección de al menos un punto a considerar de la línea central.  En último lugar se ha conseguido una gran mejora en tiempo de ejecución gracias a la metodología basada en TD. Se han obtenido tiempos de ejecución significativamente inferiores a los que se conseguían con un algoritmo tradicional de esqueletonización, basado en adelgazamiento, que se utilizaba hasta ahora en aplicaciones orientadas al análisis de hueso trabecular.  Finalmente, con los resultados de este proyecto, se conseguirá obtener estructuras esqueletonizadas del hueso trabecular, siendo la metodología extensible a otras áreas anatómicas y por tanto, aportando una mejora relevante en el campo de la extracción de biomarcadores de imagen a partir del postproceso de imágenes médicas. Mª Isabel Mompeán Escartí Página 73 7. Futuras líneas de investigación y recomendaciones Durante el presente proyecto, gracias a la investigación realizada sobre el estado del arte, y a través de la aplicación de las soluciones desarrolladas a casos prácticos, se han identificado diferentes temas que podrían resultar de gran interés para la materia. A continuación se destacan las más relevantes. 1) Distinción de placas y barras del hueso trabecular para representar el esqueleto como superficies y curvas. El algoritmo presentado en este proyecto extrae un esqueleto que representa la línea central de la estructura del hueso trabecular. Ello es necesario para caracterizar las propiedades morfológicas y poder cuantificar parámetros tales como el grosor trabecular, la separación media entre trabéculas o el índice trabecular. Sin embargo, un modelo de la estructura trabecular que represente además, las placas de hueso como superficies 2D en el esqueleto y las barras (uniones entre placas) como curvas 1D en el mismo, sería ideal para realizar un estudio en profundidad de la caracterización del hueso trabecular. Esta representación del esqueleto solo es aplicable en imágenes 3D. La pérdida de masa ósea es un factor importante para evaluar la resistencia del hueso, pero no es determinante para establecer el riesgo de fractura. Por este motivo, es esencial detectar la naturaleza de la perdida de densidad ósea por su deterioro microestructural. Los cambios estructurales se producen por la reducción de la densidad de las placas de hueso, convirtiéndose en barras o incluso llegando a crear desconexiones. Así pues, el estudio de la relación placas vs. barras representado en el esqueleto como superficies vs. curvas, sería de gran ayuda para el seguimiento de enfermedades que afectan directamente a la estructura del hueso trabecular y las respuestas a sus tratamientos. 2) Conversión del esqueleto en EF basados en barras. Esta idea consiste en representar el esqueleto del hueso trabecular como una red simple de elementos de barras rectas. Para ello, es necesario transformar un conjunto de vóxeles alineados formando la pared de una trabécula en un único elemento que será una barra recta con las propiedades geométricas de la trabécula que representa. Su representación simplificada y la disminución del número de elementos, podrían reducir significativamente el tiempo de computación. Con este modelo de representación, se podría tratar y modificar por separado cada elemento permitiendo nuevos estudios topológicos. A su vez, podría simplificar la evaluación de las simulaciones con el fin de obtener una mejor comprensión del comportamiento mecánico del hueso trabecular. Mª Isabel Mompeán Escartí Página 74 3) Desarrollo y aplicación sobre otros tejidos del cuerpo humano Gracias al avance en nuevas técnicas de obtención de imagen médica, que permiten la adquisición in vivo con resoluciones espaciales elevadas, se puede visualizar nuevas imágenes de tejidos del cuerpo humano de gran complejidad microscópica con un gran detalle microestructural. La gran complejidad de estas nuevas imágenes justifica la investigación de herramientas software, entre ellas, los procesos de esqueletonización que permitan la descripción de la imagen para el diagnóstico y seguimiento de situaciones patológicas. A continuación, se comentan algunas de las regiones más complejas de la anatomía humana y la posible aplicación del proceso de esqueletonización para la extracción y cuantificación de parámetros morfológicos. El árbol de las vías respiratorias junto con los bronquios, presentan una estructura ramificada, por ello la extracción de la línea central representaría un modelo simplificado de su estructura. Con un proceso de segmentación del esqueleto, se podrían clasificar las ramas principales, los bronquios y la ramificación de los bronquiolos, para extraer características morfológicas y estudiar algunas de las patologías que afectan al sistema respiratorio. Entre las aplicaciones del tejido vascular, cabe destacar el estudio de imágenes 3D vasculares coronarias. El desarrollo de un método capaz de proporcionar una reconstrucción exacta y completa de la red vascular sería interesante en simulaciones de la hemodinámica coronaria, la autorregulación o el suministro de oxígeno al tejido. Además, se podría identificar las subredes de la malla vascular mediante la detección automática de todas las raíces y el seguimiento de la ramificación de los vasos que forman cada subred. Por otro lado, también se podría aplicar un método de segmentación y esqueletonización en imágenes 3D de la arteria carótida. La aplicación de un proceso de detección y segmentación de la arteria carótida interna, tras un algoritmo de esqueletonización basado en TD, podría ser útil en estudios que evaluaran el grado de estenosis, por medio de la distancia de la línea central a la pared del vaso, para cuantificar el riesgo de accidente cerebro-vascular en imágenes in vivo obtenidas mediante angiografía por RM o TC. De forma similar, se puede proceder para analizar imágenes 3D de la aorta y ramas iliacas, representando un modelo de los vasos como la línea central a los contornos activos, calculando el diámetro de estos y detectando la bifurcación de las ramas iliacas. Página 75 Referencias Abdel-Hami GH, Yang YH. Multiresolution Skeletonization: An Electrostatic Field-Based Approach. Conference Proceedings IEEE Image Processing 1994. Austin, TX. 1994; 1:949-953. Alberich-Bayarri A. In vivo morphometric and mechanical characterization of trabecular bone from high resolution magnetic resonance imaging [tesis]. Valencia, España: Universidad Politécnica de Valencia; 2010. Amenta N, Choi S, Kol-luri RK. The power crust, unions of balls, and the medial axis transform. Computational Geometry: Theory and Applications. 2000; 19: 127-153. Assmus A. Early History of X Rays. Beam Line. 1995; 25(2): 10-24. Attali D, Boissonnat JD, Edelsbrunner H. Stability and Computation of the Medial Axes. A State of the Art Report. En: Möller T, Hamann B Russell RD. Editors. Mathematical Foundations of Scientific Visualization, Computer Graphics, and Massive Data Exploration. Berlin: SpringerVerlag; 2009. 109-126. Bertrand G. A Parallel thinning algorithm for medial surfaces. Pattern Recognition Letter. 1995; 16 (9): 979-986 Bitter I, Kaufman AE, Sato M. Penalized-distance volumetric skeleton algorithm. IEEE Transactions on Visualization and Computer Graphics. 2001; 7(3): 195-206. Blum H. A transformation for extracting new descriptors of shape. En: Whaten-Dunn W, editor. Models for Perception of Speech and Visual Form. Cambridge, MA: MIT Press; 1967: 362-380. Borgefors G, Nyström I, Sanniti di Baja G. Computing Skeletons in Three Dimensions. Pattern Recognition. 1999; 32(7):1225-1236. Bouix S, Siddiqi K, Tannenbaum A. Flux driven automatic centerline extraction. Medical Image Analysis. 2005; 9(3):209-21. Breu H, Gil J, Kirkpatrick D, Werman M. Linear Time Euclidean Distance Transform Algorithms. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1995; 17(5): 529-533. Brunner D, Brunnett G. Divergence Analysis of Discrete 2-d Shapes. Journal of WSCG 2010: 2010: 57-63. Buckland-Wright JC. Quantitative radiography of osteoarthritis. Annals of the Rheumatic Diseases. 1994; 53: 268-275. Cárdenes R, Bogunovic H, Frangi AF. Fast 3d centerline computation for tubular structures by front collapsing and fast marching . Conference Proceedings IEEE 17th International Conference on Image Processing, Septiembre 2010; Hong Kong. 2010: 4109-4112. Cornea ND, Silver D, Min P. Curve-skeleton properties, applications, and algorithms. IEEE Transactions on Visualization and Computer Graphics. 2007; 13(3): 530-548. Página 76 Costa Carballo CM. Computacion Medica: Pasado, Present y Futuro. Revista General de Información y Documentación, Servicio de Publicaciones Universidad Complutense, Madrid. 1997; 7(1): 63-92. Dequeker J, Mokassa L, Aerssens J, Boonen S. Bone density and local growth factors in generalized osteoarthritis. Microsc Res Tech. 1997; 37: 358-71. Desco M, Vaquro JJ. Más de un siglo de imagen médica. Arbor CLXXVII. 2004; 698: 337-364. Dimitrov P, Damon JN, Siddiqi K. Flux invariants for shape. Conference Proceedings Computer Vision and Pattern Recognition 2003. 2003: 835-841. Feldkamp LA, Goldstein SA, Parfitt AM, Jesion G, Kleerekoper M. The direct examination of three-dimensional bone architecture in vitro by computed tomography. Journal of Bone and Mineral Research 1989; 4: 3-11. Fernández-Tresguerres Hernández-Gil I, Alobera Gracia MA, Canto Pingarrón M, Blanco Jerez L. Physiological bases of bone regeneration I. Histology and physiology of bone tissue. Med Oral Patol Oral Cir Bucal. 2006; 11(1): E47-51. Ferreirós Domínguez J. Impacto de la imagen por resonancia magnética (IRM) en la práctica médica. An. R. Acad. Nac. Farm. 2004; 70: 33-39. Gagvani N, Silver D. Animating Volumetric Models. Academic Press Professional. 2001; 63(6): 443-458. Gagvani N, Silver D. Parameter Controlled Skeletonization of Three Dimensional Objects. Technical Report CAIP-TR-216, Computer Aids For Industrial Productivity Center, Rutgers University, Piscataway, NJ, 1997. Gomberg BR, Saha PK, Song HK, Hwang, SN, Wehrli FW. Topological Analysis of Trabecular Bone MR Images. IEEE Transactions on Medical Imaging. 2000; 19(3): 166-174. Gomberg BR, Saha PK, Wehrli FW. Topology-based orientationanalysis of trabecular bone networks. Med Phys. 2003; 30: 158-168. Gong W, Bertrand G. A simple parallel 3D thinning algorithm, in: International Conference on Pattern Recognition, 1990; 188-190. González Salamea CG, Gustavo Molina P. La Informática Médica y los Sistemas de Información. 2003 [citado 17 Septiembre 2011]. Disponible en: http://medicinadefamiliares.cl/trabajos.htm Grigorishin T, Yang YH. Skeletonization: An Electrostatic Field-Based Approach. Pattern Analysis and Applications. 1998; 1: 163-177. Hassouna MS, Farag AA. Variational curve skeletons using gradient vector flow. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2009; 31(12): 2257–2274. Página 77 Hildebrand T, Laib A, Muller R, Dequeker J, Ruegsegger P. Direct three-dimensional morphometric analysis of human cancellous bone: microstructural data from spine, femur, iliac crest, and calcaneus. Journal of Bone and Mineral Research. 1999; 14(7): 1167-1174. Hsieh J. Computed tomography: principles, design, artifacts, and recent advances. Bellingham, Washington: SPIE Press; 2003. Kass M, Witkin A, Terzopoulos D. Snakes: active contours models. International Journal of Computer Vision. 1987; 1: 321-331. Kong TY, Rosenfeld A. Digital topology: Introduction and survey. Computer Vision, Graphics and Image Processing. 1989; 48(3): 357-393. Laib A, Newitt DC, Lu Y, Majumdar S. New model-independent measuresof trabecular bone structure applied to in vivo high resolution MR images.Osteoporosis International 2002; 13: 130-136. Lam L, Lee SW, Suen CY. Thinning Methodologies-A Comprehensive Survey. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1992; 14(9): 869-885. Latecki LJ, Li QN, Bai X, Liu WY. Skeletonization using SSM of the Distance Transform. Conference IEEE Image Processing 2007; 2007: 349-352. Lieutier A. Any Open Bounded Subset of R Has the Same Homotopy Type than Its Medial Axis. Computer-Aided Design. 2004; 36(11): 1029-1046. Lohou C, Bertrand G. A 3D 12-Subiteration Thinning Algorithm Based on P-Simple Points. Discrete Applied Math. 2004; 139: 171-195. Ma CM, Sonka M. A Fully Parallel 3D Thinning Algorithm and Its Applications. Computer Vision and Image Understanding. 1996; 64(3): 420-433. Ma CM, Wan SY, Lee JD. Three-Dimensional Topology Preserving Reduction on the 4-Subfields. IEEE Transactions Pattern Analysis and Machine Intelligence. 2002; 24(12): 1594. Ma CM, Wan SY. A Medial-Surface Oriented 3D Two-Subfield Thinning Algorithm. Pattern Recognition Letters. 2001; 22: 1439. Ma WC, Wu FC, Ouhyoung M. Skeleton Extraction of 3D Objects with Radial Basis Functions. En: Leif Kobbelt, Peter Schröder, Hugues Hoppe, editors. Conference Proceedings IEEE Shape Modeling and Applications 2003; Seoul, Korea. 2003: 207-215 Majumdar S, Genant HK, Grampp S, Newitt DC, Truong VH, Lin JC, et al. Correlation of trabecular bone structure with age, bone, mineral density, and osteoporotic status: in vivo studies in the distal radius using high-resolution magnetic resonance imaging. Journal of Bone and Mineral Research 1997; 12: 111-118. Manzanera A, Bernard T, Preteux F, Longuet B. A Unified Mathematical Framework for a Compact and Fully Parallel n-D Skeletonization Procedure. En: Latecki LJ, Melter RA, Mount DM, Wu AY. SPIE Proceedings Vision Geometry VIII 1999; Denver, USA. 1999; 3811: 57-68. Página 78 Morgenthaler DG. Three dimensional simple points: serial erosion, parallel thinning, and skeletonization. TR-1005, Computer Vision Laboratory, Computer Science Center, Univ. of Maryland, College Park. 1981. Mosekilde L. Consequences of the remodelling process for vertebraltrabecular bone structure: a scanning electron microscopy study (uncoupling ofunloaded structures). Bone and Mineral 1990; 10: 13-35. Naf M, Székely G, Kikinis R, Shenton M, Kubler G. 3D Voronoi skeletons and their usage for the characterization and recognition of 3D organ shape. Comput. Vision Graphics Image Process. 1997; 66(2): 147-161. Natural Medicine Seminars. Human Anatomy & Physiology Chapter 6 Bones & Skeletal Tissue. [citado 13 Marzo 2012]. Disponible en: http://www.naturalmedicineseminars.com NystroKm I, Borgefors G. Synthesising objects and scenes using the reverse distance transformation in 2D and 3D, in: Braccini C, DeFloriani L, Vernazza G, editors. IEEE Proceedings International Conference on Image Analysis and Processing 1995; Washington, DC. 1995: 441446. Ogniewicz RL, Kubler O. Hierarchic Voronoi Skeletons. Pattern Recognition. 1995; 28: 343-359. Otsu N. A threshold selection method from gray-level histogram. IEEE Transactions Systems, Man and Cybernetics. 1979; 9(1): 62-66. Palagyi K, Kuba A. A parallel 3D 12-subiteration thinning algorithm. Graphical Models and Image Processing. 1999; 61(4): 199-221. Palágyi K. A 3D fully parallel surface-thinning algorithm. Theoretical Computer Science. 2008; 406(1-2): 119-135. Perroton L, Aktouf Z. Parallelization of an Efficient 3D Thinning Algorithm. En: Shi H, Coffield PC, editors. SPIE Proceedings Parallel and Distributed Methods for Image Processing 1997; San Diego, California. 1997; 3166: 128-134. Pothuaud L, Van Rietbergen B, Charlot C, Ozhinsky E, Majumdar S. A new computational efficient approach for trabecular bone analysis using beam models generated with skeletonized graph technique. Computer Methods in Biomechanics and Biomedical Engineering. 2004; 7(4): 205-213. Pothuaud L, Van Rietbergen B, Mosekilde L, Beuf O, Levitz P, Benhamou CL, et al. Combination of topological parameters and bone volume fraction better predicts the mechanical properties of trabecular bone. Journal of Biomechanics 2002; 35: 1091-1099 Pudney C. Distance-Ordered Homotopic Thinning: A Skeletonization Algorithm for 3D Digital Images. Computer Vision and Image Understanding. 1998; (72)3: 404-413. Ramírez EI, Ortiz PA, Shouwenaars FR, Ruiz CO. Modelado de hueso trabecular mediante paquetería de elemento finito basándose en estructuras de Voronoi. Ingeniería Mecánica Tecnología y Desarrollo. 2007; 2(5): 151-156. Página 79 Rumpf M, Telea A. A Continuous Skeletonization Method Based on Level Sets.En: Ebert D, Brunet P, Navazo I, editors. IEEE TCVG Symposium on Visualization 2002. Konstanz, Alemania. 2002: 151-157 Saha PK, Chaudhuri BB, Majumder D. A New Shape Preserving Parallel Thinning Algorithm for 3D Digital Images. Pattern Recognition. 1997; 30(12): 1939-1955. Saha PK, Chaudhuri BB. 3D Digital Topology under Binary Transformation with Applications. Computer Vision and Image Understanding. 1996; 63(3): 418-429. Saha PK, Chaudhuri BB. Detection of 3-D simple points for topology preserving transformations with application to thinning. IEEE Transactions on Panern Analysis and Machine Intelligence. 1994; 16(10):1028-1032. Saha PK, Gomberg BR, Wehrli FW. Three-dimensional digital topological characterization of cancellous bone architecture. International Journal of Imaging Systems and Technology 2000; 11: 81-90. Sanniti di Baja G, Svensson S. A New Shape Descriptor for Surfaces in 3D Images. Pattern Recognition Letters. 2002; 23(6): 703-711. Science photo library. Bone Tissue. Science Photo Library Ltd. [citado 13 Marzo 2012]. Disponible en: http://www.sciencephoto.com Seeman E. Invited Review: Pathogenesis of osteoporosis. Journal of applied physiology 2003; 95: 2142-2151. SEER Training. Structure of Bone Tissue. U. S. National Institutes of Health, National Cancer Institute [citado 13 Marzo 2012]. Disponible en: http://training.seer.cancer.gov/anatomy/skeletal/tissue.html Seiler B, Levitt B. Unique surgical simulation and technology center opens at the University of Maryland Medical Center. University of Maryland Medical Center. 2006. [citado 8 Noviembre 2011]. Disponible en: http://www.umm.edu/news/releases/archives.htm She FH, Chen RH, Gao WM, Hodgson PD, Kong LX H.Y.Hong. Improved 3D Thinning Algorithms for Skeleton Extraction. Digital Image Computing: Techniques and Applications 2009; Melbourne, Australia. 2009: 14-19. Shilane P, Min P, Kazhdan M, Funkhouser T. The Princeton Shape Benchmark. Proceedings of the International Conference on Shape Modeling and Aplications 2004; Genova, Italia. 2004: 167-178. Siddiqi K, Bouix S. Hamilton-Jacobi Skeletons. International Journal of Computer Vision. 2002; 48(3): 215-231. Silva MJ, Gibson LJ. Modeling the mechanical behavior of vertebral trabecular bone: Effects of age-related changes in microstructure. Bone. 1997; 21(2):191-199. Página 80 Svensson S, Arcelli C, Sanniti di Baja G. Finding Cavities and Tunnels in 3D Complex Objects. 12th International Conference on Image Analysis and Processing 2003. Mantova, Italy. 2003: 342-347. Teichmann M, Teller S. Polygonal Approximation of Voronoi Diagrams of Triangles in Three Dimensions. Tech Rep 766, Lab of Comp. Sci. MIT. 1997. Ulrich D, Van Rietbergen B, Laib A, Ruegsegger P. The ability of three-dimensional structural indices to reflect mechanical aspects of trabecular bone. Bone 1999; 25(1): 55-60. Van Uitert R, Bitter I. Subvoxel precise skeletons of volumetric data based on fast marching methods. Medecal Physics. 2007;34(2):627-38. Vicent L. Efficient computation of various types of skeletons. Proceedings of SPIE Medical Imaging V 1991; San Jose, CA. 1991; 1445: 297-311. Vlaardingerbroek MT, den Boer JA. Magnetic Resonance Imaging. Theory and practice. 3rd ed. Berlin, Alemania: Springer; 2003. Wade L, Parent RE. Automated Generation of Control Skeletons for Use in Animation. The Visual Computer. 2002; 18(2): 97-110. Wang T, Basu A. A note on ‘A fully parallel 3D thinning algorithm and its applications’. Pattern Recognition Letters. 2007; 28(4): 501-506. Wehrli FW. Structural and functional assessment of trabecular and cortical bone by micro magnetic resonance imaging. Journal of Magnetic Resonance Imaging. 2007; 25(2): 390-409. Wolff J. Das Gesetz der Transformation der Knochen. Berlin, DE: Hirschwald, 1892. World Health Organization. Assessment of fracture risk and its application to screening for postmenopausal osteoporosis, WHO Technical Report Series. World Health Organization. Geneva. 1994; 843: 1-129, p.4. Wu FC, Ma WC, Liou P,Liang R, Ouhyoung M. Skeleton Extraction of 3D Objects with Visible Repulsive Force. En: Proceedings of Computer Graphics Workshop 2003; Hua-Lien, Taiwan. 2003.