Tomografía Computarizada

Tomograf´ıa Computarizada. Kristiansen L. Lara D. Departamento de F´ısica Universidad de Santiago de Chile (USACH) Avenida Ecuador 3493, Estaci´on Central, Santiago de Chile. 30 de diciembre de 2016

1.

Introducci´ on y prop´ osito

da (ver figura 2), de la cual nos referiremos por su terminolog´ıa en ingl´es: CT scan.

En una radiograf´ıa de rayos x, tal como se ve en la figura 1, la estructura tridimensional del Consecuentemente, el objetivo del presente cuerpo es representada por una imagen bidi- art´ıculo es poder: mensional. Durante el procesado de la imagen, Proveer de una introducci´on cualitativa una de las dimensiones de la imagen se pierde. para el lector no especializado en el tePor ende, todos los planos del paciente que son ma, y as´ı poder facilitar una apertura a paralelos al film de rayos x son superpuestos, textos m´as t´ecnicos y especializados. uno encima de otro. Comprender la f´ısica detr´as del proceso de creaci´on de rayos x. Entender como la creaci´on de rayos x origina ciertos fen´omenos u ´tiles para la formaci´on de una imagen. Familiarizarse con la matem´atica detr´ as de la formaci´on de la imagen, enmarcado en una geometr´ıa de haces paralelos. El art´ıculo se dividir´a en dos partes: la primera parte contendr´a la f´ısica de la producci´ on de rayos x y las caracter´ısticas y ventajas de la formaci´on de una imagen utilizando ciertos fen´omenos que ocurren. Ya en la segunda parte, se abordar´a la matem´atica detr´as de este proceso, en comuni´on con la geometr´ıa de haces paralelos, llegando al concepto de Transformada de Radon y el Fourier Slice Theorem. En base a estas dos partes se formar´a una conclusi´on, recapitulando las ideas principales, para la simplicidad del lector.

Figura 1: Rayos x del pecho. Notemos la superposici´on de las partes del cuerpo, una encima de otra. Entendiendo que esto presenta una clara desventaja a la hora de realizar un correcto diagn´ostico de enfermedades del paciente, se ha propuesto la creaci´ on de la tomograf´ıa de rayos x, o llamada tambi´en tomograf´ıa computariza-

1

Figura 2: Tomograf´ıa computarizada de rayos x [2] de cuarta generaci´on.

2.

F´ısica de la CT scan

de tungsteno (Z = 74) con U = 100 [kV ], entonces η = 0, 0081 %, lo cual es menor al 0, 1 % En esta secci´ on, explicar´e de manera m´as cua- de eficiencia. La eficiencia es algo dif´ıcil de melitativa que cuantitativa la creaci´ on de rayos x jorar, ya que es directamente proporcional a U y los efectos que ella produce en un tejido. y Z, y en ambos casos, ambas cantidades son dif´ıcil de aumentar experimentalmente.

2.1.

Creaci´ on de rayos x

Ahora bien, un rayo x con longitud de onda en el rango de 10 [nm] (124 [eV ]) a 0, 1 [nm] (12, 4 [keV ]) se dice que son rayos x suaves, debido a que no pueden penetrar l´aminas gruesas de materiales. Estos rayos no se utilizan en radiolog´ıa por obvias razones. Los utilizados son aquellos que van de 0, 1 [nm] (12, 4 [keV ]) a hc E = hν = 0, 01 [nm] (123 [keV ]). Pero la pregunta que λ nos surge es la siguiente: ¿por qu´e no se utiAhora bien, si llevamos electrones a alta velocilizan rayos x con menor longitud de onda si dad, entonces se producir´ an rayos x. La energ´ıa penetran m´ as? cin´etica se transformar´ a en radiaci´ on electromagn´etica. En este caso, la eficiencia η de la Rayos x utilizados producci´on de rayos x vendr´ıa dada por la reLongitud de onda [nm] Energ´ıa [keV ] laci´on [1] 0, 1 a 0, 01 12, 4 a 123 η = aU z Cuadro 1: Rayos x utilizados, debido a su ma−9 −1 , donde a = 1, 1 × 10 [V ] es una constante, yor profundidad. U el potencial de aceleraci´ on medido en [V ] y Z el n´ umero at´ omico del elemento con que se tra- La respuesta es porque simplemente es dif´ıcil baje. Para asentar ideas, si tenemos un objetivo producirlos1 . La luz es una onda electromagn´etica, eso es bien sabido. No obstante, debido a los avances de la f´ısica moderna, se ha descubierto que la luz viene cuantizada en peque˜ nos paquetes llamados fotones.

1 Aunque si soy m´ as preciso con la respuesta, la diferencia radica en que con las m´ aquinas que producen rayos x directamente acelerando electrones en un metal de n´ umero at´ omico alto (como tungsteno) puedo controlar mejor

2

chos tipos de colisiones. La mayor´ıa de estos encuentros involucran peque˜ nas transferencias de energ´ıa desde un electr´on a altas velocidades a electrones que son golpeados fuera de los ´ atomos, originando ionizaci´on en los ´atomos obFigura 3: Espectro de energ´ıa para distintos ra- jetivos. Este tipo de interacci´on no producen yos. rayos x, sino electrones secundarios y eventualmente calor. Es m´as, el 99 % de la energ´ıa que Entonces, los fotones de rayos x son produci- entra se convierte en calor. dos cuando una sustancia es bombardeada por electrones a alta velocidad. Cuando estos elec- Por tanto, si bien aqu´ı se producen tres intertrones interact´ uan con el target, ocurren mu- acciones, me centrar´e en tres.

Figura 4: Interacciones del electr´ on con el target y su relaci´on con el espectro de rayos x. (a) corresponde a la radiaci´ on de frenado generada cuando electrones a altas velocidades son desacelerados por el campo el´ectrico del n´ ucleo del targe. (b) es la radiaci´on caracter´ıstica producida cuando los electrones de altas velocidades interact´ uan con las capas el´ectricas del targe, y (c) es cuando el electr´ on a altas velocidades golpear directamente al n´ ucleo, convirtiendo toda su energ´ıa cin´etica en energ´ıa de rayos x. Se utiliz´o un tubo de rayos x a potencial de 120 [kV ] con filtraci´ on adicional para remover los fotones de bajas energ´ıas [3]. las variables. Los electrones son acelerados por una diferencia de potencial (como se vio reci´en) en un tubo de rayos x. En este caso, puedo controlar la energ´ıa m´ axima del espectro de rayos x resultante, y podemos controlar esto tal como poner un interruptor en on u off cuando se producen los rayos. No obstante, con los rayos gamma, esto no es tan f´ acil, ya que la radiaci´ on proviene del n´ ucleo de a ´tomos radioactivos y no se puede simplemente poner en on u off. S´ olo puedo apantallar la radiaci´ on.

3

2.1.1.

particular rango de energ´ıa en funci´on de su distancia al n´ ucleo.

Radiaci´ on de frenado

O llamado tambi´en Bremsstrahlung, del alem´an ’frenado’. Corresponde a la parte (a) Las capas electr´onicas son numeradas correlatide la figura 4. vamente, partiendo de la m´as cercana al n´ ucleo, y se identifican mediante letras: Imaginemos que vamos en una competencia corriendo muy, muy r´ apido y hay un tumulto de Capa K, n = 1. Es la m´as interior, pregente que va aumentando su n´ umero a medisente en todos los elementos qu´ımicos (2 da que vamos llegando al final de la carrera. electrones m´ax.) Cuando reci´en nos topemos con la gente nos Capa L , n = 2 (8 electrones m´ax.) costar´a poco o nada seguir corriendo, pero a medida que vayamos llegando al n´ ucleo (mayor Capa M , n = 3 (18 electrones m´ax.)” [4] tumulto de gente) nos costar´ a cada vez m´as poder seguir avanzando, gastando mucha m´as Capa N´ umero cu´antico Electrones energ´ıa para poder seguirnos moviendo. Este es principal n m´aximos electr´onica el esp´ıritu de la radiaci´ on de frenado, y como K 1 2 podemos notar en el eje x de la figura 4(a), la L 2 8 energ´ıa va aumentando a medida que los elecM 3 18 trones interact´ uan con el campo el´ectrico del n´ ucleo, desaceler´ andose. Cuadro 2: Capas electr´onicas, con su nombre, n´ umero cu´antico principal y electrones m´ axiEn este caso, la intensidad de radiaci´on I mos por capa. I∝

Z 2 z 4 e6 m2

(1)

, donde Z es el n´ umero at´ omico, z es el valor cuantizado de la carga el´ectrica, e la carga del electr´on tal que Q = ze, donde Q es la carga, y m la masa de la part´ıcula cargada. Esta ecuaci´ on indica que un electr´ on es 3 millones de veces m´ as eficiente al generar radiaci´on de frenado que una part´ıcula masiva como protones o part´ıculas alfa. Esta es la raz´ on pr´actica Figura 5: Capas electr´onicas, o en su terminopor la que se utilizan electrones de altas velo- log´ıa en ingl´es: shell. cidades. Al chocar el electr´on con un electr´on de una Notemos adem´ as que I aumenta r´ apidamente capa interior, ´este u ´ltimo se libera. Por tanto, con el n´ umero at´ omico Z. cuando el agujero por el electr´on faltante en la capa es llenado por un electr´on de una ca2.1.2. Colisi´ on con las capas electr´ oni- pa exterior, se emite radiaci´on. En el modelo de B¨ohr - recordemos -, los electrones ocupan ´ orbicas tas con niveles cu´anticos de energ´ıa espec´ıficos. “Una capa electr´ onica, capa de electrones o cu- Por ejemplo, las energ´ıas de las capas K, L, M , bierta de electrones puede pensarse como el y N del tungsteno son 70, 11, 3, y 0, 5 [keV ], conjunto de ´ orbitas seguidas por un grupo de respectivamente. Entonces, cuando un electr´ on electrones alrededor del n´ ucleo de un ´atomo. de una capa K es liberado y el agujero es lleCada capa puede contener un cierto n´ umero nado por un electr´on de la capa L, se genera m´aximo de electrones, y est´ a asociada con un un fot´on de rayos x de 59 [keV ]. De manera 4

2.2.1.

similar, cuando un electr´ on de la capa M se mueve a la capa K, se produce un fot´on con energ´ıa 67 [keV ]. Esto es ilustrado en los peaks caracter´ısticos de la figura 4(b). Notemos que cada elemento en la tabla peri´ odica tiene una u ´nica energ´ıa de uni´ on de capas, por lo que las energ´ıas caracter´ısticas de rayos x ser´an u ´nicas por cada ´ atomo [3].

2.1.3.

Efecto fotoel´ ectrico

El fot´on de rayos x, o fotoelectr´on, incidente dona toda su energ´ıa para liberar un electr´ on de una capa del ´atomo. Por tanto, el fot´on deja de existir. Ahora bien, de lo expuesto anteriormente, el agujero dejado es llenado por un electr´on de una capa exterior, generando radiaci´on; produciendo un ion positivo. Para tejidos (o parecidos) [3], la energ´ıa de uni´on de un electr´on en la capa K es muy peque˜ na (500 [eV ] aprox), de forma que el fotoelectr´on se queda con casi toda esta energ´ıa y, dado que es muy peque˜ na, no alcanza a viajar muy lejos (menos que las dimensiones de una c´elula humana) antes de que decaiga. Es m´ as, casi todos los rayos x producidos en pacientes por efecto fotoel´ectrico son reabsorbidos.

Colisi´ on directa con el n´ ucleo

El tercer fen´ omeno para la producci´ on de rayos x ocurre cuando el electr´ on a altas velocidades choca directamente con el n´ ucleo y su energ´ıa aparece completamente como radiaci´on de frenado. La energ´ıa de rayos x producida por esta interacci´on representa el l´ımite superior del espectro de rayos x. La probabilidad de que ocurra esta colisi´ on es casi nula, por lo que la mag- Es posible demostrar [5] que la probabilidad de nitud ser´a muy baja en el espectro. interacci´on P es

2.2.

Interacci´ on de rayos x con la ma, o bien teria

Pphotoelectric ∝ E −3

(2)

Pphotoelectric ∝ Z −3 (3) La energ´ıa t´ıpica que se ocupa para fotones de rayos x, en aplicaciones de CT, es aproximada¿Y qu´e significa esto? Que consecuentemente, mente 20 [keV ] a 140 [keV ]. tejidos con peque˜ nas diferencias en su n´ umero En la interacci´ on de los rayos x con la materia at´omico producir´an grandes diferencias en las aparecen tres fen´ omenos sumamente interesan- probabilidades de efectos fotoel´ectricos. Esto, tes, y muy conocidos por todos los f´ısicos: efecto en consecuencia, resulta en diferentes tasas de fotoel´ectrico, efecto Compton, y dispersi´on de absorci´on de los fotones de rayos x y lleva a grandes contrastes entre diferentes tejidos. Rayleigh.

Figura 6: Ilustraci´ on de interacciones con efecto fotoel´ectrico.

5

Figura 7: Dibujo esquem´ atico del efecto Compton. El fot´on de rayos x choca y libera un electr´ on y un fot´on. 2.2.2.

Efecto Compton

te. Dado que una peque˜ na porci´on de la energ´ıa del fot´on es absorbida, la radiaci´on absorbida El efecto Compton consiste en el aumento de por el paciente es mucho menor que aquella obla longitud de onda de un fot´ on cuando este tenida en el efecto fotoel´ectrico. choca con un electr´ on libre y pierde parte de su energ´ıa. ¿Y qu´e hay de la probabilidad de que ocurra el efecto Compton? Bueno, la probabilidad deh ∆λ = (1 − cos θ) pende de la densidad de electrones del mateme c umero at´omico Z. La falta de deNotemos que la variaci´ on que se produce en la rial, no del n´ pendencia con respecto al n´ umero at´omico Z, longitud de onda s´ olo depende del ´angulo de dispersi´on, donde h es la constante de Planck, provee poca o nada de informaci´on con respecme la masa del electr´ on, c la velocidad de la luz to al contraste entre diferentes tejidos, ya que la diferencia de densidades de electrones entre y θ el ´angulo de dispersi´ on. diferentes tejidos es muy peque˜ na, por tanto, Esta es una interacci´ on muy importante al nivel debemos minimizar esto. de los tejidos. En esta interacci´ on, la energ´ıa del fot´on de rayos x incidente es considerablemente 2.2.3. Dispersi´ on de Rayleigh mayor que la energ´ıa de uni´ on del electr´on. Un ´ltima interacci´on corresponde a la interacfot´on incidente golpea un electr´ on y lo libera La u del ´atomo. El fot´ on incidente de rayos x es de- ci´on menos importante para los investigadores flectado con una p´erdida parcial de su energ´ıa del ´area: dispersi´on de Rayleigh, o scattering inicial, como muestra la figura 7. Por tanto, una coherente. interacci´on de Compton produce un ion positi- En este caso, no hay energ´ıa que se conviervo, un electr´ on y un fot´ on dispersado. ta en energ´ıa cin´etica, as´ı que no hay ionizaci´on. Dado que no hay energ´ıa que se convierta en energ´ıa cin´etica, no nos interesa mucho. No obstante, avances recientes permitir´ıan utilizar la dispersi´on de Rayleigh para caracterizaci´ on de huesos [6] [7].

Notemos que θ puede ir de 0 a 2π, por lo que la informaci´ on que nos puede proporcionar el efecto Compton es muy baja con respecto al lugar de donde se llevaba a cabo la interacci´on, debido a la gran posibilidad que tiene el fotoelectr´on de moverse luego del choque.

Si deseamos entender la importancia que tieEn t´erminos de lo que nos interesa: las colisio- ne cada una de estas interacciones, entonces nes no son u ´nicas, por lo que hay muchas coli- conviene realizar un gr´afico del porcentaje de siones antes de que el fot´ on abandone al pacien- interacci´on versus la energ´ıa para cada efecto. 6

Figura 8: Gr´ afica del porcentaje de interacci´on versus la energ´ıa para el efecto fotoel´ectrico, efecto Compton y dispersi´ on de Rayleigh. La gr´afica se ha realizado en agua [8]. puede ser entendido por el hecho de que m´ as energ´ıa puede ser transferida usando efecto fotoel´ectrico que usando efecto Compton.

Notemos de la figura 8 que el porcentaje de interacci´on del efecto fotoel´ectrico decrece a medida que aumenta la energ´ıa. Por otro lado, el porcentaje de interacci´ on del efecto Compton aumenta mientras hay m´ as energ´ıa.

Podemos dividir la gr´afica en tres secciones. La primera zona cubre el espectro de energ´ıa hasta los 50 [keV ]. En esta zona, la absorci´on via efecto fotoel´ectrico domina. En la siguiente zona, representada por la regi´on oscura (entre 50 y 90 [keV ]), tanto el efecto fotoel´ectrico como Compton son importantes. En la tercera zona, es el remanente entre 90 y 150 [keV ]. Aqu´ı domina el efecto Compton.

Ahora bien, m´ as importante a´ un, conviene realizar una gr´ afica de energ´ıa transferida como una funci´ on de la energ´ıa de los fotones de rayos x. Como se ve de la figura 9, la mayor parte de la energ´ıa de los rayos x a bajas energ´ıas se transfiere via efecto fotoel´ectrico. Este fen´omeno

Figura 9: Porcentaje de energ´ıa transferida versus energ´ıa. 7

La figura 10 es la gr´afica de µ versus la energ´ıa, medida en [keV ]. La gr´afica fue hecha para yodo, hueso, agua y tejido blando.

Varias investigaciones y desarrollos en energ´ıa dual est´an basadas en las diferentes zona donde domina cada efecto. Para una lectura m´as profunda de esto, ver cap´ıtulo 12 de la referencia 3.

De las conclusiones que podemos obtener para la curva:

Estos tres efectos realizan algo as´ı como un efecto red, ya que algunos de los fotones son absorbidos o dispersados. En otras palabras, son los fotones de rayos x son atenuados cuando pasan por un material. La atenuaci´on, como bien sabemos, puede ser expresada por un decaimiento exponencial (para el caso monocrom´atico). Esta ley tambi´en es bien conocida por los f´ısicos: corresponde a la ley de BeerLambert.

2.3.

El agua y los tejidos blandos tienen un ´ındice de atenuaci´on sumamente similar. Por tanto, el agua se utiliza para calibraci´on de CT para as´ı asegurar precisi´ on de los n´ umeros de CT para tejidos blandos. Segundo, el yodo tiene un coeficiente de atenuaci´on alto; en consecuencia, se utiliza como agente de contraste, por lo que podemos - por ejemplo - ver los vasos sangu´ıneos m´as opacos en una angiograf´ıa. En ellas, el contraste basado en yodo es seguido por CT scans para investigar la integridad de los vasos e identificar la posible presencia de stenosis. Una inspecci´on m´as de cerca de la atenuaci´on de la curva del yodo muestra un repentino cambio en el coeficiente de atenuaci´ on a los 33, 2 [keV ], el cual es la energ´ıa de ligaz´on2 de los electrones de la capa K. A la energ´ıa de rayos x un poquito menor a 33, 2 [keV ], el coeficiente de atenuaci´ on del yodo es 32, 31 [cm−1 ], mientras que la atenuaci´on salta a 176, 69 [cm−1 ] para las energ´ıas de rayos x por poco menores a 33, 2 [keV ]. Este fen´omeno es llamado “K-edge”3 para el yodo. Los K-edge son diferentes para distintos materiales. Por ejemplo, el K-edge es 69, 5 [keV ] para el tungsteno.

Ley de Beer-Lambert

Podemos expresar este decaimiento del que hemos hablado como I = I0 e−(τ +σ+σr )L

(4)

, donde I e I0 son las intensidades incidentes y transmitidas, respectivamente. L es el espesor del material y τ, σ, σr son los coeficientes para el efecto fotoel´ectrico, Compton y dispersi´on de Rayleigh del material, respectivamente. Si defino µ = τ + σ + σr , entonces: I = I0 e−µL

(5)

, donde µ es el coeficiente de atenuaci´on lineal del material. La ecuaci´ on anterior es la llamada ley de Beer-Lambert.

Tercero, la curva correspondiente a los huesos tambi´en tiene una discontinuidad, pero en 4.0[keV ] para los huesos. Esto se debe al K-edge propio del calcio.

Notemos que µ ser´ a una funci´ on de los fotones de rayos x incidentes, y que queremos que µ sea lo m´as alto posible. 2 3

La energ´ıa que se requiere para desarmar un sistema en partes separadas. Algo as´ı como peaks en la capa K.

8

Figura 10: Coeficientes de atenuaci´on lineales para distintos materiales. Proyecci´ on como integral de l´ınea Consideremos ahora el caso de un objeto no uniforme (i.e., un objeto hecho de distintos maEntonces tenemos la ley de Beer-Lambert (5). teriales, cada uno con distintos coeficientes de Supongamos que los rayos x son monoenerg´etiatenuaci´on lineales) pero con el mismo espesor. cos y que el espesor del cuerpo es ∆x. As´ı, (5) Es decir, dividimos la torta en partes iguales. puede ser expresada como 2.3.1.

I = I0 e−µ∆x

(6)

Figura 11: Ilustraci´ on de la atenuaci´on de un material para un rayo x monocrom´atico. (a) La atenuaci´ on de rayos x monoenerg´eticos en objetos uniformes viene dada por la ley de BeerLambert. (b) Cualquier cuerpo no uniforme puede ser dividido entre m´ ultiples elementos. Dentro de cada elemento podemos asumir la existencia de coeficientes de atenuaci´on. La ley de BeerLambert puede seguir siendo aplicada, pero un poco m´as distinta que antes. 9

3.1. I = I0 e−µ1 ∆x e−µ2 ∆x e−µ3 ∆x . . . e−µn ∆x Pn  = I0 e− i=1 µi ∆x ln ! n X ln I = ln I0 − µi ∆x

Preliminares

Consideremos la imagen:

i=1

 g ≡ − ln

I I0

 =

n X

µi ∆x

i=1

La cantidad g es llamada proyecci´ on, y la interpretaci´on f´ısica de esto ser´ a vista en la secci´on matem´atica de este art´ıculo. Por el momento, qued´emonos con la definici´ on de que la proyecci´on g viene dada por la relaci´ on  g = − ln

I I0

 =

n X

µi ∆x

(7)

i=1

¿Pero qu´e pasa si ∆x → 0, es decir, si el espesor se hace diferencial? Bueno, entonces claramenFigura 12: Arriba: Rayos x paralelos incidiendo te sobre un cuerpo. Abajo: Perfil de absorci´ on.   Z I g = − ln = µ(x)dx (8) I0 L De la figura (12), notemos que si los rayos x El problema de esto es que no conocemos µ. vienen de manera paralela uno con respecto Es mas, si lo conocieramos tendr´ıamos la ate- a otros, entonces si yo realizo una especie de nuaci´on en cada parte del cuerpo, lo que nos transformada de Fourier, es decir, si me paso permitir´ıa reconstruir la imagen. M´as a´ un: si a otro espacio, podr´ıa tener algo as´ı como el tenemos µ entonces (suponiendo que pudi´era- perfil de absorci´on. Este perfil de absorci´ on es mos resolver la integral), ser´ıamos capaces de la transformada de Radon. Observemos el siobtener la proyecci´ on. Es decir, si obtenemos guiente ejemplo realizado en MATLAB para un proyecci´on, resolvemos el problema. c´ırculo dibujado en PAINT. Otra manera de ver la proyecci´ on es que es la suma continua de todas los coeficientes de atenuaci´on en una sola l´ınea del objeto.

3.

Matem´ atica de la CT scan

Todo parte en 1917 cuando Johann Radon, en marco de la geometr´ıa integral, elabora lo que llamamos hoy d´ıa como transformada de Radon. La verdad es que estudiar esto desde un pun- Figura 13: C´ırculo blanco dibujado en paint. to matem´ atico es sumamente dif´ıcil, por lo que dar´e una pincelada desde el punto de vista f´ısi- Luego, los perfiles de absorci´on (o m´as bien co. transformadas de Radon) se ven como: 10

observemos la figura (17).

0°

Entonces, hagamos el experimento con dos c´ırculos. 90°

Figura 14: Transformada de Radon a 0◦ y 90◦ respectivamente. Podemos pensar el perfil de absorci´ on de la siguiente manera: imaginemos que tenemos a un mago frente a nosotros y detr´ as de ´el mantiene oculta nuestra carta en su mano. Si nos quedamos frente a ´el (0◦ ) no veremos la carta. Si nos movemos hacia su costado (90◦ ) ya podremos ver que tiene algo oculto. Si nos movemos y quedamos tras el mago (180◦ ), entonces claramente podremos ver la identidad de nuestra carta, la forma en la que la tiene sujeta e incluso si esta tiene arrugas o est´ a estropeada. La misma esencia se puede aplicar a la transformada de Radon. Si tenemos un c´ırculo frente a nosotros, el perfil de absorci´ on ser´a el mismo tanto a 0 como 90◦ , ya que es un c´ırculo; es igual en cualquier parte que se le mire. Para un an´alisis m´ as profundo de manera de asociarlo directamente con el concepto de tomograf´ıa,

Figura 15: Dos c´ırculos. 0°

90°

Figura 16: Transformada de Radon a 0◦ y 90◦ . El problema de la superposici´on es evidente: a 90◦ se genera la superposici´on de los dos c´ırculos, lo que genera que veamos en el perfil de absorci´on un mont´ıculo.

Figura 17: Concepto general de una tomograf´ıa. (a) Cuando se ve un objeto semitransparente desde un ´ angulo, algunas esferas internas del objeto se superponen. Dado que no conocemos a priori lo opaco de cada esfera, es imposible para nosotros estimar la forma, intensidad y n´ umero de esferas. (b) Despu´es de que el objeto es rotado, las sombras producidas por las esferas ya no se superponen, y cada esfera puede ser vista de manera individual. Cuando se rota una mayor cantidad de veces, la forma, intensidad y n´ umero de esferas pueden ser estimados precisamente.

11

Ahora bien, si realizo la transformada de Radon inversa de la figura (14) (el c´ırculo solo) a 90◦ , obtendr´ıamos una figura de la forma:

Figura 20: N´ umero de proyecciones en una imagen con geometr´ıa de haces paralelos. Escaner de primera generaci´on. ¿Y si aumento dram´aticamente el n´ umero de Figura 18: Transformada de Radon inversa a proyecciones? Es decir, ¿si aumento dram´ atica90◦ . mente el n´ umero de ´angulos en que proyecto? Bueno, la imagen tambi´en cambia dram´aticamente. Uno creer´ıa que esto no nos dice nada, pero eso no es cierto. La verdad es que esto fue hecho s´olo a un ´ angulo fijo. Si realizo la transformada inversa con ´ angulo θ de 0◦ a 179◦ con salto de ◦ 20 , obtengo:

Figura 21: Reconstrucci´on de la imagen usando transformada de Radon inversa, pero ahora con un barrido de ´angulos de 0◦ a 179◦ con salto de 0, 1◦ . Uno estar´ıa tentado a decir a que la figura 21 representa fielmente la imagen 13, no obstante, observemos un peque˜ no defecto que ocurre: (si miramos muy bien) hay unas l´ıneas perpendiculares a la figura que se cruzan en todo el plano. Estos peque˜ nos defectos pueden ser desastrosos a la hora de detectar alguna enfermedad en el paciente. Y m´as a´ un, el problema Figura 19: Reconstrucci´ on de la imagen usando es a´ un m´as grave: en el medio (de la imagen 19) transformada de Radon inversa, con a´ngulo de hay una especie de halo que ilumina la imagen, 0◦ a 179◦ , pero con 20◦ de salto. impidiendo que se vea completamente.

12

Figura 22: Esquema cl´ asico que se muestra en textos especializados del tema tales como [3]. ρ1 , ρ2 , ρ3 y ρ4 son rectas que etiquetan a cada rayo x pasando por el cuerpo situado en el origen del plano xy. Si sumamos cada coeficiente de atenuaci´on de cada recta, entonces obtenemos el perfil de absorci´ on que se muestra en la parte de arriba de la figura en el eje trasladado x0 . Observemos que cada punto del perfil de absorci´on - al que se llega via transformada de Radon - y que no es m´ as que la proyecci´ on de cada punto, le corresponde una recta ρ a ´angulo θk . Un punto cualquiera ρj a ´ angulo θk del perfil de absorci´on ser´a la proyecci´on g(ρj , θk ) y en particular, si seguimos la l´ınea ρ3 , obtenemos g(ρ3 , θk ).

3.2.

Transformada de Radon

Consideremos el siguiente dibujo ilustrativo:

tronco y queremos llegar al final de este para pisar tierra firme. Si deseo caminar sobre el tronco - y no tener que pasarme a otros -, entonces - valga la redundancia - lo u ´nico que debo hacer es simplemente seguir caminando derecho y ni siquiera pensar en la idea de pasarme a un segundo tronco. S´olo seguir mi camino.

Bueno, la idea anterior es la misma que sumar sobre la l´ınea de rayos x y finalmente obtener la proyecci´on. EsR decir, si yo
logro caminar por todo el rayo x L µ(x)dx , mi premio al final Figura 23: Recta ρ a ´ angulo θ del camino ser´a obtener la proyecci´on del rayo En este caso, la recta perpendicular ρ puede ser x (g (ρj , θk )). Si sumo sobre todos los rayos x, entonces tendr´e todas las proyecciones de todos escrita en t´ermino de sus componentes: los rayos. ρ = x cos θ + y sin θ Matem´aticamente, lo anterior viene dado por la Ahora observemos el dibujo 22 para mostrar delta de Dirac δ(x − a). Todos los puntos ser´ an mejor el punto al que quiero llegar. nulos a excepci´on del punto a, que es donde la Imaginemos que estamos caminando sobre un distribuci´on se hace infinita. Con esto en mente, definir´e la transformada de Radon como: 13

Z

∞

Z

∞

I(x, y) δ(x cos θ + y sin θ − ρ) dx dy

g(ρ, θ) = −∞

(Transformada de Radon)

(9)

−∞

, donde I(x, y) es la imagen (siendo m´as espec´ıfico: funci´ on) que depende de x e y. Notemos de la definici´ on que la delta garantiza que siempre estemos dentro de la recta, ya que en la primera parte de la delta s´ olo estamos escribiendo ρ de otra manera simplemente: en funci´on de su ´ angulo.

3.3.

Fourier slice theorem

Ahora determinaremos la relaci´on entre la transformada de Fourier de cada una de las proyecciones y la imagen original que estamos tratando de reconstruir.

Recordemos que la transformada de Fourier viene definida por: Si bien se asumen ciertos conocimientos preZ ∞ vios, nunca est´ a de sobra mencionar la propief (t)e−iωt dt (10) F [f (t)] = dad m´as importantes de la delta: −∞ Z

Consideremos un ´angulo fijo de proyecci´ on θ y tomemos la transformada de Fourier 1D con respecto a ρ.

∞

f (x)δ(x − a)dx = f (a) −∞

Por ende - repito -, la delta ser´ a nula siempre que estemos fuera de la recta, es decir, tal que x cos θ +y sin θ 6= ρ. De esta forma, se garantiza que estemos recorriendo el rayo x incidente.

En este caso, ρ act´ ua como el “dominio temporal”, por lo que nos estar´ıamos pasando al “espacio de frecuencias” ω, todo esto a θ fijo. En la imagen (25) podemos ver este cambio que genera la transformada de Fourier.

En el mundo real, las integrales son sumas. Ahora bien, si visualizamos g(ρ, θ) como una imagen, entonces esto recibe el nombre de sinograma.

Figura 25: Cambio de espacios originado por la transformada de Fourier. As´ı: Z

∞

G(ω, θ) =

g(ρ, θ)e−i2πωρ dρ

(11)

−∞

Figura 24: Sinograma para la figura de dos c´ırculos. En el eje de las abscisas va el ´angulo θ y en el eje de las ordenadas el nuevo eje x0 .

Recordando ahora la transformada de Radon (ecuaci´on (9)), si reemplazo g(ρ, θ) escrito como la transformada de Radon (9) y la reemplazo en la expresi´on anterior de la transformada

14

conectar algo que conocemos - la proyecci´ on de de Fourier (11), obtenemos: la imagen - con algo que no conocemos - la imaZ ∞Z ∞ Z ∞ gen que deseamos obtener. ¿Y qu´e conecta esI(x, y) G(ω, θ) = −∞ −∞ −∞  (12) tas dos cosas? La transformada de Fourier. La transformada de Fourier act´ ua como un puente −i2πωρ ×δ(x cos θ + y sin θ − ρ)dxdy e dρ entre lo que conocemos y lo que no. Luego: Z ∞Z ∞ I(x, y) G(ω, θ) = −∞ −∞ Z ∞ × δ(x cos θ + y sin θ − ρ)e−i2πωρ −∞  × dρ dxdy

(13)

Recordando las propiedades de la delta, evaluamos ρ como x cos θ + y sin θ. Por tanto: Z ∞Z ∞ G(ω, θ) = I(x, y)e−i2πω(x cos θ+y sin θ) −∞

−∞

× dxdy (14) Sea u = ω cos θ y v = ω sin θ. Entonces: Z ∞Z ∞ G(ω, θ) = I(x, y)e−i2π(ux+vy) dxdy −∞

−∞

(15) La ecuaci´ on anterior es claramente la transformada de Fourier 2D sobre la l´ınea u = ω cos θ y v = ω sin θ. Esto es lo que se conoce como el Fourier slice theorem o projection-slice theorem. Es m´ as, resulta interesante ver que si aplicamos la tranformada de Fourier inversa podemos obtener la imagen I(x, y) en t´erminos de la proyecci´ on G(ω, θ).

Figura 26: Gr´afico v versus u de la transformada de Fourier 2D. Si realizamos un gr´afico v vs u, que ser´ıan las componentes de la transformada de Fourier 2D, entonces obtendr´ıamos la imagen y sus proyecciones, que ser´ıan las l´ıneas rectas dentro de la imagen, en la figura 26. No obstante, como se ve en la figura, se producir´ıa el ya mencionado halo en el medio. La raz´on de esto es que la parte central de la transformada de Fourier 2D est´a sobremuestreada, generando que haya demasiado “peso” en la parte de al medio. Para poder “bajar” este “peso”, viene a mano el algoritmo llamado retrproyecci´ on filtrada, o en su nombre en ingl´es: filtered backprojection, ya que las bajas frecuencias son el problema.

Resumiendo entonces: la transformada de Fourier en 1D de la proyecci´ on de la imagen, es la transformada de Fourier en 2D de la imagen. 3.4.

Filtered backprojection

Repito de otra forma, aunque quiz´ as sea lo mismo: La transformada de Fourier (en una di- Si renombro G(ω, θ) = F (u, v), entonces, al mensi´on) de la proyecci´ on, es lo mismo que la aplicar la transformada de Fourier inversa a transformada de Fourier en (dos dimensiones) (15), obtenemos: de la imagen en s´ı. Z ∞Z ∞ Entonces, hemos desarrollado el formalismo paF (u, v)ei2π(ux+vy) dxdy (16) ra reconstruir una imagen ya que hemos podido I(x, y) = −∞

4

Algo as´ı como “rdrdθ”.

15

−∞

Si cambiamos a coordenadas polares, el Jacobiano dudv = ωdωdθ4 , entonces: Z 2π Z ∞ F (ω cos θ, ω sin θ) I(x, y) = (17) 0 −∞ i2π(x cos θ+y sin θ)ω ×e ωdωdθ Notemos que F (ω cos θ, ω sin θ) = F (u, v) = G(ω, θ). El cambio s´ olo lo hice para que se pudiera ver mejor el cambio a coordenadas polares.

2. Multiplicar cada transformada de Fourier (G(ω, θk )) por |ω|.

3. Calcular la transformada de Fourier inversa del paso anterior.

4. Integrar (en el mundo real: sumar) sobre todos los ´angulos para obtener I(x, y).

Entonces: 2π

Z

Z

∞

G(ω, θ)ei2π(x cos θ+y sin θ)ω

I(x, y) = −∞

0

× ωdωdθ π

Z

Z

∞

|ω| G(ω, θ)ei2π(x cos θ+y sin θ)ω

= −∞

0

× dωdθ (18) , consideramos valor absoluto por el hecho que G(ω, θ + π) = G(−ω, θ). As´ı: Z π Z

∞

|ω| G(ω, θ)

I(x, y) = 0

×e

−∞ i2πωρ

dω

ρ=x cos θ+y sin θ

 dθ

Figura 27: Multiplicaci´on de la funci´on filtro |ω| por una funci´on ventana h(ω). El resultado (19) es una ventana de Hamming (“Hamming window”).

La raz´on de a˜ nadir el filtro es poder eliminar las bajas frecuencias y tambi´en preveer posibles errores num´ericos. Me explico: si tenemos una funci´on peri´odica y realizo la transformaci´on de Fourier, es posible que el programa que est´e ocupando (en este caso, MATLAB) arroje un error. ¿Por qu´e? Porque la funci´on, al termiEsto es lo que se llama retroproyecci´ on filtrada. nar su dominio, vuelve a repetirse nuevamente, “Retro”, porque aplicamos la transformada de pero el problema es que en este punto de reFourier inversa, y “filtrada” porque est´a filtra- petici´on la funci´on puede que no sea peri´odica, da por la funci´ on |ω|. arrojando errores num´ericos. Notemos que si |w| no aparece en la integral, tendr´ıamos la transformada de Fourier inversa de la proyecci´ on G(ω, θ). Entonces, la ecuaci´on anterior es la transformada de Fourier inversa de G(ω, θ), pero multiplicada por una funci´on filtro |ω|.

Resumiendo entonces, si deseo determinar la Con el Hamming se soluciona el problema del imagen I(x, y), la forma del slice del paciente halo en el medio. y del sector que deseamos estudiar, debemos realizar el siguiente algoritmo: Ahora vayamos a un ejemplo concreto: conside1. Determinar la transformada de Fourier remos una imagen simulada de un corte cef´ alico 1D de cada proyecci´ on. transversal de la cabeza de un paciente. 16

4.

Figura 28: Simulaci´ on de un corte cef´alico transversal de la cabeza. Aplicando entonces lo anteriormente aprendido podemos obtener la reconstrucci´ on de la imagen sin ninguna clase de error - a excepci´on de ciertos problemas de resoluci´ on - eliminando las extra˜ nas l´ıneas perpendiculares de la figura 23, gracias al Hamming. Obteniendo as´ı:

Figura 29: Reconstrucci´ on final de la imagen. No obstante, se realiz´ o una interpolaci´on para dar una mayor densidad de puntos y obtener una mejor resoluci´ on. Problema finalizado.

Resumen

Un electr´on movi´endose a altas velocidades es la piedra angular de la creaci´on de rayos. La creaci´on de los mismos, involucra la radiaci´ on de frenado, la colisi´on con las capas electr´onicas y la colisi´on directa del electr´on con el n´ ucleo. Al interactuar los rayos x con la materia, se originan tres fen´omenos: efecto fotoel´ectrico, efecto Compton y dispersi´on de Rayleigh. Dentro de la ciencia m´edica, resulta m´as importante el efecto fotoel´ectrico y efecto Compton, aunque tampoco debe descartarse el potencial que la dispersi´on de Rayleigh puede tener para futuros avances. No obstante, es interesante mencionar que los tres efectos tienen algo en com´ un: un efecto red. Dicho efecto red viene modelado por la ley de Beer-Lambert, de la cual, para un cuerpo con espesor infinitesimal, es posible encontrar una expresi´on logar´ıtmica negativa y una expresi´on asociada a la integral del coeficiente de atenuaci´on lineal. Estas dos expresiones son equivalentes a la proyecci´on de la imagen. Para determinar la proyecci´on de la imagen en una sola l´ınea - y en una geometr´ıa de haces paralelos - es necesario definir una herramienta matem´atica llamada transformada de Radon. La transformada de Radon est´ a intimamente relacionada con la transformada de Fourier: si se realiza la transformada de Fourier unidimensional de la proyecci´on, se obtiene una identidad equivalente a la transformada de Fourier bidimensional de la imagen misma. Este teorema recibe el nombre de Fourier slice theorem. En virtud de lo anterior, si realizamos la transformada de Fourier inversa, es posible determinar una expresi´on directa para la imagen en funci´on de su proyecci´on. No obstante, para evitar dificultades del sobremuestreo de la imagen con respecto a las partes de baja frecuencia, es necesario filtrar la transformada a˜ nadiendo una funci´on filtro, generando lo que se conoce como ventana de Hamming. Con la ventana de Hamming es posible eliminar el problema del halo en el centro y reconstruir fidedignamente - al integrar en θ - la imagen deseada.

17

Figura 30: Tipos de geometr´ıa utilizados en los esc´aneres.

5.

Lo que falt´ o

Este trabajo fue s´ olo un esbozo cualitativo de lo que es la f´ısica y la matem´ atica de la tomograf´ıa computarizada. Faltaron mencionar muchas cosas, como por ejemplo la linealidad de la transformada de Radon o mencionar que la transformada inversa viene dada por el teorema de Plancherel.

el lector desea profundizar en esto, existe un excelente software hecho por Kevin M Rosenberg y puede ser descargado de forma gratuita www.ctsim.org para ilustrar mejor todo esto.

No obstante, no se han abordado problem´aticas m´as importantes, tales como: 1. La geometr´ıa con la que se trabaj´o. Notemos que la geometr´ıa utilizada fue dada por una geometr´ıa de haces paralelos. Existen distintas clases de geometr´ıa: fan beam y cone beam (tipo ventilador y de cono. Ver figura 30). Sin embargo, este problema nos lleva al siguiente punto sin mencionar: 2. La generaci´ on de tomograf´ıa computarizada. Existen siete clases de generaciones, pero en la actualidad se emplean en su mayor parte m´ aquinas de tercera y cuarta generaci´ on [2]. Para dar una idea, la imagen 20 representa a una esc´aner de primera generaci´ on, la de figura 31 a uno de segunda, la figura 2 a uno de tercera y la figura 32 a uno de cuarta generaci´ on. Notemos que en los esc´aneres de tercera y cuarta generaci´ on, la geometr´ıa no es de haces paralelos. Ahora bien, si 18

Figura 31: Escaner de segunda generaci´ on.

3. Hsieh, Jian, “Computed Tomography, Principle, Design, Artifacts and Recent Advances”, second edition. 4. Link de Wikipedia: Capa electr´onica. 5. C. L. Morgan, “Basic Principles of Computed Tomography”, University Park Press, Baltimore (1983). 6. D. L. Batchelar, W. Dabrowski, and I. A. Cunningham, “Tomographic imaging of bone composition using coherently scattered x rays,” Proc. SPIE 3977, 353?361 (2000).

Figura 32: Escaner de cuarta generaci´on.

6.

Referencias 1. Krestel A C and Slanney M 1988 “Principles of Computerized Tomographic Imaging (New York: IEEE)”. 2. Michael, Greg. “X-ray computed tomography”. The Centre for Medical, Health and Environmental Physics, Queensland University of Technology, GPO Box 2434, Brisbane 4074, Australia.

19

7. M. S. Westmore, A. Fenster, and I. A. Cunningham, “Tomographic imaging of the angular-dependent coherent-scatter cross section,” Med. Phys. 24 (1), 3?10 (1997). 8. H. E. Johns and J. R. Cunningham, The Physics of Radiology, Charles C. Thomas Publisher, Ltd., Springfield, IL (1983).