Reconstrucción de alta resolución en tomografía axial computerizada : cálculo de una matriz del sistema polar por el método de Joseph

Reconstrucci´ on de alta resoluci´ on en tomograf ´ ıa axial computerizada C´alculo de una matriz del sistema polar por el m´etodo de Joseph Proyecto Final de Carrera Presentado por: Amadeo Iborra Carreres Dirigido por: Dra. Mar´ıa Jos´e Rodr´ıguez ´ Alvarez Dr. Antonio Soriano Asensi Escuela T´ecnica Superior de Ingenier´ıa Infom´atica Universidad Polit´ecnica de Valencia 27 de septiembre de 2010 ´ Indice general 1. Introducci´on 5 1.1. Presentaci´on del problema . . . . . . . . . . . . . . . . . . . . 5 1.2. Objetivo del proyecto . . . . . . . . . . . . . . . . . . . . . . 7 1.3. Descripci´on de la soluci´on . . . . . . . . . . . . . . . . . . . . 10 1.4. Descripci´on de la estructura de la memoria . . . . . . . . . . 11 2. El m´etodo de Joseph 13 2.1. Planteamiento te´orico . . . . . . . . . . . . . . . . . . . . . . 13 2.2. Implementaci´on.......................... 14 3. Antecedentes, metodolog´ıa y materiales empleados 17 3.1. M´etodos de reconstrucci´on . . . . . . . . . . . . . . . . . . . . 17 3.1.1. Hasta el presente . . . . . . . . . . . . . . . . . . . . . 17 3.1.2. Investigaci´on actual . . . . . . . . . . . . . . . . . . . 18 3.2. Resoluci´on del sistema de ecuaciones mediante el MLEM . . . 19 3.3. Tratamiento de la matriz del sistema . . . . . . . . . . . . . . 20 3.3.1. Mallado.......................... 20 3.3.2. C´alculo de los elementos de la matriz del sistema . . . 20 4. Resultados 23 4.1. Complejidad temporal . . . . . . . . . . . . . . . . . . . . . . 23 4.1.1. Influencia en el c´alculo de la matriz A......... 23 4.1.2. Influencia en las iteraciones del MLEM . . . . . . . . 25 4.2. Calidad de la imagen reconstruida . . . . . . . . . . . . . . . 25 5. Conclusiones 33 5.1. L´ıneasfuturas........................... 34 3 4 Cap´ıtulo 1 Introducci´on 1.1. Presentaci´on del problema Lo que pretendemos con la tomograf´ıa axial computerizada (TAC) es obtener una imagen en tres dimensiones del interior de un objeto sin destruirlo. Se trata de obtener una imagen en tres dimensiones a partir de m´ultiples proyecciones en dos dimensiones. La reconstrucci´on de im´agenes mediante TAC puede abordarse mediante m´etodos algebraicos. Para ello se plantea el sistema lineal de ecuaciones que describe el proceso de adquisici´on. La resoluci´on de dicho sistema de ecuaciones es el proceso de reconstrucci´on. AX =B(1.1) Donde Xrepresenta la imagen a reconstruir, Brepresenta las mediciones obtenidas del TAC y Arepresenta la matriz del sistema: la geometr´ıa del esc´aner. Arecoge las contribuciones de cada p´ıxel a cada rayo en cada rotaci´on. La imagen resultante representa el ´area que el esc´aner ha barrido en todas sus rotaciones. Este ´area se conoce como field of view (FOV) y se subdivide en p´ıxeles (conocidos como v´oxels). Considerando que cada uno de los p´ıxeles en que se subdivide el FOV es una de las inc´ognitas a resolver mediante el sistema de ecuaciones 1.1, cada ecuaci´on representa la medida de un p´ıxel del detector en cada una de las proyecciones. En funci´on de la resoluci´on del detector empleado, el n´umero de proyecciones y la resoluci´on deseada de la imagen reconstruida, el sistema de ecuaciones descrito en 1.1 puede alcanzar tama˜nos considerables. As´ı, en el caso de una adquisici´on de Nrotaciones con un detector con Sp´ıxeles, con el que se pretende reconstruir un FOV subdividido en nxnp´ıxeles, tenemos nxninc´ognitas y NxS 5 ecuaciones. Este sistema es un sistema disperso y muy grande (128x128 p´ıxeles y 300 proyecciones son datos habituales de este tipo de im´agenes). Bajo demanda de alta resoluci´on alcanza tama˜nos del orden de varios GB. La matriz del sistema plantea, pues, el problema principal al tratarse de una matriz que alcanza tama˜nos muy grandes. De la forma de los p´ıxeles en que se subdivide el FOV depende el aprovechamiento de las simetr´ıas que presenta el funcionamiento natural de TAC [1]. As´ı, un mallado de p´ıxeles cuadrados, aprovecha tan solo cuatro simetr´ıas, mientras que un mallado polar, permite aprovechar tantas simetr´ıas como posiciones de rotaci´on tiene el esc´aner, como muestra la figura 1.1. (a) Simetr´ıa Cartesiana (b) Simetr´ıa Polar Figura 1.1: La zona sombreada representa la menor ´area que resulta sim´etrica al resto en un TAC de dos dimensiones y Nrotaciones. Con este dise˜no polar del pixelado del FOV la matriz del sistema se reduce de forma que s´olo es necesario calcular un trozo de una parte de la matriz (A1) para tener la matriz completa. El resto se obtiene por rotaci´on. Si se expresan las simetr´ıas del sistema de forma matricial, el vector de proyecciones se puede descomponer en un conjunto de Nproyecciones (Bi), donde cada Biact´ua sobre todos los p´ıxeles del FOV. A su vez, Ase puede descomponer en Nsubmatrices Ai, que afectan a Bicomponentes distintas y cada Aise puede obtener a partir de rotar la submatriz A1. En consecuencia, Ase puede expresar como una matriz de bloques en la que el primer bloque (A1) corresponde al bloque calculado y el resto son 6 rotaciones de este, es decir, A=        A1 A2 A3 . . . AN        =        A1 A1R A1R2 . . . A1RN−1        donde las diferentes potencias Rirepresentan la permutaci´on o rotaci´on de los p´ıxeles necesaria para obtener la submatriz deseada a partir de la matriz A1, obteniendo al final la matriz completa. Una vez elegida la estructura de Ase deben hallar los valores de cada uno de los elementos de la matriz, que representan la contribuci´on de cada p´ıxel del FOV a cada uno de los rayos del esc´aner en cada rotaci´on. Para, a continuaci´on, resolver el sistema de ecuaciones 1.1 y obtener la imagen del objeto bajo estudio (X). 1.2. Objetivo del proyecto Este proyecto se desarrolla en el marco de la investigaci´on que, el grupo de Matem´aticas Multidisciplinares para la Industria, la Ciencia y los Servicios P´ublicos (MATMICS) en el Departamento de Matem´atica Aplicada, est´a llevando a cabo para el desarrollo de varios m´etodos de reconstrucci´on de im´agenes para un TAC de alta resoluci´on que se est´a desarrollando ´ıntegramente en la Comunidad Valenciana. Nuestro grupo de trabajo tiene dos objetivos complementarios para resolver este problema: Obtener una matriz que describa el sistema que pretendemos estudiar, de forma que esta matriz se pueda almacenar f´acilmente en la memoria RAM de un ordenador y a su vez que se pueda acceder a los elementos de la matriz de forma sencilla para resolver el sistema de ecuaciones 1.1. Resolver el sistema de ecuaciones 1.1 con calidad y rapidez. Este trabajo desarrolla el primer objetivo de nuestro grupo (construir la matriz de sistema). Se propone el desarrollo del c´alculo de los valores de la matriz del sistema para el c´alculo de integrales de l´ınea por el m´etodo de Joseph [2] (ver figura 1.2). Este m´etodo debe calcular la intersecci´on entre los rayos y un mallado polar compuesto por sectores como se muestra en la figura 1.3, i.e., el m´etodo de Joseph debe resolver el problema ilustrado en la figura 1.4. 7 Figura 1.2: Representaci´on gr´afica del c´alculo de integrales de l´ınea por el m´etodo de Joseph en un mallado cartesiano tradicional. Figura 1.3: Mallado por sectores polares para la matriz del sistema. La imagen representa un sector del FOV. 8 Figura 1.4: Representaci´on gr´afica del calculo a realizar. Cada rayo kcubre una porci´on de ´area aij del ´area total del p´ıxel ij. 9 16 Cap´ıtulo 3 Antecedentes, metodolog´ıa y materiales empleados Durante este cap´ıtulo revisaremos el trabajo, t´ecnicas y algoritmos que han servido de base para la elecci´on de algunas de las l´ıneas de investigaci´on en las que ha trabajado nuestro grupo y las cuales han hecho posible la realizaci´on de este proyecto. 3.1. M´etodos de reconstrucci´on 3.1.1. Hasta el presente Durante la primera etapa de la implantaci´on del TAC se resolvi´o el problema de la reconstrucci´on mediante m´etodos directos. Estos m´etodos contemplan el sistema de ecuaciones 1.1 y lo resuelven, bien invirtiendo la matriz del sistema o descomponi´endola en una matriz triangular y otra con buenas propiedades y usando sustituci´on regresiva. El principal problema de los m´etodos directos es el gran tama˜no que alcanza la matriz del sistema bajo demanda de resoluciones aceptables. Este hecho los hace inabordables en la pr´actica incluso con el crecimiento de la potencia de las m´aquinas actuales. Estos primeros m´etodos directos dieron paso a m´etodos iterativos y anal´ıticos. Los m´etodos iterativos usan la matriz del sistema (A) y las mediciones del esc´aner (B) para converger a una imagen (X) que da soluci´on al problema. Estos m´etodos iterativos se basan en una funci´on de minimizaci´on (cada m´etodo propone una diferente) y resuelven el problema como un problema de optimizaci´on. 17 Los m´etodos anal´ıticos se basan en una operaci´on llamada backprojection que permite analizar de qu´e p´ıxeles provienen los diferentes valores medidos por el esc´aner. Este tipo de m´etodos requieren muchas rotaciones del esc´aner (del orden de cientos) para evitar el artefacto estrella que ocurre cuando se realizan pocas proyecciones. Por supuesto se buscan m´etodos en los que la p´erdida de calidad en la imagen sea menor y se establece un compromiso entre la eficiencia (en t´erminos de tiempo hasta la convergencia y uso de espacio) y la calidad del resultado (en t´erminos de la calidad de la imagen, como resoluci´on, contraste, ausencia de artefactos, etc). 3.1.2. Investigaci´on actual El hecho de que los m´etodos iterativos sean realizables en la pr´actica con el fin de conseguir im´agenes de calidad aceptable no solo les ha hecho imponerse sobre los m´etodos directos, si no que ha potenciado la mejora y propuesta de nuevos m´etodos iterativos que cada vez alcanzan unos resultados mejores. En la actualidad, de entre los m´etodos anal´ıticos se impone el Filtered Backprojection [15] (FBP) y m´etodos basados en ´este, al resultar los que mejor eficiencia y calidad han demostrado en sus resultados. Las actuales l´ıneas de investigaci´on abiertas en nuestro grupo son las m´as prometedoras de un amplio espectro de t´ecnicas exploradas por el grupo en los pasados a˜nos. Parte de esta exploraci´on ha constituido la tesis doctoral de C. Mora Mora [3], presentada en 2008, los resultados de la cual, nos sirven de base. La exploraci´on del grupo en estos ´ultimos a˜nos concluye que el uso de m´etodos iterativos y matrices del sistema que aprovechen al m´aximo las simetr´ıas inherentes al esc´aner resulta el camino m´as prometedor. Adem´as presta especial atenci´on a un m´etodo en concreto: Maximum Likelihood Expectation Maximization [7, 8] (MLEM), que pertenece a la familia de m´etodos iterativos y dentro de ´esta a una subfamilia de m´etodos estad´ısticos. Los m´etodos iterativos estad´ısticos poseen todas las ventajas de los m´etodos iterativos y adem´as intentan modelar el ruido mediante funciones estad´ısticas. Por ejemplo el MLEM modela el ruido como una distribuci´on de Poisson y obtiene muy buenos resultados en cuanto a eficiencia y calidad de la imagen. Este m´etodo ha sido implementado por el grupo y es el m´etodo para el que se implementa el calculo de integrales por el m´etodo de Joseph, que constituye este proyecto. 18 3.2. Resoluci´on del sistema de ecuaciones mediante el MLEM Los m´etodos iterativos intentan resolver el sistema de ecuaciones 1.1 mediante una serie de estimaciones. Durante una iteraci´on dada k, estos algoritmos proyectan una imagen Xkcon la informaci´on de A, la geometr´ıa del esc´aner, y la comparan con B, la adquisici´on del esc´aner. Calculan una correcci´on de esta comparaci´on y la aplican a Xkobteniendo Xk+1. La diferencia entre cada m´etodo iterativo radica en la forma de comparaci´on y de correcci´on. El proceso iterativo empieza generando una imagen X0arbitraria (por ejemplo, una imagen uniforme, inicializada a cero, si el proceso de correcci´on consiste en una suma, o inicializada a uno, si el proceso de correcci´on consiste en una multiplicaci´on, como es el caso del MLEM). El algoritmo MLEM tiene la siguiente forma: Imagenk+1 =Imagenk×CalculoCorrecci´on MedidaEsc´aner P royecci´on(Imagenk) El t´ermino P royecci´on(Imagenk) se calcula mediante Pm j=1 aijxk jque deber´ıa ser igual a bi. Mientras que el t´ermino MedidaEsc´aner es, de hecho, bi. Mediante el cociente, el MLEM realiza la comparaci´on de su estimaci´on de los p´ıxeles cruzados por el rayo icon la adquisici´on del esc´aner para el rayo i. Por simplicidad, durante el siguiente paso de la explicaci´on nos referiremos a estos coeficientes de comparaci´on como τi(un coeficiente por cada rayo i). La funci´on CalculoCorrecci´on() se calcula mediante Pn i=1 τiaij. Este proceso es conocido como backprojection. El resultado de la backprojection debe dividirse por Pn i=1 aij de forma que quede normalizado. Este c´alculo del factor de correcci´on se efect´ua para cada p´ıxel jde la imagen. Componiendo los t´erminos de la explicaci´on llegamos a la expresi´on del algoritmo MLEM: xk+1 j=xk j 1 n X i=1 aij n X i=1 bi m X j0=1 aij0xk j0 aij que se eval´ua para cada p´ıxel jde la imagen durante cada iteraci´on k. 19 3.3. Tratamiento de la matriz del sistema 3.3.1. Mallado Este proyecto gira en torno a la matriz de sistema, parte fundamental para la reconstrucci´on de la imagen. La partici´on que se realice del FOV se ver´a reflejada en la matriz del sistema. Se han investigado muchos tipos de mallado posibles y en concreto nuestro grupo ha investigado a fondo mallados cuadrados, sectoriales (con radio constante y relaci´on de aspecto constante) y circulares [1, 3]. El mallado cuadrado consiste en dividir el FOV en p´ıxeles cuadrados de igual altura y anchura. Pero a˜nade problemas como el escaso aprovechamiento de las simetr´ıas inherentes en el funcionamiento del TAC. Otro de los inconvenientes es que el ´area pixelada no se ajusta adecuadamente a la forma del FOV (que es circular) por lo que existe un ´area que forma parte del pixelado y que no forma parte del FOV en todas las posiciones de rotaci´on del esc´aner. En el mallado polar de radio constante el espaciado radial de los p´ıxeles se fija a un tama˜no concreto. Si observamos la figura 1.3 ∆rpermanecer´ıa constante y ∆kvariar´ıa por cada anillo de forma que los p´ıxeles fuesen lo m´as grandes posible. El mallado sectorial de relaci´on de aspecto constante mantiene la relaci´on entre la altura y anchura del sector conforme se va alejando del centro del FOV. Si observamos la figura 1.3 la relaci´on entre ∆ry ∆kpermanece a la unidad. A diferencia de el mallado polar de radio constante, los anillos polares no tienen siempre la misma anchura. El mallado circular de relaci´on de aspecto constante considera circunferencias inscritas en sectores de relaci´on de aspecto constante y ofrece la ventaja de un c´alculo del ´area de intersecci´on entre el p´ıxel y el rayo un tanto m´as simple, al requerir tan solo el c´alculo del ´area de un sector del circulo. No obstante un mallado circular no cubre toda el ´area del FOV. Entre los c´ırculos quedan huecos que constituyen un problema en cuanto a la calidad de la imagen. 3.3.2. C´alculo de los elementos de la matriz del sistema Se han propuesto formas y m´etodos para el c´alculo de los valores de la matriz del sistema que representan las intersecciones entre los p´ıxeles del FOV y los rayos en cada rotaci´on del esc´aner. En concreto, nuestro grupo ha investigado a fondo entre otros el m´etodo de Siddon [16] (v´ease figura 3.1), el m´etodo de Joseph (v´ease figura 3.2) y el m´etodo de los cubos [17, 18, 3] 20 (v´ease figura 3.3). De los cuales, los resultados m´as prometedores se han obtenido mediante el uso del m´etodo de Joseph. Por ello, resulta interesante adaptar el m´etodo de Joseph para el c´alculo de las intersecciones de los rayos en los p´ıxeles de un mallado sectorial de relaci´on de aspecto constante. Figura 3.1: Representaci´on gr´afica del c´alculo de integrales de l´ınea por el m´etodo de Siddon en un mallado cartesiano tradicional. Figura 3.2: Representaci´on gr´afica del c´alculo de integrales de l´ınea por el m´etodo de Joseph en un mallado cartesiano tradicional. 21 Figura 3.3: Representaci´on gr´afica del c´alculo de integrales de l´ınea por el m´etodo de cubos en un mallado cartesiano tradicional. 22 Cap´ıtulo 4 Resultados Resulta de especial inter´es comparar el m´etodo de Joseph contra el m´etodo ya implementado por el grupo para el c´alculo de las intersecciones de los rayos con los sectores polares, como ya hemos comentado anteriormente. Esta comparaci´on se realiza en t´erminos de eficiencia temporal y calidad de la imagen obtenida. Tras la implementaci´on del m´etodo de Joseph se comprueba que ´este como el m´etodo ya implementado por el grupo anteriormente, tienen una complejidad espacial constante (O(1)) lo cual resulta ideal para la tarea en cuesti´on y deja la eficiencia espacial fuera de discusi´on. 4.1. Complejidad temporal Con el fin de determinar si existe alguna diferencia significativa entre el tiempo de generaci´on de la matriz del sistema al usar uno u otro m´etodo, se han llevado a cabo diferentes ejecuciones en una misma m´aquina con carga del sistema (ajena a la propia ejecuci´on que nos ocupa) nula. Dicha m´aquina posee un procesador modelo PC Intel Core 2 Quad funcionando a una frecuencia de 2.50 GHz y 3.5 GB de memoria RAM. 4.1.1. Influencia en el c´alculo de la matriz A Tras la ejecuci´on de la generaci´on de matrices para la reconstrucci´on de im´agenes con diferentes resoluciones y n´umero de proyecciones se han obtenido los resultados mostrados el los siguientes gr´aficos. La pareja de gr´aficos de la figura 4.1 muestra la comparativa entre el m´etodo implementado anteriormente y el m´etodo de Joseph respecto a los tiempos de generaci´on de la matriz del sistema (A). Podemos observar que desde el c´alculo con una resoluci´on de 0.5 mm 23 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 150 200 250 300 350 400 450 500 550 600 Tiempo de construcción de A (s) Número de proyecciones Pixelado de 0.5 mm (resolución 158x158) método Joseph cálculo de área de sectores 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 150 200 250 300 350 400 450 500 550 600 Tiempo de construcción de A (s) Número de proyecciones Pixelado de 0.25 mm (resolución 316x316) método Joseph cálculo de área de sectores Figura 4.1: Gr´afico comparativo del tiempo de reconstrucci´on de la matriz del sistema (A) de ambos m´etodos. (pixelado de 158x158) hasta el c´alculo con una resoluci´on de 0.25 mm (pixelado de 316x316) la tendencia en el aumento de del tiempo de c´alculo 24 respecto al aumento de proyecciones se mantiene constante. Adem´as tambi´en se mantiene la diferencia de tiempos entre ambos m´etodos. Podemos decir que el m´etodo de Joseph presenta una ligera ventaja temporal durante la generaci´on de la matriz. En todos los casos el m´etodo de Joseph tarda, aproximadamente, un 30 % menos del tiempo invertido por el m´etodo implementado anteriormente por el grupo, pero con una resoluci´on de 316x316 p´ıxeles esta diferencia representa a penas unas d´ecimas de segundo. Por lo anterior, aun habiendo diferencias en los tiempos de c´alculo de ambos m´etodos, estas diferencias no tienen mayor importancia. Al menos por el momento, pues es posible que, con el fin de obtener m´as resoluci´on o trasladar el problema a la reconstrucci´on en tres dimensiones, esta ventaja temporal pueda tener alg´un inter´es, pero no en este momento. 4.1.2. Influencia en las iteraciones del MLEM La siguiente pareja de gr´aficos (figura 4.2) muestra la comparativa entre los tiempos invertidos por el algoritmo MLEM para completar una iteraci´on en funci´on del n´umero de proyecciones y de la resoluci´on de la imagen a reconstruir, usando la matriz del sistema (A) calculada mediante ambos m´etodos. Se debe tener en cuenta que los tiempos representados en los gr´aficos resultan del c´alculo de una sola iteraci´on. As´ı pues el algoritmo MLEM realiza entre diez y treinta iteraciones para converger a la imagen reconstruida, por tanto, el tiempo total de la reconstrucci´on ser´a del orden de entre diez y treinta veces mayor que los tiempos reflejados en los gr´aficos. No obstante no se observan apenas diferencias entre ambos m´etodos. El tiempo que el algoritmo MLEM invierte en realizar una proyecci´on es, por tanto, el mismo, independientemente de que se haya usado cualquiera de los dos m´etodos en cuesti´on. As´ı pues, desde el punto de vista temporal, no existe diferencia apreciable entre ambos m´etodos. 4.2. Calidad de la imagen reconstruida Durante esta secci´on se muestran im´agenes reconstruidas mediante el algoritmo MLEM a partir de mediciones obtenidas con el m´odulo CT de Albira (correspondiente al phantom de la figura 4.3). Para el estudio de la calidad de las im´agenes se han seleccionado diferentes regiones de inter´es (region of interest), a las que a partir de ahora nos referiremos como ROI, se˜naladas en la figura 4.4 mediante ´ındices (del 1 al 5) que adem´as nos permiten saber de qu´e material est´an hechas estas ROI. Adem´as en la figura 4.4 se ha a˜nadido un detalle del phantom real para mejorar la comprensi´on de la forma de las im´agenes que se reconstruyen. 25 32 Cap´ıtulo 5 Conclusiones Los ensayos realizados han sido en circunstancias tales que pretenden modelar escenarios generales. Se ha usado un phantom 5 materiales de distintas densidades y se ha evitado el efecto del ruido. Para ambos m´etodos, se ha comparado el tiempo de generaci´on de la matriz del sistema y el tiempo de cada iteraci´on de la reconstrucci´on en funci´on del n´umero de proyecciones y resoluci´on. Tambi´en se ha evaluado la convergencia del valor medio de p´ıxel reconstruido en funci´on del n´umero de proyecciones y resoluci´on para ambos m´etodos. Tras el an´alisis de los resultados que ofrece el m´etodo de Joseph frente al m´etodo anterior desarrollado por el grupo podemos ver que ambos m´etodos tienen el mismo coste temporal, tanto en la generaci´on de la matriz del sistema como en la influencia en la iteraci´on del MLEM, e influyen de igual manera en la calidad de las im´agenes, tanto en cuanto a la pronta estabilizaci´on del valor medio de los p´ıxeles como en cuanto a la desviaci´on est´andar. No es posible, entonces, discernir cual de los dos m´etodos es preferible al otro. En una circunstancia en particular, es posible que uno de los dos m´etodos aventaje al otro. Una vez obtenida la medici´on del esc´aner, si el usuario no queda satisfecho con un m´etodo elegido, puede interesarle volver a reconstruir la imagen usando el otro m´etodo. De esta forma, no se descarta uno de los dos m´etodos, sino que se propone que ambos m´etodos formen parte del software de reconstrucci´on que se est´a desarrollando por nuestro grupo, ofreciendo la ventaja de que la reconstrucci´on de im´agenes mediante matriz del sistema de sectores polares posea dos alternativas de generaci´on a elecci´on del usuario. 33 5.1. L´ıneas futuras Paralelamente al enfoque de este trabajo, contin´uan abiertas varias l´ıneas de investigaci´on. Entre ellas la de resolver el sistema 1.1 de forma directa mediante la descomposici´on de Aen QR, de forma que A=QR y el nuevo sistema a resolver es RX =QTB, siendo Runa matriz triangular y Quna matriz ortogonal. Este proceso de descomposici´on constituye la mayor parte del tiempo en la reconstrucci´on de forma directa. Puesto que si no se var´ıa la resoluci´on de la imagen a obtener o la geometr´ıa del esc´aner la matriz Ano cambia, tampoco cambian las matrices QyRpor lo que estas matrices pueden ser precalculadas y eliminadas del coste de la reconstrucci´on. Nuestro grupo plantea en [19] un algoritmo eficiente para realizar esta descomposici´on mediante rotaciones de Givens de forma que el coste del algoritmo est´a asociado con el n´umero de elementos no nulos de A. En el caso de la reconstrucci´on de im´agenes TAC, Aposee alrededor del 1 % de elementos no nulos, lo cual permite la descomposici´on QR de matrices del sistema para resoluciones realistas y reduce el proceso de reconstrucci´on de im´agenes a una multiplicaci´on matricial (QTB) y a sustituci´on regresiva (debido a que Res triangular), dos procesos altamente eficientes y r´apidos. 34 Bibliograf´ıa [1] Mar´ıa Jos´e Rodr´ıguez ´ Alvarez, Filomeno S´anchez, Antonio Soriano, Amadeo Iborra, and Cibeles Mora. Exploiting Symmetries for Weight Matrix design in CT Imaging. Mathematical and Computer Modelling, 2010. (Aceptado). [2] Peter M. Joseph. An Improved Algorithm for Reprojecting Rays through Pixel Images. IEEE Transactions on Medical Imaging, 1(3):192– 196, 1982. [3] Cibeles Mora. M´etodos de reconstrucci´on volum´etrica algebraica de im´agenes topogr´aficas. Aplicaci´on a un TAC de peque˜nos animales y a un Simulador-TAC. PhD thesis, Universidad Polit´ecnica de Valencia, Mayo 2008. [4] Yair Censor and Gabor T. Herman. On some optimization techniques in image reconstruction from projections. Applied Numerical Mathematics, 3(5):365–391, 1987. [5] Michel Defrise, Fr´ed´eric Noo, and Hiroyuki Kudo. A solution to the long-object problem in helical cone-beam tomography. Physics in Medicine and Biology, 45(3):623, 2000. [6] Xuan Liu, C. Comtat, C. Michel, P. Kinahan, M. Defrise, and D. Townsend. Comparison of 3-D reconstruction with 3D-OSEM and with FORE+OSEM for PET. IEEE Transactions on Medical Imaging, 20(8):804–814, 2001. [7] L. A. Shepp and Y. Vardi. Maximum Likelihood Reconstruction for Emission Tomography. IEEE Transactions on Medical Imaging, 1(2):113–122, 1982. [8] K. Lange and R. Carson. EM reconstruction algorithms for emission and transmission tomography. Journal of Computer Assisted Tomography, 8(2):306–316, 1984. 35 [9] G.D. Tourassi, C.E. Jr. Floyd, M.T. Munley, J.E. Bowsher, and R.E. Coleman. Improved lesion detection in SPECT using MLEM reconstruction. IEEE Transactions on Nuclear Science, 38(2):780–783, 1991. [10] Y.K. Ng, I. Orlic, S.C. Liew, K.K. Loh, S.M. Tang, T. Osipowicz, and F. Watt. A PIXE micro-tomography experiment using MLEM algorithm. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 130(1-4):109–112, 1997. [11] R.G. Wells, P.H. Simkin, P.F. Judy, M.A. King, P.H. Pretorius, H.C. Gifford, and P. Schneider. Maximizing the detection and localization of Ga-67 tumors in thoracic SPECT MLEM (OSEM) reconstructions. Nuclear Science Symposium Conference Record, 1998. IEEE, 2:1367– 1371, 1998. [12] Reinhard M¨oller. A systolic implementation of the MLEM reconstruction algorithm for positron emission tomography images. Parallel Computing, 25(7):905–920, 1999. [13] J.E. Ortu˜no, J.J. Vaquero, G. Kontaxakis, M. Desco, and A. Santos. Dise˜no de un tom´ografo pet para peque˜nos animales con geometr´ıa octagonal: estudios preliminares. Proceedings CASEIB 2003 - XXI Annual Congress of the Spanish Society of Biomedical Engineering, pages 183–186, 2003. [14] N. Bissantz, B.A. Mair, and A. Munk. A multi-scale stopping criterion for MLEM reconstructions in PET. Nuclear Science Symposium Conference Record, 2006. IEEE, pages 3376–3379, 2006. [15] L. A. Feldkamp, L. C. Davis, and J. W. Kress. Practical cone-beam algorithm. Journal of the Optical Society of America, 1(6):612–619, 1984. [16] R. Siddon. Fast calculation of the exact radiological path length for a three dimensional CT array. Medical Physics, 12:252–255, 1985. [17] H. Turbell. Cone Beam Reconstruction using Filtered Backprojection. PhD thesis, Link¨opings University, Department of Electrical Engineering, 2001. [18] T. M. Benson and J. Gregor. Framework for iterative cone-beam micro CT reconstruction. IEEE Transactions on Nuclear Science, 52(5):1335–1340, 2005. [19] Mar´ıa Jos´e Rodr´ıguez ´ Alvarez, Filomeno S´anchez, Antonio Soriano, and Amadeo Iborra. Sparse Givens resolution of large system of linear equa36 tions: Applications to image reconstruction. Mathematical and Computer Modelling, 52(7-8):1258 – 1264, 2010. Mathematical Models in Medicine, Business and Engineering 2009. 37