Visualizacion-de-volumenes-usando-superficies-rastreadas-con-normales-suaves-obtenidas-por-art

Universidad Nacional Autónoma de México Universidad Nacional Autónoma de México\nBiblioteca Central\nDirección General de Bibliotecas de la UNAM 2011\nFátima Vicenta Quintal Flores\nTodos los Derechos Reservados Para mi familia y mis queridos GoMaPe, por todos los años compartidos Agradecimientos\n\nEl presente trabajo representa el esfuerzo de varias personas y organizaciones que directa o indirectamente colaboraron en su desarrollo. Agradezco al Dr. Edgar Garduño por haberme confiado este trabajo y por la dirección del mismo, por su paciencia al leer, opinar y corregir este trabajo pero sobre todo por los conocimientos que me ha compartido. Agradezco también al Dr. Jorge Márquez, Dr. Ernesto Bribiesca, Dr. Frnándo Arámbula y al Dr. Jorge Urrutia, los sinodales que me apoyaron en todo momento y aportaron valiosos comentarios y correcciones en este trabajo.\n\nGracias al Posgrado en Ciencia e Ingeniería de la Computación, al IIMAS y al CONACYT por el apoyo que recibí y me permitió cursar la maestría. También gracias al DGAPA proyecto IN101108 y al Departamento de Ciencias de la Computación por el equipo y las instalaciones que fueron fundamentales para el desarrollo de este trabajo. También agradezco a Lulú, Diana, Montse, Rosy y Liz por ayudarme con trámites, equipo y demás molestias que les ocasioné.\n\nAgradezco muy especialmente a toda mi familia por su apoyo y comprensión, en especial a mis padres y mi abuelita por hacerse cargo de mis responsabilidades mientras estoy lejos de casa. Gracias a Dios por cuidarlos. Gracias también a mis amigos paty y rodrigo por su compañía y apoyo incondicional, a irem por formar parte de nuestro equipo y a jorge por compartirme sus conocimientos de graficación y comics, pero sobre todo por su amistad. Por último pero no menos importante gracias a mis amigos de Mérida por acordarse y estar pendientes de mí. Resumen\n\nEn la actualidad existen varias tecnologías que producen imágenes 3D conocidas también como volúmenes. En este trabajo nos interesan aquellas imágenes 3D producidas a partir de proyecciones por medio de algoritmos basados en expansión de series. Estos métodos asumen que la imagen 3D es una función que se aproxima por medio de una combinación lineal de funciones base. Cuando las funciones base son “suaves” y continuas, los métodos producen una función continua, que posteriormente es discretizada para producir la imagen 3D. A partir la imagen 3D se pueden obtener mallas que representan superficies implícitas. Estas mallas son desplegadas usando técnicas de graficación por computadora, en particular de visualización por superficie y utilizan las normales para generar imágenes 2D que preservan las pistas de profundidad de la malla. Estas técnicas requieren de un algoritmo para determinar precisamente la superficie de los objetos. Uno de los métodos más populares es el llamado Marching Cubes [34]. Otro método es el algoritmo de rastreo de superficies propuesto por Artzy et al [2], en rejillas cúbicas simples. Algunas de las ventajas del Algoritmo de Artzy frente al de Marching Cubes es que no recorre todos los vóxeles y que produce superficies que son el equivalente discreto de una superficie de Jordan. Sin embargo, produce mallas compuestas por rectángulos, lo cual resulta en una apariencia “cuboid”.\n\nEn esta tesis proponemos usar el hecho de que para imágenes 3D producidas por medio de una combinación de funciones base suaves se pueden tener las normales de dicha aproximación suave e incorporarlas a la malla producida por el Algoritmo de Artzy para generar una representación visiblemente más suave. Índice general\n\nAgradecimientos\n\nResumen\n\n1. Introducción\n\n2. Reconstrucción a partir de Proyecciones\n2.1. Adquisición de Datos .................................. 5\n2.2. Proyecciones .......................................... 6\n2.3. Métodos de reconstrucción ............................. 7\n2.4. Técnicas de Reconstrucción Algebraica (ART) ......... 9\n2.5. Funciones Kaiser-Bessel Generalizadas (blobs) ....... 13\n2.6. Parámetros para Reconstrucción ....................... 16\n\n3. Visualización por Superficies\n3.1. Algoritmo de Extracción de Superficie de Lorensen y Cline .... 25\n3.2. Algoritmo de Seguimiento de Superficies de Artzy et al .... 28\n3.3. Visualización de Mallas ................................. 34\n3.3.1. Iluminación ........................................... 36\n3.3.2. Sombreado ............................................ 42\n3.3.3. Mapeo de Relieves (Bump Mapping) .................. 44\n3.4. Normales de una Combinación lineal de Blobs ........... 45\n\n4. Metodología propuesta para el Cálculo de Normales\n4.1. Discretización del volumen obtenido por ART ............ 48 4.2. Adquisición de la malla poligonal y cálculo de las normales ................ 50\n4.3. Experimentos ...................................................... 54\n4.3.1. Maniquíes (Phantoms) ................................... 54\n4.3.2. Datos reales ............................................... 61\n4.3.3. Conclusiones ............................................ 67\n4.3.4. Trabajo Futuro .......................................... 67\n\nBibliografía ............................................................. 69 Índice de figuras\n\n2.1. Esquema general de los procesos involucrados en imaginología biomédica 3D. .... 5\n2.2. Dos tipos de proyecciones. (a) proyección en paralelo y (b) proyección en cono ........ 8\n2.3. Esquema que muestra el proceso de ART para una función en 2D y píxeles. ....... 11\n2.4. Esquema que muestra el funcionamiento básico del método de Kaczmarz [3]. .......... 13\n2.5. Diferentes funciones base para (2.4) de [18]. .................... 14\n2.6. (a) Espectro normalizado de un blob. (b) Perfil de un blob con diferentes valores de α. .. 15\n3.1. Diferentes tipos de Rejilla ............................... 22\n3.2. Cubo lógico del algoritmo Marching Cubes. Los vértices con puntos rojos representan el conjunto I_a y los demás al conjunto E_a. A cada vértice se le asigna el valor 0 o 1 dependiendo si su valor de densidad es mayor o menor que τ. El conjunto de ceros y unos forman un byte que indica el índice en un arreglo que contiene la configuración correspondiente de triángulos. 26\n3.3. Las 15 configuraciones base de Marching Cubes. ............... 27\n3.4. Se muestra la ω-adyacencia del vóxel azul con sus 6 vóxeles adyacentes en rosa. ..... 29\n3.5. Se muestra la δ-adyacencia del vóxel azul con sus 18 vóxeles adyacentes en rosa. ..... 30\n3.6. Configuración de direcciones válidas en el Algoritmo de Artzy. .... 31 3.7. Direcciones válidas en el Algoritmo de Artzy.\n3.8. Proceso de renderizado en graficación por computador consiste en pasar de la representación en polígonos hasta el despliegue en pantalla. . . . 35\n3.9. Si cada vértice posee un valor diferente de color, el proceso de rasterización se encarga de interpolar los colores y asignarlos a los píxeles correspondientes. . . . . . . . . . . . . . . . 35\n3.10. Ley de reflexión. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36\n3.11. (a) Reflexión Especular, (b) Reflexión Difusa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37\n3.12. La visualización de una escena depende de la posición del observador, las propiedades del objeto así como del tipo y posición de las fuentes de luz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38\n3.13. Iluminación difusa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39\n3.14. Iluminación difusa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40\n3.15. Iluminación especular. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40\n3.16. Ilustración del modelo de iluminación de Phong [45], mostrando los componentes de iluminación ambiental, difusa, especular y el efecto de la suma de todos ellos. . . . . . . . . . . . . . . . . . . 42\n3.17. Sombrado Plano. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43\n3.18. Una malla más un mapa de alturas producen el Bump Mapping. . . . . . 44\n4.1. Proceso de discretización para el punto pi y el coeficiente cj. . . . . . . 49\n4.2. Esquinas para el vóxel ki. . . . . . . . . . . . . . . . . . . . . . . . . . 51\n4.3. Normal nq para el vértice qi. . . . . . . . . . . . . . . . . . . . . . . . . 52\n4.4. Algunas rebanadas de las proyecciones de una esfera vistas con la biblioteca Ximpp. . . . . . . . . . . . . . . . . . 55 4.6. Ejemplo de vértice con normal cero. La superficie pasa sobre el centro del vóxel pero no sobre todas sus esquinas, lo que puede producir normales con valor cero. . . . . . . . . . 57\n4.7. Algunas rebanadas de las proyecciones del maniquí 2 vistas con la biblioteca Ximpp. . . . . . . . . . . . . . . . 58\n4.8. Resultados del maniquí mostrado con Δv = 100. La Figura (a) muestra la malla obtenida por el Algoritmo de Artzy original, (b) muestra el resultado de aplicar las normales con el método propuesto y (c) muestra el resultado de la malla obtenida con Marching Cubes. . . . . . . . . . . 59\n4.9. Resultados del maniquí 2 mostrado con Δv = 256. La Figura (a) muestra la malla obtenida por el Algoritmo de Artzy original, (b) muestra el resultado de aplicar las normales con el método propuesto y (c) muestra el resultado de la malla obtenida con Marching Cubes. . . . . . . . . . . 60\n4.10. Algunas rebanadas de las proyecciones del maniquí 3 vistas con la biblioteca Ximpp. . . . . . . . . . . . . . . . 62\n4.11. Resultados del maniquí 3 mostrado con Δv = 100. La Figura (a) muestra la malla obtenida por el Algoritmo de Artzy original, (b) muestra el resultado de aplicar las normales con el método propuesto y (c) muestra el resultado de la malla obtenida con Marching Cubes. . . . . . . . . . . 63\n4.12. Algunas rebanadas de las proyecciones de la molécula DnaB-DnaC vistas con la biblioteca Ximpp. . . . . . . . . . 64\n4.13. Resultados de la molécula DnaB-DnaC mostrada con Δv = 256. La Figura (a) muestra la malla obtenida por el Algoritmo de Artzy original, (b) muestra el resultado de aplicar las normales con el método propuesto y (c) muestra el resultado de la malla obtenida con Marching Cubes. . . . . . . . . 65\n4.14. Resultados de la molécula DnaB-DnaC mostrada con Δv = 320. La Figura (a) muestra la malla obtenida por el Algoritmo de Artzy original, (b) muestra el resultado de aplicar las normales con el método propuesto y (c) muestra el resultado de la malla obtenida con Marching Cubes. . . . . . . . . 66 Capítulo 1\nIntroducción\n\nSe considera que el sentido más desarrollado de los humanos es la vista, por lo que la visualización es la manera más fácil de comunicar información, especialmente cuando es compleja o es dada en grandes cantidades. Esta es una de las razones por las cuales la grafificación por computador se ha convertido en el medio que permite aprovechar la abstracción matemática y modelado para comunicar información visualmente [23]. La visualización científica es una rama interdisciplinaria de la ciencia, que permite la transformación de datos científicos y abstractos en imágenes. Un proceso típico en visualización trata de encontrar, para una serie de datos, una representación visual adecuada que permita analizarlos o evaluarlos [11].\n\nUsualmente, cuando se obtienen los datos, estos incluyen varios parámetros, que en ocasiones poseen un alto grado de dimensionalidad. El almacenarlos en sistemas de administración de bases de datos puede ocasionar que sólo sea posible verlos por partes o pequeñas porciones, lo que puede ser un inconveniente a la hora de analizarlos pues podría ser complicado a simple vista distinguir relaciones entre ellos. A diferencia de estos sistemas, la visualización facilita el proceso de análisis y comprensión de los datos (algo relacionado a la geometría) a una mayor escala; lo que permite reconocer patrones de comportamiento o sus relaciones (algo relacionado a la topología) [23]. Por estos motivos se han desarrollado técnicas para transformar la información en imágenes, las cuales son la base de la visualización científica. En este trabajo nos interesan ciertos tipos de visualización científica que permiten obtener representaciones visuales de un tipo específico de funciones conocidas como funciones de densidad. Existen herramientas computacionales que generan discretizaciones de este tipo de funciones por lo que es posible transformar en imágenes datos provenientes de aparatos de medición como telescopios, microscopios, satélites y, en el caso de la biomedicina, de aparatos que usan técnicas de imaginología médica tales como tomografías (CT), resonancia magnética (MRI), radiografía, ecografía, o microscopía electrónica tomográfica.\n\nLos aparatos de medición o adquisición de datos muestran la función de densidad en tres dimensiones y, como se verá más adelante, en la práctica producen un conjunto de vóxeles que representan una discretización de la función de densidad, la cual se puede visualizar por medio de varias técnicas de graficación por computadora. Estas técnicas consisten en generar una imagen bidimensional a partir de proyectar un modelo tridimensional creado por medio de programas de cómputo [43]. Existen dos técnicas principales que permiten visualizar las funciones de densidad: la visualización directa de volúmenes (volume rendering) y la visualización por superficies (surface rendering). La primera, como su nombre lo indica, permite visualizar directamente o proyectar toda la función discretizada, mientras que la segunda, como se verá más adelante, visualiza o proyecta, solamente la frontera del objeto de interés.\n\nExisten varios métodos para realizar visualización por superficie, los cuales proporcionan diferentes niveles de detalle. Algunos suelen producir resultados más reales, pero también suelen ser más lentos debido a que tanto el hardware como el software gráfico han sido optimizados para producir, graficar, proyectar y iluminar polígonos. Por esa razón las técnicas de visualización por superficie usando proyección de polígonos son usados ampliamente en visualización científica y son las que interesan en este trabajo. Existen varios algoritmos que permiten realizar visualización por superficie y difieren básicamente en la forma en la que determinan los límites o frontera del objeto de interés. Uno de los más populares es el llamado Marching Cubes, publicado originalmente en 1987 por Lorensen y Cline [34]. Otro algoritmo para obtener la superficie es el propuesto originalmente por Artzy, Frieder y Herman [2] (nombrado a lo largo de ese trabajo como Algoritmo de Artzy et al) pero que a diferencia del Marching Cubes produce superficies visiblemente menos suaves; sin embargo, posee propiedades matemáticas más adecuadas como se verá más adelante.\n\nUna vez generada la malla que representa la superficie del objeto ésta se puede visualizar usando técnicas estándar de graficación por computadora. Como explicaremos más adelante, las normales a esta superficie juegan un papel muy importante en la generación de la imagen final. Una forma de alterar la apariencia final de la representación de los objetos consiste en modificar \"apropiadamente\" las normales a las mallas [15].\n\nEn este trabajo proponemos un método que permite producir representaciones más suaves de las mallas generadas por el Algoritmo de Artzy por medio de la incorporación de normales \"suaves\". Estas normales son obtenidas del modelo \"suave\" utilizado por algunas técnicas que aproximan la función de densidad original.\n\nExisten varios métodos que permiten obtener una estimación de la función de densidad cuando ésta no se puede medir directamente. Estos métodos reconstruyen (o estiman) la función de densidad a través de sus proyecciones (estas técnicas son usadas en microscopía, astronomía e imaginología médica).\n\nExiste una familia de algoritmos de reconstrucción que aproximan la función de densidad por medio de una combinación lineal de funciones base. Cuando estas funciones son suaves, el resultado es una función de densidad suave. Posteriormente para su manipulación y almacenamiento en computadoras, esta aproximación es discretizada.\n\nDebido a que se conoce el método de reconstrucción utilizado, se posee entonces, información de la función de densidad continua, que combinada con técnicas de graficación tradicionales deben permitir obtener una apariencia suave de la superficie a pesar de usar la malla producida por el Algoritmo de Artzy et al.\n\nLa organización de este trabajo es la siguiente. En el Capítulo 2 se describe la técnica de reconstrucción que se usa en este trabajo para obtener el modelo suave. En el Capítulo 3, se abordan tanto los algoritmos Marching Cubes y el de Artzy, como las técnicas de graficación por computadora que permitirán obtener una visualización suave de la superficie. El Capítulo 4 contiene la metodología propuesta para suavizar la superficie obtenida con el algoritmo de Artzy y los resultados obtenidos, así como las conclusiones y propuestas para trabajos futuros. Capítulo 2\nReconstrucción a partir de Proyecciones\n\nLa imaginología biomédica es el campo que se ocupa del desarrollo y uso de dispositivos de imágenes, y técnicas para obtener imágenes anatómicas internas [4]. Algunas técnicas como la tomografía computarizada (CT), resonancia magnética (MRI), radiografía, ecografía (también llamada ultrasonido), imaginología molecular y microscopía electrónica desempeñan un papel cada vez más importante en la práctica clínica y la investigación básica de ciencias de la vida; en parte porque permiten visualizar y analizar la estructura de objetos que pueden ser muy complicados o difíciles de observar de otra forma y en muchos casos con la ventaja de no alterar, al menos no de una manera drástica, su naturaleza.\n\nEl esquema general de un proceso de imaginología biomédica en 3D se muestra en la Figura 2.1. En cada paso se pueden emplear diferentes técnicas y la elección e implementación de cada una depende de la aplicación biomédica para la cual se está empleando. Para el desarrollo de este trabajo nos interesa en particular la visualización así como el método empleado para la reconstrucción. En este capítulo se describe brevemente cómo se obtienen los datos que se usan en la reconstrucción. Incluye también una descripción básica del algoritmo que permite obtener el muestreo de una función tridimensional a partir de sus proyecciones; la cual, será empleada en el proceso de obtención de superficies visualmente suaves. Figura 2.1: Esquema general de los procesos involucrados en imaginología biomédica 3D.\n\n2.1. Adquisición de Datos\n\nLos aparatos de medición utilizados en imaginología biomédica funcionan capturando, por medio de detectores, la energía emitida por alguna fuente después de que ésta ha interactuado con el objeto biológico (o espécimen). En general, algunos instrumentos utilizan radiación electromagnética cuya energía es proporcional a su frecuencia; entre mayor su frecuencia, mayor su energía. El diseño de los instrumentos se puede categorizar en tres esquemas básicos: transmisión, reflexión y emisión.\n\nAlgunas técnicas como la ecografía que emplea ondas de presión, son ejemplo de esquemas de reflexión ya que en su modalidad de operación principal se emiten ondas en el intervalo ultrasonico en dirección del espécimen a estudiar y, pasado un lapso corto de tiempo, un arreglo de detectores recibe la energía reflejada por los diferentes tejidos del cuerpo; la energía regresada está relacionada a las diferentes interfaces entre tejidos de distintas densidades.\n\nLos equipos que hacen uso de radioisótopos representan ejemplos típicos de los esquemas de emisión. En estos equipos se introduce un radioisótopo ligado a un compuesto específico para cierta función biológica. Una vez que el compuesto es ingerido, los fotones de alta energía y las partículas subatómicas producidas por el proceso de decaimiento nuclear interactúan con la materia de los tejidos. Dicha interacción es detectada por sensores alrededor del cuerpo. Otro ejemplo de este esquema es el MRI (imagen por resonancia magnética) usando radiación electromagnética (EM).\n\nFinalmente, instrumentos como los radiografos, los tomógrafos computarizados (también conocidos como TAC) y los fluoroscopios, son ejemplos clásicos del esquema de transmisión. En este esquema el objeto a examinar se coloca entre una fuente de energía (por ejemplo, un emisor de rayos x) y un arreglo de detectores. La energía viaja en dirección del objeto e interactúa con el tejido biológico. Los detectores capturan el resultado acumulado de dicha interacción.\n\nLa radiación utilizada por varios de estos instrumentos posee un gran contenido energético y por ello, las partículas o fotones que son parte de esta radiación tienden a viajar en trayectorias rectilíneas. 2.2. Proyecciones\n\nPor todo lo anterior, los instrumentos detectan la acumulación de la radiación con la materia a lo largo de las trayectorias rectilíneas. De manera más formal, este proceso se puede expresar como:\n\n[Pν] (~σ,~x) = ∫R1 ν(~x + τ~σ)dτ, (2.1)\n\ndonde ν representa la función de densidad, ~x es un punto en el espacio que se representa por medio de un vector y ~σ es un vector de dirección (vector normalizado). Entonces, el operador P representa la integral de línea recta sobre la trayectoria definida por la recta ~x + τ~σ, con τ ∈ R. La transformada (2.1) es mejor conocida como la Transformada del Rayo [30].\n\nEn este trabajo también representaremos a los vectores ~x ∈ Rn, alternativamente k̄ ∈ Zn, por medio de una tupla de n números reales (o n-tuple) representada por (x1, x2, ..., xn), alternativamente (k1, k2, ..., kn).\n\nUn conjunto de integrales de línea se conoce como proyección. Una notación alternativa para expresar una proyección es por medio de la Transformada de Radón 8\n\nJ = 1 - p donde p = No. de proyecciones\nl = -q...q donde q = No. de integral de línea\n[P V ] | (σ̄ 1 , x̄𝑞 )\n\n(a)\n\nS\n\nC\n\nB\n\n- q\n\nj = 1 ... p\n\nl\n\n[-pV ] | (σ̄ 1 , x̄𝑞 )\n\nC\n\nS\n\nB\n\n(b)\n\nFig. 2.2: Dos tipos de proyecciones. (a) proyección en paralelo y (b) proyección en cono. 9\n\nPor otra parte, los métodos basados en expansión de series asumen que cualquier función de densidad puede aproximarse por medio de la combinación lineal de funciones base como sigue\n\nν (x̄) ≈ ∑\n\nj = 1\n\nJ\n\nc\n\nj b\n\nj (x̄ - p̄j ) ,\n\n(2.4)\n\ndonde x̄ es un vector de posición; es decir, un punto en el espacio, {b\nj } es el conjunto de funciones base, cada una ponderada por el correspondiente coeficiente c\nj . El conjunto {p̄\nj } representa los puntos donde se encuentran los orígenes de las funciones base {t\nj }. En esta aproximación, las funciones base son conocidas de forma anticipada, mientras que los valores de los coeficientes c\nj son determinados por el algoritmo de reconstrucción. En este trabajo estamos interesados en este tipo de algoritmos, en particular en los métodos conocidos como Técnicas de Reconstrucción Algebraica (ART por sus siglas en inglés) [27].\n\n2.4. Técnicas de Reconstrucción Algebraica (ART)\n\nSustituyendo (2.4) en (2.1) se tiene\n\n[P V ] | (σ̄ , x̄ ) ≈ ∫\n\nR\n\n∑\n\nj = 1\n\nJ\n\nc\n\nj b\n\nj (x̄ + τσ̄ - p̄j ) dτ ≈ ∑\n\nj = 1\n\nJ\n\nc\n\nj ∫\n\nR\n\nb\n\nj (x̄ + τσ̄ - p\nj ) dτ,\n\nque resulta en la aproximación\n\n[P V ] | (σ̄ , x̄ ) ≈ ∑\n\nj = 1\n\nJ\n\nc\n\nj [P b\nj ] | (σ̄ , x̄ ) .\n\n(2.5)\n\nEs decir, dado que se puede definir una proyección para todo (σ̄, x̄), entonces la ecuación (2.5) indica que la transformada de rayo de la función de densidad ν (x̄) es similar a la suma de las transformadas de rayo de las funciones {b\nj } sopesas por los coeficientes { 10\n\nse caracterizan por medio de (¯σ\nm , x̄\nm ), para 1 ≤ m ≤ M , entonces tenemos:\n\ny\nm ≈ ∑\n\nj = 1\n\nJ\n\nc\n\nj [ x\nm ] ≈ c\n\nj ℓ\n\nmj ,\n\n(2.6)\n\ndonde y\nm es la m-ésima medición y ℓ\nmj representa la integral de línea de la función base j en la dirección de la línea m, es decir, representa la intersección de una función base y una línea en particular. Entonces para cada línea M se tiene\n\ny\n1 =\n\n...\n\n+ c\n\nj ℓ\n\n2J\n\n.\n\n.\n\n.\n\ny = c\n\n1 ℓ\n\nM 1 + c\n\n2 ℓ\n\nM 2 + .. + ℓ\n\nM J .\n\nClaro está que se puede expresar el problema por medio de un sistema de ecuaciones lineales ȳ = Lc.\n\nEn este momento hacemos una pausa para ilustrar esta aproximación en una función de densidad bidimensional. Para este ejemplo usaremos la Figura 2.3 como referencia. Por simplicidad usaremos píxeles cuadrados como funciones base, definidos como\n\nb\nj (x̄) = \n\n{ \n\n1, si − 1/2 ≤ |x̄| < 1/2 ,\n\n0, e otro caso,\n\n(2.7)\n\ndonde x̄ ∈ R2 .\n\nPara este ejemplo usamos 100 funciones base (J = 100) para representar la función. Cada medición y\nm se puede expresar como m = c\n\n1 ℓ\n\n1 m + c\n\n2 ℓ\n\n2 m2 + · · · + c\n\nN ℓ\n\nJ ;\n\nen particular para la medición y\nm + 1 de la Figura 2.3, la suma queda como y\nm + 1 =\n\nc\n\n51 ℓ\n\nm + 1 , 51 + c\n\n52 ℓ\n\nm + 1 , 52 + c\n\n61 ℓ\n\nm + 1 , 61 + c\n\n62 ℓ\n\nm + 1 , 62 + c\n\n63 ℓ\n\nm + 1 , 63 + c\n\n7 ℓ\n\nm 1 , 73 + c\n\n7\n\nm + 1 , 98\nn\n.\n\nComo se puede observar en la Figura 2.3 los elementos c\nj para y\nm + 1 corresponden a las celdas verdes. Cada una de las mediciones ℓ\nm,j se puede calcular fácilmente. ya que es la intersección del píxel j y la línea (o rayo) (\\bar{x}_m, \\bar{x}_m). Dichos segmentos de línea pueden ser precalentados rápidamente. Podemos notar que esta explicación puede extenderse sin perder generalidad al caso de ℝ^3 y vóxeles cúbicos.\n\nReceptor\n\nEmisor\n\nFigura 2.3: Esquema que muestra el proceso de ART para una función en 2D y píxeles.\n\nEn la práctica el sistema ȳ = Lc está sobredeterminado debido a la necesidad de obtener muchas mediciones para compensar el ruido (aunque recientemente se estudiaron métodos donde el sistema está subdeterminado [13, 17]). Además, si los valores de M y J fueran pequeños entonces se podría utilizar cualquier método de inversión de matrices para resolver ȳ = Lc; sin embargo, en la práctica estos valores suelen ser grandes, por ejemplo, en una rejilla de tamaño 256 × 2 × 6 (considerado moderado), el valor de J es igual a 65,536. También en la práctica las líneas m son vistas como delgadas bandas de ancho σ que en muchas ocasiones son aproximadamente iguales al ancho de los píxeles, así considerando que el ancho de las líneas debe ser aproximadamente el ancho de los píxeles, se requieren entonces J rayos, es decir, M = J. Así, se tienen M × J ecuaciones, o sea, 65, 536 × 65, 5 4.2 × 10^9 y más aún, si se requieren cálculos de doble precisión cada elemento de la matriz ocuparía 8 bytes por lo que una sola matriz requeriría de 8 × 4.2 gigabytes de almacenamiento [3].\n\nEntonces, debido a lo anterior y principalmente a que en el sistema ȳ = L la matriz es dispersa (la mayor parte de elementos son cero), este sistema se puede resolver por medio de la generalización del método de Kaczmarz [31], el cual es un método iterativo definido como sigue:\n\ny_m − \\sum_{i=1}^{J} e_{m,i} c_{i}^{(k)} \\left/ \\sum_{i=1}^{J} (e_{m,i})^2 \\right.,\n\n(2.8)\ndonde λ^{(k)} es un parámetro de relajación que determina la cantidad de actualización para la iteración k.\n\nUna imagen formada por el conjunto de coeficientes (o valores de densidad) de los píxeles (en el caso de dos dimensiones) representados por {c_j} puede verse como un punto en el espacio J-dimensional, por lo tanto cada ecuación y_m puede verse como un hiperplano. De esta forma si la solución al sistema de ecuaciones ȳ = Lc es única, ésta será el punto de intersección de dichos planos. En general el método propuesto por Kaczmarz consiste en dar un punto inicial a partir del cual se realizan una serie de proyecciones perpendiculares, es decir, este punto inicial (denotado por c^{(0)}_1, c^{(0)}_2, ..., c^{(0)}_J) y representado por el vector \\vec{c}^{(0)} en la Figura 2.4 se proyecta a la primera línea (y_1 en la Figura 2.4), el punto resultante \\vec{c}^{(1)} se proyecta a la segunda línea (y_2). Cuando se llega a la última línea, se proyecta nuevamente a la primera y así sucesivamente; por lo que, si la solución es única, el método converge al punto de intersección de los hipersuperficies.\n\nAlgunos instrumentos de adquisición de datos tales como PET [29], producen la reconstrucción usando la metodología anterior. En microscopía electrónica se ha demostrado que usar el algoritmo ART produce mejores resultados que cuando se usa alguno de los basados en transformadas [18]. Las funciones base usadas en el ejemplo anterior no son suaves, sin embargo, muchos objetos se pueden aproximar empleando funciones base como las que se verán a continuación.\n\n(a) Proyectar los puntos sobre los hiperpianos.\n\n(b) Si existe la solución el método converge a ese punto.\n\nFigura 2.4: Esquema que muestra el funcionamiento básico del método de Kaczmarz [3].\n\n2.5. Funciones Kaiser-Bessel Generalizadas (blobs)\n\nLas funciones base de (2.7) son no-continuas (lo cual resulta en una función ν(x̄) discontinua) y aunque ellas han sido utilizadas en aplicaciones prácticas, para algunas aplicaciones biomédicas se usan funciones base continuas (tal como en TEM [46]), con las cuales se han producido mejores resultados que con las funciones de (2.7).\n\nEn este trabajo usamos funciones base radialmente simétricas que son continuas y con una transición suave en el intervalo [1, 0], definidas como b(\n donde es la distancia radial desde el centro de la función, es la función de Bessel de orden, a es el radio de la función (o soporte, es decir donde la función no vale cero) y α es el parámetro que controla su forma. En el resto de este trabajo nos referiremos a la función de (2.9) como blob. Entonces los parámetros m (un entero no negativo), a y α (números reales no negativos) controlan la suavidad y forma del blob por lo que la elección de estos valores es importante ya que afectan los resultados del algoritmo de reconstrucción. Estas funciones fueron propuestas para reconstrucción por Lewitt [33] y son una generalización de las funciones Kaiser-Bessel usadas en procesamiento digital de señales.\n\n[PV]\n\n(a) Cubos como funciones base. \n(b) Blobs como funciones base.\n\nFigura 2.5: Diferentes funciones base para (2.4) de [18].\n\nEn el ejemplo visto anteriormente de reconstrucción ART para el caso de dos dimensiones (ver Figura 2.3) las funciones base fueron celdas cuadradas, para el caso de tres dimensiones como se puede observar en la Figura 2.5 (a), las celdas resultan ser celdas cúbicas; sin embargo, la función reconstruída resulta no tener una transición suave debido a los cambios bruscos de las celdas. La Figura 2.5 (b) muestra el caso en el que se usan los blobs como funciones base y debido a su característica de suavidad permiten obtener reconstrucciones más suaves que al usar celdas cúbicas.\n\nEstas funciones son ampliamente usadas en procesamiento de imágenes pues poseen soporte compacto al ser acotadas en el intervalo [1,0] y al poseer, en la práctica, ancho de banda limitado tal como puede observarse en la Figura 2.6.\n\nFigura 2.6: (a) Espectro normalizado de un blob. (b) Perfil de un blob con diferentes valores de α.\n\nComo se mencionó anteriormente el parámetro m del blob, es un entero positivo que determina el orden o número de derivadas que posee en el punto a. El parámetro α se relaciona con la tasa de cambio de la función por lo que en visualización generalmente es asociado con la suavidad de la función. Cuando m = 2 el blob es una función suave que tiene derivada en todo el soporte y como aún no se ha demostrado que valores de m mayores a 2 sean más efectivas, entonces tal y como se hace en [19], usaremos dicho valor. En la Figura 2.6(b) podemos observar que el parámetro α afecta la forma del blob, por lo que, mientras más pequeño sea este valor la función decrece más lentamente. 2.6. Parámetros para Reconstrucción\n\nCuando se usa la función (2.9) para la aproximación de (2.4) es necesario seleccionar apropiadamente la distribución {pj} de los centros de los blobs y los parámetros α y a de los mismos. En [36] se propone un criterio para seleccionar estos parámetros. A continuación hacemos un resumen de este criterio. Primero, definimos la transformada de Fourier de una función g en Rn como\n\nĝ(ξ) = (2π)−n/2 ∫Rn g(x)e−i〈x,ξ〉dx, (2.10)\npara todo ξ ∈ Rn. Una relación importante entre la convolución de dos funciones y sus transformadas de Fourier está dada por el siguiente teorema.\n\nTeorema 1 (Convolución). Sean f y g dos funciones en Rn, entonces\n\nf * g = (2π)n/2 f̂ × ĝ y f × g = (2π)n/2 ĝ × f̂. (2.11)\n\nLas distribuciones de puntos que son de interés para este trabajo son las siguientes. El primer conjunto es el más usado en procesamiento de imágenes y está definido por\n\nGA = {Δk̄|k̄ ∈ Z3}, (2.12)\n\ndonde Δ es un número real positivo que se conoce como la distancia de muestreo (ver Figura 3.1 (a)). Es fácil ver que para cualquier punto k̄ ∈ GA su vecindario de Voronoi es un vóxel cúbico definido por\n\nvóxel (k̄) = {x̄ ∈ R3|Δ \n\n(k̄i − 1/2) Δ \n\n(k̄i + 1/2), para 1 ≤ i ≤ 3}. (2.13)\n\nEs fácil ver que los voxeles de (2.12) son cubos y por ello el conjunto GA se conoce como rejilla cúbica simple (sc). Cabe mencionar que para referirnos a un vóxel cúbico es suficiente con referirnos al punto k̄ que le da origen.\n\nOtro conjunto de gran interés en procesamiento de imágenes es definido por Para que esta ecuación aproxime de la mejor forma la transformada de Fourier de una función constante, la cual es un pulso en el origen, es apropiado seleccionar la función del blob de tal forma que b(2, α, a; R) sea cero en las posiciones de FΔ que se encuentren más cercanas al origen; en otras palabras ξ = √πmΔ . La transformada de Fourier de un blob está definida como sigue[33]\nb(2, α, a; R) = α3α2 I2(α)\n{\n I z 2 (√a2−(aR)2) (√a2−(aR)2) 7/2, si aR > α,\n I z 1 (√(aR)2−α2) (√(aR)2−α2) 7/2, si aR > α,\n}\n(2.20)\ndonde I z 2 nunca se hace cero y I z 1 (√ se hace cero por primera vez cuando 2 = 6.987832. Por lo tanto, es fácil ver que se tiene la siguiente relación\nα = √2π (a Δ )2 − 6.9879322. A este conjunto se le conoce como la rejilla cúbica centrada en el cuerpo (bcc) y el vecindario de Voronoi asociado a cada punto de esta rejilla es un octaedro truncado (ver Figura 3.1(b)).\n\nFinalmente, nos interesa el conjunto\n\n= {Δl|k|k ∈ Z3 y k1 + k2 + k3 ≡ 0 (mod 2)} .\n\nA este conjunto de puntos de le conoce como rejilla cúbica centrada en la cara (fcc) y el vecindario de Voronoi asociado a cada uno de sus puntos es un rombo dodecaedro (ver Figura 3.1(c)).\n\nUn tren de pulsos sobre una rejilla Γ se define como IIIΓ = ∑γ∈Γ δxγ. Por lo tanto podemos tener trenes de pulsos sobre las rejillas GΔ, BΔ y FΔ definidas como IIIGΔ, IIIBΔ y IIIFΔ.\n\nEn [5, 18] se ha demostrado que las transformadas de Fourier de estos trenes de pulsos se definen como sigue\n\nIIIGA = (√2π Δ )3 IIIGx3, \nIIIBA = 1√2 (√π Δ )3 IIIFx3,\nIIIFA = √2 (√π Δ )3 IIIBx3.\n\nSe ha demostrado en [37, 14] que la rejilla bcc es la más \"eficiente\" para mostrar el espacio Euclideano 3D cuando una función de ancho de banda limitado y espectro radialmente simétrico. Para reconstrucción, Matej y Lewitt [36] demostraron que la rejilla más deseable para {pj} es la bcc. Entonces, con lo definido anteriormente queda seleccionar los parámetros α, a y Δ. Para la combinación lineal de blobs con coeficientes cj = 1 para 1 ≤ j ≤ J se debe tener una aproximación de una función Para seleccionar un solo conjunto de parámetros α y a se realizó el siguiente procedimiento. Se selecciona una aproximación ν(x̄) = c1b(2, α, m; |x̄ − p̄1|) ∗ c2b(2, α, m; |x̄ − p̄2|) para c1,2 = 1 y una separación |p̄1 − p̄2| = 2Δ. Para dicha aproximación se busca la isosuperficie para ν(x̄) = 12 que sea apenas convexa. Para hallar dicha superficie se utiliza los parámetros obtenidos para aproximar las normales de la esfera de densidad igual a uno. Bajo dichas circunstancias se halló que los parámetros que obtienen una superficie apenas convexa son los siguientes α = 13.36, a = 2.4, m = 2 y Δ = 1√2. Por lo tanto, estos son los valores que usaremos en este trabajo. Capítulo 3\n\nVisualización por Superficies\n\nUna vez que se obtiene la aproximación de una función de densidad ν por medio de (2.4) con las funciones base definidas por (2.9) es posible visualizarla por dos métodos básicos. El primero utiliza el método raycasting para hallar el conjunto\n\nSτ = {x̄|ν(x̄) = τ},\n\n(3.1)\ndonde τ es un valor real. La hipótesis detrás de esta ecuación es que existe un umbral (o isovalor) τ para el que los valores de ν mayores a τ definen los puntos dentro de los objetos de interés, mientras que valores de ν menores a τ definen los puntos en el exterior de los objetos de interés. Por lo tanto, el conjunto de (3.1) define puntos en la frontera o superficie de los objetos de interés. El conjunto Sτ también es conocido como iso-superficie o superficie implícita y representa puntos de valor constante en un espacio volumétrico, es decir, es un conjunto de nivel de la función continua ν cuyo dominio es el espacio tridimensional. En el caso de imaginología biomédica representa regiones de las imágenes en donde hay un valor particular de densidad dado por τ. La selección de umbral o isovalor en la ecuación (3.1) es en sí un proceso tradicional de segmentación conocido como umbralización (o thresholding).\n\nEn la práctica, el método de raycasting se puede usar para visualizar la aproximación (3.1) usando 2.4 directamente. El método consiste, básicamente, en \"lanzar\" varios rayos desde el plano de proyección (típicamente la pantalla de la computadora) hacia el arreglo de puntos {p̄j} que tienen asociados los valores del conjunto {c̄j}. Estos rayos pueden ser perpendiculares al plano de proyección (perspectiva ortogonal) o en la dirección de un punto fijo en el espacio (proyección en perspectiva). Para cada rayo se encuentran las M funciones que interceptan dicho rayo. Para dicha línea se encuentra el punto más cercano al plano de proyección tal que ∑M m=1 cmbm(x̄ - p̄m) - τ = 0. Existen varios métodos que realizan este proceso de manera eficiente como en [15] y en [44]. Cabe mencionar que en grafi cación por computadora el modelado de objetos por medio de (2.4) se conoce como modelo blobby, propuesto por [8] y se han sugerido varias funciones base tales como las propuestas por [36].\n\nUn segundo método para visualizar la función continua ν consiste en, primero, discretizarla de la siguiente forma\n\nVΓ(k̄) = ν(x̄) × IIIΓ,\n\n(3.2)\ndonde IIIΓ representa el tren de pulsos distribuido o arreglado de una forma específica [18]. Arreglos típicos son los conjuntos 2.12, 2.14 y 2.15 que se describieron en el capítulo anterior.\n\nSin embargo, la forma más común de realizar la discretización de la función ν de (2.4) es\n\nVGA(k̄) = ν(x̄) × IIIGA.\n\n(3.3)\nEl uso de las rejillas (2.14) y (2.15) es poco común a pesar de poseer mejores propiedades para el procesamiento de imágenes y de ser más óptimas que la definida en (2.12) [14]. En este trabajo usaremos la discretización de ν definida por (2.12). Después de un proceso de adquisición de imágenes y, generalmente de un proceso de conversión de datos y filtrado, se obtiene una imagen 3D. A partir de aquí se pueden seguir varios caminos para visualizar un objeto pero generalmente se clasifi can en dos grupos. El primer grupo asume que la imagen 3D VGA de (2.4) se puede entender como un conjunto de imágenes 2D de regiones contiguas apiladas. Por lo tanto, se puede representar el objeto en la imagen 3D a partir de los contornos del objeto presentes en cada imagen 2D. El segundo grupo representa a los objetos en\n\n(a) Rejilla SC. (b) Rejilla BCC. (c) Rejilla FCC.\n\nFigura 3.1: Diferentes tipos de Rejilla\n\nVGA, tomando en cuenta toda la imagen 3D. Entre los métodos del primer grupo se encuentra el algoritmo presentado por E. Keppel [32] donde se reconstruye a partir de una pila de contornos de las imágenes transversales en dos dimensiones. Varios autores han hecho alguna modificación de este algoritmo como los presentados en [9] y [16] pero básicamente consisten en definir primero los contornos de los objetos de interés en cada una de las imágenes ya sea de forma interactiva o utilizando algún método de detección de bordes. Posteriormente los contornos de imágenes contiguas se conectan para poder formar una estructura tridimensional, esto generalmente se realiza por medio de triangulación. Así, el paso principal de estos métodos es la técnica empleada para conectar los contornos, pero en la práctica éstos pueden ser formas extremadamente complejas, tanto, que de una imagen a la siguiente los cambios pueden variar ampliamente. Por esa razón, la conexión de los contornos generalmente se basa en heurísticas donde hay que definir cuál se aplica a cada caso, por lo que esos métodos casi no son empleados actualmente [24].\n\nEn el caso de los métodos que emplean toda la imagen 3D, el siguiente paso en el proceso de visualización consiste en identificar los diferentes objetos representados en el volumen de tal forma que se puedan seleccionar aquellos vóxeles que representan los objetos que se desea visualizar. Este paso se lleva a cabo por medio de un proceso conocido como segmentación. La forma más simple consiste en binarizar los datos con un umbral (o isovalor) de intensidad para distinguir, por ejemplo, los huesos de otros tejidos (como los músculos) en el caso de tomografías computarizadas. Este método de umbralización asume que se tiene una distribución binomial en los valores de densidad de la función VGA . Esta situación no siempre se cumple y por lo tanto existen técnicas más sofisticadas, así que la investigación en esta área es abundante.\n\nPara los métodos que trabajan directamente en toda la imagen VGA existen dos métodos básicos para realizar su visualización.. El primero de ellos se conoce como visualización directa de volúmenes (volume rendering). El segundo es la visualización indirecta o por superficie (surface rendering); haciendo énfasis en que este puede ser ejecutado por métodos estándares de graficación por computadora. Los métodos basados en superficie y los basados en vóxeles tienen sus méritos; sin embargo, la decisión de cuál emplear en una aplicación especifíca depende de los recursos computacionales disponibles así como de los objetivos de la visualización [24].\n\nLa visualización directa de volúmenes consiste en asignar un valor de opacidad a cada vóxel de VGA (proceso conocido como clasificación o función de transferencia [48]). Posteriormente, todos los vóxeles son proyectados y sus opacidades se acumulan a lo largo de la dirección de proyección. Cabe mencionar que existen otros esquemas además del de acumulación; tales como seleccionar el valor máximo, Maximum Intensity Projection (MIP) [22, 38].\n\nLa visualización por superficies se basa en (3.1); sin embargo, sobre la función discretizada VGA . Por lo tanto, los diversos métodos que emplean esta técnica encuentran el conjunto Fr que aproxima al conjunto Sr de (3.1). Existen muchos métodos en la literatura para hallar el conjunto Fr pero básicamente pueden ser clasificados en dos categorías: basados en descomposición de vóxeles y los basados en crecimiento de regiones de la superficie [1]. Dentro de la primera categoría el método más conocido es el Marching Cubes propuesto por Lorensen y Cline en 1987 [34], más adelante se describe con un poco más de detalle pero básicamente consiste en definir un cubo formado por ocho vóxeles contiguos y determinar si las esquinas de este cubo (centros de los ocho vóxeles) se encuentran dentro o fuera de la superficie. De este modo, dependiendo de la configuración del cubo, se realiza una triangulación generando así por lo menos un triángulo por cada vóxel que atraviesa la superficie. Como la superficie queda definida por pequeños triángulos entonces generalmente se requiere de algoritmos de postprocesado y simplificación. También se han realizado variantes al Marching Cubes o modificaciones para resolver ambigüedades o prevenir agujeros. En [39] se presenta un resumen de los diferentes trabajos realizados a partir de Marching Cubes.\n\nCon respecto a los algoritmos basados en crecimiento de regiones de la superficie para encontrar Fr , se encuentra el propuesto por Stoddart, Illingworth y Windeat [28]. Este algoritmo es conocido como Marching Triangles y produce una malla a partir de una nube de puntos. Básicamente consiste en encontrar un triángulo perteneciente a la superficie y ir añadiendo vértices que formen triángulos con las aristas pertenecientes a la malla y que cumplan con los requisitos de una triangulación de Delaunay. Comparado con los métodos basados en vóxeles este algoritmo puede generar mallas con mayor calidad y mejor equip espacidas; sin embargo, para superficies con grandes curvaturas puede generar un gran número de triángulos por lo cual se han propuesto algoritmos como el presentado en [1].\n\nSin embargo, el Marching Cubes, se ha establecido como el estándar de visualización para superficies [39]. Este algoritmo produce un conjunto Mr que aproxima una isosuperficie Sr . Por otra parte, propuesto por Artzy, Frieder y Herman [2], es un algoritmo que produce un conjunto Ar que aproxima al conjunto Sr , burdamente (crea representaciones cuboides) pero posee ciertas propiedades interesantes que se describirán más adelante.\n\nAsí, una vez que se obtiene el conjunto Fr (ya sea Mr o Ar ) por cualquiera de los métodos antes mencionados, la representación final se produce usando técnicas estándares de graficación por computadora [12, 15]. Como veremos más adelante, los vectores normales a la superficie juegan un papel muy importante en la generación de la representación final. Por ejemplo, las normales asociadas al conjunto Ar producido por la implementación original del Algoritmo de Artzy et al contribuyen a la apariencia cuboide de las representaciones obtenidas por éste método. Sin embargo, en este trabajo debido a sus propiedades matemáticas estamos interesados en usar el Algoritmo de Artzy para producir el conjunto Fr , pero suavizaremos su representación final por medio de asignar normales \"apropiadas\" a los vértices de la malla Ar . Sabe