Métodos iterativos de reconstrucción tomográfica en SPECT Carles Maria FaIcon FaIcon
ADVERTIMENT. La consulta d’aquesta tesi queda condicionada a l’acceptació de les següents condicions d'ús: La difusió d’aquesta tesi per mitjà del servei TDX (www.tesisenxarxa.net) ha estat autoritzada pels titulars dels drets de propietat intel·lectual únicament per a usos privats emmarcats en activitats d’investigació i docència. No s’autoritza la seva reproducció amb finalitats de lucre ni la seva difusió i posada a disposició des d’un lloc aliè al servei TDX. No s’autoritza la presentació del seu contingut en una finestra o marc aliè a TDX (framing). Aquesta reserva de drets afecta tant al resum de presentació de la tesi com als seus continguts. En la utilització o cita de parts de la tesi és obligat indicar el nom de la persona autora.
ADVERTENCIA. La consulta de esta tesis queda condicionada a la aceptación de las siguientes condiciones de uso: La difusión de esta tesis por medio del servicio TDR (www.tesisenred.net) ha sido autorizada por los titulares de los derechos de propiedad intelectual únicamente para usos privados enmarcados en actividades de investigación y docencia. No se autoriza su reproducción con finalidades de lucro ni su difusión y puesta a disposición desde un sitio ajeno al servicio TDR. No se autoriza la presentación de su contenido en una ventana o marco ajeno a TDR (framing). Esta reserva de derechos afecta tanto al resumen de presentación de la tesis como a sus contenidos. En la utilización o cita de partes de la tesis es obligado indicar el nombre de la persona autora.
WARNING. On having consulted this thesis you’re accepting the following use conditions: Spreading this thesis by the TDX (www.tesisenxarxa.net) service has been authorized by the titular of the intellectual property rights only for private uses placed in investigation and teaching activities. Reproduction with lucrative aims is not authorized neither its spreading and availability from a site foreign to the TDX service. Introducing its content in a window or frame foreign to the TDX service is not authorized (framing). This rights affect to the presentation summary of the thesis as well as to its contents. In the using or citation of parts of the thesis it’s obliged to indicate the name of the author.
Departament de Física Aplicada i Electrónica Universitat de Barcelona Programa de Micro y Optoelectrónica Física. Bienio 1993-1995
Métodos iterativos de reconstrucción tomográfica en SPECT
Caries Maria Falcon Falcon
Barcelona, septiembre de 1999. Memoria presentada para optar al título de Doctor en Ciencias Físicas. Realizada bajo la dirección de: Dr. Doménec Ros i Puig Dr. Ignasi Juvells i Prades
Agradecimientos: Quisiera expresar mi agradecimiento al doctor Doménec Ros, del Laboratori de Biofísica i Bioenginyeria de la Facultat de Medicina y al doctor Ignasi Juvells, del Laboratori d'Óptica de la Facultat de Física, directores de esta tesis, por su muy generosa ayuda en la realización, discusión y corrección de este trabajo.
También quisiera expresar mi agradecimiento al doctor Javier Pavía, por su constante colaboración en la parte experimental y en la discusión de los resultados, al doctor David Fuster, por su colaboración en la aplicación médica del capítulo 7 y al doctor Jorge Setoain Quinquer, jefe del Servei de Medicina Nuclear del Hospital Clínic i Provincial de Barcelona y al resto de personal de este servicio por las facilidades dadas en la utilización de sus instalaciones.
Así mismo, quisiera expresar mi agradecimiento a la doctora Mar Rotger por las numerosas horas dedicadas a la resolución de problemas informáticos y al resto de del Laboratori de Biofísica i Bioenginyeria, el doctor Daniel Navajas, el doctor Ramón Farré, Deborah, Santi, Miguel Ángel, Nuria, los Jordis, Marina, Ester y Albert por las constantes muestras de apoyo recibidas.
Finalmente, quisiera expresar mi agradecimiento a mi familia y amigos. Sin su apoyo incondicional este trabajo no podría haberse llevado a cabo.
Esta tesis ha sido parcialmente financiada la Comisión Interministerial de Ciencia
y
Tecnología
(CICYT)
SAF96/0062 y SAF99/0137
mediante
los
proyectos
SAL91/314,
índice
1. Introducción 1.1. Reconstrucción tomográfica de imágenes en SPECT
1 1
1.2. Objetivos de la tesis
2
1.3. Estructura de la tesis
3
2. Planteamiento del problema 2.1. Descripción del proceso de obtención de las proyecciones
5 5
2.1.1. La emisión 2.1.2. La interacción de la radiación con el medio material
5 8
2.1.3. La detección 2.2. Planteamiento matemático de la reconstrucción tomográfica
10 14
2.2.1. Descripción analítica 2.2.2. Discretización Discretización de la imagen Discretización de las proyecciones Discretización de las ecuaciones 2.2.3. Particularización al caso de SPECT 2.2.4. Interpretación estadística del proceso de formación de la proyección 2.2.5. Problemas matemáticos de la reconstrucción 2.2.6. Reconstrucción: métodos analíticos y métodos iterativos 3. Metodología 3.1. Adquisición de proyecciones experimentales 3.2. Simulación de un estudio 3.2.1. Cálculo de la matriz de transición 3.2.2. Tratamiento de la dispersión
16 18 19 20 22 24 27 28 30 35 36 37 38 38
3.2.3. Inclusión de la PSF 3.2.4. Inclusión de la atenuación
40 43
3.2.5. Empaquetamiento de la matriz de pesos 3.2.6. Simulación del ruido sobre las proyecciones
45 45
3.2.7. Cálculo de la matriz de pesos 3.2.8. Generación de las proyecciones
46 47
3.2.9. Estudios simulados 3.3. El proceso de reconstrucción 3.4. La evaluación de imágenes mediante funciones de mérito
48 52 53
3.4.1. Figuras de mérito globales 3.4.2. Figuras de méritos sobre regiones de interés
55 55
3.4.3. Figuras de mérito estadísticas
57
4. La retroproyección
filtrada
59
4.1. Formulación matemática 4.1.1. Filtros 4.1.2. Teorema de la sección central 4.1.3. Reconstrucción tomográfica por retroproyección 4.1.4. Compensación de las degradaciones Corrección del ruido mediante filtros paso-baja
59 60 60 62 63 64
Corrección de la PSF. Filtro de Metz Corrección de la atenuación. El método de Chang 4.1.5. Implementación de la FBP 4.1.6. El algoritmo IFBP 4.1.7. Rango de validez de la FBP. Valoración de las aproximaciones realizadas en su implementación 4.2. Descripción de las pruebas realizadas 4.3. Resultados obtenidos 4.4. Discusión de los resultados 5. Métodos de reconstrucción algebraicos
67 69 72 74 75 76 78 89 91
5.1. Formulación matemática 5.1.1. Método de las proyecciones de Kaczmarz 5.1.2. Método de proyecciones sobre conjuntos convexos (POCS^ 5.1.3. Limitaciones del método ART. Parámetro de relajación
91 92 95 100
5.1.4. El algoritmo utilizado 5.2. Descripción de las pruebas realizadas
103 105
5.3. Resultados obtenidos 5.4. Discusión de los resultados 6. Métodos de reconstrucción estadísticos
106 114 117
6.1. Formulación matemática 117 6.1.1. Estimadores de máxima verosimilitud 117 6.1.2. Aplicación de los estimadores de máxima verosimilitud a la reconstrucción tomográfica 6.1.3. Métodos de aceleración del proceso Parámetro de aceleración
119 123 124
Subconjuntos Ordenados 6.1.4. Rango de validez de los métodos estadísticos 6.1.5. La validación cruzada como criterio de detención del proceso iterativo
127 128 130
6.1.6. Algoritmo MLE
132
6.2. Estudio de la CVR
134
6.2.1. Pruebas realizadas
135
6.2.2. Resultados obtenidos
135
6.2.3. Discusión de los resultados 6.3. Estudio del parámetro de aceleración 6.3.1. Pruebas realizadas 6.3.2. Resultados obtenidos 6.3.3. Discusión de los resultados 6.4. Subconjuntos ordenados 6.4.1. Pruebas realizadas 6.4.2. Resultados obtenidos 6.4.3. Discusión de los resultados 7. Comparación de resultados 7.1. Comparación de los resultados de los capítulos 4, 5 y 6 7.1.1. Recopilación de resultados para el modelo de Jaszczak 7.1.2. Discusión de los resultados 7.2. Ejemplo de aplicación 7.2.1. Descripción de las pruebas realizadas 7.2.2. Resultados obtenidos 7.2.3. Discusión de los resultados
140 141 141 142 150 151 151 152 160 163 164 164 167 170 170 172 177
8. Conclusiones 8.1. Capítulo 3: Simulación de proyecciones
179 179
8.2. Capítulo 4: IFBP 8.3. Capítulo 5: ART 8.4. Capítulo 6: MLE 8.5. Capítulo 7: Comparación de métodos 8.6. Publicación de los resultados de esta tesis 8.7. Desarrollo futuro de las investigaciones Apéndice A. Siglas, abreviaturas y nombres de variables
179 180 181 184 185 185 187
Bibliografía.
Capítulo 1
Introducción
pág. 1
1. Introducción 1.1. Reconstrucción tomográfica de imágenes en SPECT La Medicina Nuclear es una especialidad médica incluida e n el ámbito del diagnóstico por la imagen en q u e s e evalúa el funcionamiento d e un d e t e r m i n a d o órgano o tejido mediante el seguimiento exterior d e un f á r m a c o m a r c a d o c o n un isótopo radiactivo. La intensidad d e emisión de los fotones g a m m a provenientes del radiofármaco es proporcional a su concentración y, por tanto, indica una mayor o m e n o r actividad fisiológica. A s i g n a n d o una escala d e colores o grises al nivel de emisión en los diferentes puntos d e un órgano o tejido se obtiene una imagen q u e refleja su estado. La información extraída mediante las técnicas d e la Medicina Nuclear es un c o m p l e m e n t o cada v e z m á s importante para obtener un diagnóstico
médico
preciso.
Puede
encontrarse
un sumario
d e sus
diferentes aplicaciones en [DOM-80], [MAI-83], [MET-86] o [GEL-87]. En [FRE-84]
[PHY-87]
puede
encontrarse,
además,
una breve
reseña
histórica d e las diversas técnicas y aplicaciones.
Las técnicas e m p l e a d a s en Medicina Nuclear p u e d e n clasificarse en dos g r u p o s : planares y tomográficas (ECT: Emission
Computed
Tomography;
T o m o g r a f í a C o m p u t e r i z a d a de Emisión). En éstas últimas, se establece la distribución tridimensional del f á r m a c o a partir d e la detección d e la radiación
en
diferentes
direcciones
por
medio
de
técnicas
de
reconstrucción t o m o g r á f i c a \ A su vez, dentro de E C T pueden distinguirse dos técnicas: S P E C T Tomography-
o SPET
{Single
Photon
Emission
[Computed]
T o m o g r a f í a [Computerizada] por Emisión d e Fotón Ünico)
^ La idea original de la reconstrucción tomográfica parte de un trabajo de 1917 del matemático austríaco J. Radon [RAD-86]. En éste, se demuestra que puede determinarse exactamente un objeto bidimensional (tridimensional) a partir de sus proyecciones unidimensionales (bidimensionales) en conjunto infinito de direcciones.
Capítulo 1
Introducción
pág. 2
e n el q u e el isótopo radiactivo posee un núcleo en estado excitado q u e al reducir s u nivel de energía emite un fotón-;K y P E T (Positrón Tomography-
Emission
T o m o g r a f í a por Emisión de Positrones), e n la que el isótopo
emite un positrón q u e se aniquila, de forma casi inmediata, con electrón,
resultando
simultáneamente
dos
fotones-y opuestos
que
son
un
detectados
[TER-80]. A u n q u e P E T obtiene imágenes de
mayor
calidad, la instalación requerida en esta técnica es sensiblemente m á s compleja y cara que la utilizada en S P E C T , por lo q u e su uso es m u c h o menos común.
Debido a la existencia de una gran cantidad de factores degradantes, q u e se pueden agrupar s e g ú n su origen, e n : físicos (la emisión radiactiva, la interacción
de
los fotones
emitidos
con
el cuerpo
del paciente,
la
detección), intrínsecos al planteamiento del problema (la digitalización d e la
imagen
y
las
proyecciones)
y
inherentes
a
cada
método
reconstrucción tomográfica empleado, la calidad d e las imágenes
de en
S P E C T no es buena, obteniéndose imágenes de poca resolución y con una presencia importante de ruido.
1.2. Objetivos de la tesis El creciente uso de
ECT, ha llevado a un Importante esfuerzo
de
investigación
finales
de
desde
de
los
años
70
sobre
técnicas
reconstrucción y corrección de degradaciones que permitan mejorar la calidad d e las i m á g e n e s obtenidas. El a u m e n t o de las prestaciones de los equipos informáticos y las mejoras técnicas que han experimentado los equipos de g a m m a g r a f í a han permitido desarrollar y potenciar el uso de nuevos m é t o d o s d e reconstrucción m á s precisos. En este contexto y gracias
a
la
colaboración
entre
el
Laboratori
de
Biofísica
i
Capítulo 1
Introducción
Bioenginyería
del Departamant
el Laboratori Barcelona,
de Fisiología
d'Óptica
de la Facultat
y el Servei
de Medicina
de Barcelona
pág. 3
I de la Facultat
de
Medicina,
de Física, a m b o s de la Uníversitat Nuclear
del Hospital
Clínic i
de
Provincial
y e n m a r c a d a en los proyectos S A L 9 1 / 0 3 1 4 y S A F 9 6 / 0 0 6 2
de la C I C Y T (Comisión Interministerial de Ciencia y Tecnología), se ha abierto una línea d e investigación sobre reconstrucción tomográfica de imágenes
en
SPECT.
Esta tesis
describe
y
presenta
los
trabajos
realizados e n este ámbito. Los objetivos son: •
Implementación de una simulación numérica de proyecciones.
•
Implementación d e los métodos d e reconstrucción más básicos.
•
Implementación
de
un
método
de
evaluación
objetiva
de
las
reconstrucciones obtenidas. •
Estudio de la variación de la calidad de la imagen e n función de los parámetros d e los que d e p e n d e cada método de reconstrucción a fin de determinar sus características intrínsecas y la idoneidad de su uso.
•
Desarrollo d e una metodología q u e permita determinar e n cada caso concreto con q u é método y en q u é condiciones se obtienen mejores resultados.
Posteriormente, una vez introducidos los métodos de reconstrucción, se concretan m á s detalladamente ios objetivos particulares de c a d a uno de los estudios realizados.
1.3. Estructura de la tesis Esta tesis consta de 8 capítulos. Tras este primer capítulo de introducción, se prosigue en el capítulo 2 con una descripción detallada de la física q u e interviene e n S P E C T , así c o m o del planteamiento m a t e m á t i c o de la reconstrucción tomográfica. A continuación, en el capítulo 3 s e e x p o n e la metodología
de
los
estudios
realizados,
obtención
de
datos
(experimentales y por simulación) y la evaluación de resultados. En los tres capítulos siguientes, capítulos 4 , 5 y 6, se desarrollan y analizan tres m é t o d o s d e reconstrucción tomográfica. En c a d a caso se describe la
Capítulo 1
formulación
Introducción
matemática
y
la
pág. 4
implementación
del
método
de
reconstrucción y se realizan diversos estudios para establecer la relación entre
las características
de
la imagen
obtenida c o n
los
parámetros
variables existentes. El capítulo 7, consta d e dos partes. En la primera s e establecen las s e m e j a n z a s y diferencias entre los m é t o d o s desarrollados en los capítulos 4, 5 y 6 y se contrastan c o n los resultados existentes e n la bibliografía. En la s e g u n d a , se evalúan los resultados obtenidos e n un determinado
estudio
clínico,
como
muestra
de
aplicación
de
la
metodología desarrollada. Finalmente, e n el capítulo 8 se e x p o n e n las conclusiones de la tesis.
Capítulo 2
Planteamiento del problema
pág. 5
2. Planteamiento del problema La reconstrucción d e un objeto a partir de s u s proyecciones es un problema general c o n muchos ámbitos de aplicación. En este capítulo se realiza
una descripción
del problema
general
d e la
reconstrucción
tomográfica d e i m á g e n e s y su particularización a imágenes d e S P E C T . En primer lugar, s e estudia el proceso de formación d e las proyecciones, tanto
los
fenómenos
físicos
que
intervienen,
como
los
factores
d e g r a d a n t e s q u e afectan su calidad. Posteriormente, se introduce el formalismo m a t e m á t i c o asociado al problema general de la reconstrucción tomográfica y s u aplicación a S P E C T . T a m b i é n e s introducida la notación y
nomenclatura
q u e serán
utilizadas en adelante.
El
planteamiento
m a t e m á t i c o p u e d e encontrarse expuesto en m u c h o s trabajos, c o m o por
ejemplo [HER-80], [BAR-81], [KAK-84] o [VIE-88].
2.1. Descripción
del
proceso
de
obtención
de
las
proyecciones Para poder modelizar las proyecciones deben ser estudiados todos los m e c a n i s m o s q u e intervienen e n la formación d e las mismas, d e s d e la emisión hasta la detección c o n g a m m a c á m a r a s . Conocer el origen y efecto d e las distintas degradaciones es básico para poder realizar su corrección.
Un
sumario
detallado
de
todos
estos
efectos
puede
encontrarse, por ejemplo, en [JAS-80], [BAR-81], [TSU-94] o [ROS-95a]. A continuación, se describe el proceso d e obtención de las proyecciones: e m i s i ó n , transporte d e fotones en el interior d e l cuerpo y detección, así c o m o d e las degradaciones q u e se producen e n c a d a paso.
2.1.1. La emisión La emisión radiactiva es un proceso estocástico. La variable q u e indica el n ú m e r o d e emisiones producidas durante un intervalo predeterminado de tiempo, si dicho inten/alo de t i e m p o es mucho menor q u e el periodo de
Capítulo 2
Planteamiento del problema
pág. 6
semidesintegración del isótopo, es una variable aleatoria con distribución de Poisson P{Á,x),
cuya expresión matemática es:
(ec2-1)
P(X-x) = ^ : — ^
x!
d o n d e el parámetro Á d e la distribución es una constante positiva, propia del proceso, q u e corresponde al valor medio de desintegraciones emitidas si se hicieran infinitas observaciones [TUB-63]. Este parámetro d e p e n d e del periodo d e semidesintegración del isótopo radiactivo utilizado, de la concentración del m i s m o y de la duración de la adquisición, es decir, del intervalo de tiempo durante el cual se contabiliza la radiación.
Esta
distribución expresa la probabilidad de q u e se produzcan un determinado número x (variable entera) de emisiones radiactivas durante el tiempo fijado de adquisición.
La desviación típica de una distribución d e Poisson P[Á,X)
es igual a ^¡Á .
Esto significa q u e a u n q u e la dispersión crezca con el número medio d e fotones emitido, el error relativo (desviación típica respecto el valor medio) es porcentualmente m á s importante cuanto menor sea el número de fotones emitidos en el intervalo predeterminado de tiempo. Para reducir la presencia
de
ruido,
sería
conveniente
disponer
de
un
número
de
emisiones grande. Esto, sin embargo, no es generalmente factible, pues implicaria a u m e n t a r la dosis suministrada de radiofármaco o prolongar el tiempo de adquisición. La primera solución está descartada de a n t e m a n o por razones de protección radiológica del paciente. Alargar en e x c e s o el tiempo de las exploraciones, ya de por sí largo, d e b e
descartarse,
también, por d o s motivos principalmente. En primer lugar, por la relación exploraciones/día q u e se pueden realizar en cualquier centro de Medicina Nuclear. En s e g u n d o lugar, una adquisición muy prolongada a u m e n t a considerablemente el riesgo de obtener unas proyecciones inutilizables en
Capítulo 2
la
Planteamiento
reconstrucción
por
del problema
movimiento
pág. 7
involuntario
del
paciente.
C o n s e c u e n t e m e n t e , la presencia d e ruido es uno d e los factores m á s importantes de degradación d e las imágenes.
Otros aspecto a tener en cuenta es q u e la probabilidad de emisión no se mantiene constante e n el tiempo. En consecuencia, si no se corrige este efecto
sobre
las proyecciones,
el número
d e cuentas
obtenido
en
diferentes ángulos d e proyección no será equiparable y s e generan artefactos en la i m a g e n reconstruida. Hay dos circunstancias q u e alteran el ritmo p r o m e d i o decaimiento adquisición
d e emisión d e fotones. En primer
de la actividad es
prolongada
semidesintegración
radiactiva con el tiempo (más d e
del radiofármaco
30 minutos) es corta
lugar, está el (decay).
y
el
Si la
periodo
(para el isótopo
de más
habitual, el ^^"^Tc, e s d e 6 horas), la reconstrucción queda sensiblemente afectada. La corrección d e este efecto es muy sencilla. Basta multiplicar la proyección en c a d a ángulo por un factor de normalización q u e tenga en cuenta
el decaimiento
de
la actividad.
En segundo
lugar,
ha
de
considerarse la cinética de evolución del radiofármaco en el interior del órgano e n estudio [LIN-91]. Si el tiempo d e adquisición es prolongado, la concentración
de
radiofármaco
puede
modificarse
por
razones
metabólicas. La corrección d e este efecto es difícil, pues s u p o n e conocer de a n t e m a n o la biocinética del fármaco. No obstante, la variación de la emisión e n los estudios usuales es despreciable y raramente e s tenida en cuenta.
El último efecto degradante a tener en cuenta e n la emisión es la variación
de
la
posición
de
los
puntos
emisores
por
movimiento
involuntario del paciente [COO-92] [BOT-93]. Si esto ocurre, los datos obtenidos e n las proyecciones antes del movimiento del paciente no son compatibles
con los posteriores.
Dependiendo
d e la naturaleza
del
movimiento, aparecen artefactos d e mayor o m e n o r importancia en la imagen reconstruida.
Capítulo 2
Planteamiento del problema
pág. 8
2.1.2. La interacción de la radiación con el medio material El transporte d e los fotones por el interior del cuerpo, d e s d e el punto emisor hasta el exterior, q u e d a afectada por la presencia d e materia. Un porcentaje d e la radiación emitida interacciona c o n los á t o m o s q u e c o m p o n e n los tejidos. Está estimado [ROS-95a] q u e para los fotones d e 140
K e V (energía
c o n q u e emite
generalizado e n Medicina
el
isótopo
^^"^Tc,
de
uso m u y
Nuclear, con el que se han realizado las
diferentes adquisiciones experimentales para este trabajo) la interacción m á s c o m ú n c o n la materia es el efecto C o m p t o n . T a m b i é n se produce la absorción fotoeléctrica, y, en menor medida, dispersión elástica. Esta interacción radiación-materia tiene dos efectos sobre las proyecciones: •
Parte d e los fotones emitidos en una z o n a determinada del tejido, q u e debieran llegar a un punto determinado del detector, no lo h a c e n porque s o n absorbidos o desviados. Este f e n ó m e n o es d e n o m i n a d o atenuación
del
medio
[DUT-80].
La
atenuación
produce
una
disminución del n ú m e r o de fotones detectados en relación a los emitidos. Esta disminución d e p e n d e d e la distancia q u e recorre el fotón dentro del c u e r p o y de la naturaleza de los diferentes tejidos q u e atraviesa. Cuanto mayor sea el grueso d e tejido atravesado por un fotón, m á s probable e s s u interacción c o n los átomos q u e c o m p o n e n el tejido y, c o n s e c u e n t e m e n t e , su atenuación. A s u v e z , el tejido ó s e o tiene un coeficiente d e atenuación m a y o r que los tejidos blandos. Por tanto, la atenuación es no isótropa. A m o d o d e ejemplo, para poder estimar la importancia de este efecto, en un estudio cerebral c o n fotones d e 140 KeV, s e calcula q u e entre el 20 y el 2 5 % de los fotones emitidos s o n a t e n u a d o s [TSU-94]. Para el mismo tipo d e radiación, s e necesita unos 5 c m d e grosor d e un tejido blando, con índice d e atenuación similar al del agua, para atenuar la mitad d e fotones q u e lo atraviesan [JAS-80]. La no corrección d e la atenuación hace q u e las imágenes
tomográficas
reconstruidas
presenten
una
distribución
Capítulo 2
Planteamiento
del problema
pág. 9
incorrecta d e la actividad, siendo mayor en la periferia d e la m i s m a y m e n o r en el centro [GEL-88]. Un porcentaje d e los fotones q u e no deberían llegar a un determinado punto del detector lo hacen porque s o n desviados por el c h o q u e con un á t o m o o electrón (dispersión o scattering
Compton) (Figura 2-1).
Los fotones dispersados, al llegar al detector con una dirección diferente d e la correspondiente a s u punto emisor, aportan información falsa sobre la localización del radiofármaco. Los fotones q u e han sufrido algún tipo d e dispersión puede ser entre el 20 y el 4 0 % d e los fotones detectados [ROS-95a]. A s í , los fotones detectados s o n la s u m a de los q u e no h a n sufrido dispersión alguna (fotones primarios) y d e los q u e han sido dispersados. Analizando el espectro d e energía de los fotones q u e inciden en el detector, los fotones primarios p u e d e n se distinguidos d e los dispersados inelásticamente, pues estos últimos tienen
una energía
menor.
Por ejemplo,
un fotón
de
140 KeV
dispersado un ángulo de 52° sale con una energía de 126 K e V [ROS95a]. A tal efecto, las g a m m a c á m a r a s p o s e e n un m e c a n i s m o que permite calcular
la energía d e la radiación
incidente, a u n q u e su
resolución es baja. Por ejemplo, la resolución (FWHM)^ en la medición de la energía e n un cristal d e centelleo de y o d u r o sódico (Nal) para radiación d e 140 K e V es del orden del 1 3 . 2 % [LJU-90a]. Esto implica q u e , para captar una cantidad significativa d e los fotones primarios, la g a m m a c á m a r a ha d e contar todos los fotones que llegan con una energía dentro d e un intervalo (ventana de energía) con centro en la energía de emisión del isótopo y un margen d e desviación d e ± 1 3 . 2 % de e s e valor. Los fotones q u e han sufrido una dispersión inelástica con pérdida de energía superior al 1 3 . 2 % p u e d e n distinguirse d e los primarios y no s o n tenidos en cuenta en la detección. En cambio, los
^ FWHM: Full Width at Half Máximum. Intervalo de las x alrededor del máximo en el cual la función toma un valor mayor o igual que la mitad del valor máximo.
Capítulo 2
Planteamiento
pág. 10
del problema
fotones d i s p e r s a d o s q u e han mantenido su longitud d e onda, o la han variado m í n i m a m e n t e , pero han c a m b i a d o su dirección d e propagación son indistinguibles d e los primarios. Reducir el intervalo d e energía, significa eliminar parte del efecto d e la dispersión pero a costa d e reducir el n ú m e r o d e fotones primarios detectados c o n el consiguiente a u m e n t o d e ruido. La dispersión t a m b i é n d e p e n d e d e la profundidad y de la naturaleza d e l tejido [ROS-90], así como d e la geometría d e l cuerpo
[FRE-91].
La
no
corrección
de
la
dispersión
produce
reconstrucciones borrosas [TSU-94].
paciente
detector
Figura 2-1. Parte de los fotones emitidos desde A que deberían llegar a C no lo hacen por ser absorbidos por el medio o desviados de su trayectoria (atenuación). Por otra parte, fotones emitidos en B pueden llegar a C si son desviados por un choque con un átomo (dispersión).
2.1.3. La detección Para considerar la d e g r a d a c i ó n que se produce en la detección hay q u e tener e n cuenta la estructura y funcionamiento d e una g a m m a c á m a r a . Una descripción detallada puede encontrarse en
[DUT-80], [JAS-80],
[BAR-81], [ZAR-92] o [CRO-95]. La g a m m a c á m a r a , o cámara de Anger, está c o m p u e s t a por un cabezal giratorio e n el que están los detectores d e radiación y un ordenador q u e procesa los datos obtenidos. El m e c a n i s m o
Planteamiento
Capítulo 2
pág. 11
del problema
d e detección está compuesto por un cristal de centelleo, generalmente Nal con impurezas de Tallo. El fotón incidente e n el cristal v a cediendo parte de su energía a los electrones del cristal c o n que interacciona. La energía q u e a b s o r b e n los electrones hace q u e suban a un nivel de energía superior. A l desexcitarse, al cabo d e unas millonésimas
de
s e g u n d o , estos electrones emiten un fotón luminoso. De esta m a n e r a , el cristal d e centelleo e s un transductor que convierte un fotón y en un conjunto d e fotones d e luz visible. Cuanto m á s energético sea el fotón incidente,
más
intensa
será
la
señal
luminosa
de
respuesta
(más
electrones serán excitados y, por tanto, más fotones serán emitidos). Junto
al
cristal
de
centelleo
está
acoplada
una
batería
de
fotomultiplicadores conectados a un procesador. A partir d e la señal recibida por los diferentes fotomultiplicadores se puede deducir la posición del impacto del f o t ó n . Así se contabiliza j u n t a m e n t e la cantidad d e fotones detectados (número d e cuentas), su energía (a través de la intensidad de la señal luminosa) y su posición e n el plano de proyección. En este proceso de detección se produce una pequeña adición de ruido sobre la señal (ruido electrónico) debido a: •
La eficiencia del detector. No todos los fotones que inciden e n el detector son detectados. A l g u n o s no llegan a producir señal porque no interaccionan con el cristal de centelleo.
Colimador
Fotón incidente
III Cristal de centelleo Fotomultiplicadores Colector de la señal
b)
Figura 2.2. a) Estructura de una gammacámara. b) Disposición transversal de los fotomultiplicadores
Capítulo 2
•
Planteamiento
del problema
pág. 12
D e s p u é s d e realizar una detección, el detector tiene un tiempo muerto e n el q u e es inoperante (de 1 a 5|as, según la g a m m a c á m a r a
[TSU-
94]). Cualquier fotón q u e incidiera en e s e intervalo no sería detectado. •
El localizador del impacto sobre el cristal de centelleo tiene
una
resolución limitada d e unos 3 m m [TSU-94].
Otro factor f u n d a m e n t a l a tener en cuenta es la digitalización del plano de proyección. El plano se cuadricula en unidades de superficie (píxeles). El t a m a ñ o de los mismos p u e d e ser predefinido por el técnico que realiza la adquisición. Los fotones q u e inciden en cada uno d e los píxeles son acumulados en el correspondiente contador, produciéndose cierta pérdida de información al no tenerse en cuenta la posición exacta del impacto. Del t a m a ñ o d e píxel escogido en las proyecciones d e p e n d e la resolución de la reconstrucción (objeto de mínimo t a m a ñ o q u e puede ser distinguido en la imagen). No obstante, no puede definirse arbitrariamente pequeño, y a que, al reducirlo, también se reduce el n ú m e r o de cuentas por píxel, y, por tanto, se a u m e n t a el ruido. Se ha de buscar una solución de c o m p r o m i s o entre ruido y resolución: para evitar que el ruido sea m u y importante, con la consiguiente pérdida de fiabilidad de los datos, se ha de
perder
resolución y contabilizar conjuntamente los fotones provenientes d e una región amplia de tejido, es decir usar t a m a ñ o s de píxel grandes.
Por otra parte, para determinar la posición del emisor hay que conocer, a d e m á s d e la posición del impacto, la dirección de incidencia del fotón. C o m o la emisión radiactiva es isótropa en el espacio y la radiación g a m m a no p u e d e ser focalizada con ningún tipo de lente, para poder garantizar
que
los fotones
incidentes
lleguen
con
una
determinada
dirección se utilizan colimadores. Los colimadores utilizados constan de una lámina de plomo de unos 4 ó 5 c m de espesor en la que s e han realizado una serie d e agujeros transversales por los q u e los fotones llegan ai cristal de centelleo. De esta manera, solamente la radiación q u e llega al cabezal en la dirección que definen estos agujeros es capaz de
Capítulo 2
Planteamiento
pág. 13
del problema
atravesar el colimador y ser detectada (en realidad, hay un pequeño porcentaje d e la radiación q u e atraviesa las paredes de los colimadores, fenómeno
conocido
como
penetración
septal,
y, a
su vez, existe
dispersión d e los fotones con los átomos q u e c o m p o n e n el colimador [MAT-57] [BAR-81] [DEV-90]). Esto permite asegurar q u e los fotones detectados en un punto del cabezal han incidido en la dirección q u e define el agujero del colimador. Si los agujeros del colimador fueran de calibre infinitesimal, para cada ángulo, la proyección de un punto emisor sería un punto, aquel pixel d e las proyecciones que 'apunta' al pixel emisor d e la i m a g e n . C o m o el calibre d e los agujeros del colimador es finito,
hay fotones
q u e inciden
e n el detector
con un ángulo
no
necesariamente perpendicular al m i s m o . El resultado de esto, junto con la penetración septal y dispersión d e fotones dentro d e l colimador, es q u e la imagen d e un punto emisor es una m a n c h a o PSF (Point Function)
Spread
(Figura 2-3). Cuanto m á s alejado esté el punto emisor del
detector tanto m a y o r será la P S F [MET-80]. La presencia del colimador, por otra parte totalmente necesaria, produce otro efecto: reduce en un factor
del orden
d e 1000 los fotones
detectados
[KNO-83]
consiguiente a u m e n t o del ruido.
posición del detector
colimador
PSF
punto emisor Número de cuentas
letectores
Figura 2-3. Esquema del origen de la PSF
c o n el
Capítulo 2
Planteamiento
del problema
pág. 14
S e g ú n los agujeros del colimador sean paralelos o estén focalizados sobre una recta o sobre un punto, se hablará de colimadores paralelos, fan-beam
o cone-beam
respectivamente. En este trabajo, s e
utilizan
s o l a m e n t e colimadores de agujeros paralelos.
Existe una clara relación entre el colimador, el t a m a ñ o de píxel e n las proyecciones, el ruido y la resolución. A fin de reducir el nivel de ruido debe
obtenerse
el
mayor
número
de
cuentas
posibles
proyecciones. Para ello, puede trabajarse con colimadores de
en
las
mayor
eficiencia y utilizar un t a m a ñ o d e píxel m a y o r en las proyecciones. A m b a s actuaciones p r o d u c e n una pérdida de resolución (pérdida de detalle) e n las
proyecciones
que
conlleva
una
pérdida
de
resolución
en
las
reconstrucciones. Determinar la relación óptima entre ruido y resolución d e p e n d e de la finalidad d e la imagen obtenida [TSU-84].
Finalmente, un mal reglaje de las g a m m a c á m a r a s p u e d e producir, a su vez, d e g r a d a c i ó n de las imágenes. Los efectos m á s c o m u n e s son el m a p a de sensibilidad no uniforme de los detectores, falta de linealidad y mal ajuste del centro d e rotación. No obstante, estos dos efectos son fácilmente detectables e n un control de calidad y su corrección es sencilla
[JAS-80] [ENG-86] [TSU-94].
2.2. Planteamiento
matemático
de
la
reconstrucción
tomográfica El proceso de obtención de las proyecciones de un objeto p u e d e ser representado c o m o un sistema lineal, en el que la entrada es el objeto y la salida es el conjunto de proyecciones. La linealidad del sistema es fácil de asumir; un objeto q u e t e n g a el doble de actividad g e n e r a proyecciones que t a m b i é n tiene el d o b l e de cuentas por píxel (en realidad, esto será cierto mientras el t i e m p o muerto del detector sea despreciable frente el tiempo
medio
entre
dos
impactos).
Por tanto,
matemáticamente
la
Capítulo 2
Planteamiento
del problema
pág. 15
reconstrucción tomográfica es un problema d e inversión d e sistemas lineales. El m a r c o matemático para definir la reconstrucción tomográfica es el del análisis funcional. La imagen a reconstruir es una función q u e se ha de determinar a partir de otra función, las proyecciones, y un operador q u e las relaciona. A u n q u e el formalismo general d e la reconstrucción tomográfica se desarrollará e n dicho marco, t a m b i é n se incluyen otras visiones del problema, como la algebraica, q u e aparece al implementar n u m é r i c a m e n t e el problema, y la estadística.
Las reconstrucciones tomográficas q u e se llevan a cabo en S P E C T s o n de
imágenes
tridimensionales
a
partir
de
sus
proyecciones
bidimensionales. Sin embargo, si cada sección (corte tomográfico) de la imagen
tridimensional
se
puede
tratar
independientemente
de
las
restantes, el problema queda reducido a varias reconstrucciones (tantas c o m o secciones se t o m e n de la imagen) bidimensionales a partir de proyecciones unidimensionales. En la realidad suponer las diferentes secciones independientes es una aproximación, ya q u e , tanto la P S F c o m o la dispersión d e fotones interrelacionan a varias secciones entre sí. No obstante, esta aproximación e s utilizada
habitualmente.
En caso
contrario, la reconstrucción es d e m u c h a mayor complejidad y requiere un esfuerzo d e c o m p u t a c i ó n notablemente superior. En este trabajo todas las reconstrucciones realizadas son d e imágenes bidimensionales. Por tanto, t o d o el formalismo y las siguientes definiciones están referidos a la reconstrucción d e imágenes bidimensionales a partir de sus proyecciones unidimensionales.
S u generalización
a
reconstrucción
de
imágenes
tridimensionales a partir d e proyecciones bidimensionales e s inmediata [KAK-84] [HER-80]. El proceso de reconstrucción es abordado
tanto
d e s d e el punto d e vista analítico como numérico. Posteriormente, se particulariza el formalismo a la reconstrucción tomográfica en S P E C T .
Capítulo 2
Planteamiento del problema
pág. 16
2.2.1. Descripción analítica Se define u n a i m a g e n c o m o una función bidimensional q(x,y) positiva, finita y s o l a m e n t e distinta d e cero e n u n a región acotada del plano (soporte acotado) [GON-77]. En el caso d e la Medicina Nuclear, la función q(x,y) representa la concentración d e radiofármaco e n cada punto d e l corte t o m o g r á f i c o e n cuestión.
Esta función será s i e m p r e d e cuadrado integrable. Por tanto, p u e d e ser interpretada c o m o un e l e m e n t o del espacio ¿2 (espacio d e Hilbert - lineal, n o r m a d o - d e las funciones d e cuadrado integrable) [LOU-83] [MAD-91].
Se define c o m o rayo cualquier recta q u e atraviese el soporte d e la imagen. La ecuación d e l rayo puede ser determinada e n función d e d o s variables. U s u a l m e n t e , s e utilizan el ángulo 0 que f o r m a la perpendicular al rayo c o n el e j e d e las abscisas y la distancia t al origen d e c o o r d e n a d a s (Figura 2-4). Entonces, la ecuación del rayo e s :
X-COS0 + y-sind
=t
Se d e n o m i n a rayo-suma a la integral d e q(x,y) a lo largo d e un rayo. S u expresión para el rayo definido por 0 y por t e s :
Po{t)= lj{x,y)-ds
(ec2-2)
d o n d e ds e s el diferencial d e camino. Sustituyendo ds por su expresión sobre el rayo:
Pg{t)= r r q{x,y)-S{x-cos0
+ y-sine-t)-dx-dy
(ec 2-3)
Capítulo 2
Planteamiento
Dáa. 17
del problema
soporte función
de
la
Figura 2.4 Esquema del rayo que atraviesa el soporte de una función
La proyección Pg{t) e n la dirección definida por 0
es la función d e /
obtenida al considerar el conjunto de todos los rayos-suma en esa dirección (Figura 2-5). A cada valor d e t se le hace corresponder el rayosuma
Pg{t).
La proyección d e q(x,y) es la función bidimensional
P{0,t)
obtenida al
considerar el conjunto de las proyecciones Pg[t) para cada ángulo 0. El funcional q u e asocia la función c o n su proyección, entendida en estos términos, s e d e n o m i n a la transformada de R a d o n d e la función [KAK-84]. S e escribe:
q{x,y)^^P{Qj) (ec 2-4)
P{Q,t) =
^(q{x,y))
Planteamiento
Capítulo 2
pág. 18
del problema
Figura 2-5. Generación de las proyecciones en una dirección dada.
2.2.2. Discretización Las expresiones del apartado anterior están f o r m u l a d a s para
objetos
analíticos. D a d o q u e e n Medicina Nuclear se trabaja con i m á g e n e s y proyecciones digitalizadas, a continuación se escribirán para
objetos
discretizados. E n primer lugar, se realiza un proceso d e discretización (sampling)
d e la imagen y de las proyecciones y, posteriormente, se
transcribe la expresión analítica del proceso de proyección (ec 2-3) e n su análoga discreta. El proceso d e discretización, q u e es necesario para poder realizar los cálculos en la reconstrucción, produce efectos sobre la calidad d e la reconstrucción. Un tutorial s o b r e discretización y pérdida d e información p u e d e encontrarse e n [JER-77]. Otros m u c h o s textos tratan con m a y o r o m e n o r reconstrucción
profundidad c ó m o afecta la discretización
tomográfica,
[VIE-88] [FAR-90] [ROS-95a].
por ejemplo
[SMI-73]
[BAR-81]
a la
[LOU-83]
Capítulo 2
Planteamiento
Discretización
del problema
pág. 19
de la imagen
Para digitalizar u n a imagen, su soporte se acota mediante un rectángulo. Este rectángulo e s dividido en unidades básicas d e superficie, cada una d e las cuales recibe el nombre d e píxel (del inglés picture volume
element-
e n el caso tridimensional).
element;
vóxel -
Generalmente, tanto el
soporte d e la imagen c o m o los píxeles son cuadrados. A cada píxel se le asigna el valor medio d e la imagen e n su interior. Así, una imagen digitalizada es una matriz de números reales. Por c o m o d i d a d , en v e z d e trabajar
con matrices y tener q u e manejar d o s índices, se
utilizan
vectores, q u e precisan un solo índice. El vector imagen s e construye mediante la sucesión d e filas de la matriz imagen. En lo sucesivo, cuando se hable d e una imagen, se sobreentenderá q u e está digitalizada sobre un matriz cuadrada y expresada mediante un vector. La variable q u e nos da el t a m a ñ o del píxel (longitud d e sus lados) será tp. A l número de filas y c o l u m n a s d e la i m a g e n digitalizada s e le representa con la variable M. La imagen q u e d a , entonces, a l m a c e n a d a e n un vector d e MM c o m p o n e n t e s ( M f i l a s d e M c o m p o n e n t e s cada una). A l producirse la digitalización d e la imagen se ha d e tener en cuenta q u e se está acotando sus c o m p o n e n t e s espectrales. Una función discretizada uniformemente está limitada en el espacio
d e las frecuencias.
En efecto,
no habrá
componentes
de
frecuencia mayor a la frecuencia d e Nyquist, correspondiente a una sinusoide d e 2 píxeles d e periodo [PRA-76] [GON-77]:
1 w
Por
tanto,
hay q u e asegurar
frecuencias superiores a fenómeno
de aliasing
wmax,
=
q u e en la i m a g e n
original
no
haya
ya q u e , en caso contrario, aparecería el
[GON-77].
En la práctica, el ruido sobre las
Capítulo 2
Planteamiento
del problema
pág. 20
proyecciones limita la cantidad d e píxeles utilizados (ver 2.1.1). El valor usual d e M e s d e 6 4 o 128, dependiendo del tipo de imagen^.
Discretización de las proyecciones La primera discretización a q u e es sometida la proyección e s en la variable angular. Las proyecciones sólo s o n adquiridas sobre un n ú m e r o finito de ángulos. La variable ángulos
d e proyección
proyección
necesarios
indicará, a partir de ahora, el n ú m e r o d e
utilizados. para
Sobre
reconstruir
el número
de ángulos
de
la imagen correctamente, e n
[KAK-84] s e recomienda el uso d e tantas proyecciones c o m o bins haya por proyección para evitar el efecto del aliasing.
En el mismo sentido, en
un reciente trabajo [CAO-96] s e realiza un estudio d e los artefactos q u e aparecen e n la imagen d e b i d o s al f e n ó m e n o del aliasing
en función del
n ú m e r o d e ángulos d e proyección utilizados. Se concluye q u e utilizar un n ú m e r o impar d e á n g u l o s mejora los resultados, sobre todo si no s e considera el efecto d e la atenuación. N o obstante, el número usual d e proyecciones utilizado es d e 6 0 o 120 por s e r divisores d e 360.
La proyección para c a d a uno d e los ángulos está, a s u vez, discretizada e n una sucesión de valores reales. La recta sobre la q u e se proyecta s e divide en una serie de s e g m e n t o s de igual longitud llamados bins (píxeles en
caso
de
proyecciones
bidimensionales).
El valor
que se
hace
corresponder a cada bin d e la proyección e s la integral d e la proyección entre s u s extremos. Para el biny-ésimo:
PÁj)=tPMdt
(ec2-5)
^ Se emplea una potencia de 2 para poder utilizar el algoritmo FFT optimizado (Fast Fourier Transform) en la transformada de Fourier de la imagen.
Planteamiento
Capítulo 2
del problema
pág. 21
Si se c o m b i n a la fórmula anterior c o n la expresión integral d e Pg(t) (ec 22) resulta:
Pe{j)-
l^^q(x,y)-ds-dt
(ec2-6)
Así, c u a n d o la proyección esté digitalizada, se entenderá por rayo, no la recta q u e cruza el soporte de la i m a g e n , sino la banda de anchura tb=ti+i-t¡ q u e lo cruza. La variable t es sustituida por un índice q u e indica el orden de los bins. El valor q u e se le asigna al bin y-ésimo d e la proyección es la integral d e la función q(x,y) sobre la superficie q u e delimitan los dos e x t r e m o s del bin (Figura 2-6). En lo sucesivo, el t a m a ñ o del bin será designado
con
la
variable
tb.
En
reconstrucción
tomográfica
con
colimadores paralelos, s e hace coincidir el n ú m e r o de bins de cada proyección c o n el n ú m e r o de las filas y columnas d e la imagen expresada matricialmente, así c o m o sus respectivos t a m a ñ o s {tb = tp). Por lo tanto, la proyección entendida c o m o función bidimensional de las variables 6 y t q u e d a reducida a TV filas d e M valores reales, e s decir, un vector d e MxTV c o m p o n e n t e s . La imagen formada por las
filas d e las M c o m p o n e n t e s
de las proyecciones dispuestas matricialmente se la denomina sinograma.
Figura 2-6. La contribución a un bin de la proyección
Capítulo 2
Planteamiento del problema
pág. 22
Otro factor a tener e n cuenta es q u e las proyecciones e n los diferentes ángulos d e proyección s e t o m a n de igual longitud ( m i s m o n ú m e r o d e bins de igual t a m a ñ o ) . P u e d e ocurrir a v e c e s que la proyección para
un
determinado ángulo s o b r e p a s e los límites del segmento donde se realiza la digitalización. En e s e caso, la truncación de las proyecciones p u e d e comportar artefactos e n la reconstrucción [MUL-96].
Discretización de las ecuaciones C u a n d o s e discretiza espacialmente la imagen, la integral doble de la expresión (ec 2-6) se p u e d e d e s c o m p o n e r en una serie de integrales sobre los píxeles q u e recorre el rayo. C o m o en cada píxel la función a integrar e s uniforme, el valor de la integral sobre cada píxel es el valor de la función por la integral d e la superficie intersección del píxel y el rayo, es decir, el área de la intersección rayo-píxel.
Pe{j)='
[
*rayo
q{x,y)-ds-dt
MxM
= Zíf.
, . n
.Q{x,y)-dx-dy
(ec 2-7)
MxM JJptxei-t^
^fayo-j
d o n d e q) es el valor p r o m e d i o de la función en el píxel. La integral doble de este último término es un valor numérico A}, que d e p e n d e del píxel y del bin considerados. C o r r e s p o n d e al área de la intersección del pixel con el rayo.
Usualmente, este factor se normaliza dividiéndolo por el área del píxel. En este caso, e n vez de multiplicar este factor por el valor promedio de la función e n el píxel, se multiplica por la actividad total del m i s m o . La expresión (ec 2-7) se reescribe:
Planteamiento
Capítulo 2
MxM
Dáa. 23
del problema
ff
dx-dv MxM t
;=1
„
/=1
De esta m a n e r a AJÍ está normalizado a la unidad. En todo caso, contribución
del
píxel
/-ésimo
de
la imagen
al bin j - é s i m o
es la de
proyecciones (Figura 2-7). Escribiendo el vector proyección c o m o
las la
sucesión de vectores proyección para los distintos ángulos y utilizando una notación vectorial [HER-80], esta última expresión queda:
(ec 2-8)
o p =A q
La matriz A f o r m a d a por el conjunto de factores AJÍ es conocida c o m o matriz d e transición o matriz de pesos. Contiene la información necesaria para la obtención d e las proyecciones a partir d e la imagen.
N
Figura 2-7. Un píxel no contribuye a un bin de la proyección con toda su actividad. Su contribución es proporcional al área de píxel inscrita dentro del rayo. Por ejemplo, los dos píxeles sombreados aportan una proporción de actividad diferente al rayo dibujado.
Capítulo 2
Planteamiento
del problema
pág. 24
2.2.3. Particularización al caso de SPECT Una v e z presentado
el formalismo
general de las
reconstrucciones
tomográficas, a continuación se desarrolla este mismo formalismo para el caso
particular
degradaciones
de
SPECT,
hace
que
las
pues
la
presencia
ecuaciones
deban
de ser
las
distintas
ligeramente
modificadas. S e han d e tener en cuenta la atenuación, la P S F , la dispersión y el ruido, y c ó m o se incorporan estas degradaciones al formalismo.
C u a n d o s e tiene en cuenta la atenuación, el valor de la proyección para un ángulo d a d o no p u e d e obtenerse mediante la integral de la (ec 2-3). Debe introducirse en la misma una función de peso q u e incorpore el efecto d e la pérdida d e valor e n la proyección debido a la atenuación para cada rayo. Puesto q u e la atenuación no es isótropa, esta función d e peso dependerá, a s u vez, d e las variables d e la proyección:
Pg (í) =
q{x,y)- w{x, y,t,0) • S[x • cosd + y • sinO - {)• dx • dy
Por ejemplo, si la atenuación es uniforme en todo el objeto, esta última expresión se escribe [LOU-83]:
Pe{f)=
q{x,y)-Qxp{--ju-lg)-S{x-cosd
+y-sin0-t)-dx-dy
(ec 2-9)
donde ju e s el coeficiente d e atenuación y h es la distancia del punto (x,y) al borde del objeto en la dirección 0. Esta expresión suele d e n o m i n a r s e transformada d e R a d o n atenuada. Si el mapa d e atenuación no es uniforme,
la exponencial
d e esta expresión debe
sustituirse
por la
exponencial inversa d e la integral de línea desde cada punto al extremo del objeto en la dirección 0 del coeficiente de atenuación punto a punto ju{x,y) ( m a p a d e atenuación del objeto) por el camino recorrido por el rayo.
Capítulo 2
-^^(0=
Planteamiento
(í{x,y)-&^p
-
pág 25
del problema
'ju{x,y)-dl
•S{x-cos0
+
y-sin6-t)-dx-dy
d o n d e la integral dentro de la exponencial se extiende desde el punto {x,y) hasta el punto correspondiente del detector. En la práctica, la atenuación producida por el aire entre el paciente y el detector es despreciada y sólo se calcula la atenuación que se produce dentro del paciente.
Si se incluye la P S F , sobre un punto de la proyección, a d e m á s de la contribución
del
rayo
perpendicular,
otros
rayos
provenientes
de
direcciones ligeramente diferentes pueden contribuir a un bin d a d o (Figura 2-8). En este caso, la función delta q u e marca la dirección del rayo en (ec 2-9) d e b e ser suprimida y la función de peso d e b e incluir la contribución d e cada punto de la imagen sobre cada punto de la proyección.
Figura 2-8. Esquema de los distintos rayos que contribuyen a un punto de la proyección cuando existe PSF.
Capítulo 2
Planteamiento del problema
pág. 26
De hecho, t o d a la información del proceso físico d e proyección puede incluirse s i m u l t á n e a m e n t e en esta función d e peso, la atenuación, la P S F y la dispersión.
(ec2-10)
Po{t)= [^[j{x,y)-w(x,y,t,0)-dx-dy
El ruido no p u e d e incluirse d e una m a n e r a determinista e n esta función d e pesos d a d a s u naturaleza aleatoria. S e incluye c o m o una función d e residuos q u e s e s u m a a la proyección del objeto. La expresión d e la proyección c o n todas las degradaciones presentes e n S P E C T q u e d a :
^e(0=
[^[j{x,y)-w{x,y,t,e)-dx-dy
+ n,{t)
(ec2-11)
donde ng(t) representa el ruido. Prescindiendo del ruido o incorporándolo en las proyecciones, la reconstrucción tomográfica s e plantea c o m o la resolución d e una ecuación integral d e Fredholm d e primera
especie
[DEM-89]. A l discretizar esta última expresión, al igual q u e se hacía en el apartado anterior, s e obtiene:
p = Aq
(ec2-12)
+n
d o n d e la matriz A está f o r m a d a por los elementos:
w(x,y,e,,tf)-dx-dy A ^ . = ^
^
con j = k-M
+h
(ec2-13)
^p
y n representa el ruido (« e s a su vez un vector d e MxN c o m p o n e n t e s , cada una d e ellas e s el ruido del correspondiente bin d e la proyección).
Capítulo 2
Planteamiento
del problema
pág. 27
En la (ec 2-12), s e v e que la proyección se c o m p o n e de una parte de señal (determinista) y otra de ruido (aleatoria) q u e , según el caso, puede ser m á s o m e n o s significativa.
En r e s u m e n , se ha obtenido una formulación de la obtención de las proyecciones a partir de la imagen. Una imagen es un e l e m e n t o de un espacio vectorial d e dimensión MxM.
La proyección se obtiene mediante
el producto de la imagen por la matriz de transición y la adición de ruido. Proyecciones obtenidas de una m i s m a imagen e n las mismas condiciones d e adquisición difieren sólo en el vector ruido.
2.2.4. Interpretación estadística del proceso de formación de la proyección. S e ha m e n c i o n a d o con anterioridad que la emisión radiactiva es un proceso aleatorio. La cantidad de fotones emitidos por un pixel de la imagen está sujeto a fluctuaciones en el tiempo. Por tanto, diferentes muestras de actividad de un m i s m o pixel, darían lugar a
diferentes
resultados. En este sentido, la variable que describe la emisión de un pixel de la imagen durante un periodo predeterminado de t i e m p o puede ser considerada c o m o una variable aleatoria de Poisson. Una imagen es un vector de variables aleatorias de Poisson independientes entre ellas. La proyección e s una realización concreta de una función lineal de dichas variables. El e l e m e n t o AJÍ de la matriz de transición (normalizado a 1) p u e d e ser interpretado c o m o la probabilidad d e q u e un fotón q u e es emitido
por el pixel i de
la imagen
incida sobre
un bin J de
las
proyecciones (probabilidad condicional de detectar un fotón e n el bin J si ha sido emitido por el pixel i) [SHE-82]. A su vez, el número d e fotones detectados por un bin también está sujeto a fluctuaciones e n el tiempo, b á s i c a m e n t e d e b i d o a la variabilidad e n el número de fotones emitidos por unidad de tiempo en cada pixel de la imagen. Por consiguiente, el valor de cada bin puede ser considerado c o m o una variable aleatoria. Tanto en el
Capítulo 2
Planteamiento
del problema
pág. 28
caso de q u e no haya degradaciones [SHE-82] [LLA-84] c o m o en el caso de q u e e s t é n presentes [HER-90] cada una d e estas variables p u e d e ser descrita mediante una variable aleatoria d e Poisson. Dos adquisiciones diferentes d e proyecciones corresponden a dos muestras distintas del mismo conjunto d e variables aleatorias.
2.2.5. Problemas matemáticos de la reconstrucción. La reconstrucción tomográfica es un problema que, debido a s u propia estructura, a la presencia d e ruido y a los procesos d e digitalización q u e se llevan a cabo, está m a l planteado (ill-posed)
[DAV-83] [LOU-83] [JOY-
84] [NAT-88] [DEM-89]. Este hecho tiene tres posibles consecuencias: •
El
problema
puede
no
tener
solución,
debido
a
la
posible
incongruencia e n los datos a causa del ruido. Este hecho e s m u y relevante e n sistemas sobredeterminados. En este, en [TIH-63] s e e x p o n e q u e una ecuación integral c o m o (ec 2-10) no tiene solución si el núcleo d e la integral tiene un grado mayor d e suavidad q u e la proyección. Esto s u c e d e cuando este núcleo incluye la P S F y las proyecciones tienen una importante presencia d e ruido. Es lógico pensar q u e este resultado es extrapolable al sistema discretizado. •
Puede
q u e exista
solución, pero q u e no s e a única. En objetos
analíticos, el hecho d e disponer solamente de un conjunto finito d e proyecciones implica q u e la reconstrucción no s e a única [SMI-73]. E n [LOU-81] [LOU-83] s e d e m u e s t r a q u e , e n este caso, existen i m á g e n e s "invisibles" e n las direcciones consideradas, es decir, d e proyección nula e n esas direcciones. Si estas Imágenes s e s u m a n a la imagen reconstruida se obtiene una nueva solución. Estas imágenes, llamadas imágenes
fantasma
(ghosf),
constituyen
el núcleo
o
/cerne/ del
operador proyección y tienen una estructura m u y oscilante [LOU-81]. No está evaluado c o m o afectan las diferentes degradaciones a este resultado. A l digitalizar las imágenes este f e n ó m e n o prevalecerá o no en función del n ú m e r o d e píxeles d e la imagen y del n ú m e r o d e
Capítulo 2
Planteamiento del problema
pág. 29
ángulos d e proyección considerados, es decir, si el sistema lineal resultante está suficientemente determinado o no. •
Puede que tenga
solución pero la imagen reconstruida
no varíe
c o n t i n u a m e n t e c o n las proyecciones. Este último hecho significa q u e una p e q u e ñ a variación en las proyecciones p u e d e dar una variación m u y g r a n d e e n la imagen (amplificación del ruido) [LLA-82] [VIE-88]. C o n todo, cada método d e reconstrucción tiene sus técnicas de regularización^ q u e nos garantizan q u e la solución tiene cierto grado d e proximidad a la solución ideal. La regularización t a m b i é n puede solventar el problema de las imágenes ghost [LOU-81].
Por
otra
parte,
circunstancias
del
la
corrección
mal
que
planteamiento
puede del
hacerse
problema
en no
algunas garantiza
estabilidad en la reconstrucción. El buen planteamiento es condición necesaria pero no suficiente para q u e pequeñas variaciones en los datos no d e n lugar a u n a solución m u y alejada de la ideal. En el planteamiento numérico, se ha d e determinar el buen o m a l condicionamiento del problema [BER-86]. La reconstrucción tomográfica es un problema mal condicionado"* [DEM-89]. El hecho d e ser un problema mal planteado y/o mal condicionado e s un factor determinante a la hora d e decidir las estrategias de resolución.
^ Se entiende por regularización de un problema mal planteado, al paso a otro problema equivalente, dependiente de un parámetro a positivo, que nos garantice que puntos próximos a ios datos están relacionados con puntos próximos a las soluciones. El grado de proximidad lo da el parámetro a. Para a = O, el problema se reduce al original [TIK74]. En ios capítulos siguientes se profundiza sobre este punto. ^ Si un problema numérico está mal condicionado, pequeñas perturbaciones en los datos producen una gran perturbación en la solución. Es decir, al invertir el problema el ruido
existente en los datos se amplifica [BAK-80] [BER-86] [GOL-89].
Capítulo 2
Planteamiento
2.2.6. Reconstrucción:
del problema
métodos
pág. 30
analíticos
y
métodos
iterativos. En los apartados anteriores, se ha realizado un desarrollo d e la obtención de las proyecciones a partir d e la imagen desde diferentes puntos d e vista. En primer lugar, s e ha considerado la imagen c o m o una función analítica y la proyección c o m o una transformación funcional de la m i s m a . Posteriormente se ha introducido la discretización d e la imagen y el sistema d e e c u a c i o n e s d e la proyección, q u e ha sido interpretado d e d o s maneras diferentes, determinista y estadísticamente.
En este punto se d e s c r i b e n , a grandes rasgos, las estrategias existentes de inversión del problema, es decir, encontrar la imagen a partir d e unas proyecciones conocidas (reconstrucción tomográfica) y las ventajas e inconvenientes d e unas y otras.
Básicamente, se clasifican los métodos e n tres grupos: •
La imagen e s considerada c o m o una función analítica y la proyección pura
(sin función
degradación).
En
de
ponderación,
tal
caso,
la
es decir,
sin incluir
reconstrucción
ninguna
tomográfica
es
equivalente a invertir la expresión de la (ec 2-4):
Este problema, resoluble mediante el teorema d e la sección central según s e verá e n el capítulo 4 , tiene solución exacta si se tiene la proyección para t o d o s los valores d e la variable angular
[SMI-73]
[KAK-84]. No obstante, en la práctica se dispone sólo de un n ú m e r o finito
de
proyecciones,
afectadas
además
por
diferentes
degradaciones, c o n lo q u e la aplicación del t e o r e m a de la sección central e n la reconstrucción d e este tipo de proyecciones d a una solución a p r o x i m a d a .
Las degradaciones
de las proyecciones s o n
Capítulo 2
Planteamiento
del problema
pág. 31
tratadas y a s e a previamente a la reconstrucción, mediante el filtrado de las proyecciones, o posteriormente, c o n algoritmos de corrección sobre la i m a g e n reconstruida. Este m é t o d o d e reconstrucción es d e n o m i n a d o retroproyección filtrada ( F B P - Filtered Filtered
Back
Projection-
si
se
Back Projection utilizan
o IFBP-
algoritmos
Iterative
iterativos
de
corrección de las degradaciones) y se estudia e n el capítulo 4 . •
Interpretación determinista. S e invierte el sistema lineal (ec 2-12) m e d i a n t e algún m é t o d o algebraico iterativo. Estos métodos reciben el n o m b r e genérico de métodos d e reconstrucción algebraicos (ARTAlgebraic
•
Reconstruction
Tecfiniques)
y se estudian en el capítulo 5.
Interpretación estadística. S e resuelve mediante un método estadístico el sistema d e ecuaciones (ec 2-12), siendo p el conjunto d e datos, A una matriz c o n o c i d a y q las incógnitas. El m é t o d o habitual es el de m á x i m a verosimilitud d e la imagen (imagen q u e tiene las mayor probabilidad d e generar las proyecciones obtenidas). S e hablará de m é t o d o s d e reconstrucción estadísticos ( M L E - Máximum Estimator).
Likelihood
S o n el objeto d e estudio del capítulo 6.
Los m é t o d o s d e reconstrucción IFBP, algebraicos y estadísticos s o n m é t o d o s iterativos y precisan el conocimiento d e la matriz d e transición para el cálculo d e las proyecciones, al contrario q u e FBP. Esto le d a a la retroproyección filtrada una gran ventaja frente ellos: una e n o r m e facilidad y rapidez d e computación. En la actualidad, es el procedimiento estándar e n los e q u i p o s d e g a m m a g r a f í a comerciales. No obstante, este m é t o d o no permite
incorporar
la
corrección
de
las
degradaciones
en
la
reconstrucción, c o n lo que las imágenes resultantes son d e baja calidad. En c a m b i o , el uso d e los métodos iterativos permite incluir en la matriz de pesos varios de los efectos degradantes d e la imagen [GUL-85] [TSU-94]. Así, al realizar las sucesivas iteraciones, estos efectos se v a n corrigiendo y se obtiene una imagen d e mayor calidad [STA-88]. En particular, los efectos q u e se p u e d e n corregir en la reconstrucción s o n la atenuación y la P S F . T a m b i é n podría incluirse la dispersión e n la matriz del sistema.
Capítulo 2
Planteamiento
del problema
pág. 32
a u n q u e no suele hacerse por la complejidad de los cálculos que requiere. A d e m á s , los m é t o d o s iterativos son m á s robustos e n los casos de pérdida parcial d e las proyecciones por movimiento involuntario del paciente u otro motivo [CAR-86] [AND-89]. El inconveniente esencial de los m é t o d o s iterativos e s su gran costo computacional, tanto en tiempo de cálculo c o m o e n cantidad de m e m o r i a . Puede calcularse con facilidad c ó m o crece la memoria
requerida
para trabajar con
la matriz de pesos en
un
ordenador. Para los valores estándar utilizados en la rutina clínica M=64 y A''=60, la memoria necesaria para la matriz de transición de cada corte tomográfico es d e 6 4 x 6 0 x 6 4 x 6 4 x 4 bytes, siendo 4 los bytes q u e ocupa un numero real e n float. Esto da un total d e m á s de 62 Mb. Si se trabajara con el doble de direcciones de proyección, su t a m a ñ o se multiplicaría por dos, pero si la imagen s e digitalizara, por ejemplo, e n matrices 128x128, el t a m a ñ o de la matriz crecería e n un factor 16 (más d e 200 Mbi). A su vez, el n ú m e r o d e operaciones aritméticas realizadas e n la reconstrucción depende directamente del número de elementos de la matriz de pesos. Esto conlleva q u e el tiempo necesario de computación en los m é t o d o s iterativos sea notablemente mayor al de la retroproyección filtrada. Para evitar los problemas d e memoria utilizada, cabría la posibilidad d e calcular los elementos d e la matriz en cada iteración, pero e s o haría m u c h o m á s lento el proceso de reconstrucción. No obstante, los elementos de esta matriz son mayoritariamente cero (matrices tipo sparse)
y se
puede
almacenar d e una m a n e r a compacta s e g ú n se describe en el siguiente capítulo. A u n q u e
las matrices así e m p a q u e t a d a s
pueden
reducir
su
t a m a ñ o e n m á s de un 7 5 % , sus dimensiones son un factor limitante m u y importante a tener e n cuenta a la hora de utilizar m é t o d o s iterativos.
Otro factor
a
determinar
es
el
número
de
iteraciones
que
deben
realizarse. S e g ú n se e x p o n e en los capítulos posteriores, los diferentes métodos iterativos de reconstrucción tomográfica suelen mostrar
una
convergencia parcial a la imagen ideal, con un acercamiento en las primeras iteraciones y un distanciamiento
posterior si el número
de
Capítulo 2
Planteamiento
del problema
pág. 33
iteraciones crece indefinidamente. Por esto, d e b e ser determinado en cada caso, preferentemente de f o r m a automática, el instante preciso en q u e el proceso iterativo debe interrumpirse para obtener los mejores resultados. Así m i s m o , es deseable poseer m e c a n i s m o s de regulación de la velocidad de convergencia del proceso iterativo, para acelerarlo o retardarlo según c o n v e n g a .
Capítulo 2
Planteamiento
del problema
pág. 34
Capítulo 3
Metodología
pág. 35
3. Metodología El f u n c i o n a m i e n t o d e los métodos d e reconstrucción puede evaluarse objetivamente analizando la s e m e j a n z a entre las características de la reconstrucción d e unas proyecciones con las d e la imagen original de q u e provienen. La obtención de una imagen y de sus proyecciones puede conseguirse de d o s maneras: a partir de un m o d e l o real, del q u e se c o n o z c a la g e o m e t r í a y la actividad d e cada z o n a , y de sus proyecciones experimentales o a partir de un modelo numérico y de un algoritmo de simulación. A u n q u e las proyecciones experimentales son preferibles a las simuladas, pues la simulación requiere algunas aproximaciones, d a d o que los m é t o d o s de reconstrucción d e b e n ser probados e n un gran n ú m e r o de condiciones para poder extraer conclusiones generales, e s
necesario
disponer d e un simulador que permita generar las proyecciones con gran facilidad. La validación de la simulación, que posibilita extrapolar los resultados de la simulación a casos reales, se consigue c o m p a r a n d o los resultados
obtenidos
mediante
las
proyecciones
experimentales
y
s i m u l a d a s d e un m i s m o modelo. Otra posibilidad q u e ofrece la simulación es
la
generación
de
proyecciones
que
no
incluyan
todas
las
degradaciones. D e esta manera, p u e d e estudiarse por separado el efecto de las diferentes degradaciones de las proyecciones sobre la imagen reconstruida. Finalmente, para valorar la calidad d e las reconstrucciones y cuantificar su s e m e j a n z a con la imagen original s e utilizan funciones de mérito ( F D M e n lo sucesivo).
El e s q u e m a de los diferentes estudios realizados e n este trabajo es: •
O b t e n c i ó n de proyecciones, tanto experimentales como simuladas por ordenador.
•
Análisis de los m é t o d o s de reconstrucción y ajuste de los parámetros de
ios q u e
dependen
mediante
la evaluación
con
FDM
reconstrucciones obtenidas c o n las proyecciones simuladas.
de
las
Capítulo 3
•
Metodología
pág. 36
Contrastación de los resultados obtenidos mediante simulación c o n los obtenidos d e las proyecciones experimentales.
En este capítulo, pues, s e e x p o n e n la metodología e m p l e a d a para la obtención
de
proyecciones,
tanto
experimentalmente
como
por
simulación, la reconstrucción y la evaluación de los resultados.
3.1. Adquisición de proyecciones experimentales El modelo físico utilizado consta de un t a n q u e de metacrilato de f o r m a cilindrica d e 2 2 c m de diámetro exterior c o n cierre hermético. En el interior del m i s m o , están dispuestos longitudinalmente 6 cilindros de metacrilato de 5, 4 , 3, 2, 1'5 y 1 c m d e diámetro respectivamente (Figura 3-1). Este t a n q u e s e llena de una disolución a c u o s a de ^^"^Tc, d e manera q u e la emisión radiactiva es uniforme e n todo el recipiente salvo en los cilindros, d o n d e es nula. Este m o d e l o será referenciado como m o d e l o de Jaszczak.
Figura 3-1 Sección digitalizada (512x512) del modelo de Jaszczak
Capítulo 3
Metodología
pág. 37
Las adquisiciones d e las proyecciones se han llevado a c a b o e n 60 direcciones equiespaciadas en 360° mediante u n a g a m m a c á m a r a A P E X S P 4 - E L S C I N T del Servicio de Medicina Nuclear del Hospital Provincial
de Barcelona.
Clínic
i
Han sido utilizados d o s colimadores distintos,
uno d e baja resolución (colimador- A ) y otro de alta resolución (colimadorB), a m b o s d e agujeros paralelos^
Para determinar los parámetros de la adquisición, s e han tenido e n cuenta los parámetros usuales utilizados e n la rutina clínica. Dado el t a m a ñ o del m o d e l o físico, las proyecciones s e han adquirido con los parámetros utilizados n o r m a l m e n t e en los estudios de cerebro. De esta m a n e r a se ha fijado
el
número
de
cuentas,
el
radio
de
giro
del
cabezal
de
la
g a m m a c á m a r a y el t a m a ñ o del píxel. La proyecciones experimentales obtenidas
han sido corregidas de decaimiento
de radiactividad.
Los
p a r á m e t r o s de adquisición están descritos en la tabla 3 - 1 .
Colimador Número de cuentas por corte tomográfico
A
B
400.000 (400KC)
200.000 (200KC)
17 cm
17cm
0.4587 cm
0.4717 cm
Radio de giro del cabezal Tamaño del píxel
Tabla 3 - 1 . Parámetros de adquisición de los estudios.
3.2. Simulación de un estudio S e g ú n se vio e n el capítulo anterior, las proyecciones de un objeto numérico se calculan con la expresión (ec. 2-12):
p = A-q + n
^ Aunque la estructura de ambos colimadores es análoga, el colimador A, por tener agujeros de mayor calibre, permite el paso de una mayor cantidad de fotones que el colimador B (mayor eficiencia), a costa de perder precisión en la dirección de incidencia (peor resolución, PSF de mayor tamaño)
Capítulo 3
Metodología
pág. 38
Para aplicar esta e x p r e s i ó n , es preciso construir la matriz de transición que
incluya
tanto
la
parte
geométrica
como
el
efecto
de
las
degradaciones. Posteriormente, a la proyección obtenida al multiplicar la imagen por esta matriz d e pesos se le a ñ a d e el ruido.
3.2.1. Cálculo de la matriz de transición El e l e m e n t o Ay de la matriz de pesos incluye el factor geométrico, la respuesta del detector, la atenuación y la dispersión. Por esto, p u e d e expresarse c o m o función d e tres factores:
donde el factor J>, tiene e n cuenta la dispersión de los fotones emitidos, atji tiene e n cuenta la atenuación y r,, es el factor geométrico de respuesta del
detector
(PSF).
Bajo
la asunción
de
que
los
tres
fenómenos,
dispersión, P S F y atenuación son independientes entre ellos, la función puede ser factorizada [LIA-92]:
Aj,=sjrat^,-r,
Cada factor s e obtiene por separado, de una forma exacta o a p r o x i m a d a , a partir del c o n o c i m i e n t o de la geometría y composición del conjunto objeto-detector. A continuación, se detalla el cálculo de cada uno de ellos.
3.2.2. Tratamiento de la dispersión Para cada bin de la proyección, los fotones detectados es la s u m a d e fotones
primarios
(que
no han sufrido dispersión
alguna) y
fotones
dispersados: p = p"''+
p'""
(ec3-1)
Capítulo 3
Metodología
pág. 39
La d e g r a d a c i ó n d e b i d a a la dispersión depende, entre otros factores, de la distancia del e m i s o r al detector y d e la geometría del objeto [FRE-91] [FRE-93]. S u inclusión en la matriz d e pesos d e manera exacta es un p r o b l e m a d e gran complejidad q u e e s c a p a a los objetivos d e e s t e trabajo. Por
tanto,
se
ha
buscar
una
forma
alternativa
de
simular
esta
degradación.
La mayoría de correcciones de la dispersión e n estudios reales están b a s a d a s e n análisis espectral d e la energía de los fotones detectados. La detección simultánea y separada d e fotones con energías dentro d e dos o m á s intervalos disjuntos, uno centrado en la energía d e emisión del isótopo y los otros e n energías menores, permite hacer una estimación de los
fotones
Compton
y
sustraerlos
del total.
Entre
la
abundante
bibliografía referente a este tema, puede destacarse [JAS-85], [KOR-88],
[YAN-88], [MAS-89], [LJU-90a], [KIN-92] y [BUV-95]. No obstante, una simulación la dispersión que incluya la evaluación de la energía de los fotones
dispersados
precisa
utilizar complejos
métodos
basados
en
procesos d e M o n t e Cario [BEC-82], [FLO-84], [LJU-89], [FRE-90] o [ROS-
90].
A l g u n o s autores h a s propuesto tratamientos alternativos d e la dispersión como
por e j e m p l o
[AXE-84],
[FLO-85]
y
[MSA-87].
Estos
métodos
s u p o n e n q u e el efecto de la dispersión de c a d a punto emisor es una función Function)
d e respuesta
sobre
la proyección
( S R F : Scatter
Response
q u e p u e d e ser modelizada. Si se realiza la aproximación d e
suponer la S R F invariante para todos los puntos d e la imagen, entonces, la obtención d e los fotones dispersados se consigue convolucionando la proyección c o n esta S R F . En este trabajo, siguiendo los resultados d e [MSA-87], para modelizar la dispersión se utiliza la S R F :
/ ( x ) = 0.035 •exp(-0.2-x)
(ec 3-2)
Capítulo 3
r^etodología
pág. 40
La expresión (ec 3-1) s e escribe:
p = pp" + pvUf
y
puede
ser
invertida
pasando
transformada de Fourier (FT- Fourier
(ec
al d o m i n i o
frecuencial
mediante
la
Transform):
FT(/>)
(ec 3-4)
FT(6 + / )
En r e s u m e n , utilizando
3-3)
la aproximación
anteriormente
citada, no
es
preciso incluir el efecto d e la dispersión e n la matriz d e transición, y a q u e puede corregirse por deconvolución sobre las proyecciones mediante esta última expresión. Otra importante ventaja de usar esta aproximación e s q u e la dispersión es m á s fácilmente simulable mediante la expresión ec 33. De esta m a n e r a , la ecuación de la generación de la proyección a partir de la i m a g e n , para proyecciones libres d e dispersión, p u e d e escribirse:
/7/"=4,(ar,,r,).^,.+«,
(ec3-5)
3.2.3. Inclusión de la PSF C u a n d o s e hace referencia a la PSF, respuesta del detector a un punto emisor, s e entiende e n sentido amplio. Es decir, no sólo se incluye el factor geométrico d e respuesta debido a la geometría del colimador, sino también t o d o s los f e n ó m e n o s degradantes asociados al detector, incluida la penetración septal, la dispersión dentro del colimador, el ruido del circuito d e localización d e impactos en el cristal de centelleo, etc. La parte geométrica
puede
modelizarse
matemáticamente
a
partir
del
conocimiento d e la estructura del colimador [MET-80], pero las restantes
Capítulo 3
Metodología
pág. 41
son d e difícil tratamiento matemático y precisan m é t o d o s d e M o n t e Cario para s u cálculo [DEV-90]. No obstante, la P S F d e un punto para un colimador d e agujeros paralelos es ajustable a una gaussiana [FOR-89] [HEB-92] c o n dispersión dependiente d e la distancia emisor-detector.
Por otra parte, y a q u e las imágenes están digitalizadas, se tiene q u e calcular la respuesta d e un pixel emisor a partir d e la P S F d e un punto emisor. La P S F d e un pixel es la integral sobre todos los puntos del pixel d e la P S F de un punto. Si se realiza la aproximación d e suponer q u e todos
los
puntos
del
pixel
tienen
igual
función
respuesta
( i n d e p e n d i e n t e m e n t e d e la posición del punto dentro del pixel), la P S F del pixel
se puede
respuesta
calcular
mediante
con la proyección
la convolución
geométrica
del pixel.
d e dicha
función
Entonces,
si la
proyección del pixel es de dimensiones mucho m e n o r e s q u e la gaussiana, el resultado d e convolucionar
u n a gaussiana
c o n dicha función es
a p r o x i m a d a m e n t e una gaussiana d e mayor dispersión y el ajuste de la P S F d e un pixel c o n una gaussiana es bueno. En caso contrario, es decir si la dispersión d e la P S F es m u y pequeña e n relación al t a m a ñ o de la proyección del pixel, prevalecerá la forma d e esta proyección y el ajuste a una gaussiana será deficiente. Esto ocurre a distancias m u y cortas del colimador. En r e s u m e n , a las distancias usuales emisor-detector, la P S F d e un pixel emisor e s correctamente aproximable por una gaussiana.
Para determinar la variación d e la P S F con la distancia del punto emisor al detector, se ha realizado una medición experimental d e la P S F en función d e la distancia para los d o s colimadores d e agujeros paralelos utilizados c o n el m o d e l o físico. H a n sido adquiridas proyecciones planas d e un punto e m i s o r a distintas distancias emisor-detector. El punto emisor s e ha c o n s e g u i d o blindando un pequeño depósito q u e contenía una gota d e una disolución d e ^^"^Tc, por todas partes m e n o s por un orificio. Para c a d a proyección, s e ha ajustado numéricamente una gaussiana a la P S F obtenida.
Mediante
una regresión
se ha c o m p r o b a d o
q u e , en las
Capítulo 3
pág. 42
Metodología
distancias usuales emisor-detector (de 10 a 50 cm), la variación de la dispersión (o) d e la gaussiana e n función de la distancia (z) del punto emisor al detector es lineal, obteniéndose las siguientes relaciones:
Colimador
A :
a = 0.0275 -z + O.l
coef.corr
= 0.9998
Colimador
B:
cr = 0.0172-7 +0.2
coef.corr
= 0.9997
d o n d e tanto a c o m o z están expresados e n c m .
detector
colimador
PSF
punto emisor
Figura 3-2. La PSF recoge aproximadamente el mismo número de cuentas con independencia de la distancia al detector (abarca un mismo ángulo sólido), aunque su dispersión aumenta con dicha distancia.
La PSF indica c ó m o se distribuye la energía emitida por cada punto de la imagen e n los diferentes píxeles de la proyección. Una
observación
importante es q u e el n ú m e r o de fotones detectados e s independiente d e la distancia. Esto es d e b i d o a que, dentro del rango d e distancias e n q u e se trabaja, el ángulo sólido de radiación detectada s i e m p r e e s el m i s m o (Figura 3-2). C o m o el área de la P S F e s proporcional a la radiación detectada, t o d a s las P S F han de tener el área normalizada a un m i s m o valor. Q u e d a indeterminado un factor de escala c o m ú n a todas las P S F . Puede s u p o n e r s e las gaussianas de área unidad. Cualquier otro valor, daría una reconstrucción igual pero reescalada.
Capítulo 3
Metodología
pág. 43
PSF del pixel /
>i
y
y+ly+2y+3...
Figura 3-3. Determinación del factor geométrico. Contribución del pixel / de la imagen al bin J las proyecciones
Una vez calculada la PSF, el factor geométrico R,, e s el área d e la P S F q u e c o r r e s p o n d e al pixel j de las proyecciones (Figura 3-3). Es la parte de radiación emitida e n el pixel i de la imagen (sobre un total de 1) q u e llega al bin j d e las proyecciones.
3.2.4. Inclusión de la atenuación La presencia d e un medio atenuante entre el emisor y el detector provoca q u e sólo un porcentaje de la radiación emitida llegue al receptor, el resto e s absorbido por el medio. Este porcentaje no d e p e n d e de la intensidad de e m i s i ó n . D e p e n d e sólo del grosor y naturaleza del m e d i o atenuante q u e la radiación atraviesa. La intensidad / q u e atraviesa
un
medio
a t e n u a n t e d e grosor d es proporcional a la intensidad emitida Ig. El factor de proporcionalidad at, es positivo, menor que uno y varía d e manera e x p o n e n c i a l inversa con el grosor [TUB-63] :
d o n d e ju e s el coeficiente de atenuación del m e d i o , que d e p e n d e de la naturaleza del m i s m o y de la energía de la radiación y s e determina e x p e r i m e n t a l m e n t e e n cada caso, a u n q u e para algunos c o m p u e s t o s está tabulado.
Por e j e m p l o , el coeficiente
de atenuación
del
agua
para
radiaciones de 140 KeV (energía d e los fotones y del ^^"^Tc) es de 0.15
cm-^ [TUB-63], [HAR-84].
Capítulo 3
l\/letodología
Si la radiación atraviesa diferentes
pág. 44
medios
atenuantes el factor d e
reducción d e la radiación e s el producto d e todos los factores:
exp(-^n,í/,)
at = g-^"'' • e-^''':..=
(ec 3-6)
k
d o n d e juk e s el coeficiente d e atenuación del ^-ésimo medio y ¿4 es el espesor del m i s m o q u e atraviesa el rayo e n consideración.
Para poder calcular el factor atenuante atji, e s necesario conocer el m a p a de atenuación
d e la i m a g e n , es decir,
el valor
del coeficiente
de
atenuación e n cada p u n t o d e la imagen. Este m a p a d e atenuación está n o r m a l m e n t e digitalizado c o n el mismo n ú m e r o de píxeles que la i m a g e n . El
mapa
d e atenuación
puede
ser determinado
experimentalmente
reconstruyendo las i m á g e n e s d e transmisión de la radiación de una fuente externa
al
paciente
simultáneamente
tomadas
al estudio
en
diferentes
direcciones
previa
o
[BAI-87] [BAR-97] o c o n la información
procedente d e imágenes obtenidas con T C (Tomografía Computerizada) o R M N (Resonancia Magnética Nuclear). Una vez determinado el m a p a de atenuación, el factor atenuante atj¡ d e la radiación del píxel / d e la imagen q u e e s detectada e n el bin j d e las proyecciones viene d a d o por la expresión análoga a e c 3-6, moviendo el índice k s o b r e todos los tejidos atravesados:
k
Se ha d e mencionar u n a última aproximación utilizada para simplificar los cálculos. S e ha supuesto q u e la atenuación era igual para toda la P S F d e un punto. C o m o coeficiente d e atenuación para una P S F s e ha t o m a d o la de su rayo central. D a d a la poca variación de los caminos para los diferentes rayos q u e p r o d u c e n la P S F d e un píxel emisor y la uniformidad
Capítulo 3
Metodología
pág. 45
de la mayoría d e los m a p a s d e atenuación utilizados, esta aproximación tiene efectos m í n i m o s . Los distintos espesores ¿4 pueden
calcularse
fácilmente utilizando el algoritmo descrito en [SID-85].
3.2.5. Empaquetamiento de la matriz de pesos. D a d o q u e e n un bin d e las proyecciones contribuyen pocos píxeles d e la i m a g e n , los e l e m e n t o s d e la matriz d e pesos s o n mayoritariamente nulos (matrices tipo sparse).
Esto permite guardarla d e manera
compacta
utilizando el m é t o d o descrito e n [BRO-89]. Un elemento d e la matriz q u e d a d e t e r m i n a d o por un índice fila, un índice columna y un valor real. Para a l m a c e n a r u n a matriz de este tipo basta c o n guardar tres vectores, el primero q u e c o n t e n g a todos los elementos d e la matriz diferentes de cero
(en
la
práctica
mayores
que
un
valor
umbral)
tomados
s e c u e n c i a l m e n t e por filas; el s e g u n d o q u e contenga el índice c o l u m n a de los e l e m e n t o s d e l primer vector y el tercero q u e contenga el índice del primer vector e n q u e se cambia d e fila, esto e s , el índice del elemento cuyo índice fila es m a y o r q u e el índice fila del e l e m e n t o anterior.
3.2.6. Simulación del ruido sobre las proyecciones El ruido producido durante la emisión y la detección causa q u e el n ú m e r o d e fotones d e t e c t a d o s e n cada bin no sea el q u e le correspondería teóricamente d a d a s la intensidad emisora de c a d a pixel y la contribución de c a d a pixel. S e g ú n se señaló e n el capítulo anterior, cada bin d e las proyecciones p u e d e ser descrito mediante una variable d e Poisson. El parámetro d e esta distribución d e probabilidad es s u valor m e d i o q u e , por otra parte, es el valor esperado [CUA-81]. Para simular el ruido e n cada bin, s e realiza un sorteo d e una variable de Poisson de parámetro igual a su valor teórico (sin ruido, pero c o n las d e m á s degradaciones) mediante las rutinas estándar d e generación d e números aleatorios encontradas, por ejemplo, en [PRE-92]. Dado q u e la dispersión d e una distribución de
Capítulo 3
Metodología
pág. 46
Poisson d e p e n d e d e la raíz cuadrada d e su valor e s p e r a d o , la variación del valor del pixel será porcentualmente m á s importante respecto el valor esperado cuanto m e n o r s e a éste.
3.2.7. Cálculo de la matriz de pesos A fin de recapitular t o d o lo expuesto hasta ahora sobre la matriz de pesos, a continuación se describe el algoritmo implementado para su cálculo. I. Reserva de m e m o r i a , lectura de los parámetros (tamaño del pixel, número d e píxeles, á n g u l o s de proyección, tipo de colimador) y del m a p a de atenuación. II. Para c a d a bin de la proyección e n cada ángulo a) Cálculo
del
rayo
que
incide
en
el
centro
del
bin
y
de
sus
intersecciones c o n los píxeles d e la imagen. b) Para c a d a pixel d e la i m a g e n . - Cálculo d e la distancia lateral, distancia ortogonal del pixel al rayo. - Cálculo de la distancia del pixel a la proyección. - Cálculo d e la g a u s s i a n a q u e se ajusta a la P S F correspondiente a e s a distancia. - Determinación del factor geométrico a partir de la gaussiana y la distancia lateral. - Si el factor g e o m é t r i c o e s mayor q u e un valor de corte • Cálculo del factor de a t e n u a c i ó n a partir de las intersecciones del rayo con los píxeles y el m a p a de atenuación. . Producto del factor g e o m é t r i c o c o n el factor d e atenuación • Si el peso resultante es mayor q u e el valor d e corte, guardar el p e s o y su posición dentro d e la matriz. C o n ello se c o n s i g u e evitar
el
crecimiento
excesivo
de
la
matriz
con
pesos
irrelevantes. III. G u a r d a r la matriz d e pesos según la estructura q u e tiene definida.
Capítulo 3
Metodología
pág. 47
3.2.8. Generación de las proyecciones La simulación numérica se realiza mediante un programa q u e realiza los procesos descritos con anterioridad: I. Lectura d e la i m a g e n a proyectar y de la matriz d e pesos. Lectura de los p a r á m e t r o s de la simulación (número de realizaciones de ruido, n ú m e r o inicial para la serie d e números aleatorios necesaria para el sorteo de las variables de Poisson, número de cuentas de la proyección). II. Cálculo de las proyecciones
mediante el producto de la imagen por la
matriz d e pesos. III. C o n v o l u c i ó n con la S R F para la simulación d e la dispersión. IV. Cálculo del n ú m e r o de cuentas d e las proyecciones. Ajuste del n ú m e r o de c u e n t a s al d e s e a d o mediante el reescalado d e las proyecciones. V. Para cada realización de ruido: a) Modificación del valor q u e inicializa las series aleatorias. b) Simulación del ruido. Para c a d a bin de las proyecciones: - Sorteo d e una variable d e Poisson de parámetro igual al valor del bin. - Sustitución del valor del bin por el valor obtenido e n el sorteo. c) Corrección d e dispersión por deconvolución c o n la S R F . d) Grabación de las proyecciones e n un archivo. Si no se simula ruido, los pasos convolucionar
y
deconvolucionar
III y V-(c) s e anulan, puesto con
la
misma
SRF
no
altera
que la
proyección inicial. En este caso, la simulación d e la dispersión no tiene ningún efecto.
El proyector ha sido objeto de diversos controles para garantizar su correcta
implementación.
En
primer
lugar,
se
han
simulado
las
proyecciones sin degradaciones de objetos de los q u e se podían calcular las
proyecciones
superposición
analíticamente
(objetos
d e elipses- las proyecciones
definidos analíticas
mediante de
una
la
elipse
p u e d e n encontrarse e n [KAK-84]) a fin de verificar los programas de proyección
por
comparación
de
las
proyecciones
obtenidas.
Capítulo 3
Metodología
pág. 48
Posteriormente, se han proyectado sin ruido y sin atenuación i m á g e n e s tipo delta para c o m p r o b a r las d i m e n s i o n e s de la P S F resultante. La atenuación
se
proyecciones
ha con
comprobado atenuación
mediante uniforme
la
y
obtención
sin
PSF
ni
analítica
de
ruido
su
y
c o m p a r a c i ó n c o n las proyecciones obtenidas mediante el simulador en las mismas condiciones. La validez de las aproximaciones realizadas e n el simulador s e obtiene c o n la corroboración de la coincidencia d e
los
resultados obtenidos e n simulaciones numéricas con los obtenidos c o n las proyecciones experimentales.
3.2.9. Estudios simulados Han
sido
generadas
gran
cantidad
de
proyecciones
con
diferentes
parámetros d e adquisición de diversos m o d e l o s numéricos. De todos ellos se han o b t e n i d o diferentes realizaciones de ruido mediante diferentes sorteos e n c a d a uno d e los bin de las proyecciones. La estructura del simulador ha permitido generar proyecciones con ninguna, alguna o todas las d e g r a d a c i o n e s incluidas. De esta manera, se ha podido estudiar por separado el efecto d e c a d a una de estas degradaciones en la i m a g e n resultante, tanto si la degradación era corregida e n la
reconstrucción
c o m o si no. A d e m á s , s e ha podido c o m p a r a r la reconstrucción obtenida con el m o d e l o original, lo q u e ha permitido evaluar objetivamente
la
calidad d e la reconstrucción.
En los estudios numéricos simulados se ha utilizado principalmente los modelos
numéricos
reflejados en
la Figura
3-4.
El primero e s
una
digitalización del m o d e l o real disponible e n una matriz d e 64x64 píxeles. Se hace notar q u e esta imagen es binaria. Los bordes de los m á r g e n e s del objeto no tienen valores intermedios (grises) c o m o correspondería si se digitalizara la i m a g e n , sino que valen O o 1 según su valor teórico sea mayor o m e n o r que 0'5. Esto se ha hecho así para q u e las regiones d e interés no tuvieran los bordes difusos y se pudieran evaluar las diferentes
Capítulo 3
Metodología
pág. 49
F D M c o n precisión. En contrapartida, se lia perdido parecido con el m o d e l o real, s o b r e todo en los círculos de radio menor. El m a p a de atenuación de este modelo se ha definido uniforme con / / = 0 ' 1 5 c m ' \ igual q u e el m o d e l o real. El s e g u n d o modelo es la simplificación de una sección de
tórax.
Esta
formado
por
una
elipse
de
actividad
uniforme,
correspondiente a los tejidos blandos del tórax, d o s regiones elípticas sin actividad, q u e simulan los pulmones, un aro d e mayor actividad, que representa una sección del corazón, y dos círculos de actividad media que corresponden
al esófago y columna vertebral. En e s t e caso
la
atenuación no es constante en todo el objeto, sino que cada región tiene la atenuación del tejido que representa. En la Figura 3-4, j u n t o a los modelos
numéricos
numeración
con
se
muestran
sus
mapas
de
atenuación
y
la
la q u e serán referidas las distintas regiones e n
la
evaluación de las F D M .
En la Figura 3-5, s e muestran los sinogramas del modelo de Jaszczak. Los
5
primeros
corresponde
a
están una
obtenidos
adquisición
mediante
simulación
experimental
del
y
modelo
el
último
real.
El
s i n o g r a m a (a) corresponde a la proyección sin degradaciones del modelo numérico
(proyecciones
de
los
círculos
negros
superpuestas
a
la
proyección del f o n d o uniforme, q u e es sinusoidal por estar el modelo ligeramente descentrado). En (b), se incluyen
las degradaciones
de
atenuación y P S F . S e obsen/a la difuminación d e los contornos (debida a la PSF) y la disminución del efecto de los círculos negros e n algunas proyecciones (debido a la atenuación). Los tres sinogramas siguientes c o r r e s p o n d e n a tres realizaciones d e ruido d e las proyecciones
con
atenuación, P S F y dispersión con 50 Kc (1 Kc=1.000 cuentas), 2 0 0 Kc y 800 Kc respectivamente. Se observa el aumento de ruido al disminuir el n ú m e r o de cuentas, así como una buena concordancia entre el modelo real c o n el m o d e l o simulado correspondiente (200 Kc).
Capítulo 3
pág. 50
Metodología
u = 0'5 cm
a)
b)
Figura 3-4 Principales modelos numéricos utilizados en las simulaciones. a) digitalización del modelo de Jaszczak a 64x64 píxeles, su mapa de atenuación y la numeración de las regiones de interés, b) modelo numérico que simula una sección del tórax, su mapa de atenuación y la numeración de las regiones de interés.
Capítulo 3
Metodología
•
pág. 51
Figura 3-5. Aspecto del sinograma del modelo de Jaszczak. a) Proyecciones sin degradaciones, b) Proyecciones con degradación de atenuación y PSF. c) d) e) Proyecciones con atenuación, PSF y ruido (correspondiente a 50 Kc, 200 Kc y 800 Kc respectivamente), f) Proyecciones experimentales del modelo real (200 Kc).
Capítulo 3
Metodología
pág. 52
3.3. El proceso de reconstrucción En
este
punto
se
expone
los
rasgos
comunes
a
todas
las
reconstrucciones realizadas, sea cual s e a el método d e reconstrucción empleado.
Las
consideraciones
particulares
de
cada
método
están
expuestas e n los capítulos concernientes a esos m é t o d o s : •
D a d a la g r a n cantidad de variantes d e los métodos de reconstrucción, s e han estudiado los algoritmos m á s genéricos d e cada especie. S e ha procurado utilizar la m e n o r cantidad posible d e restricciones sobre los m é t o d o s para q u e los resultados f u e r a n lo más generales posibles.
•
Los estudios se han realizado sobre uno o dos m o d e l o s solamente. No se ha realizado una evaluación a f o n d o de la d e p e n d e n c i a d e los resultados c o n el m o d e l o . Los resultados, por tanto, están vinculados al modelo. No obstante, la metodología de estudio de la d e p e n d e n c i a de
la
calidad
de
la
imagen
respecto
a
los
parámetros
de
la
reconstrucción es g e n e r a l . Este procedimiento p u e d e ser e m p l e a d o para estudiar la reconstrucción óptima de cualquier modelo de q u e s e disponga. •
Todas
las
proyecciones
numéricamente
como
las
reconstruidas, reales,
han
sido
tanto
las
corregidas
simuladas mediante
deconvolución de dispersión utilizando la expresión (ec 3-4) con la S F R definida en (ec 3-2). •
Para realizar las reconstrucciones y d e m á s cálculos s e ha utilizado un o r d e n a d o r H P 9 0 0 0 / 7 2 0 workstation (57 MIPS, 17 M F L O P S ) .
Como
lenguaje de programación se ha e m p l e a d o el C s o b r e HP-UNIX. •
A la hora d e p r o g r a m a r los métodos d e reconstrucción, en todos los casos la i m a g e n s e ha asignado a vectores float c o n precisión sencilla (4 bytes por pixel). A u n q u e visualmente la precisión que ofrecen los vectores char
(1 byte por pixel) o short
int (2 bytes por pixel) e s
suficiente, e n los cálculos se ha utilizado más precisión para permitir una evolución continua de la imagen reconstruida e n las diferentes iteraciones. El uso d e precisión doble se desestimó porque a u m e n t a el
Capítulo 3
tiempo
Metodología
de cálculo, c o n s u m e
más
pág. 53
memoria y no aporta
mejoras
sensibles a la reconstrucción. •
En los m é t o d o s de reconstrucción y corrección iterativos, para poder evaluar
la
evolución
de
las
imágenes
respecto
al
número
de
iteraciones e m p l e a d a s , a cada iteración s e ha grabado la i m a g e n . Posteriormente, sobre cada una de las imágenes q u e f o r m a n estas s e c u e n c i a s se ha calculado el valor de las diferentes F D M . Esto ha permitido obtener su variación en función del número de iteraciones. •
D a d a la naturaleza aleatoria del ruido, para poder asegurar q u e los resultados son lo m á s generales posibles y desligar las conclusiones d e la realización particular de ruido, todos los cálculos han sido hechos s o b r e 10 realizaciones diferentes de ruido y se han obtenido
los
p r o m e d i o s y la dispersión. En algunos casos puntuales e n q u e la dispersión
resultaba
elevada,
se
han
llegado
a
promediar
100
realizaciones d e ruido diferente. •
El h e c h o de expresar los resultados con p r o m e d i o s no implica q u e no s e haya estudiado el comportamiento de los métodos e n realizaciones concretas de ruido. En los casos más significativos, a d e m á s
de
evaluar los promedios, se ha evaluado el comportamiento singular de las diferentes realizaciones de ruido para determinar la existencia o no de anomalías e n los resultados particulares. •
Para visualizar las imágenes resultantes, se han convertido a formato 1 byte o 2 bytes. En el primer caso se han normalizado de m a n e r a que su valor m á x i m o fuera 255 (1 byte puede expresar valores d e O a 255).
3.4. La evaluación de imágenes mediante funciones de mérito El primer problema q u e aparece cuando se quiere valorar una imagen reconstruida, es q u e difícilmente su calidad será óptima e n todos los aspectos. Por ejemplo, suele ocurrir q u e una imagen con b u e n contraste tiene una importante presencia de ruido y una imagen sin ruido suele
Capítulo 3
Metodología
pág. 54
tener los contornos m u y difuminados. Dependiendo d e la finalidad d e la imagen, una característica u otra puede ser aceptable o no. Por ejemplo, si se quiere detectar lesiones de t a m a ñ o pequeño, la presencia d e ruido será m u y molesta, c o n lo q u e será preferible imágenes m á s suaves, pero si s e quiere cuantificar niveles d e actividad, será preciso disponer d e un buen contraste. La evaluación de la imagen d e p e n d e d e l propósito para el que ha sido adquirida [HAN-90] [HER-91a].
Existen, b á s i c a m e n t e , d o s procedimientos de valoración d e los m é t o d o s de reconstrucción: los m é t o d o s subjetivos y el cálculo de funciones d e mérito sobre las i m á g e n e s reconstruidas. E n el primer caso, las i m á g e n e s obtenidas
en
las
diferentes
pruebas
se
presentan
a
un
número
significativo d e o b s e r v a d o r e s cualificados q u e valoran diferentes aspectos de la m i s m a , c o m o p u e d e ser la detectabilidad de una z o n a d e diferente actividad, el aspecto general de la i m a g e n , el ruido, etc. C o n s u s resultados
se
construyen
las d e n o m i n a d a s
curvas
ROC
(Receiver
Operating Characteristic) c o n las q u e s e determina el mejor m é t o d o o parámetros d e reconstrucción e n cada c a s o [DEB-88] [HAN-90] [GUI-91].
La
valoración
de
los resultados
puede
hacerse
también
mediante
funciones q u e cuantifiquen las características que determinan la calidad de la i m a g e n . Estas funciones son d e n o m i n a d a s figuras o funciones d e mérito (FDM). H a n sido propuestas diferentes F D M para evaluar la calidad
de las imágenes e n [DEB-88] [HER-89] [HER-91b] [MILL-92], [KIN-94] o [MUR-94]. En este trabajo se empleará únicamente e s t e método, pues e n todos los casos se c o n o c e d e a n t e m a n o la imagen ideal q u e debiera dar la reconstrucción, lo q u e permite el uso d e muy diversas F D M . Las F D M utilizadas e n este trabajo pueden clasificarse en tres categorías, las referentes a la globalidad d e la imagen, las referentes a una o varias regiones d e interés d e la imagen y las funciones estadísticas, utilizadas fundamentalmente estadísticos.
para evaluar i m á g e n e s reconstruidas c o n m é t o d o s
Capítulo 3
Metodología
pág. 55
3.4.1. Figuras de mérito globales Estas F D M cuantifican el parecido de la imagen reconstruida con el m o d e l o original. Las funciones utilizadas son: •
Coeficiente d e correlación (CC) entre la reconstrucción y la imagen original. Está definido por la expresión:
CC =
d o n d e VM significa valor medio (suma de valores dividido por el n ú m e r o de píxeles) y el superíndice
designa a la imagen de
referencia ( i m a g e n original). El coeficiente d e correlación tiene un valor m á x i m o d e 1 correspondiente al caso d e q u e a m b a s imágenes coincidan e x a c t a m e n t e . No obstante, es preciso hacer notar q u e esta función es invariante por cualquier transformación lineal q u e afecte a la imagen entera, por ejemplo, sumarle un m i s m o valor a todos los píxeles de la i m a g e n o multiplicar el todos los píxeles de la imagen por un m i s m o factor. Esto refleja, q u e C C es alto c u a n d o a m b a s imágenes s e parecen e n forma y los valores de los píxeles se
mantienen
relacionados c o n una función lineal, a u n q u e no coincidan e n valor exactamente.
3.4.2. Figuras de méritos sobre regiones de interés Estas F D M se utilizan para valorar los detalles de la imagen y la correcta relación entre diferentes partes de la misma. Para su uso, e s preciso contar c o n una m á s c a r a de la imagen que defina las diferentes regiones. Es c a s o de estudios simulados, la máscara se obtiene a partir de la imagen original. Para estudios reales, la máscara puede ajustarse e n una imagen reconstruida o definirse a través del conocimiento de la geometría del m o d e l o físico e m p l e a d o . Las F D M de este tipo utilizadas s o n :
Capítulo 3
Metodología
pág. 56
C o n t r a s t e {CON) d e una región 1 s o b r e otra región 2 que la incluye. Definido por la e x p r e s i ó n :
CON =
VM\-VM2 VMl + VM2
d o n d e VM e s el valor medio dentro d e la región correspondiente. El contraste, así definido, oscila entre O (una región contraste
nulo
equivale a decir q u e tiene el mismo valor que su periferia y, por tanto no e s distinguible) y el valor teórico correspondiente a los valores exactos del modelo. Para regiones d o n d e el valor teórico e s O, el contraste m á x i m o e s 1 . Coeficiente d e variación (CV) de una región 1, definido por
CF = ^ ^ - 1 0 0 VMl
d o n d e a^ e s la desviación estándar del valor de los píxeles de la z o n a 1. Para regiones q u e tienen un valor uniforme, C V es cero. La relación señal-ruido {SNR: Signa! to Noise Ratid)
de una región 1
sobre otra región 2 q u e la incluye, definida por:
SNR =
VMl - VM2 G2
Es decir, S N R c o m p a r a las diferencias de valores medios de dos z o n a s c o n el ruido d e la zona envolvente. Una S N R baja indica q u e una región no e s diferenciable de las fluctuaciones de la zona q u e la incluye. Por tanto, s e deseará q u e el valor de esta F D M sea lo m á x i m o posible.
Capítulo 3
¡[Metodología
pág. 57
3.4.3. Figuras de mérito estadísticas Otra m a n e r a d e evaluar imágenes es realizar algún tipo de test de fiabilidad estadística. Estas F D M s e utilizan, s o b r e todo, para evaluar imágenes
reconstruidas con m é t o d o s estadísticos
[VEK-87]
[COA-91]
[HUD-94], a u n q u e igualmente p u e d e n ser e m p l e a d a s en otros casos. Las F D M estadísticas utilizadas en este trabajo s o n : •
Verosimilitud L(Xi,...,Xn)
directa
{DL:
Direct
Likelihood)
La
verosimilitud
d e u n a muestra d e una variable aleatoria es el valor d e la
función d e d e n s i d a d cuando las variables aleatorias tienen el valor q u e indica la muestra [CUA-81] [CUA-96]. Para el caso de reconstrucción tomográfica d e imágenes, s e g ú n se explica c o n detalle e n el capítulo 6, las proyecciones son una realización de las variables aleatorias. Por tanto, la verosimilitud es la probabilidad las proyecciones condicionada al valor de los píxeles de la imagen reconstruida. Da una medida de la probabilidad d e q u e una imagen reconstruida hubiera g e n e r a d o las proyecciones d e q u e se dispone. Por tanto, es deseable q u e adquiera el valor m á s alto posible (valores entre O y 1). Considerando cada bin d e la proyección pj c o m o una variable aleatoria d e Poisson,
.\fx.\/
L(q) = PÍp\q) = Yle
-
'MXM /=i
F r e c u e n t e m e n t e , en lugar de la verosimilitud, s e calcula s u logaritmo, p u e s s u evaluación resulta m á s sencilla y, p o r ser el logaritmo una función m o n ó t o n a creciente, el valor en q u e alcanza su m á x i m o coincide
con
el
de
la verosimilitud.
Algunos
autores,
incluyen
directamente el logaritmo e n la definición d e la verosimilitud [HUD-94]. En este trabajo se utiliza esta definición. Entonces, su valor oscila entre -oo y O y s u expresión es:
Capítulo 3
Metodología
MxNf
./=i
V
pág. 58
MxM
MxM
'=1
'=1
Verosimilitud c r u z a d a {CL: Cross Likelihood).
C u a n d o se dispone d e
dos conjuntos de proyecciones, se puede calcular la
probabilidad
condicionada de obtener el s e g u n d o conjunto de proyecciones c o n la imagen reconstruida a partir del primer conjunto de proyecciones. A n á l o g a m e n t e a DL, s e calcula su logaritmo. Su expresión es:
MxNf
./=1
V
MxM
MxM
/=1
/=!
d o n d e el superíndice 'B' hace referencia a que se trata del s e g u n d o conjunto d e proyecciones.
Capítulo 4
La retroproyección
filtrada
pág. 59
4. La retroproyección filtrada La
retroproyección
filtrada
(FBP) es
el
método
de
reconstrucción
tomográfica q u e , por s u sencillez y facilidad de aplicación, es utilizado d e m a n e r a estándar e n los equipos convencionales d e g a m m a g r a f í a . S u aplicación a E C T está descrita por primera vez e n [SHE-74].
En el planteamiento d e la FBP, la reconstrucción d e la imagen a partir de las proyecciones s e consigue por inversión de la transformada d e R a d o n por m e d i o d e un algoritmo basado en el teorema d e la sección central. La fórmula d e inversión resultante es teóricamente exacta para imágenes y proyecciones analíticas. No obstante, la aplicación de este m é t o d o d e reconstrucción
en
casos
reales
supone
asumir
una
serie
de
aproximaciones. El uso d e dichas aproximaciones y la presencia d e las d e g r a d a c i o n e s hace q u e , en general, las imágenes obtenidas sean d e baja calidad. Existen, no obstante, diferentes m é t o d o s d e corrección o c o m p e n s a c i ó n d e las degradaciones, ya sea sobre las proyecciones antes de la reconstrucción o sobre la imagen una v e z reconstruida. En este capítulo, se e x p o n e la formulación matemática e impiementación d e la FBP y el estudio d e la calidad d e la imagen en función de los parámetros del filtro e m p l e a d o y el número d e iteraciones d e un método iterativo de corrección d e la atenuación.
4.1. Formulación matemática. En primer lugar, se recuperan algunos resultados del Análisis Funcional q u e serán precisos e n el desarrollo d e la F B P . Pueden
encontrarse
r e s ú m e n e s de esta teoría, por ejemplo, en [GON-77], [BAR-81], [LEW-83] o [HER-87]. A continuación, se expone la teoría q u e permite construir el algoritmo
de
matemática
la de
degradaciones.
F B P y, finalmente, las
distintas
se muestra
técnicas
de
la
impiementación
compensación
de
las
Capítulo 4
La retroproyección filtrada
pág. 60
4.1.1. Filtros. Se entiende p o r filtrado d e u n a función f(x) mediante la función d e filtro lineal e s p a c i a i m e n t e invariante g(x) a la función obtenida d e convolucionar la función original c o n la función de filtro
/ ( x ) - / ( x ) * g ( x ) = £ / ( « ) • g{x - a) da
S e g ú n el t e o r e m a de convolución, esto e s equivalente a la modificación d e s u e s p e c t r o d e frecuencias c o n la transformada d e Fourier d e la función d e filtro, y a q u e
F T ( / ( x ) * g ( x ) ) = F T ( / ( x ) ) • F T ( g ( x ) ) = F{U) • G{U)
es decir, la t r a n s f o r m a d a d e Fourier d e u n a convolución e s el producto d e t r a n s f o r m a d a s d e Fourier d e las d o s funciones. Por tanto, se p u e d e filtrar una función mediante el paso al espacio d e las frecuencias:
/ ( x ) =/(x)*g(x) = FT-'(F(W) • G M )
A su vez, el filtrado sucesivo por d o s o más filtros p u e d e
hacerse
s i m u l t á n e a m e n t e por la propiedad asociativa de la convolución:
f{x)*h{x)
= FT-'{F{u)-G{u)-H{u))
(ec- 4-1)
4.1.2. Teorema de la sección central El t e o r e m a d e la s e c c i ó n central relaciona la transformada d e Fourier unidimensional de las proyecciones c o n la transformada bidimensional d e l objeto. E n efecto, s e a q(x,y) u n a i m a g e n y Q(u,v) s u transformada d e Fourier. S e tienen las siguientes relaciones:
Capítulo 4
La retroproyección
Q(u,v) =
filtrada
pág. 61
q{x,y) • exp{^-27ri(ux + vy)^dxdy (ec-4-2)
qix,y)
=
Q{ii,v) • exp{2;ri{ux + vy)^dudv
Sea Sofw) la transformada de Fourier unidimensional de la proyección en la dirección &.
Sg{w)= r Pg{t)-exp{-27riwt)
dt
(ec-4-3)
A partir d e estas relaciones, se puede formular el siguiente t e o r e m a :
Teorema
de la sección
central:
S e a Q(w,0) el valor d e Q(u,v) a lo largo d e
la línea q u e forma una ángulo 6'con el eje u. Entonces:
(ec-4-4)
Q{w,0) = S,{w)
La demostración para el ángulo d e proyección O e s inmediata. Basta con expresar Q(u,0) e n función de la proyección:
Q(u,0) -
q{x,y) • exp{-27vi
ux)dxdy
•00 J-CO ICO
q{x,y)dy =
• exp(-2;r/ia:)í&
Pa{t)exp{-27ri ut)dt = Sq{U)
ya q u e la integral entre corchetes d e la segunda igualdad corresponde al valor d e la proyección para un ángulo igual a cero. En la última igualdad se ha sustituido la variable m u d a x por la variable t para q u e la expresión coincidiera con (ec-4-3).
Capítulo 4
La retroproyección
pág. 62
filtrada
La generalización a un ángulo cualquiera, puede hacerse mediante un cambio d e las c o o r d e n a d a s cartesianas a otras giradas un ángulo 0.
'
señó»'
COS0
^-sen^
COS0J
yy)
Entonces:
P,{t)=
[
q{s,t)ds
con lo q u e
Sg{w) = £ P^{t)-exp{-27riwt) q{t,s)ds
dt
-expi^-lniwt)
q{t,s) • exp(-2;ri
dt
wt) dt ds
ya q u e 5 = 0 c o r r e s p o n d e a la recta q u e pasa por el origen e n la dirección 0.
4.1.3. Reconstrucción tomográfica por retroproyección. A partir del t e o r e m a d e la sección central puede implementarse un algoritmo d e reconstrucción tomográfica. C o n este propósito, s e escribe la igualdad
(ec-4-2)
expresando
las
frecuencias e n c o o r d e n a d a s polares
li^^y)
=
variables
{u,v)
del
espacio
{w,0):
Q{^,'V) • exp[2;^/(wx + vy)) dudv Q(W,
0) •
exp(2;r / w • (x eos 0 + y sen 0)^
wdwd0
de
Capítulo 4
La retroproyección
pág. 63
filtrada
d o n d e la w q u e multiplica a los diferenciales e n la última
expresión
proviene d e l j a c o b i a n o del cambio d e variables. Si se cambian los límites de
integración
Q{-w,9)
de
las variables
= Q{w,0+7r),
polares,
teniendo
en
cuenta
que
y utiliza el hecho q u e í = xcos6' + >^sen6', esta última
expresión p u e d e escribirse c o m o :
q(x,y)=
Q[w,0)-exp(2;riw-1) Jo J-cc
\
/
V
w dwdO /
Haciendo uso del t e o r e m a de la sección central
qix,y)
=
Sg{w) • expilTii w-t)w
dw de (ec-4-5)
Pg{t)d9=
[ Pg{xcos9 +
ysene)d0
d o n d e Pg{t) corresponde a la expresión entre corchetes y es el filtrado d e la proyección m e d i a n t e un filtro de r a m p a \ Esta última igualdad, indica q u e p u e d e obtenerse la imagen a partir de la integración d e la proyección filtrada. A l proceso d e sumar para todos los ángulos y sobre todos los puntos del rayo el valor de la proyección se le d e n o m i n a retroproyección.
4.1.4. Compensación de las degradaciones. Las fórmulas hasta ahora deducidas, son exactas para
proyecciones
ideales. A l aplicarlas a proyecciones de S P E C T , debido a las diferentes degradaciones, el resultado e s d e baja calidad. A continuación, se describe c ó m o p u e d e mejorarse la calidad la reconstrucción resultante por medio d e técnicas d e compensación d e las degradaciones.
^ Se conoce como filtro de rampa al filtro cuya descomposición espectral es el valor absoluto de la recta identidad H(u)=\u\.
Capítulo 4
La retroproyección
filtrada
Corrección del ruido mediante filtros
pág. 64
paso-baja.
En ECT, se ite q u e el ruido es blanco, es decir, que tiene una d e s c o m p o s i c i ó n espectral uniforme^ [BAR-81] [ROS-95a]. Dado q u e e n las imágenes con q u e s e trabaja en S P E C T las c o m p o n e n t e s de baja y media frecuencia tienen m u c h o m á s peso q u e las de alta frecuencia, es e n la alta frecuencia d o n d e el ruido predomina sobre la señal. El uso del filtro de rampa amplifica este f e n ó m e n o al potenciar las altas frecuencias, provocando la existencia de un ruido e n forma de picoteado sobre las reconstrucciones.
La
principal estrategia
para
limitar el ruido en
la
reconstrucción e s filtrar las proyecciones con un filtro que reduzca las c o m p o n e n t e s de alta frecuencia (filtro paso-baja). La inevitable reducción de las altas frecuencias d e la imagen producida en el filtrado diminuye la definición contornos) frecuencias
de
la
imagen
[GON-77]. menos
(suavizado
Cuanto
mayor
y
fenómeno
sea
de
la reducción
Gibbs de
importancia tendrá el ruido, pero más
en
las
los altas
suavizada
resultará la reconstrucción. En cada caso d e b e ser determinado el filtrado que proporcione la relación deseada entre el nivel de ruido y la suavidad de la i m a g e n .
Diversos autores han propuesto modificaciones del filtro de rampa para tratar el p r o b l e m a del ruido. Generalmente, estos filtros pueden ajustarse a cada situación mediante uno o m á s parámetros, de forma q u e pueda decidirse las características del filtrado. S e ha de hacer notar q u e la utilización del filtro de rampa modificado equivale a un doble filtrado, e n primer lugar por el filtro d e rampa, necesario para la F B P , y en s e g u n d o lugar, con un filtro paso-baja (ec 4-1). Es por ello q u e t o d o s los filtros q u e aparecen a continuación están multiplicados por el filtro de rampa, y
^ Suponer que el ruido es blanco es aproximadamente equivalente a suponer que el ruido en cada bin es independiente (no está correlacionado) respecto el de los otros bins [BAR-81].
La retroproyección
Capítulo 4
filtrada
pág. 65
ocasionalmente por otros factores, y representados con dicho filtro de referencia. En un abuso de lenguaje, se designan los filtros resultantes con el n o m b r e del filtro paso-baja utilizado. En [PUC-97], hay un extenso sumario
de
los
diferentes
filtros
utilizados
en
las
principales
g a m m a c á m a r a s del mercado. Los filtros más usuales t a m b i é n
están
descritos e n [GON-77] [HER-87] [APE-88] o [TSU-94]. Son: •
Filtro d e S h e p p - L o g a n , definido e n el espacio d e las frecuencias por la expresión: í
w H(W) = w • sinc 2w..
^
d o n d e la función sinc(x) está definida por:
sinc(x) =
sen{^ • x) n-x
y yvmax es la frecuencia máxima (frecuencia d e Nyquist). En la Figura 4-1 se puede observar su perfil frente al del filtro de rampa. El filtro de S h e p p - L o g a n tiene una resolución similar a la del filtro de rampa pero reduce la amplificación del ruido e n las altas frecuencias.
Filtro de Shepp-Logan
H(SN)
Figura 4 - 1 . Filtro de Shepp-Logan (línea continua) obtenido al multiplicar el filtro de rampa (línea discontinua) con un filtro paso-baja (línea puntos). Este último filtro está representado a una escala diferente (valor máximo igual a 1)
La retroproyección
Capítulo 4
pág. 66
filtrada
Filtro d e Butterworth, definido por la expresión:
w
i El perfil d e
este
filtro está
1+
w b-w...
controlado
por tres
parámetros,
el
parámetro d e resolución a ( 0 < a < 0 . 9 ) , la frecuencia de corte 6 y el orden del filtro c. El primero de los parámetros sirve para potenciar las frecuencias
de
la
proyección
que
resultan
disminuidas
en
la
adquisición a causa de la respuesta del detector. Dado q u e el efecto de este parámetro es de difícil predicción, suele t o m a r el valor 0. Los otros dos son responsables de la suavidad del filtro. Su variación da una
flexibilidad
propiedades.
En
muy la
grande
para
crear
Figura 4-2 s e ve
filtros
con
la variación
diferentes
del filtro
de
Butterworth en función de sus parámetros.
Filtro de Butterworth. a=0 b=0.2 c=1,3, 6, 9
Filtro de Butterworth. a=0 b=0.2, 0.3, 0.4, 0.5 c=6
H(w,
a)
b)
Figura 4-2. Variación del filtro de Butterworth al variar sus parámetros: a) variación respecto el parámetro c, b) variación respecto el parámetro b.
La retroproyección
Capítulo 4
pág. 67
filtrada
Filtro d e H a n n i n g , definido por la expresión:
H[W)=
w\-{\ + a-\wfj
0.5 + 0.5- eos
W
1
v
^
J )
Este filtro d e p e n d e de tres parámetros a n á l o g o s a los del filtro de Butterworth. Su variación con dichos parámetros queda reflejada en la Figura 4-3.
Filtro de Hanning. a=0 b=1.5 c=1,2, 3, 5
Filtro de Hanning. a=0 b=1,1.5, 2, 2.5 c=3
H(w)
a)
b)
Figura 4-3. Variación del filtro de Hanning en relación a sus parámetros, a) variación respecto el parámetro b. b) variación respecto el parámetro c.
Corrección de la PSF. Filtro de Metz La corrección de la PSF dependiente del punto emisor (distancia emisordetector) no es abordable en el algoritmo de la F B P . Para tratar la PSF, se ha d e prescindir d e su dependencia con la distancia y suponer que es la m i s m a para c a d a punto emisor (se toma una P S F promedio, e n general la correspondiente al centro de la imagen). A u n q u e la c o m p e n s a c i ó n del efecto de la P S F p u e d e realizarse por filtrado d e la imagen reconstruida.
La retroproyección
Capítulo 4
pág. 68
filtrada
el método m á s c o m ú n es la deconvolución sobre las proyecciones d e esta PSF promedio.
El filtro de Metz está definido por la expresión:
H{w) = w
S{w)
con
S{w)=C-e
siendo C y k dos constantes que determinan la P S F promedio. Es la composición d e tres filtros: el filtro de r a m p a , un filtro paso-baja para reducir el ruido y el filtro inverso de la P S F promedio (Figura 4-4a). La Figura 4-4b muestra la variación del filtro d e Metz con el exponente c.
Filtro de Metz. a=1 b=0.01 c=2, 4, 6 ,8
Obtención del filtro de Metz.
H(w)
H(w)
a)
b)
Figura 4-4. a) Obtención del filtro de Metz (línea gruesa) como producto de tres filtros: el filtro de rampa (línea discontinua), un filtro paso-baja (línea fina) y la inversa de la PSF (línea de puntos). Los dos últimos filtros están escalados a un valor máximo de 2. b) Variación del filtro de Metz con el parámetro c
Capítulo 4
La retroproyección
filtrada
pág. 69
Corrección de la atenuación. El método de Chang. El efecto
de
la
atenuación
es
la disminución
de
cuentas
en
las
proyecciones. La retroproyección de dichas proyecciones produce una i m a g e n q u e , a su vez, contiene un menor n ú m e r o de cuentas q u e la imagen real. El déficit de cuentas no tendría m a y o r trascendencia si fuera uniforme en toda la reconstrucción. Sin embargo, está más a c e n t u a d o en los píxeles centrales de la imagen q u e en los periféricos [GEL-88]. La rectificación de este f e n ó m e n o , importante para la correcta cuantificación e interpretación de las imágenes, ha sido objeto d e numerosos estudios. Los diferentes m é t o d o s de corrección desarrollados pueden clasificarse básicamente
en
tres
grupos:
pre-corrección
(corrección
sobre
las
proyecciones m e d i a n t e filtrado) [GLI-91] [KIN-91] [MAZ-93], corrección en la reconstrucción mediante la modificación del algoritmo FBP ( W F B P Weighted
Filtered
Back
Projection)
u otros m é t o d o s [BEL-79] [CEN-79]
[GUL-81] [GUL-85] [MAN-88a] [MAN-88b] [SIN-88] [ZEN-91] [PAN-96] [WEL-97] y post-corrección (corrección en la i m a g e n reconstruida) [CHA-
78] [CHA-79] [WAL-81] [MAN-87] [TSU-89] [PAN-93].
De entre los m é t o d o s de post-corrección destaca el método desarrollado por C h a n g [CHA-78] [CHA-79], q u e es el m é t o d o empleado en este trabajo. La estrategia utilizada en este método, para paliar el déficit de cuentas en la reconstrucción es la normalización del número cuentas en cada píxel con un factor que c o m p u t e , en promedio, la atenuación de los rayos emitidos d e s d e dicho píxel. Para la deducción de la fórmula de corrección en el método de C h a n g , en primer lugar, se proyecta y retroproyecta una imagen tipo 5 en un medio con atenuación uniforme, sin tener en consideración otras degradaciones. La imagen puede expresarse mediante:
q{x,y)=AS{x-XQ,y-y^)
Capítulo 4
La retroproyección
Siendo A s u intensidad y
{xo,yo)
pág. 70
filtrada
s u ubicación. Si se sustituye esta i m a g e n
en la expresión (ec 2-9):
Pg{t) = AS ( x o cos6' + S Q n 6 - t ) - Q x p { - j U - l g )
Su t r a n s f o r m a d a d e Fourier es :
Sg{w)=A-exp{2m
• ( x g e o s 0 + yo senO))• exp(-/u-lg)
La proyección para el ángulo 0 filtrada c o n el filtro de rampa vale:
Pg(t) =
donde
""" A • exp(-ITTÍW• ( x q eos6* +
la integración
está
limitada
sen(9-1))• e x p ( - / / / ^ ) w d w
por los valores
máximos
d e la
frecuencia. La retroproyección de esta expresión e s :
q(x, y ) =
A • exp{2n:iw{{x - Xg )cos 0 + {y- yo )sen 0)) • exp(- julg ) wdwd0
que e n el punto
{xo,yo)
es:
^(^o'J^o) =
['"'"
Á •exp(-///^)wJw d0
Se define el factor d e corrección para el punto
{xo,yo)
c o m o el valor q u e s e
obtendría e n el caso d e q u e no hubiera atenuación (//=0) dividido por el calculado c u a n d o hay atenuación:
Capítulo 4
La retroproyección
c(Xo,>'o)=
filtrada
pág. 71
^
exp(-/.
I,)d0
q u e c u a n d o se utiliza funciones digitalizadas es:
i-texp(-/./í)
donde
es la distancia d e s d e el pixel k hasta el e x t r e m o d e la imagen en
el ángulo diy e\ sumatorio se extiende a las A/^proyecciones.
Una v e z deducida para una imagen puntual, la corrección de una imagen cualquiera se consigue calculando estos coeficientes
para todos
los
puntos d e la i m a g e n (una imagen cualquiera es la s u m a d e imágenes puntuales).
La
imagen
es
multiplicada,
pixel
a
pixel, por
el
factor
correctivo correspondiente. La compensación e n una imagen no puntual m e d i a n t e este procedimiento es aproximada y a q u e en un bin de la proyección intervienen diferentes píxeles de la i m a g e n , cada uno con una atenuación
distinta.
Si el mapa
de
atenuación
no es
uniforme,
la
exponencial d e esta última expresión debe sustituirse por:
=
• , -ív
^
(ec-4-6)
;=i /=yt
d o n d e e n el s e g u n d o sumatorio sólo s e suma para los píxeles incluidos en el rayo q u e parte del pixel emisor y llega al detector. /^^ es la longitud d e la intersección
del
rayo con el pixel y-ésimo
(Figura 4-5)
calcularse fácilmente con el algoritmo descrito e n [SID-85].
y
puede
Capítulo 4
La retroproyección
filtrada
pág. 72
El resultado d e s p u é s d e multiplicar la imagen por estos factores
de
corrección mejora los resultados de la reconstrucción, aunque, a v e c e s , no
suficientemente.
Esto
ocurre
generalmente
cuando
el
mapa
de
atenuación no es uniforme en todo el objeto. Entonces, es posible aplicar este m é t o d o iterativamente. Para ello s e d e b e disponer de un proyector que nos permita obtener las proyecciones reconstruida. originales.
Estas
proyecciones
La diferencia
entre
resultantes de la
calculadas
ambas
se
es tomada
contrastan como
imagen con
una
las
nueva
proyección q u e se reconstruye y se corrige al igual q u e la imagen original (filtro de Metz y Chang). Esta segunda reconstrucción (imagen diferencia) se s u m a a la primera obteniéndose una segunda i m a g e n . El proceso puede repetirse tantas v e c e s c o m o sea preciso para la obtención d e una imagen d e calidad.
píxel j-ésimo
/ /
rayo central del píxel Á;-ésimo al detector en el ángulo de proyección Q,
Figura 4-5. Intersección de un rayo con un píxel
4.1.5. Impiementación de la FBP. Según
la
expresión
(ec
4-5),
la
imagen
puede
obtenerse
por
retroproyección de la proyección filtrada mediante un filtro de rampa. Parte de
las d e g r a d a c i o n e s
pueden
ser corregidas de una
manera
aproximada, y a sea a priori sobre las proyecciones, mediante un filtrado.
Capítulo 4
La retroproyección
o, a posteríori,
filtrada
pág. 73
s o b r e la imagen. En este trabajo s e utilizará un filtro de
Metz para c o m p e n s a r el ruido y la P S F y el m é t o d o de C h a n g iterativo para c o m p e n s a r la atenuación. Estos son los m é t o d o s de corrección m á s c o m u n e s y c o n mejores resultados. Las proyecciones s e s u p o n d r á n corregidas de a n t e m a n o de dispersión. En caso contrario, t a m b i é n podría incluirse esta corrección en el filtro.
En la implementación numérica del algoritmo han d e tenerse e n cuenta una serie d e consideraciones: •
Las
funciones
están
discretizadas
en
una
serie
de
puntos
equidistantes. Por tanto las expresiones integrales d e b e n sustituirse por sumatorios y la transformada d e Fourier por la transformada de Fourier
discretizada
Transform,
(se e m p l e a
el algoritmo
FFT: Fast
Fourier
descrito, por ejemplo e n [GON-77] o e n [PRE-92]).
El uso d e la trasformada d e Fourier discreta para calcular convoluciones plantea
un problema. La transformada d e Fourier discreta de la
imagen
discretizada
no
coincide
con
la
discretización
de
la
t r a n s f o r m a d a continua de la imagen [BAR-81]. Es preciso orlar la i m a g e n y las proyecciones
unidimensionales
con ceros
(en este
trabajo se han añadido tantos ceros como píxeis de la imagen), antes d e proceder a s u transformación de Fourier, sino la hipótesis de periocidad d e la T F produce la aparición de solapamientos. •
La retroproyección en el caso discreto presenta, t a m b i é n , alguna particularidad. El valor del bin y-ésimo d e las proyecciones
debe
repartirse por t o d o s los píxeles d e la imagen q u e contribuyen al rayos u m a q u e incide e n e s e pixel. El hecho q u e el rayo no pase por los centros d e los píxeles implica la utilización d e una interpolación para repartir el valor d e la proyección a los píxeles adyacentes al rayo. La interpolación utilizada es la lineal.
Capítulo 4
La retroproyección
filtrada
pág. 74
4.1.6. El algoritmo IFBP. A continuación, a m o d o d e compendio d e lo expuesto hasta el m o m e n t o , se describe
el algoritmo
de
la FBP
con
corrección
iterativa
de
la
atenuación por el m é t o d o d e C h a n g . Los pasos de los q u e consta s o n :
I. Inicialización a) La i m a g e n se inicializa a cero.
q=0.
b) Mediante el m a p a d e atenuación del objeto se calcula los factores de corrección de c a d a pixel d e la imagen (imagen d e corrección). II. Reconstrucción mediante F B P . a) Se filtran las proyecciones. El filtrado se realiza con el paso al espacio d e frecuencias, producto por la transformada del filtro y regreso
mediante
la transformada
inversa
a
las
componentes
espaciales. b) C a d a bin de la proyección filtrada se retroproyecta, es decir, s e distribuyen las cuentas del mismo a lo largo d e todo el rayo q u e incide e n ese bin. Para ello: - S e calcula el rayo q u e incide en el centro del bin. - S e recorren t o d o s los píxeles de la imagen y s e calcula la distancia d d e su centro al centro del rayo. - Si e s t a distancia e s m e n o r que 1 (medida e n píxeles) se le a ñ a d e al pixel de la imagen el valor del bin multiplicado por el factor 1-d c) Si e s preciso, se g r a b a la imagen. III. Corrección de la atenuación a) Se multiplica la imagen por la i m a g e n de corrección calculada a partir del m a p a d e atenuación (corrección Chang) b) Si s e d e s e a iterar el proceso, se repite el siguiente proceso tantas v e c e s c o m o sea necesario: - S e proyecta la i m a g e n actual - S e obtiene la diferencia de la proyección original (sin filtrar) con la proyección calculada.
Capítulo 4
La retroproyección
filtrada
pág. 75
- La proyección diferencia se reconstruye c o m o la original (mediante filtrado, retroproyección y corrección d e Chang) y se s u m a a la primera, obteniéndose la segunda imagen.
4.1.7. Rango de validez de la FBP. Valoración de las aproximaciones realizadas en su implementación. En la exposición del desarrollo de la FBP, se ha mencionado un conjunto d e supuestos y aproximaciones q u e condicionan el uso d e e s t e método. En primer lugar, s e g ú n se explica e n 4.1.3, debido a las simetrías de las c o o r d e n a d a s polares se puede realizar la retroproyección tanto a partir de la integración entre O y 27t como a partir de la integración entre O y 7t de las proyecciones filtradas. Sin e m b a r g o , si se d i s p o n e de proyecciones en un conjunto de ángulos distinto de estos dos, al aplicar la FBP aparecen imágenes
muy
proyecciones
defectuosas.
han
de
Esto
desestimarse
ocurre por
cuando
movimiento
parte
de
las
involuntario
del
paciente durante la adquisición. La F B P es m u y sensible a la pérdida de proyecciones en alguno de los ángulos. T a m b i é n es importante q u e no haya truncaciones e n las proyecciones. La reconstrucción mediante FBP de proyecciones truncadas produce artefactos e n la imagen reconstruida
[OGA-84] [MÜL-96].
Otro factor a tener e n cuenta es q u e , aunque inicialmente el ruido de cada píxel
es
independiente,
el
proceso
de
reconstrucción
tiende
a
correlacionar el ruido de píxeles adyacentes. Esto significa q u e si, por ejemplo, un píxel determinado tiene un valor superior a su valor esperado, los píxeles vecinos también tenderán a tener un valor superior a su valor e s p e r a d o . La corrección del ruido mediante un filtro paso-baja reduce la varianza total del ruido pero a u m e n t a este efecto. La consecuencia sobre la imagen es q u e desaparece el ruido llamado d e "sal y pimienta" (ruido diferente
píxel
a
píxel)
pero
pueden
aparecer
artefactos
o
grumos
(cúmulos de varios píxeles con un valor superior o inferior a su valor
Capítulo 4
La retroproyección
filtrada
pág. 76
esperado) [ROS-95a], c u y a presencia e s m e n o s deseable, pues p u e d e n inducir a error sobre la existencia de z o n a s de mayor o menor actividad en un estudio.
Finalmente, la corrección de la atenuación con el método de C h a n g , presupone disponer del m a p a de atenuación de a n t e m a n o , ya sea c o n imágenes d e transmisión (TC o gammagráficas) o c o n el conocimiento a priori
de la g e o m e t r í a y composición del objeto. A d e m á s , si se quiere
aplicar iterativamente, e s preciso implementar un proyector q u e incluya las degradaciones. C o n todo ello, la principal ventaja d e la FBP respecto los m é t o d o s iterativos (facilidad de cálculo) desaparece, pues aplicar esta corrección repercute e n un incremento notable del costo computacional.
4.2. Descripción de las pruebas realizadas Para realizar este estudio han sido obtenidas proyecciones del m o d e l o experimental
descrito
en
el
apartado
3 . 1 . con
los
parámetros
de
adquisición: Colimador B de agujeros paralelos (ver 3.2.3.) 6 0 proyecciones equiespaciadas en 360°. 6 4 bins d e 0.4717 c m en las proyecciones. Radio d e giro de la g a m m a c á m a r a de 17 c m . 2 0 0 Kc d e promedio por corte tomográfico. En total han sido o b t e n i d a s 10 proyecciones independientes del m o d e l o . T a m b i é n has sido g e n e r a d o s numéricamente 3 conjuntos de proyecciones de 100 Kc, 2 0 0 Kc y 800 Kc respectivamente. Estos conjuntos
han
constado d e 10 proyecciones, cada una con una realización diferente d e ruido y los m i s m o s parámetros geométricos que el modelo real. Estas proyecciones incluyen la simulación d e atenuación, PSF, dispersión y ruido.
Las
proyecciones
han
sido
corregidas
de
dispersión
por
deconvolución de la S R F correspondiente y reconstruidas mediante F B P utilizando
un
filtro
de
Metz,
con
los
valores
de
la
gaussiana
Capítulo 4
La retroproyección
filtrada
pág. 77
correspondientes a la PSF del centro de la imagen, y corrección de atenuación por el m é t o d o de C h a n g iterativo. S e han reconstruido las proyecciones para c a d a valor del parámetro c del filtro de Metz (ver 4.1.4) c o m p r e n d i d o entre 0.5 y 5 de 0.5 e n 0.5. Las i m á g e n e s reconstruidas han sido e v a l u a d a s c o n las figuras de mérito: CC, CV, CON y SNR. CON
y SNR
Las F D M
se han evaluado sobre todos los circuios (la graficación
c o r r e s p o n d e al círculo 1, dado que es el que m á s píxeles tiene y, por tanto, estadísticamente los resultados ofrecen m e n o r dispersión) y
CV
sobre la zona 7 del modelo (fondo del modelo) (las zonas del modelo están definidas e n 3.2.9).
Igualmente, d a d o q u e el método de Chang estaba pensado originalmente para i m á g e n e s c o n un mapa de atenuación uniforme, para poder evaluar sus prestaciones e n imágenes con atenuación no uniforme han sido g e n e r a d a s un conjunto de proyecciones del m o d e l o de tórax (ver 3.2.9) con 2 0 0 kc. De este modelo no ha sido posible tener un modelo real, por lo q u e
los resultados
obtenidos
por simulación
no han
podido
ser
corroborados experimentalmente. Los parámetros de la simulación han sido los m á s usuales e n las pruebas clínicas: C o l i m a d o r A de agujeros paralelos (ver 3.2.3.) 6 0 proyecciones equiespaciadas e n 360°. 6 4 bins d e 0.566 c m en las proyecciones. Radio d e giro de la g a m m a c á m a r a de 20.4 c m . 2 0 0 Kc de promedio por corte tomográfico.
Para evaluar e s t e modelo han sido calculadas las mismas F D M . En las gráficas están representadas CON
y SNR
de la zona 3
(miocardio)
respecto a su interior y de la zona 1 (pulmón) respecto al f o n d o del m o d e l o (las z o n a s del modelo están definidas e n 3.2.9). El número de iteraciones efectuadas ha sido de 15.
La retroproyección
Capítulo 4
pág. 78
filtrada
4.3. Resultados obtenidos En este capítulo se ha creído conveniente mostrar bastantes i m á g e n e s para q u e el lector pueda establecer la relación que hay entre los valores de las F D M y el a s p e c t o de la i m a g e n . En capítulos posteriores
la
presencia d e i m á g e n e s e n los resultados e s mucho menor.
En
la
Figura
4-6,
se
puede
observar
la
evolución
de
la
imagen
reconstruida con las iteraciones. Las imágenes sucesivas (dispuestas por filas) s o n el resultado obtenido en las diferentes iteraciones del proceso de reconstrucción para el modelo simulado de 2 0 0 K c y utilizando un filtro de Metz con e x p o n e n t e = 1 . La primera imagen (arriba a la izquierda) corresponde a F B P sin corrección de atenuación. Es m á s oscura q u e las d e m á s puesto q u e al no estar corregida de atenuación, el n ú m e r o total d e cuentas e s claramente inferior. Si se modificara la escala de grises para observar c o n mayor nitidez esa imagen, se produciría la saturación d e n u m e r o s o s píxeles de las imágenes restantes.
. ?' •
Figura 4-6. Evolución de la imagen con el número de iteraciones para la reconstrucción mediante IFBP con filtro de Metz (exponente=1) del modelo de Jaszczak simulado con 200KC.
Capítulo 4
La retroproyección
filtrada
pág. 79
Cualitativamente, e n la secuencia s e aprecia q u e las imágenes s o n m á s s u a v e s al principio y, a medida q u e se itera, el ruido va amplificándose paulatinamente. T a m b i é n se observa q u e las primeras i m á g e n e s s o n m e n o s contrastadas q u e las últimas. Cuantitativamente, sobre cada una de las i m á g e n e s s e evalúan las diferentes F D M . Los resultados expuestos e n las gráficas siguientes son el promedio sobre 10 realizaciones d e ruido de estos valores.
La Figura 4-7 muestra la variación del coeficiente d e correlación (CC), el coeficiente de variación (CV) y el contraste {CONl) y relación señal-ruido {SNRlf
del círculo d e mayor t a m a ñ o en función del número de iteraciones
e n el m o d e l o d e J a s z c z a k simulado d e 200Kc (4 gráficas primeras) y para el m o d e l o real equivalente (4 gráficas restantes). Las diferentes líneas d e las gráficas corresponden a diferentes valores del e x p o n e n t e e n el filtro de Metz utilizado e n la reconstrucción
La línea continua representa los
valores obtenidos c o n el exponente igual a 0'5, la línea de puntos y rayas para el e x p o n e n t e 1 , las líneas de rayas son de los exponentes 2, 3, y 4 respectivamente (la longitud d e las rayas es m e n o r cuanto m a y o r es el exponente) y la línea de puntos corresponde al exponente 5. En la iteración 1 s e indica los valores obtenidos mediante F B P sin corrección de atenuación.
Las iteraciones
siguientes
se refieren
a
las
sucesivas
iteraciones del m é t o d o de C h a n g . Estos criterios s e han mantenido en t o d a s las gráficas d e F D M en función del número d e iteraciones d e este capítulo.
Respecto la variación del valor d e las F D M con el número de iteraciones, se observa q u e no todas tienen un comportamiento igual. Mientras CC, CV y SNRl
alcanzan un valor óptimo tras el q u e e m p e o r a n , CONl
sigue
^ El número 1 en CONl y SNRI hace referencia a la zona 1 del modelo (círculo negro de mayor tamaño- ver 3.1). En lo sucesivo se seguirá este criterio.
Capítulo 4
La retroproyección
A
2
2
6
8
10 12 14 16
pág. 80
filtrada
o
2
4
6
8
10 12 14 16
núnnero de iteraciones
rimero de iteraciones
4
4
6
8
10 12 14 16
6
8
10 12 14 16
njrero de ¡teradones
número de ¡teradones
4
4
6
8
ID 12 14 16
número de ¡teradones
6
8
10 12 14 16
número de ¡teradones
tu
o
2
4
6
8
10 12 14 16
número de ¡teradones
0
2
4 6
10 12 14 16
número de ¡teradones
Figura 4-7. l\/lodelo simulado de Jaszczak 200 Kc (cuatro gráficas superiores) y modelo real equivalente (cuatro gráficas inferiores). Representación del CC, CV, CONl y SNRl en función del número de iteraciones. La iteración 1 corresponde a FBP sin corrección de atenuación y las siguientes a las sucesivas iteraciones de IFBP. La línea continua corresponde a un exponente 0'5 y las demás a un exponente 1, 2, 3, 4 y 5 respectivamente.
Capítulo 4
La retroproyección
filtrada
pág. 81
a u m e n t a n d o hasta el final. La degradación d e CC, CV y SNRl
c o n el
n ú m e r o d e iteraciones es m á s acentuada cuanto mayor es el exponente del filtro. El n ú m e r o d e iteraciones precisas para alcanzar el valor m á x i m o e s m e n o r cuanto m a y o r es el exponente. Se observa, también, q u e para e x p o n e n t e s altos, el mejor valor obtenido e n las diferentes F D M es inferior al obtenido con exponentes m e n o r e s . Esto indica q u e globalmente la imagen e s d e peor calidad. Por otra parte, si el valor del e x p o n e n t e es bajo, la imagen e s m e n o s ruidosa [CV tiene un valor menor) y un mejor valor d e CC, pero el contraste es menor. Los resultados obtenidos c o n el m o d e l o real t i e n e n gran concordancia respecto los obtenidos c o n el m o d e l o simulado.
En la Figura 4-8, las tres primeras filas muestran las imágenes obtenidas para
el
modelo
d e Jaszczak
simulado
c o n 200Kc.
Las tres
filas
c o r r e s p o n d e n a las reconstrucciones realizadas c o n tres valores distintos del e x p o n e n t e
del filtro de Metz, 0.5, 1, y 3 respectivamente.
Las
c o l u m n a s c o r r e s p o n d e n correlativamente a la imagen obtenida c o n FBP, a la i m a g e n e n la iteración d o n d e CC es máxima y a la imagen obtenida a las 16 iteraciones. Las tres últimas filas muestran las mismas imágenes para el modelo real. En a m b a s series puede observarse, en concordancia c o n las gráficas, q u e el contraste a u m e n t a con las iteraciones y con el e x p o n e n t e al igual q u e la presencia d e ruido sobre la reconstrucción. En la primera c o l u m n a puede apreciarse una distribución no uniforme de la actividad e n la i m a g e n . La periferia del cilindro es m á s activa q u e el centro. Esto está d e acuerdo c o n la no corrección d e la atenuación sobre la imagen [GEL-88].
En la Figura 4 - 9 , s e puede c o m p r o b a r cómo varían los resultados c o n el n ú m e r o d e cuentas. Las gráficas corresponden al modelo d e Jaszczak simulado, las 4 primeras con lOOKc y las 4 restantes con 800Kc. En estas gráficas s e p o n e d e manifiesto q u e cuanto mayor es el n ú m e r o de
Capítulo 4
La retroproyección
filtrada
pág- 82
Figura 4-8. Imágenes del modelo simulado 200Kc (tres filas primeras) y del real (tres restantes). Para cada modelo, las tres filas corresponden sucesivamente al exponente 0.5, 1 y 3 del filtro de Metz. Las 3 columnas son, correlativamente, la imagen de FBP, la imagen en la iteración donde CC es máximo y en la iteración 16.
Capítulo 4
La retroproyección
filtrada
pág. 83
cuentas y, por consiguiente, menor la importancia del ruido, el n ú m e r o d e iteraciones q u e p u e d e n hacerse sin q u e la imagen s e degrade es mayor, a l c a n z á n d o s e valores de las F D M mejores. La m e n o r presencia de ruido permite un mejor ajuste de la imagen a las proyecciones. Otro hecho destacable es q u e el uso d e exponentes altos está limitado por el número d e cuentas. S e g ú n s e observa e n la gráfica d e CC para 100Kc, sólo con los
d o s valores
m á s bajos
del exponente
se
encuentran
valores
aceptables. Para valores superiores del exponente las imágenes tienen valores peores d e CC y CON!.
Para 800Kc, el uso d e e x p o n e n t e s
m a y o r e s d a b u e n o s resultados a un n ú m e r o m e n o r d e iteraciones.
Por otra parte, para los valores inferiores del e x p o n e n t e (0.5 y 1), CONI parece
ser independiente
del nivel
de
ruido.
Las líneas
continua
(exponente 0.5) y d e punto-raya (exponente 1) d e las gráficas d e 100 Kc y 800 Kc s o n prácticamente superponibles (también c o n las m i s m a s líneas del m o d e l o de 2 0 0 Kc - Figura 4-7 -). No ocurre lo mismo c o n SNR], claramente
dependiente
del n ú m e r o
de cuentas.
Para
exponentes
superiores, CONI alcanza peores valores con 100 Kc.
En la Figura 4 - 1 0 s e muestran las imágenes d e los modelos d e Jaszczak s i m u l a d o c o n 100 K c (las tres filas superiores) y 800 Kc (las tres filas inferiores), en las m i s m a s condiciones de la Figura 4-8. En ellas se pone de manifiesto tanto la mayor calidad d e las i m á g e n e s como las mejores prestaciones u s a n d o exponentes altos para 800 K c q u e para 100 Kc. La diferencia de calidad entre las imágenes de 100 Kc y 8 0 0 Kc con e x p o n e n t e 3 es m u c h o mayor q u e la diferencia d e calidad observada para el e x p o n e n t e 0.5.
Capítulo 4
La retroproyección
1.00
1 ' 1 ' 1 • 1 • 1 '
1
1
filtrada
pág- 84
1
-
-
ü ü 0.80
'
-
- .• 0.70
m
X
\ .\ \
\
• , , . i\ \ . 4
>
6
X 8
, . 10
1 • 12 14
16
4
número de iteraciones
6
8
10
12
14
16
14
16
número de iteraciones
CC z
4
6
8
10
12
14
16
o
2
4
número de iteraciones
6
8
10
12
número de iteraciones
1
1
1
I
I
1
1 '
O O
I
2
4
6
8
10
12
14
I
I 10
I
16
número de iteraciones
I
1
12
1 14 16
número de iteraciones
1
-
a.
O ü
I
/
^
//•vi
V
I
^
^ \
z OT 3
1 .
. 1 1
2
4
6
8
10
12
número de iteraciones
14
16
D
2
4
6
1 8
1 1 10 12
I 14
I 16
número de iteraciones
Figura 4-9. Modelo de Jaszczak simulado con 100 Kc (cuatro gráficas superiores) y 800 Kc (cuatro gráficas inferiores). Representación del CC, CV, CONl y SNRl en función del número de iteraciones. La iteración 1 corresponde a FBP sin corrección de atenuación y las siguientes a las sucesivas iteraciones de IFBP. La línea continua corresponde a un exponente 0'5 y las demás a un exponente 1, 2, 3, 4 y 5 respectivamente.
Capítulo 4
La retroproyección
pág. 85
filtrada
fe
1,
/
-
'
Tí
.
.
" R -
Figura 4-10. Imágenes del modelo simulado con 100 Kc (tres primeras filas) y 800 Kc (tres restantes). Para cada modelo, las tres filas corresponden sucesivamente al exponente 0.5, 1 y 3 del filtro de Metz. Las 3 columnas son, correlativamente, la imagen de FBP, la imagen en la iteración donde CC es máximo y en la iteración 16.
La retroproyección
Capitulo 4
pág. 86
filtrada
En la Figura 4 - 1 1 , esta representado el mejor valor alcanzado por las diferentes F D M e n f u n c i ó n del exponente para el m o d e l o de J a s z c z a k simulado (200 Kc). S e observa q u e para los tres valores primeros del exponente, los resultados s o n parecidos, aunque algo mejores con el exponente 1 para SNRl.
E n cambio, los mejores resultados obtenidos con
exponentes superiores s o n peores. Esto significa q u e , globaimente para exponentes altos la reconstrucción e s d e peor calidad. S o l a m e n t e el m á x i m o d e CONl s e mantiene independiente del valor del exponente. No obstante, c a b e destacar q u e , de acuerdo con la Figura 4 - 9 , mientras el valor m á x i m o d e CONl s e alcanza en la última iteración para e x p o n e n t e s bajos, para e x p o n e n t e s altos se alcanza con F B P sin corrección d e atenuación ( i m á g e n e s c o n distribución n o uniforme d e cuentas entre el centro y la periferia d e la imagen).
Mínimos de CV en función del exponente del filtro de Metz
Máximos ce en función del exponente del filtro de Metz
40.0 > 35.0 "30.0
4
5
Máximos CONI en función del exponente del filtro de Metz
5 O O ^
1.0, 0.9 0.8 0.7 0.?
w •g O .i
0. 1 0.. •ra 0.. E 0.: 0. o.c
3
4
Máximos SNR1 en función del exponente del filtro de Metz o: Z
0
2
exponente
exponente
4.0 3.0 2.0
E 1.0
2
3
4
exponente
3
4
5
exponente
Figura 4-11 Mejores valores de las FDM en función del exponente del filtro utilizado en la reconstrucción. Valores correspondientes ai modelo de Jaszczak simulado (200 Kc).
Capítulo 4
La retroproyección
pág. 87
filtrada
La Figura 4-12 m u e s t r a los resultados obtenidos c o n el m o d e l o d e tórax (con m a p a de a t e n u a c i ó n no uniforme) c o n 2 0 0 Kc. Las F D M m o s t r a d a s s o n ce,
SNR d e las z o n a s 1 (pulmón) y 3 (corazón) (ver Figura
CV, CONy
3-4). El contraste teórico d e la z o n a 1 e s 1 y el d e la z o n a 3 e s 0.5.
1.00
> O
2
4
6
8
10
12
14
4
16
6
8
10
12
14
16
1
1
número de iteraciones
número de iteraciones
' 1 • 1' 1
1
1
1
•
8
m
4
W
3
. 1
10
12
14
16
O
número de iteraciones
2
1
1
4
6
,
1
8
1
.
10
1
.
1
12
,
1
14
16
14
16
número de iteraciones
ce 2
O ü
4
6
8
10
12
número de iteraciones
14
16
4
6
8
10
12
número de iteraciones
Figura 4-12 Modelo de tórax simulado 200 Kc. Representación del CC, CV, C0N3, SNR3, CON] y SNR] en función del número de iteraciones. La iteración 1 corresponde a FBP sin corrección de atenuación y las siguientes a las sucesivas iteraciones de IFBP. La línea continua corresponde a un exponente 0'5 y las demás a un exponente 1, 2, 3, 4 y 5 respectivamente.
Capítulo 4
La retroproyección
filtrada
pág. 88
A u n q u e el connportamiento es m u y parecido al observado en las gráficas del modelo d e Jaszczak, sí existe una diferencia remarcable. En ningún caso los mejores resultados para CC, SNRl
y SNR3 s e obtienen en la
iteración 2 (correspondiente a F B P con corrección d e C h a n g no iterativa), siempre s o n obtenidos a partir d e la iteración 3. T a m b i é n es destacable los valores bajos obtenidos por CONl en F B P .
En la Figura 4 - 1 3 p u e d e observarse las imágenes d e este modelo para tres e x p o n e n t e s diferentes (0.5, 1 y 3 dispuestos correlativamente e n las filas) y para la F B P (primera columna), m á x i m o de CC (segunda columna) e iteración 16 (tercera columna). P u e d e comprobarse la poca definición de los p u l m o n e s en la primera columna, y el aumento del contraste y del ruido c o n las iteraciones y el exponente.
Figura 4-13. Imágenes del modelo de tórax simulado (200 Kc). Las tres filas corresponden respectivamente al exponente 0.5, 1 y 3 del filtro de Metz. Las 3 columnas son, correlativamente, la imagen de FBP (iteración 1), la imagen en la iteración donde CC es máximo y en la iteración 15.
Capítulo 4
La retroproyección
filtrada
pág. 89
4.4. Discusión de ios resultados La primera observación q u e p u e d e liacerse es la gran similitud de los resultados obtenidos con las proyecciones s i m u l a d a s y de los resultados obtenidos con las proyecciones experimentales (Figura 4-8). Esta similitud valida la simulación implementada y las aproximaciones realizadas e n la m i s m a . Otro
h e c h o destacable e s que
la corrección
iterativa de
la
atenuación mejora los resultados d e FBP, pues, a d e m á s de corregir la distribución no uniforme de cuentas sobre el fondo, mejora los valores de las F D M . Sólo CON
obtiene en algún caso su mejor valor con FBP,
a u n q u e d a d a la incorrecta distribución de cuentas sobre el f o n d o de la i m a g e n e n estos c a s o s , esta mejora no es muy significativa.
En
cuanto
al
método
iterativo
de
reconstrucción,
se
observa
una
d e p e n d e n c i a entre e l número d e iteraciones a realizar, el e x p o n e n t e del filtro d e Metz e m p l e a d o y el ruido sobre las proyecciones. Para niveles de cuentas altos, si el exponente es bajo, s e necesitan más iteraciones para obtener resultados análogos a los obtenidos con e x p o n e n t e s m a y o r e s . En este
sentido
el
exponente
puede
ser
interpretado
como
un
factor
regulador d e la velocidad de convergencia del proceso iterativo, c o m o un factor d e aceleración del m i s m o . No obstante, e n general, el valor del e x p o n e n t e tiene un m á x i m o tras el cual la calidad de la imagen obtenida es claramente inferior. Este valor umbral d e p e n d e del número d e cuentas, siendo mayor c u a n t a s más cuentas tengan las proyecciones. En este sentido, s i e m p r e s e podría aumentar el exponente para conseguir la mejor imagen a la s e g u n d a iteración ( F B P con Chang sin iterar), pues en este caso no sería necesaria la construcción del proyector con lo que el m é t o d o sería d e m u c h a más fácil aplicación. Sin e m b a r g o , s e g ú n lo visto e n el modelo d e tórax, cuando el m a p a de atenuación no e s uniforme, a u n q u e se modifique el exponente, nunca la mejor imagen se obtiene a la iteración 2. Esto indica que, c u a n d o el m a p a de atenuación es
no
uniforme, no basta c o n una aplicación del método d e Chang para corregir la atenuación y e s preciso iterar.
Capítulo 4
La retroproyección
filtrada
pág. 90