Regularizacion-por-variacion-total-y-sensado-compresivo-en-tomografa

Vista previa del material en texto

UNIVERSIDAD NACIONAL AUTÓNOMA
DE MÉXICO
FACULTAD DE CIENCIAS
Regularización por Variación Total y Sensado
Compresivo en Tomografía
T E S I S
QUE PARA OBTENER EL TÍTULO DE:
Físico
PRESENTA:
Oscar Hurtado González
Director de Tesis:
Dr. Arnulfo Martínez Dávalos
Ciudad Universitaria, Cd. de México, 2018
 
UNAM – Dirección General de Bibliotecas 
Tesis Digitales 
Restricciones de uso 
 
DERECHOS RESERVADOS © 
PROHIBIDA SU REPRODUCCIÓN TOTAL O PARCIAL 
 
Todo el material contenido en esta tesis esta protegido por la Ley Federal 
del Derecho de Autor (LFDA) de los Estados Unidos Mexicanos (México). 
El uso de imágenes, fragmentos de videos, y demás material que sea 
objeto de protección de los derechos de autor, será exclusivamente para 
fines educativos e informativos y deberá citar la fuente donde la obtuvo 
mencionando el autor o autores. Cualquier uso distinto como el lucro, 
reproducción, edición o modificación, será perseguido y sancionado por el 
respectivo titular de los Derechos de Autor. 
 
 
 
1. Datos del Alumno:
Hurtado
González
Oscar
55 19 44 43
Universidad Nacional Autónoma de México
Facultad de Ciencias
Física
309011848
2. Datos del tutor:
Dr.
Arnulfo
Martínez
Dávalos
3. Datos del Sinodal 1:
PhD.
Mercedes
Rodríguez
Villafuerte
4. Datos del Sinodal 2:
Dr.
Rodrigo
Alfonso
Martín
Salas
5. Datos del Sinodal 3:
Dr.
Iván
Miguel
Rosado
Méndez
6. Datos del Sinodal 4:
M. en C.
Edgar
Calva
Coraza
7. Datos del trabajo escrito:
Regularización por Variación Total y Sensado Compresivo en Tomografía
124 p.
2018
2
Agradecimientos
Agradezco inicialmente a la Universidad Nacional Autónoma de México que me
abrió sus puertas como la máxima casa de estudios para poder recibir mi educación
universitaria.
En segundo lugar agradezco a todos los profesores que impartieron su conoci-
miento en las distintas aulas de la Facultad de Ciencias e hicieron de mi carrera un
camino emocionante y desafiante.
Quiero agradecer también al Dr. Arnulfo Martínez Dávalos por haberme permiti-
do realizar mi tesis bajo su tutela y por su infinita paciencia y enseñanza en en estos
últimos años. Asimismo, a mis sinodales, la Dra. Mercedes Rodríguez Villafuerte,
el Dr. Rodrigo Alfonso Matín Salas, al Dr. Iván Miguel Rosado Méndez, y al M.
en C. Edgar Calva Coraza, por su tiempo y disposición para leer y revisar mi tesis,
aportando valiosas observaciones y correcciones para terminar de pulir este trabajo.
Agradezco también al IFUNAM por haberme brindado el espacio necesario para
desarrollar mi proyecto de tesis.
Asimismo, agradezco a los proyectos de la UNAM PAPIIT IN108615, PAPIIT
IN110616, y al proyecto CONACYT Problemas Nacionales 2015-01-612 por los apo-
yos otorgados al Laboratorio de Imágenes Biomédicas del IFUNAM que hicieron
posible la realización de esta tesis.
3
Dedicatoria
A papá y mamá, a quienes amo con todo mi ser. Gracias por todo su amor, apoyo
y cuidado durante toda mi vida y en particular por la paciencia y consejos que me
dieron durante mi carrera universitaria. Gracias papá por ser mi héroe ejemplo para
mi de siempre salir adelante y no rendirme nunca. Gracias mamá por mi heroína y
ser increíblemente diligente al cuidarme y siempre estar pendiente de la familia.
A mis abuelos, Mami y Tito, por siempre estar ahí para mí y consentirme sin
medida. Los amo también. Gracias tío Luis por tu ayuda y cuidado siempre.
Agradezco a mis tíos Victor, Ruth, Rosina, Elsa y a mi abuelo Armando por sus
buenos deseos y cariño. Gracias por los buenos momentos y la comida en todas las
reuniones familiares que durante los periodos de clase, hacían los fines de semana
tiempos de recarga.
A mis mejores amigos, Tony, Isaac Cordova, Edgar, Isaac Barrientos, David y
Sergio. Fueron fundamentales cada uno en distintas etapas de mi vida en los últimos
6 años. Altos y bajos, risas y estrés, desveladas y finales de semestre, comidas, juegos,
domingos de NFL, gimnasio, entrenamientos, béisbol; tantas historias que recordaré
por siempre. Los amo.
A mis modelos a seguir por su ética profesional y su corazón en el deporte, Tom
Brady y Roger Federer; gracias por ser una inspiración y motivación cada temporada
y torneo.
A Dios, gracias.
4
Índice general
1. El problema tomográfico 9
1.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2. Transformada de Radon. . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3. Tomografía discreta . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2. Reconstrucción Tomográfica, Sensado Compresivo y Variación To-
tal 19
2.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2. Reconstrucción iterativa . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.1. Algoritmo ART . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2.2. Algoritmos SIRT . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2.3. Algoritmo MLEM . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2.4. Convergencia de los algoritmos . . . . . . . . . . . . . . . . . 23
2.3. Métodos de regularización . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3.1. Regularización en tomografía . . . . . . . . . . . . . . . . . . 28
2.4. Sensado compresivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4.1. Incoherencia y sensado de señales escasas . . . . . . . . . . . . 30
2.4.2. Fortaleza del sensado compresivo . . . . . . . . . . . . . . . . 33
2.5. Regularización por Variación Total. . . . . . . . . . . . . . . . . . . . 34
3. Desarrollo e implementación de algoritmos 38
3.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.2. Ruido en maniquís y perfiles de intensidad. . . . . . . . . . . . . . . . 39
3.3. Matriz de proyección y sinogramas generados . . . . . . . . . . . . . . 40
3.4. Reconstrucciones con SIRT y MLEM . . . . . . . . . . . . . . . . . . 42
5
3.4.1. Breve análisis de SIRT y MLEM . . . . . . . . . . . . . . . . . 44
3.5. Algoritmos implementados: SIRT+TV, MLEM+TV . . . . . . . . . . 48
3.5.1. Minimización de la variación total: planteamiento del problema 48
3.5.2. Algoritmo SIRT+TV . . . . . . . . . . . . . . . . . . . . . . . 50
3.5.3. Algoritmo MLEM+TV . . . . . . . . . . . . . . . . . . . . . . 51
3.5.4. Variación Total Anisotrópica. . . . . . . . . . . . . . . . . . . 52
3.6. Elección de parámetros . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.6.1. Metodología para elección de parámetros . . . . . . . . . . . . 62
3.6.2. Discusión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.6.3. Parámetros elegidos . . . . . . . . . . . . . . . . . . . . . . . . 69
3.7. Reconstrucciones SIRT+TV . . . . . . . . . . . . . . . . . . . . . . . 70
3.8. Reconstrucciones MLEM+TV . . . . . . . . . . . . . . . . . . . . . . 73
4. Aplicación a datos experimentales 76
4.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.2. Tomografía óptica de luminiscencia estimulada por rayos X . . . . . . 76
4.2.1. Simulación numérica. . . . . . . . . . . . . . . . . . . . . . . . 77
4.2.2. Descripción de los maniquís utilizados . . . . . . . . . . . . . 79
4.2.3. Datos utilizados, imágenes de referencia y elección de parámetros 81
4.2.4. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.2.5. Discusión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
4.3. Tomografía computarizada: Nueces y raíces . . . . . . . . . . . . . . . 112
4.3.1. Raíz de Loto . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.3.2. Nueces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
5. Conclusiones y perspectivas 116
6
Resumen
Las técnicas de reconstrucción tomográfica que se usan comúnmente, se basan en
el concepto de la transformada de Radon (Filtered Back Projection, FBP) ó en la so-
lución de un sistema lineal de ecuaciones que se plantea a partir de la discretización
del problema tomográfico. En ambos casos se tienen limitacionesdebido a varios
factores. El factor más común es la presencia de ruido durante el proceso de adquisi-
ción de datos. Otra limitante es la cantidad de ángulos de proyección que se utilizan
para irradiar el cuerpo o el intervalo angular a través del cual se posiciona el haz de
radiación. Cuando se tiene uno o una combinación de estos factores limitantes, la
calidad de las imágenes reconstruidas disminuye drásticamente en comparación a la
calidad que se obtiene en la ausencia de éstos.
Las técnicas de regularización son utilizadas para resolver problemas que están
“mal planteados” como en estos casos. Por ejemplo en presencia de ruido, estas
técnicas se encargan de disminuir la propagación del ruido durante el proceso de
reconstrucción. Cuando se tienen pocas proyecciones o intervalos angulares restrin-
gidos, las técnicas de regularización pretenden evitar o disminuir la aparición de
“artefactos” en las imágenes.
El objetivo de esta tesis fue estudiar y aplicar la técnica de regularización por
Variación Total (TV) bajo el contexto de la teoría de sensado compresivo a problemas
de reconstrucción tomográfica en presencia de ruido y con pocas proyecciones.
Para esto, la tesis fue estructurada y dividida en 5 capítulos. El primer capí-
tulo introduce los conceptos básicos del problema tomográfico, la transformada de
Radon y el planteamiento algebraico del problema. El capítulo dos introduce los
métodos iterativos para reconstrucción tomográfica con énfasis en SIRT y MLEM
(Simultaneous Iterative Reconstruction Techniques y Maximum Likelihood Expec-
tation Maximization, respectivamente); se abordan las técnicas de reconstrucción en
el contexto de las condiciones de Hadamard para problemas “bien y mal planteados”;
la teoría básica del sensado compresivo y los resultados más relevantes de tal teoría
para esta tesis; y finalmente el capítulo concluye con la técnica de regularización por
Variación Total en el contexto de la teoría del sensado compresivo y su aplicación a
la reconstrucción tomográfica.
El tercer capítulo muestra el efecto del ruido y las pocas proyecciones en recons-
trucciones hechas con SIRT y MLEM y se introducen los criterios que se utilizarán
en el resto de la tesis para evaluar la calidad de las imágenes reconstruidas; des-
pués se llega al corazón de la tesis introduciendo los dos algortimos implementados
que son versiones de SIRT y MLEM pero usando la regularización por variación
7
total: “SIRT+TV” y “MLEM+TV”. Al haber introducido estos algoritmos, el ca-
pítulo continúa con la discusión de los parámetros que utilizan estos algoritmos y
la metodología que se implementó para poder elegir el conjunto de parámetros que
logra reconstruir las imágenes de mejor calidad. Las últimas partes de este capítulo
muestran las reconstrucciones hechas de los maniquís sintéticos Shepp-Logan y el
maniquí Forbild con ruido añadido y utilizando distintos números de proyecciones
empleando la metodología para la elección de parámetros discutida previamente. Los
resultados de estas secciones mostraron que introducir la regularización por varia-
ción total logra controlar de manera evidente la propagación del ruido y la aparición
de artefactos y superan claramente la calidad de las reconstrucciones hechas solo
con SIRT o MLEM.
El capítulo cuatro, muestra la aplicación de los dos algoritmos con regularización
por variación total implementados a dos técnicas experimentales diferentes: tomogra-
fía computarizada (CT) y tomografía óptica de luminiscencia estimulada por rayos
X (TORX). Para la TORX los datos utilizados fueron obtenidos a partir de una
simulación numérica de un arreglo experimental. Se detalla el algoritmo que simula
tal arreglo, se discute brevemente la estructura del maniquí utilizado y las resolucio-
nes de muestreo utilizadas. Después se muestran todos los resultados obtenidos para
las distintas resoluciones; se muestran las imágenes reconstruidas, tablas y gráficas
que comparan la calidad de las reconstrucciones y se observó nuevamente que añadir
la regularización por variación total logra mejorar visiblemente la calidad respecto
a los algoritmos SIRT y MLEM sin regularización. El capítulo termina mostrando
las reconstrucciones de imágenes tomográficas de un tomógrafo de rayos X de la
FIPS (Finnish Inverse Problems Society); los datos son de libre distribución en su
sitio web y son tomografías de una nuez y una raíz de loto. En este caso también se
mostró que la regularización mejora la calidad de las reconstrucciones.
El último capítulo presenta las conclusiones y las perspectivas a futuro que pue-
den ayudar a continuar el estudio abarcado en esta tesis.
8
Capítulo 1
El problema tomográfico
1.1. Introducción
Una de las principales aportaciones de la Física a la Medicina en último siglo es
la obtención de imágenes tomográficas. La palabra tomografía proviene del vocablo
griego tomos que significa rebanada, y del griego graphein que significa escribir
[Suetens, 2009].
Una imagen tomográfica es una representación gráfica de la sección transversal
de un objeto reconstruida a partir de sus proyecciones. Una proyección se refiere
a la información que se obtiene a partir de la energía transmitida de un haz de
rayos X incidente en una dirección particular. Los rayos X que pasan a través del
objeto se atenúan principalmente mediante cuatro mecanismos que son absorción
fotoeléctrica, dispersión Compton, dispersión de Rayleigh y producción de pares. La
presencia de estos fenómenos atenuantes depende de las propiedades materiales del
objeto. Las propiedades atenuantes del objeto se describen con el coeficiente lineal
de atenuación que se representa por medio de una función
µ = µ(~r, E)
en donde ~r es el vector de posición de un elemento material del objeto y E es la
energía del rayo incidente.
Si el haz incidente es monoenergético entonces el coeficiente de atenuación se
reduce a µ = µ(~r), y debido a que las imágenes tomográficas son secciones transver-
sales del objeto, se puede expresar como una función de dos variables µ = µ(x, y)
en un sistema de referencia XY como se muestra en la figura 1.1
9
x
y
µ = µ(x, y)
L
`
x
y
µ = µ(x, y)
I0
I
Figura 1: Top view.
1
(a) (b)
Figura 1.1: (a) Sección transversal de un objeto representada por el coeficiente lineal de
atenuación, también se muestra un haz colimado y la línea de trayectoria. (b) Distintas
geometrías de haz de rayos X en CT, haz en paralelo (izquierda) y haz en abanico (derecha).
En la figura 1.1 también puede observarse un haz colimado de rayos X mono-
energético con intensidad incidente I0 que al atravesar el objeto se atenúa y sale con
una intensidad I. La relación entre estas dos intensidades está dada por la ley de
Beer-Lambert expresada como la integral de línea
I = I0exp
(
−
∫
`
µ(x, y)dl
)
(1.1)
en donde ` es el segmento de línea a lo largo del objeto, que representa la trayectoria
del rayo transmitido.
Dicha integral de línea corresponde a una proyección a lo largo de la línea ` .
Cuando se tiene un haz de rayos X se tienen varias proyecciones determinadas por
el conjunto de trayectorias de los rayos del haz que cubren el campo de vista (FOV,
del inglés field of view). En tomografía computarizada (CT, del inglés computed to-
mography) se usan principalmente dos tipos de geometrías para los haces utilizados,
el haz de rayos paralelos y haz de rayos en abanico (fan-beam) como se observa en la
figura 1.1(b). El haz de rayos X se rota al rededor del objeto para medir proyecciones
a muchos ángulos. Para obtener las proyecciones del haz de rayos X se utiliza un
sistema de adquisición de datos que incluye el uso de fotodiodos o contadores de
fotones para medir la intensidad de los rayos X transmitidos.
Las bases matemáticas para reconstruir una sección transversal o de manera más
precisa, la obtención de la función µ(x, y) a partir de sus proyecciones, fueron desa-
rrolladas por Johann Radon en 1917. Sin embargo el primer escáner para reconstruir
imágenes tomográficas fue desarrollado por Godfrey N.Hounsfield en 1972 y se basó
en algunos algoritmos desarrollados por Allan Cormack. Ambos recibieron el premio
Nobel en medicina en 1979.
10
1.2. Transformada de Radon.
Para estudiar la teoría matemática desarrollada por Radon en 1917 [Radon,
1986], considérese la geometría de haz de rayos paralelos mostrada en la figura 1.2.
I
FO
V
r
s y
✓
I✓(r)
I✓(r)
r
p✓(r)
r
),
y)
µ(x, y)
x
I0
Figura 2: Top view.
2
Figura 1.2: Geometría de rayos paralelos (Imagen tomada de [Suetens, 2009]).
Los rayos X forman un ángulo θ respecto al eje y y con este ángulo se define
una rotación de coordenadas. El sistema de coordenadas RS permite calcular las
integrales para cada rayo del haz. Los rayos son paralelos al eje s y se supone que
están equiespaciados respecto a la coordenada r. La intensidad de los rayos incidentes
es I y se puede calcular la atenuación de cada uno en función de la coordenada r
usando la integral de línea dada por la ecuación (1.1). Para ello se utilizan las
matrices de rotación dadas por (1.2).[
x
y
]
=
[
cos θ −sin θ
sin θ cos θ
] [
r
s
]
[
r
s
]
=
[
cos θ sin θ
−sin θ cos θ
] [
x
y
] (1.2)
Así para un ángulo fijo θ, se puede calcular la intesidad atenuada de cada uno de
los rayos del haz en función de su posición en la coordenada r. Utilizando (1.1) se
obtiene el perfil de intensidad como función de r, Iθ(r)
Iθ(r) = I exp
(∫
`r,θ
µ(r · cos θ − s · sin θ, r · sin θ + s · cos θ)ds
)
(1.3)
En donde la recta `r,θ es la recta que representa la trayectoria del haz a distancia r
del origen del sistema de coordenadas incidente a un ángulo θ. Tomando el logaritmo
11
se define lo que es el perfil de atenuación, pθ(r)
pθ(r) = −ln
Iθ(r)
I
=
(∫
`r,θ
µ(r · cos θ − s · sin θ, r · sin θ + s · cos θ)ds
)
(1.4)
Así, el perfil de atenuación corresponde a la proyección del coeficiente de atenua-
ción a lo largo de `r,θ.
Se pueden tomar distintos ángulos θ en el intervalo [0, 2π] para obtener más pro-
yecciones y construir una función P (r, θ), llamada sinograma. La figura (1.3) muestra
el maniquí Shepp-Logan, el cual es una imagen de prueba estandarizada creada por
Larry Shepp y Benjamin Logan [Shepp and Logan, 1974], y sirve como modelo de
una cabeza humana en el desarrollo y prueba de algoritmos de reconstrucción de
imágenes . La transformación de la función µ(x, y) a su sinograma corresponde a la
transformada de Radon.
R{µ}(r, θ) ≡ P (r, θ) =
(∫
`r,θ
µ(r · cos θ − s · sin θ, r · sin θ + s · cos θ)ds
)
(1.5)
50 100 150 200 250
50
100
150
200
250
0
10
20
30
40
50
60
70
80
90
100
(a)
50 100 150
θ
50
100
150
200
250
r
0
10
20
30
40
50
60
70
80
90
(b)
Figura 1.3: (a) Maniquí Shepp-Logan; se muestran las dimensiones en pixeles en los ejes y
la escala de grises correspondiente. (b) Transformada de Radon o sinograma del maniquí
Shepp-Logan.
Obtener proyecciones de manera experimental es medir las intensidades ate-
nuadas del haz de rayos X a distintos ángulos y a diferentes posiciones, es decir,
experimentalmente se mide el sinograma o transformada de Radon. Así la recons-
trucción de una imagen tomográfica es obtener el coeficiente de atenuación µ(x, y)
a partir de su sinograma R{µ}(r, θ), es decir que reconstruir la sección transversal
del objeto equivale a resolver el problema inverso a la transformada de Radon. Por
lo tanto, resolver el problema tomográfico consiste en encontrar (si es que existe) la
“transformada inversa de Radon”.
µ(x, y) = R−1 (P (r, θ))
12
La aproximación más simple a la transformada inversa de Radon consiste integrar
con respecto a θ el sinograma P (r, θ) y asignar los valores a las coordenadas (x, y)
en función de la coordenada r utilizando el cambio de variable r = xcos θ + ysin θ
dada por (1.2)
BP (x, y) =
∫ π
0
P (r, θ)dθ =
∫ π
0
P (xcos θ + ysin θ, θ)dθ (1.6)
La ecuación (1.6) se conoce como retroproyección (Backprojection). La retropro-
yección es una primera aproximación que no sirve bien pues con ella se obtienen
imágenes borrosas como se observa en la figura 1.4(a).
50 100 150 200 250
50
100
150
200
250 25
30
35
40
45
50
55
60
65
70
(a)
50 100 150 200 250
50
100
150
200
250
0
0.2
0.4
0.6
0.8
1
(b)
Figura 1.4: (a) Retroproyección. Implementada con la función iradon de MATLAB usando
el filtro ’none’, es decir la ec. (1.6) sin filtro.(b) Filtered Back Projection (FBP), imple-
mentada con la función iradon de MATLAB.
La transformada inversaR−1 se obtiene por medio del teorema de la rebanada
central . Este teorema relaciona la transformada de Radon de la función µ(x, y) con
las transformadas de Fourier. El teorema establece las siguientes ecuaciones y una
demostración puede verse en [Suetens, 2009].
Sea F{µ}(kx, ky) la transformada en 2D de la función µ(x, y).
F{µ}(kx, ky) =
∫ ∞
−∞
∫ ∞
−∞
µ(x, y)e−2πi(kxx+kyy)dxdy
y sea Pθ(k) la transformada de Fourier de pθ(r).
Pθ(k) =
∫ ∞
−∞
pθ(r)e
−2πikrdr
Haciendo θ variable, el teorema afirma la siguiente igualdad:
P(k, θ) = F{µ}(kx, ky); ∴ F−1{P(r, θ)} = µ(x, y) (1.7)
en donde
kx = kcos θ
ky = ksin θ
k =
√
k2x + k
2
y
(1.8)
13
Entonces si se mide el sinograma P (r, θ) y se calcula su transformada de Fourier en
1D en la variable r utilizando el cambio de variable dado por (1.8) y la igualdad
dada por (1.7), se puede calcular la transformada inversa de Fourier en 2D y así
obtener µ(x, y)
F−1{P(r, θ)} = µ(x, y)
Para calcular la transformada inversa de Fourier de la ecuación (1.7) se necesita usar
la versión polar de la transformada inversa de Fourier en 2D dada por
µ(x, y) =
∫ π
0
∫ ∞
−∞
P(k, θ)|k|ei2πkrdkdθ, (1.9)
con r = xcos θ + ysin θ. Definiendo
P∗(k, θ) = P(k, θ)|k|
y
P ∗(k, θ) =
∫ ∞
−∞
P∗(k, θ)ei2πkrdk
la ecuación (1.9) se convierte en
µ(x, y) =
∫ π
0
P ∗(k, θ)dθ. (1.10)
A la ecuación (1.10) se le llama retroproyección filtrada (FBP del inglés filtered
back projection). Esta ecuación es como la retroproyección dada por 1.6 pero se
introduce un filtro en el dominio de frecuencias (es decir, se multiplica la función a
retroproyectar por el filtro “ramp” |k|, en el dominio de Fourier).
La retroproyección filtrada satisface el teorema de la rebanada central y por ende
puede ser identificada como la transformada inversa de Radón. Un ejemplo de FBP
puede observarse en la figura 1.4 (b).
1.3. Tomografía discreta
Todo lo discutido en la sección anterior ha supuesto implícitamente el uso de
variables continuas para el uso de las transformadas integrales presentadas. Sin
embargo en la práctica debido a la precisión finita de los números, la aritmética
de las computadoras, y a que tampoco se mide un sinograma continuo sino un
conjunto finito de puntos P (rn, θm) (con n ym en un conjunto de índices), se necesita
discretizar el sinograma por una malla de n×m puntos. Tal discretización puede ser
definda por la geometría de la instrumentación utilizada para obtener las mediciones.
La geometría de la instrumentación incluye el tipo de haz (paralelo o de abanico), la
cantidad de ángulos θm en los que se hacen incidir el haz de rayos X, el número de
detectores utilizados en el arreglo que mide la intensidad de los rayos transmitidos
(esto define el número de rayos que tiene el haz y por tanto define las coordenadas
rn) y la resolución de la imagen que se desea reconstruir (número de pixeles que
tendrá el tamaño de las imágenes). Si tal discretización es lo suficientemente fina,
con una cantidad grande de rayos y muchas proyecciones entonces comúnmente se
14
integran tales mediciones con la versión discreta de la FBP (1.10). Por simplicidad,
la siguiente deducción ilustra el planteamiento del problema para un haz paralelo
de rayos X.
Definir la resolución de la imagen o el número de pixeles equivale a discretizar el
coeficiente de atenuación µ(x, y) en una malla de N ×N elementos pero enumerán-
dolos con un solo índice j = 1, 2, ...N2, (µj) como se muestra en la figura 1.5(a). Esta
discretización está definida por los N detectores utilizados en la instrumentación y
denotados por las letras b1,b2, ..., bN . Estos detectores también definen el número de
rayos que forman el haz, que de ahora en adelante se denotará con Nrays, es decir
se tiene que Nrays = N , por lo que estas letras se utilizarán de manera indistinta en
lo siguiente. Esta figura corresponde al ángulo θ1 = 0. Así, las intensidades medi-
das por los detectores son I1, I2, ..., IN para el ángulo θ1. La intensidad de los rayos
incidentes es I0 para todo rayo en el haz.
Al tener un conjunto discreto de coeficientes de atenuación las integrales de línea
pasan a ser sumas y la transformada de Radon se puede ver como un sistema de
ecuaciones lineales como se muestra a continuación [Kak and Slaney, 2001].
Para el rayo con intensidad atenuada Il con l = 1, 2, ...N se tiene que
Il = I0exp
− l·N∑
i=[(l−1)·N ]+1
d · µi
 .
en donde d es la longitud de un pixel como se observa en la figura 1.5(a). Debido a
que se conocen las intensidades I0 e I1, se puede denotar a la suma con la misma
etiqueta con la que se denotó a los detectores del arreglo. Despejando la suma se
puede hacer
bl ≡ −ln
(
Il
I0
)
=
l·N∑
i=(l−1)·N+1
d · µi (1.11)
Para obtener más proyecciones se necesita hacer un barrido a distintos ángulos
en el intervalo [0, π], pues por simetría los valores de las sumas obtenidas para
proyecciones con ángulos en [π, 2π] son iguales a las sumas en el intervalo anterior.
Así, tomando un conjunto finito de Nθ ángulos en el intervalo [0, π], se pueden medir
más proyecciones como se muestra en la figura 1.5(b).
Por la simetría mencionada (θ1 = θNθ , para los ángulos θk con k = 1, ..., Nθ−1, se
pueden medir más proyecciones y obtener más sumas como las de la ecuación (1.11).
Sin embargo hay que hacer algunas modificaciones pues al cambiar el ángulo, la
distancia que recorre un rayo ya no es igual a d. Por ejemplo, como se puede observar
en la figura para el último rayo con intensidad atenuada IN(k+1), la distancia que
recorre al atravesar el primer pixel que toca es di,j. Las nuevas sumas bkN+l se
enumeran de forma consecutiva a las sumas correspondientes al ángulo anterior. Por
ejemplo la última suma para θ1 corresponde a bN , por lo que las sumas para θ2 se
enumerarán de la forma bN+1, bN+2, ..., b2N .
La suma bkN+l esta dada por la siguiente ecuación. Haciendo i = kN + l
bi =
∑
j∈IN
di,j · µj, (1.12)
15
µ1
µ3
µ2
µN
µ
N+1
µ
N+2
µ
N+3
µ
2N
µj�1
µj
µj+1
µ
N2
I 0
I 0
I 0
I 1
I 2
I N
Ar
reg
lo
de
de
tec
tor
es
µ1
µ3
µ2
µN
µ
N+1
µ
N+2
µ
N+3
µ
2N
µj�1
µj
µj+1
µ
N2
I0 I0 I0
I1 I2 IN
Arreglo de detectores
✓k
b1 b2 bN
b kN
+
1
b kN
+
2
b N(
k+
1)
} d
}
d
Figura 3: Top view.
3
(a)
µ1
µ3
µ2
µN
µ
N+1
µ
N+2
µ
N+3
µ
2N
µj�1
µj
µj+1
µ
N2
I 0
I 0
I 0
I kN
+
1
I kN
+
2
I N(
k+
1)
Ar
reg
lo
de
de
tec
tor
es
d
i,j
µ1
µ3
µ2
µN
µ
N+1
µ
N+2
µ
N+3
µ
2N
µj�1
µj
µj+1
µ
N2
I0 I0 I0
I1 I2 IN
Arreglo de detectores
✓k
b1 b2 bN
b kN
+
1
b kN
+
2
b N(
k+
1)
} d
}
d
}
Figura 3: Top view.
3
(b)
Figura 1.5: (a)Discretización de la función µ(x, y). Las proyecciones se miden con el arreglo
de detectores bl. (b)Rotación a un ángulo θk de el haz de rayos paralelos. Las proyecciones
se miden con el arreglo de detectores en la nueva posición definida por ángulo dado.
16
en donde IN es un conjunto de índices que corresponde a los subíndices de los pixeles
que atraviesa el rayo que llega al el i-ésimo detector. Entonces teniendo Nrays y Nθ
ángulos se tienen en totalM = Nrays ×Nθ proyecciones. Es decir podemos escribir
~b ≡ (b1, b2, ..., bM),
de igual forma
~x ≡ (µ1, µ2, ..., µN2 ).
Viendo a estos vectores como vectores columna, podemos expresar todas las sumas
como un sistema de ecuaciones de la forma
A~x = ~b, (1.13)
en donde las entradas de la matriz A ∈ M
M×N2
corresponden a las distancias que
atraviesan los rayos para contribuir a alguna proyección. Anteriormente se utilizo la
letra di,j para denotar a tales distancias pero de ahora en adelante a los coeficientes
de la matriz A simplemente se denotarán por ai.j. A la matriz A se le llama matriz
de proyección o simplemente matriz del sistema. Los coeficientes ai,j también se
pueden interpretar como la contribución del j-ésimo pixel µj a la i-ésima medida bi.
Si se normalizan las distancias entonces estos coeficientes pueden tomar valores entre
0 y 1 y también se puede interpretar como la probabilidad de que el j-ésimo pixel
contribuya a la i-ésima proyección (interpretación que es útil, por ejemplo, para las
imágenes obtenidas con ténicas PET y SPECT).
De la figura 1.5(b) se puede apreciar que calcular los coeficientes ai,j es un asunto
geométrico que tiene que tomar en cuenta las dimensiones de los pixeles, la distancia
de separación de los rayos y los ángulos para los cuales se miden las proyecciones.
Calcular tales coeficientes puede ser muy complicado y largo de resolver y exis-
ten distintos algoritmos que permiten calcular la matriz de proyección. Uno de los
algoritmos más utilizados es el algoritmo de Siddon [Siddon, 1985].
La idea básica de tal algoritmo es considerar a los pixeles como las áreas de
intersección entre conjuntos ortogonales de líneas paralelas equiespaciadas. Se para-
metrizan cada una de las líneas que representan las trayectorias de los rayos proyec-
tados con una ecuación de la forma ~L(α) = ~p1 + α(~p2 − ~p1), con ~p1 y ~p2 siendo dos
puntos que definen cada una de las trayectorias y α ∈ [0, 1]. Entonces se calculan los
puntos de intersección con cada una de las líneas que definen el conjunto de pixeles
para cada rayo. El algoritmo trabaja de forma inicial para cada trayectoria calcu-
lando solo la intersección de “entrada” del rayo con la primera línea del conjunto y
la intersección de “salida” con la última línea del conjunto y, utilizando la geometría
ortogonal y equiespaciada, calcula las demás intersecciones intermedias de forma
recursiva. La longitud ai,j buscada está dada en términos del parámetro α con el
que se parametrizan las trayectorias. Para una explicación detallada del algoritmo
veáse [Siddon, 1985]. Mejoras al algoritmo para hacer el calculo más rápido se han
propuesto por ejemplo en [Jacobs et al., 1998]. En esta tesis se utilizaron los algo-
ritmos implementados en [Hansen and Saxild-Hansen, 2012] incluídos en el paquete
AIR TOOLS para MATLAB.
Así, el problema tomográfico se traduce a resolver el sistema lineal de ecuacio-
nes (1.13) para encontrar los valores de los coeficientes lineales de atenuación que
17
corresponden a las entradas del vector ~x. Es decir, dadas las proyecciones ~b (sino-
grama) medidas con el arreglo de detectores que rodean el objeto y la geometría
de la instrumentación que define la matriz de proyección A. En el capítulo 2, se
describen los métodos numéricos que se emplean para resolver el sistema planteado
en la ecuación 1.13.
18
Capítulo 2
Reconstrucción Tomográfica, Sensado
Compresivo y Variación Total
2.1. Introducción
En este capítulo se introducen los métodos iterativos para la reconstrucción to-
mográfica, los métodos de regularización, la teoría básica del Sensado Compresivo
y finalmente la Regularización por Variación Total. En la primera sección de este
capítulo se detallan los métodos iterativos de reconstrucción tomográfica y algunas
de sus características principales. En la segunda sección, se detallan los métodos de
regularización en el contexto de las condiciones de Hadamard para sistemas de ecua-
ciones lineales. Después, en la tercera sección, se introduce de manera condensada
la teoría de Sensado Compresivo resumiendo los resultados más relevantes para el
objetivo de esta tesis. Finalmente en la cuarta sección se habla de la Regularización
por Variación Total, su relación con la teoría del sensado compresivo y su aplicación
a la reconstrucción tomográfica.
2.2. Reconstrucción iterativa
Si se decide resolver el sistema de ecuaciones en vez de integrar el sinograma
con la FBP, se tienen que tomar en cuenta varias propiedades de los sistemas de
ecuaciones lineales para elegir un método de solución.Por ejemplo, encontrar las
soluciones a (1.13) resulta ser computacionalmente difícil si las dimensiones de la
matriz A son muy grandes y por ello los métodos típicos de solución como eliminación
Gaussiana o encontrar la matriz inversa A−1 suelen ser muy costosos en cuanto a
tiempo de cómputo. También se puede tener que el sistema esté subdeterminado,
es decir tener menos ecuaciones que incógnitas ( N2 >M), en cuyo caso la matriz
inversa no existe y por tanto tener un subespacio solución de donde se puede elegir
alguna solución ~x, por lo que no se llega a una solución única. Una alternativa para
resolver estos sistemas lineales de ecuaciones son los métodos iterativos.
Los métodos iterativos son algoritmos que, a partir de una estimación inicial
19
~x0, aplican una fórmula recursiva de la forma ~xk+1 = f(~xk) para generar una suce-
sión de vectores ~x1, ~x2, ..., ~xk, ... que son soluciones aproximadas que convergen a la
solución buscada ~x. Es decir ~xk → ~x cuando k →∞ [Loehr, 2014a].
La ventaja de los algoritmos iterativos es que con las computadoras actuales se
pueden escribir rutinas en diferentes lenguajes de programación de manera relativa-
mente sencilla y debido a la velocidad de procesamiento, estos algoritmos resultan
ser muy rápidos en cuanto a tiempo de cómputo en comparación con los métodos de
eliminación Gaussiana o descomposición de matrices como se menciona en [Beister
et al., 2012].
Existen varios tipos de algoritmos iterativos que suelen aparecer en los libros
de álgebra lineal numérica (véase [Loehr, 2014b] capítulo 10) como el algoritmo
de Richardson, el algoritmo de Jacobi y uno de los más conocidos es el algoritmo
Gauss-Seidel. Las condiciones de convergencia de estos métodos han sido estudiadas
de manera detallada y en esa misma referencia pueden verse las condiciones para la
convergencia del algoritmo Gauss-Seidel.
Por ejemplo, el algoritmo de Richarson está dado por la ecuación
~xk+1 = ~xk + (~b− A~xk) (2.1)
Los algoritmos de interés en el desarrollo de esta tesis que se implementan común-
mente en la reconstrucción de tomografías se describen brevemente a continuación.
2.2.1. Algoritmo ART
En el caso de tomografía computarizada, el primer algoritmo que se implementó
para resolver el sistema (1.13), es el algoritmo ART (del inglés Algebraic reconstruc-
tion techniques). El algoritmo ART fue introducido en [Gordon et al., 1970] y se
basa esencialmente en el algoritmo de Kaczmarz [Kaczmarz, 1937].
El algoritmo se puede escribir de la siguiente forma:
~x+ λk(bi − 〈ai, ~x〉)
ai
‖ai‖22
→ ~x, (2.2)
en donde a λk se le conoce como parámetro de relajación, ‖ · ‖ representa la
norma `2 y ai es el i-ésimo renglón de la matriz A (que es un vector de la misma
dimensión de ~x). En este algoritmo la k-ésima iteración ~xk consiste en un “barrido”
a través de las columnas de A, es decir que para obtener cada iteración ~xk, se debe
seguir la regla dada por la la ecuación (2.2) variando i desde 1 aM.
En pseudocódigo se vería de la siguiente forma
x_0=aproximacion inicial
for k=1,2,...,iterNUM
--->for i=1,2,...M
--->x_k=x_{k-1}+Lam_k (b_i- a_i*x_{k-1})(a_i/||a_i||^2)
--->return x_k
--->end
20
end
x_final=x_iterNUM
El número M corresponde aM y el iterNUM es el número de veces que se aplicará
de forma recursiva la ecuación (2.2) para obtener la secuencia ~x1, ...~xiterNUM . x* es
la solución aproximada a (1.13) ya que en la práctica no se puede hacer k → ∞.
también es común dejar el parámetro de relajación fijo para toda iteración, es decir
se hace λk = k y Kaczmarz utiliza λ = 1. Una interpretación geométrica de la
convergencia del algoritmo puede estudiarse en [Kak and Slaney, 2001] capítulo 7.
2.2.2. Algoritmos SIRT
Los algoritmos SIRT (del inglés Simultaneous Iterative Reconstrution Techni-
ques) son mejoras al algoritmo ART, pues como se discute en [Gilbert, 1972] y
en [Andersen and Kak, 1984], el algoritmo ART presenta algunas limitaciones da-
das por la naturaleza de las hipótesis que se utilizan en [Gordon et al., 1970]. Los
algoritmos SIRT llevan el nombre de “simultáneos” porque a diferencia de ART, en
donde se tiene que iterar renglón por renglón de la matriz A, todas las ecuaciones
son utilizadas al mismo tiempo en una sola iteración por medio de operaciones entre
matrices y vectores. La forma general de estos algoritmos es
~xk+1 = ~xk + λkTA
TM(~b− A~xk), (2.3)
en donde λk es nuevamente el parámetro de relajación, M y T son matrices que
determinan el algortimo en particular que se utiliza. La característica principal de
las matrices M y T es que son matrices simétricas positivas definidas. En [Hansen
and Saxild-Hansen, 2012] se mencionan algunos tipos de algoritmos SIRT como
son el método de Landweber, el método de Cimmino y SART (del inglés
Simultaneous Algebraic Reconstruction Technique).
Para el interés de esta tesis en particular, se utilizó el algoritmo SART (se eligió
tal algoritmo debido a que es el más comúnmente implementado en las referencias
consultadas durante el desarrollo de este trabajo) que está dado por
~xk+1 = ~xk + λkCA
TR(~b− A~xk), (2.4)
en donde las matrices C ∈ M
N2×N2
y R ∈ MM×M son matrices diagonales cuyos
elementos están dados por
cj,j =
1∑
i ai,j
; ri,i =
1∑
j ai,j
En diversa literatura se intercambian de manera indistinta los términos SIRT y
SART para denotar al algoritmo SART. De aquí en adelante para el desarrollo de
esta tesis se utilizará solo el término SIRT para denotar a la ecuación (2.4) a menos
que se quiera hacer una distinción explícita entre algoritmos.
21
2.2.3. Algoritmo MLEM
El alrgoritmo MLEM (del inglés Maximum Likelihood Expectation Maximiza-
tion) es aplicado para tomografía por transmisión y para tomografía por emisión
en [Lange and Carson, 1984], y es un algoritmo que se basa en la naturaleza estadís-
tica de los fotones que son emitidos y detectados en la medición de las proyecciones.
Debido a que la estadística de los fotones es distinta en ambas modalidades (emi-
sión ó transmisión), la deducción de este algoritmo cambia un poco entre las dos
modalidades pero suele usarse en primera aproximación para ambas, la deducción
que considera la estadística de emisión. La breve discusión que aquí se presenta es
una simplificación del artículo de Lange y Carson como la presenta [Bruyant, 2002].
Para la estadística de emisión de fotones (como en el caso de Tomografía por
Emisión de Positrones, PET), el vector ~x que aparece en 1.13 representa las desinte-
graciones radiactivas (PET) o emisiones que ocurren en el cuerpo que está emitiendo
los fotones, así el elemento xi representa las desintegraciones que ocurren en el pixel
j-ésimo de la imagen. Los elementos de matriz ai,j representan la probabilidad de
que el i-ésimo detector mida un fotón emitido por el pixel j y ~b es el sinograma me-
dido. El objetivo del algoritmo MLEM es encontrar una solución “general” como el
mejor candidato estimado a ~x: el número promedio de emisiones o desintegraciones
radiactivas ~x en la imagen que pueden producir el sinograma dado con la mayor
probabilidad . Esto se logra utilizando la ley de Poisson que permite predecir la
probabilidad de un número determinado de conteo de fotones a partir del número
promedio de desintegraciones. Cada paso del algoritmo consiste en dos pasos básicos:
el paso de expectación (paso E), en donde se calcula la probabilidad de que cualquier
imagen sea reconstruida en función a las mediciones, y en el paso de maximización
(paso M), la imagen que tiene la mayor probabilidad de dar los datos medidos es
encontrada. El algoritmo se reduce en aplicar la fórmula recursiva dada por 2.5.
x
(k+1)
j =
x
(k)
j∑M
i=1Ai,j
M∑
i=1
bi∑N2
j′=1Ai,j′x
(k)
j′
Ai,j (2.5)
Para deducir esta ecuación considérese el número promedio de desintegraciones
en el pixel j, xj. El número promedio de fotones emitidos del pixel j y detectados
por el detector i, es ai,jxj. El número promedio de fotones detectados por el detector
i es la suma de los números promedios de fotones emitidos por cada pixel:
bi =
m∑
j=iai,jxj.
El número de fotones emitidos de los m pixeles y detectados por el i-ésimo detector
es una variable de Poisson. Así, la probabilidad de detectar bi fotones está dada por
P (bi) =
e
−bi b
bi
i
bi!
Todas las variables de Poisson se consideran estadísticamente independientes y la
probabilidad condicional P (b|x) de observar el vector ~x cuando el mapa de emisión
22
es ~x es el producto de las probabilidades individuales P (bi). Esto da la función de
probabilidad L(~x):
L(~x) = P (b|x) =
n∏
i=1
P (bi) =
n∏
i=1
e
−bi b
bi
i
bi!
.
Para hallar el máximo valor para la probabilidad L(~x) se utiliza la derivada. Para
maximizar la expectación (paso E), usualmente se considera l(~x) = ln(L(~x)). Y
utilizando las propiedades del logaritmo natural, la ecuación anterior se convierte
en
l(~x) =
n∑
i=1
(
−bi + biln(bi)− ln(bi!)
)
, (2.6)
sustituyendo:
l(~x) =
n∑
i=1
(
−
m∑
j=i
ai,jxj + biln(
m∑
j=i
ai,jxj)− ln(bi!)
)
. (2.7)
Esta ecuación conocida como función de probabilidad, es fundamental para el al-
goritmo MLEM pues permite calcular la probabilidad de observar un conjunto de
datos de proyecciones a partir de cualquier imagen ~x. Y el paso de maximización
(paso M) es para encontrar la imagen con mayor probabilidad que es considerada
la solución óptima. Se ha demostrado en trabajos posteriores a [Lange and Carson,
1984], que l(~x) tiene uno y solo un máximo el cual se obtiene al derivar e igualar a
cero:
∂l(~x)
∂~x
= −
n∑
i=1
ai,j +
n∑
i=1
bi∑m
j′=1 ai,j′xj′
ai,j = 0. (2.8)
Multiplicando por xj y despejando se llega a la ecuación que da el esquema iterativo
de la ecuación 2.5:
xj =
xj∑n
i=1 ai,j
n∑
i=1
bi∑m
j′=1 ai,j′xj′
ai,j. (2.9)
2.2.4. Convergencia de los algoritmos
Una de las preguntas que surge naturalmente al trabajar con métodos iterativos
es cuántas iteraciones se deben hacer para que la solución k-ésima sea una “buena”
aproximación a la solución exacta ~x. Un criterio comúnmente utilizado es utilizar
niveles de tolerancia tales como
‖~xk+1 − ~xk‖2 < �x,
en donde �x define un nivel de tolerancia para la norma de la diferencia entre itera-
ciones consecutivas o bien definiendo el término de fidelidad como ‖A~xk −~b‖2,
entonces
‖A~xk −~b‖2 < �f , (2.10)
siendo �f la tolerancia para el término de fidelidad.
23
Otro criterio utilizado involucra a alguna imagen de referencia para comparar
la imagen reconstruida. Si ~xref representa la imagen de referencia vectorizada, se
puede definir el error cuadrático medio (MSE, del inglés mean squared error)
como:
MSEk =
1
MN
‖~xref − ~xk‖22, (2.11)
en donde MN es la dimension de la imagen (M ×N pixles) y el subíndice k denota
el error cuadrático medio de la k-ésima iteración. Luego, puede usarse esta medida
y así automatizar el número de iteraciones hasta que se satisfaga algún nivel de
tolerancia dado.
MSEk < �mse (2.12)
Existen teoremas que demuestran que el algoritmo ART converge a la solución bus-
cada ~x (veáse [Eggermont et al., 1981]); una de las condiciones de convergencia es
que λk ∈ (0, 2). También diversos autores han estudiado la convergencia de los al-
goritmos SIRT ( [Censor and Elfving, 2002], [Censor et al., 2008], [Jiang and Wang,
2003a], [Jiang and Wang, 2003b], [Qu et al., 2009], [Van der Sluis and Van der Vorst,
1990]); y en [Hansen and Saxild-Hansen, 2012] se hace referencia a un teorema que
resume los resultados de estos autores y que es muy práctico para los intereses de
esta tesis.
dTeorema 1 : Los iterados ~xk de (2.3) convergen a una solución ~x∗ de
argmin
x
‖A~x−~b‖2M
si y sólo si
0 < � ≤ λk ≤
2
ρs(TATMA)
− �, (2.13)
en donde ‖A~x −~b‖2M = (A~x − b)TM(A~x − b), y � es una constante fija pequeña y
arbitraria. Si además ~x0 está en el rango (o imagen) de TAT entonces la solución ~x∗
es única con norma ‖ · ‖−1T . (Aquí ρs(Q) representa el radio espectral de una matriz
Q definido como el eigenvalor positivo más grande).c
Así, los algortimos SIRT resuleven el problema de optimización conocido como
mínimos cuadrados ponderados (en inglés weighted least squares). En el caso
particular de SART, T = C yM = R resulta ser que ρs(CATRA) = 1 y así, tomando
� = 0 la convergencia está garantizada siempre y cuando el parámetro de relajación
satisfaga
0 < λk ≤ 2. (2.14)
Esta condición de convergencia es importante y fue explotada fuertemente en la
implementación de algoritmos y en el análisis de resultados de esta tesis como se
explicará en capítulos siguientes.
Por otro lado la convergencia de MLEM es estudiada por [Kaufman, 1993]. En
su artículo Kaufman establece que MLEM es equivalente a aplicar el proceso de
expectación y maximización (EM) a un problema de mínimos cuadrados ponderados
parecido al de la ecuación del teorema 1. Así que MLEM también minimiza la
distancia entre A~x y ~b.
24
2.3. Métodos de regularización
Los algoritmos algebraicos presentados en la sección anterior resuelven el sistema
(1.13) logrando reconstruir imágenes tomográficas de buena calidad bajo ciertas con-
diciones. Una característica importante que se debe mencionar acerca de los sistemas
de ecuaciones lineales como éste, es la de si el problema está o no “bien planteado”
(en inglés well-posed problem). Por definición, un sistema está “bien planteado”
si satisface las condiciones de Hadamard [Bakushinskiy and Goncharsky, 1994]
dadas por:
1. La solución existe.
2. La solución es única.
3. La solución es estable con respecto a pequeñas perturbaciones en ~b, en el
sentido de que si A~x0 = ~b0 y A~x = ~b, entonces ~x→ ~x0 cuando ~b→ ~b0.
Equivalentemente el sistema (1.13) está “bien planteado” si la matriz del sistema
A es biyectiva (invertible) y el operador inverso A−1 es estable en el sentido del punto
3. Si una de las condiciones de Hadamard no se satisface se dice que el sistema está
“mal planteado” (en inglés ill-posed problem).
La primera condición puede no satisfacerse si, por ejemplo, los datos medidos en~b
contienen algún tipo de ruido. La unicidad de la solución para la segunda condición,
depende del número de ecuaciones e incógnitas y en los problemas tomográficos
generalmente el sistema está subdeterminado. La tercera condición casi siempre se
viola en la presencia de ruido pues éste se propaga rápidamente dando una solución
alejada de la deseada [Wang et al., 2011].
Cuando se tiene un problema “mal planteado” se sigue una técnica conocida como
regularización que consiste en imponer condiciones adicionales para restringir al
problema y de alguna manera estabilizarlo.
Una manera no tan rigurosa de ver cómo la tercera condición de Hadamard
produce problemas mal planteados ante la presencia de ruido y lo que son las técnicas
de regularización, es estudiar más a detalle el sistema (1.13) con la teoría espectral
de matrices simétricas.
Dado que A ∈ M
M×N2
no es necesariamente una matriz cuadrada entonces po-
demos extender el sistema con la matriz AT
ATA~x = AT~b.
Definiendo As ≡ ATA y ~y ≡ AT~b, reescribimos el sistema como
As~x = ~y.
As ∈ MN2×N2 es una matriz positiva definida (es decir, ~r · (As~r ) > 0 para todo
vector ~r 6= ~0 ∈ RN
2
) y por tanto existen N2 eigenvalores 0 ≤ λ1 ≤ λ2 ≤ ... ≤ λN2
25
y los correspondientes eigenvectores ui que forman una base ortonormal en la cual
As tiene una representación diagonal de la forma [Doicu et al., 2010]:
As =
N2∑
i=1
λiuiu
T
i .
El número de condición está definido por el cociente entre los eigenvalores con
mayor y menor valor. Suponiendo λ1 6= 0, el número de condición de As está dado
por κ ≡ λN2
λ1
. Para simplificar la notación se suele suponer que los eigenvalores están
escalados y que λN2 = 1. Así, el número de condición está dado por el eigenvalor
mas pequeño κ = 1
λ1
.
El número de condición es una medida de la estabilidad del sistema ante pequeñas
variaciones en los datos. Supongamos que ~yδ representa un conjunto de datos con
ruido que satisface
‖~yδ − ~y‖2 < δ (2.15)
Sea ~xδ la solución que satisface As~xδ = ~yδ. Entonces utilizando la representación
espectral se puede escribir
~xδ − ~x =
N
2∑
i=1
1
λi
uiu
Ti · (~yδ − ~y).
Luego, utilizando la ortogonalidad de la base se puede medir la norma de la expresión
anterior
‖~xδ−~x‖22 =
N
2∑
i=1
1
λ2i
[uTi ·(~yδ−~y)]2 ≤
N
2∑
i=1
1
λ21
[uTi ·(~yδ−~y)]2 =
1
λ21
N
2∑
i=1
[uTi ·(~yδ−~y)]2 =
1
λ21
‖~yδ−~y‖22,
en donde la desigualdad se cumple ya que, 0 ≤ λ1 ≤ λ2 ≤ ... ≤ λN2 implica que
λ−21 ≥ λ−22 ≥ ... ≥ λ−2N2
‖~xδ − ~x‖22 ≤
1
λ21
‖~yδ − ~y‖22 = κ2‖~yδ − ~y‖22 ≤ (κδ)2,
o bien
‖~xδ − ~x‖2 ≤ κδ. (2.16)
El término ‖~xδ−~x‖2 se puede pensar como la “amplificación” del ruido o simplemente
como qué tanto se alejan dos soluciones cuando existen cambios pequeños en las
mediciones (ecuación (2.15)). Entonces para sistemas de ecuaciones con números de
condición grandes no necesariamente se satisface la tercera condición de Hadamard
en el sentido de ~xδ → ~x cuando ~yδ → ~y, pues la cota superior dada por (2.16) da
lugar a posibles grandes diferencias en las soluciones.
Los números de condición grandes corresponden a eigenvalores pequeños y es
aquí donde entra la idea fundamental de los métodos de regularización . De ma-
nera intuitiva se podría pensar en modificar o aproximar As por una familia de
matrices cuyos eigenvalores más pequeños sean alejados del cero. Otra manera
de entender los métodos de regularización es pensarlos como la aproximación de un
problema “mal planteado” por una familia adyacente de problemas “bien planteados”.
26
Un ejemplo sencillo es hacer el cambio a la matriz As dado por
Aα = As + αId ; α > 0. (2.17)
Se puede demostrar que los eigenvalores ~Li de la matriz Aα están dados por λi + α,
con i = 1, ...N2 y los eigenvectores son los mismos de As. Si tenemos ~x = A−1s ~y y
~xα = A
−1
α ~y, entonces usando la representación espectral de As y de Aα de puede
seguir un análisis como el anterior y obtener
~x− ~xα =
N
2∑
i=1
(
1
λi
− 1
λi + α
)
uiu
T
i · ~y =
N
2∑
i=1
α
λi(λi + α)
uiu
T
i · ~y.
De la expresión anterior se puede estimar el error de aproximación de esta técnica
de regularización como
E(α) = ‖~x− ~xα‖2 ≤
α
λ1(λ1 + α)
‖~y‖2. (2.18)
Utilizando un conjunto de datos ruidoso ~yδ y siguiendo un análisis análogo se puede
observar lo siguiente.
~xδα − ~xα =
N
2∑
i=1
1
λi + α
uiu
T
i · (~yδ − ~y).
de donde el error se estima como
E(α, δ) = ‖~xδα − ~xα‖2 ≤
δ
λ1 + α
. (2.19)
Luego, utilizando la desigualdad del triángulo se llega al error total
‖~x− ~xδα‖2 ≤ E(α) + E(α, δ) =
α
λ1(λ1 + α)
‖~y‖2 +
δ
λ1 + α
. (2.20)
Sin embargo como en la práctica no podemos conocer ‖~y‖2 se puede estimar ‖~y‖2 ≤
‖~yδ‖2 + δ y así
‖x− ~xδα‖2 ≤ E(α) + E(α, δ) =
α
λ1(λ1 + α)
(‖~yδ‖2 + δ) +
δ
λ1 + α
. (2.21)
Al parámetro “α” se le conoce como parámetro de regularización y como puede
observarse en la ecuación (2.21) la idea de introducir el parámetro de regularización
es acotar el error entre las soluciones correspondientes a mediciones muy cercanas.
Se puede ver que para un δ fijo y cuando α → 0, se regresa a la condición dada
por (2.16), es decir el primer término es decreciente y el segundo término crece.
Entonces debe existir algún α∗ que minimice el error total. La elección del parámetro
de regularización dependerá del nivel de ruido expresado en términos de δ y de los
eigenvalores de As. Las reglas de elección de parámetros son reglas que buscan escoger
los mejores parámetros de regularización en función del nivel de ruido δ y de los datos
ruidosos ~yδ y estas reglas tienen que ver con la definición formal de lo que son los
métodos de regularización. Una discución precisa y detallada de la definición de los
métodos de regularización puede verse en [Engl et al., 1996] y en [Burger, 2007].
27
2.3.1. Regularización en tomografía
Existen varios casos prácticos bajo los cuales la matriz de proyección en el pro-
blema tomográfico está “mal condicionada”. Un caso es cuando se tienen pocas pro-
yecciones, esto es, cuando se hacen incidir los rayos X pero con pocos ángulos; otro
caso es cuando se presentan las llamadas proyecciones incompletas, que es cuando se
toman las proyecciones en intervalos de ángulos restringidos por ejemplo un conjunto
finito de ángulos θm en el intervalo [π6 ,
5π
6
] en vez de el intervalo completo [0, π] (en
estos dos casos se tiene que la matriz esta subdeterminada y tipicamente se tienen
muchas menos ecuaciones que incógnitas). La presencia de ruido como se mencionó
anteriormente es una fuente inmediata de problemas “mal planteados” y si esto se
combina con los casos de pocas proyecciones y proyecciones incompletas se tiene una
receta directa para “plantear mal” la matriz de proyección.
Es útil, como se verá más adelante, definir un parámetro α
A
que mida el cociente
entre el número de renglones y de columnas de la matriz A.
αA =
M
N2
(2.22)
Con este parámetro es fácil ver que si α
A
= 1 entonces se tiene la misma cantidad
de ecuaciones que incógnitas por lo que en principio se puede calcular la inversa de
la matriz A pero en vez de eso en estos casos se utiliza comúnmente la FBP para
integrar el sinograma y obtener la reconstrucción debido a que es menos costoso en
cuanto a tiempo de cómputo. Si α
A
≤ 1 entonces hay más ecuaciones que incógnitas
y es en este caso en donde se puede empezar a preguntar qué tanto se “mal condi-
ciona” el problema por estos casos y cómo se podría implementar alguna técnica de
regularización para mejorar las reconstrucciones obtenidas.
Lo anterior motiva fuertemente el interés particular de esta tesis en los métodos
de regularización debido a que los problemas de pocas proyecciones y proyecciones
incompletas (o tomografías de ángulo limitado) surgen de manera natural en algunas
aplicaciones prácticas.
Por ejemplo, la tomosíntesis digital de seno, utiliza un tomógrafo que por su di-
seño geométrico y la necesidad del estudio, está forzado a tener un conjunto de pro-
yecciones incompletas (ángulo limitado). También, en una tomografía convencional
de rayos X, se pueden obtener pocas proyecciones cambiando diferentes parámetros
en la instrumentación como son la cantidad de detectores o hacer un barrido con
menos ángulos, esto supondría una cantidad de radiación menor suministrada al
paciente reduciendo los riesgos de la exposición a tal radiación. También existe una
modalidad de tomografía que se parece mucho a las modalidades de PET y SPECT
llamada tomografía óptica de luminiscencia estimulada por rayos X, en
donde un material luminiscente (es decir que emite radiación debido a la excitación
y relajación de sus estados energéticos después de ser radiado por un haz de rayos
X), es la fuente que emite la radiación medida en las proyecciones. En este caso la
función µ no corresponde a los coeficientes de atenuación sino a la distribución del
material luminiscente. Este tipo de tomografía resulta ser tardado, pues se necesita
radiar el material de manera muy precisa y fina con haces colimados de rayos X que
aproximen a un solo rayo e ir midiendo con los detectores correspondientes, cada
una de las proyecciones producidas por estos rayos. Si se quieren reducir los tiempos
28
de adquisición de datos se podría optar por utilizar menos rayos y menos ángulos
de proyección, obteniendo en total muy pocas proyecciones.
La técnica de regularización que se utilizó en este trabajo fue la regularización
por variación total. Se eligió esta técnica debido a que es una técnica que ayuda
a reducir el ruido y su amplificación, y además ayuda mucho a obtener imágenes de
buena calidad cuando se tienen pocas proyecciones debido a la naturaleza matemá-
tica de esta regularización como quedará claro más adelante cuando se introduzca
la teoría del sensado compresivo. En las siguientes dos secciones se introducen
los conceptos de sensado compresivo y regularización por variación total, respecti-
vamente.
2.4. Sensado compresivo
La adquisición, almacenamiento y procesamiento de señales suele lidiar con el
problema del “muestreo” o sampling. Dada la naturaleza discreta del proceso de
adquisición de datos, surge la pregunta de cuántas muestras de la señal se necesitanmedir para obtener suficiente información que identifique la señal y pueda utilizarse
para su almacenamiento y procesamiento. Intuitivamente se puede pensar que entre
más muestras se tomen de una señal, más completa es la información que se puede
llegar a obtener y procesar de tal señal.
Por ejemplo, en el problema tomográfico a mayor cantidad de pixeles, detectores,
y proyecciones (número de rayos y de ángulos), mayor es la cantidad de información
con la que se puede reconstruir la sección transversal de un objeto y obtener una
imagen de muy alta calidad. Otro ejemplo es en la adquisición de fotografías digita-
les en donde es común escuchar que mientras más “Megapixeles” (mayor resolución)
tenga una cámara, mayor es la calidad de la fotografía adquirida. Sin embargo en
estos dos ejemplos surgen inconvenientes cuando se tienen muchas muestras de la
señal (sinograma y fotografía, respectivamente). En el caso de la tomografía como ya
se mencionó, existen casos prácticos en donde un muestreo completo no es posible
o se necesitan tomar menos muestras para reducir los tiempos de adquisición. Para
las fotografías digitales, entre más “Megapixeles”, se requiere mayor capacidad de
almacenamiento de información y ésta suele estar limitada por el tamaño de la me-
moria de almacenamiento del dispositivo. Por ello, existen algoritmos de compresión
de datos que consisten en aplicar una trasformación a la señal, expresándola en una
base conveniente en donde se seleccionen sólo los coeficientes necesarios de su ex-
pansión, se conserve información de la base original y la cantidad de bits necesarios
para su almacenamiento se logre reducir.
La cantidad mínima de muestras necesarias para adquirir una señal está dada
por el teorema Nyquist-Shannon que establece un criterio entre las frecuencias de
muestreo y las frecuencias del espectro de Fourier de la señal adquirida. Por ejemplo
si se tiene una señal en 1D dependiente del tiempo f(t) entonces un muestreo es un
conjunto finito de mediciones fi = f(ti) con i = 1, .., n en donde ∆t = ti+1 − ti es
constante para todo i. Aquí la frecuencia de muestreo esta dada por
Fs =
1
∆t
.
29
El teorema para este tipo de señales afirma que, si las frecuencias de la transformada
de Fourier son de banda limitada (esto es que las frecuencias están dentro de un
intervalo acotado) y si la frecuencia de muestreo satisface Fs ≥ 2kmax, entonces el
muestreo fi define de manera única la señal continua f(t) (con kmax la frecuencia
del espectro de Fourier más grande). Este teorema se puede generalizar para señales
de varias variables como imágenes y aplicarse también a el muestreo de proyecciones
en el caso de tomografía.
La idea intuitiva del sensado compresivo puede entenderse de manera intro-
ductoria pensando en los algoritmos de compresión de imágenes en el sentido de
muestrear la imagen de manera comprimida, es decir, comprimir la información
directamente en el proceso de adquisición; o bien utilizando el teorema de Nyquist-
Shannon, la teoría de sensado compresivo logra reconstruir de manera única la señal
muestreada aunque la frecuencia de muestreo no cumpla Fs ≥ 2kmax. La teoría del
sensado compresivo se introdujo recientemente en los artículos de Donoho [Donoho,
2006], y en los trabajos de Romberg y Candès [Candes and Tao, 2005], [Candes
et al., 2006]. Para una revisión breve e introductoria a la teoría pueden consultarse
los artículos [Candès and Wakin, 2008] y [Romberg, 2008].
2.4.1. Incoherencia y sensado de señales escasas
El sensado compresivo se basa en dos principios: la escasez (en inglés sparsity)
de la señal, y la incoherencia que tiene que ver con la modalidad de sensado.
Para revisar de manera introductoria la teoría del sensado compresivo de aquí
en adelante se trabajará con señales discretas ~f ∈ RN .
Sea una señal ~x ∈ RN y sea una base ortonormal Ψ con elementos ~ψi con i =
1, 2, .., N . Entonces la representación de ~x en Ψ está dada por
~x =
N∑
i=1
xi ~ψi,
en donde xi = 〈~x, ~ψi〉.
Se dice que ~x es una señal escasa si la mayoría de los coeficientes xi en la base
Ψ son cero; de manera más precisa, si a lo más S coeficientes xi con S << N son
distintos de cero, se dice entonces que ~x es S-parse en la representación de la base
Ψ.
Este es el principio que permite comprimir la información guardando sólo los
coeficientes que son distintos de cero, sin embargo la compresión de datos requiere
que se conozcan los N coeficientes y las posiciones de los S coeficientes distintos
de cero. La idea del sensado compresivo es explotar esta propiedad de escasez en el
proceso de adquisición de datos.
Para hablar de incoherencia se necesita expresar en términos matemáticos lo
que es el sensado de la señal ~x. El sensado de la señal (ó mediciones de la señal ~x)
es el vector ~y cuyas componentes están dadas por los productos interiores
yk = 〈~x, ~φk〉, k = 1, ...,m
30
en donde ~φk ∈ RN es un funcional que representa la modalidad de sensado. Estos
valores se ordenan como entradas de un vector ~y y éste representa las mediciones
de la señal ~x. Por ejemplo si son deltas de Dirac, el vector de mediciones ~y, es un
vector que representa los valores de ~x en el dominio temporal o en el dominio espa-
cial; si son funcionales indicadores de pixeles entonces ~y, son los datos típicamente
obtenidos por una cámara digital; si son sinusoides el vector de mediciones repre-
senta los coeficientes de Fourier (como es el caso de la resonancia magnética); en el
caso tomográfico discutido anteriormente el vector de mediciones corresponde a ~b
(proyecciones) y ~x es la imagen tomográfica (señal) que se desea reconstruir.
La idea importante para el sensado compresivo es cómo reconstruir la señal ~x
a partir de m mediciones cuando m << N , ¿qué tan pequeño puede ser m para
obtener una buena reconstrucción de la señal?
Se pueden organizar los funcionales ~φk de tal forma que formen una base orto-
normal Φ y si cada vector se arregla como un reglón de una matriz Φ ∈ Mm×N , el
problema de sensado se puede escribir como
~y = Φ~x.
Con esta terminología se puede introducir la coherencia entre la base de sensado
Φ y la base en donde la señal es escasa, Ψ, como el número dado por
µ(Φ,Ψ) =
√
N max
1≤k,j≤m,N
|〈φk, ψj〉|. (2.23)
Esto significa que la coherencia mide la más grande correlación entre cualesquie-
ra dos elementos de las Φ y Ψ, además µ(Φ,Ψ) ∈ [1,
√
N ]. Entre más pequeño sea
el número µ(Φ,Ψ), las bases son más incoherentes; a mayores µ(Φ,Ψ) más cohe-
rencia entre las bases. El sensado compresivo está interesado en pares de bases con
baja coherencia. Sensar una señal de manera incoherente es construir una base con
funcionales ~φk que sea incoherente con la base en la señal tiene una representación
escasa.
Una modalidad de sensado comúnmente utilizada por esta teoría, es emplear
matrices de sensado aleatorio Φrand. Se ha demostrado que las matrices aleatorias
son incoherentes con cualquier base fija Ψ. La coherencia entre matrices aleatorias
y bases fijas es aproximadamente
√
2logN . Es por esto que comúnmente el sensado
compresivo se asocia con un sensado aleatorio pero esto no es una condición com-
pletamente necesaria para la teoría. Basta con encontrar bases en las que la matriz
de sensado sea incoherente con alguna base en la que la señal es escaza.
Cuando se tienen m mediciones de tal forma que al menos m < N entonces se
tiene un submuestreo y tal como se discute en [Candes et al., 2006], utilizando un
planteamiento basado en algoritmos de optimización convexa, se puede reconstruir
la señal escasa ~x resolviendo el problema
~x = argmin
~x∗
‖~x∗‖
`1
sujeto a ~y = Φ~x∗, (2.24)
en donde ‖ · ‖ representa la norma `1 dada por
∑N
i=1 |xi|. La elección de la norma `1
sobre otras normas tiene que ver con argumentos geométricos cuya discusión puede
verse de manera introductoria en [Romberg, 2008].
31
x⇤(1)
r h
H = {~x⇤/�~x⇤ = ~y} ~x
⇤
`2
x⇤(2)
H = {~x⇤/�~x⇤ = ~y}
~x⇤0
(a) (b) (c)
Figura 1: Top view.
1
Figura 2.1: (a) Bola de radio r en la norma `1: ~x∈ R2 tal que |x(1)| + |x(2)| ≤ r (b)
Optimización de (2.24) (c) Mínimos cuadrados
Intuitivamente, puede entenderse por qué la norma `1 es mejor candidata para
el problema de optimización observando la figura 2.1. La figura 2.1(a) ilustra la bola
de radio r definido con la norma `1 en R2. Esta bola es anisotrópica, es decir no es
invariante ante rotaciones. La figura (b) es el diagrama del programa que optimiza
usando la norma `1: el vector etiquetado ~x∗0 es un vector escaso (solo una de sus dos
componentes es distinta de cero) del cual se toma una sola medición. La línea H es
el conjunto de todos los vectores ~x∗ que comparten la misma medición.
La tarea del problema planteado en la ecuación (2.24) es elegir el vector que sa-
tisface esas mediciones con la norma `1 más pequeña. Para visualizar como funciona
eso, es útil imaginar una bola en la norma `1 de radio muy pequeño y gradualmente
expanderla hasta que intersecte con H. El primer punto de intersección es por defi-
nición, el vector que resuelve (2.24). La combinación de la anisotropía de la bola en
la norma `1 y la dimensión de la recta H resulta en la intersección de uno solo de
sus puntos, precisamente en donde las señales escasas son reconstruidas.
Si en comparación se utiliza la norma `2 (lo que sería un problema de mínimos
cuadrados), la bola en la norma `2 es isotrópica y usando la idea de expansión de
una bola diminuta, se ve que el punto de intersección con el conjunto H no es
necesariamente escaso (figura 2.1 (c)). En RN la diferencia se vuelve mucho mayor.
En el artículo [Candès and Romberg, 2007] se puede encontrar uno de los teore-
mas fuertes de mayor importancia para la teoría de sensado compresivo. El teorema
afirma lo siguiente.
dTeorema 2 : Sea una señal ~x ∈ RN S-sparse en alguna base Ψ. Si se tienen m
muestreos con la base Φ y además m satisface la desigualdad
m ≥ Cµ2(Φ,Ψ) · S · log(N), (2.25)
entonces la solución dada por el problema de optimización convexa (2.24) es exacta
con una probabilidad muy cercana a uno (C es una constante positiva). c
Aquí, se ve el papel que juega la coherencia de las matrices de las bases en las
que se representan las señales. Este resultado es muy importante ya que cuando se
32
tiene un sistema mal planteado, como se comentó en la sección anterior, es difícil
encontrar una solución única al problema de reconstruir la señal a partir de las
mediciones ~y. También dicha desigualdad muestra por qué comúnmente el sensado
compresivo busca utilizar bases de sensado que sean incoherentes con la base escasa
para que la cantidad de muestreos m satisfaga la desigualdad sin tener que ser tan
grande.
Con este resultado, se puede pensar que el algoritmo de minimización convexa
que resuelve (2.24) juega el papel de descompresor de datos.
2.4.2. Fortaleza del sensado compresivo
En la práctica es común que no se tengan señales escasas como se definió an-
teriormente sino aproximadamente escasas y siempre con la presencia de ruido en
el proceso de adquisición de datos. Es decir, en la práctica se presentan casos no
ideales. La teoría del sensado compresivo busca ser capaz de reconstruir este tipo de
señales no ideales. Diremos entonces que la teoría es robusta si logra reconstruir
señales aproximadamente escasas y regularizar el ruido presente en los datos. Si no
es posible reconstruir este tipo de señales ni regularizar el ruido, diremos que la
teoría no es robusta.
Una manera de estudiar la robustez de la teoría es con un concepto que es la
propiedad de isometría restringida (RIP); veáse [Candes and Tao, 2005].
Definición: Para cada entero S=1,2,... se define la constante de isometría δs
de una matriz A como el número más pequeño tal que
(1− δs)‖~x‖2`2 ≤ ‖A~x‖
2
`2
≤ (1 + δs)‖~x‖2`2 (2.26)
se satisface para todo vector S-parse ~x.
Esta propiedad es importante ya que si se satisface, entonces A aproximadamente
preserva longitudes de vectores S-sparse. Esto ayuda a demostrar que dos distintos
vectores escasos ~x1 6= ~x2 no pueden dar las mismas mediciones ~y. Supongamos que
en efecto ~x1 6= ~x2 pero A~x1 = A~x2 = ~y y que A satisface el criterio RIP con algúna
constante suficientemente pequeña para que 1− δs > 0. Entonces se tiene que
0 < (1− δs)‖~x1 − ~x2‖2`2 ≤ ‖A(~x1 − ~x2)‖
2
`2
= 0 ≤ (1 + δs)‖~x1 − ~x2‖2`2 ,
lo cual es una contradicción, por lo tanto es falso suponer que dos vectores distin-
tos S-sparse pueden dar las mismas mediciones ~y. Así, la reconstrucción que da la
optimización dada por (2.24) es única.
Señales aproximadamente escasas.
Para atacar el problema de un vector aproximadamente escaso el siguiente teo-
rema ayuda a la robustez del Sensado Compresivo.
dTeorema 3 : Sea la constante de isometría de la matriz A para vectores 2S-
sparse δ2s tal que satisfaga δ2s <
√
2 − 1; sea también ~x∗ la solución de (2.24).
33
Entonces ~x∗ satisface
‖~x∗ − ~x‖`2 ≤ C0‖~x− ~xs‖`1/
√
S y ‖~x∗ − ~x‖`1 ≤ C0‖~x− ~xs‖`1 , (2.27)
en donde C0 es una constante, y ~xs es el vector ~x con todas excepto las S componentes
más grandes forzadas a ser cero (esto es, ~xs es un vector aproximadamente escaso),
veáse [Candès et al., 2006]. c
Este teorema ayuda a la robustez de la teoría en el sentido que dado un vector
aproximadamente escaso, es posible acotar la diferencia entre la señal verdadera ~x
y la señal reconstruida ~x∗ y además, si la señal original ~x = ~xs (es decir la señal
es exactamente S-sparse), entonces la reconstrucción dada por (2.24) es única y
exacta. Además a diferencia del teorema 2, este resultado es determinista y no
habla de alguna probabilidad, lo cual lo hace más fuerte.
Señales con ruido
Para el caso de las señales con ruido se hace una pequeña modificación al pro-
blema de optimización de la siguiente forma:
~x = argmin
~x∗
‖~x∗‖
`1
sujeto a ‖A~x∗ − ~y‖`2 ≤ �, (2.28)
en donde � es una cota para la cantidad de ruido en los datos. El siguiente teorema
da una cota para la precisión de la reconstrucción dada por la solución a (2.28)
dTeorema 4 : Sea la constante de isometría de la matriz A para vectores 2S-
sparse δ2s tal que satisfaga δ2s <
√
2 − 1; sea también ~x∗ la solución de (2.28).
Entonces ~x∗ satisface
‖~x∗ − ~x‖`2 ≤ C0‖~x− ~xs‖`1/
√
S + C1�
en donde C0 y C1 son algunas constantes, y nuevamente ~xs es el vector ~x con todas
excepto las S componentes más grandes forzadas a ser cero. c
Este teorema muestra que el error de la reconstrucción está acotado por dos
términos, el primer término es el error que se tendría si no se tuviera una señal
ruidosa y el segundo es un término proporcional a la cota del ruido. Además, si
fuera una señal exactamente escasa ~x = ~xs, entonces se tendría que el error de la
solución estaría únicamente dado por ~x∗ − ~xs ≤ C1� [Romberg, 2008].
Por tanto, la teoría de sensado compresivo es robusta frente a estas dos situacio-
nes no ideales y por ello se ha intentado aplicar a distintas áreas siendo una de ellas
la reconstrucción de imágenes tomográficas.
2.5. Regularización por Variación Total.
Como se mencionó en la sección métodos de regularización, existen algunas situa-
ciones en las cuales se tienen proyecciones medidas a pocos ángulos y proyecciones
34
en intervalos restringidos, ambos casos con la presencia de ruido debido al proceso
de adquisición de los datos. Recordando el problema tomográfico dado por (1.13)
y el parámetro αA (2.22), en estos casos se tiene que M < N2 y es necesario un
método de regularización para reconstruir con mejor calidad la imagen. Al tenerM
mediciones se puede pensar en relacionar esta situación con las situaciones descritas
por la teoría de sensado compresivo para señales en RN
2
. Debido a que las imágenes
son más bien matrices de N ×N entonces se necesita plantear un poco diferente la
teoría del sensado compresivo. Como se discute en [Candes et al., 2006] la manera de
aplicar la teoría a imágenes es con la seminorma conocida como Variación Total
(TV). De manera intuitiva la variación total está relacionada con la norma `1 del
gradiente de la imagen, el cual es importantey es el vínculo más fuerte con el sen-
sado compresivo pues como se mostrará gráficamente más adelante, los gradientes
de la mayoría de las imágenes tomográficas son escasos.
Para señales en una dimensión ~x ∈ RN la variación total se define como
TV (~x) =
N−1∑
i=1
|xi+1 − xi| (2.29)
Si se define el gradiente de la señal ~x como
∇~x = (∂1, ∂2, ..., ∂N−1)
en donde ∂k = xk+1 − xk, entonces la variación total se puede escribir como
TV (~x) = ‖∇~x‖`1 (2.30)
Dos definiciones para la variación total de una señal en dos dimensiones como es el
caso de una imagen, fueron utilizadas en esta tesis (veánse [Goldstein and Osher,
2009], [Yan and Lu, 2015], [Chambolle, 2004], [Chambolle, 2005]) . La diferencia entre
ambas definiciones tiene que ver con la invariancia ante rotaciones de la imagen.
Para definir la variación total isotrópica de una imagen u ∈M
N×N , sea
∇i,ju = (∂xui,j, ∂yui,j),
el gradiente (en primera aproximación) de la imagen en el pixel i, j, en donde
∂xui,j = ui+1,j − ui,j y ∂yui,j = ui,j+1 − ui,j,
luego, la variación total isotrópica es
TV (I)(x) =
N−1∑
i=1
N−1∑
j=1
√
(∂xui,j)2 + (∂yui,j)2. (2.31)
Utilizando la misma notación, la variación total anisotrópica está dada por
TV (A)(x) =
N−1∑
i=1
N−1∑
j=1
|∂xui,j|+ |∂yui,j| (2.32)
Utilizando la variación total se han propuesto técnicas de regularización como
problemas de optimización convexa que son parte de la teoría de sensado compresivo.
35
A este tipo de regularización de le conoce como regularización por variación
total y cae dentro de las técnicas de regularización por escasez como se discute
en [Jørgensen and Sidky, 2015].
El problema de regularización por variación total es el más importante para el
problema tomográfico en el que se centra esta tesis pues la mayoría de estas imágenes
tienen gradientes escasos (por ejemplo el maniquí Shepp-Logan, figura 2.2). Debido
que la variación total se relaciona con el gradiente de la imagen, entonces se puede ver
que la variación total hereda las propiedades de escasez del gradiente y por lo tanto
se puede intentar aplicar la teoría de sensado compresivo utilizando la regularización
por variación total. Distinguiendo de casos ideales y no ideales en cuanto a presencia
de ruido se refiere, los dos problemas principales que se resolvieron en esta tesis son
( [Candes et al., 2006]):
x = argmin
x∗
TV α(x∗) sujeto a ‖A~x∗ − ~y‖`2 ≤ �, (2.33)
en donde α = I, A; ~x∗ ∈ RN2 es el vector que resulta de ordenar la imagen x ∈
M
N×N de acuerdo a la geometría dada por el problema tomográfico y a la matriz de
proyección A y � es una cota dada por la cantidad y tipo de ruido presente en las
proyecciones. Para el caso ideal en ausencia de ruido el problema se simplifica a
x = argmin
x∗
TV α(x∗) sujeto a A~x∗ = ~y, (2.34)
50 100 150 200 250
50
100
150
200
250
(a)
50 100 150 200 250
50
100
150
200
250
(b)
Figura 2.2: (a) Maniquí Shepp-Logan. (b) Gradiente del maniquí Shepp-Logan.
Como puede verse a simple vista en el gradiente del maniquí Shepp-Logan, en
donde todos los pixeles en color negro corresponden a ceros y los que son distintos
de cero se muestran en color blanco, la gran mayoría de los pixeles son cero por
lo que se puede pensar en que el gradiente tiene una representación escasa y por
tanto se puede aplicar las ideas del sensado compresivo para la reconstrucción de
las imágenes tomográficas. El maniquí tiene 256× 256 pixeles (65,536) de los cuales
27, 409 son distintos de cero mientras que el gradiente solo tiene 2, 184 elementos
distintos de cero.
Resulta ser que resolver el problema (2.33) tal y como predice la teoría de sen-
sado compresivo logra reconstruir buenas imágenes tomográficas en presencia de
36
ruido y para pocas proyecciones. Existen estrategias para resolver dicho problema
de optimización convexa, siendo los métodos de gradientes conjugados y gradientes
conjugados precondicionados los más utilizados para ello. En esta tesis se implemen-
taron dos algoritmos principales para la solución de este problema. El desarrollo de
estos algoritmos, así como su aplicación a datos sintéticos se detalla en el siguiente
capítulo y su aplicación a datos experimentales se muestra en el capítulo 4. Una
estrategia similar fue utilizada en [Sidky et al., 2006].
37
Capítulo 3
Desarrollo e implementación de
algoritmos
3.1. Introducción
Como se vio en el capítulo anterior, uno de los problemas que intenta atacar
esta tesis es el de reconstruir imágenes tomográficas cuando el sistema (1.13) está
mal planteado. También se discutió que existen varias maneras de mal plantear el
sistema y en particular, esta tesis se enfocó específicamente en el caso de pocas
proyecciones con presencia de ruido, dejando como perspectiva a futuro la apli-
cación de una estrategia similar para el problema de proyecciones incompletas en
intervalos angulares restringidos.
Recordando el sistema de ecuaciones que representa el problema inverso tomográ-
fico A~x = ~b, se tiene que A ∈M
M×N2
conM = Nrays×Nθ, siendo Nrays la cantidad
de rayos del haz definida por la cantidad de detectores y Nθ la cantidad de ángulos
en las que se miden las proyecciones como se explicó en el capítulo introductorio.
Con esto se puede ver que se tienen dos maneras de reducir la cantidad de pro-
yecciones: reduciendo el número de ángulos de proyección o reduciendo la cantidad
de rayos que se utilizan para radiar el objeto. Para reducir la cantidad de rayos
del haz, se debe reducir la cantidad de detectores utilizados en la instrumentación
como se observa en la figura 1.5, o por ejemplo en el caso de la tomografía óptica
de luminiscencia estimulada por rayos X (en el que se tiene un haz colimado que se
va “barriendo” o desplazando en distintas posiciones para cada ángulo de proyección
utilizado), reducir la cantidad de posiciones de barrido.
Utilizando los conceptos de sensado compresivo y con los algoritmos implemen-
tados que se describen en las siguientes secciones, la imagen tomográfica recons-
truida que es solución a (1.13) xr, está dada por
xr = argmin
x
TV α(x) sujeto a ‖A~x−~b‖`2 ≤ � y xi ≥ 0 ∀i = 1, .., N2 (3.1)
en donde α = I, A (indicando si se utiliza la definición isotrópica o anisotrópica de
la variación total, respectivamente), ~x ∈ RN2 es el vector que resulta de ordenar la
38
imagen x ∈M
N×N , � es una cota dada por el nivel de ruido en los datos y la última
constricción es una condición de positividad que se introduce para asegurar que todos
los coeficientes de atenuación representados en el vector ~x sean mayores o iguales a
cero, ya que un coeficiente negativo, indicaría que en lugar de atenuarse, el rayo X
que atraviesa esa zona estaría aumentando de intensidad, cosa que físicamente no
tiene sentido.
En las siguientes secciones se describen los algoritmos implementados y su apli-
cación a datos sintéticos dados por maniquís. Los maniquís utilizados fueron el ya
mostrado Shepp-Logan y el Forbild Head phantom que trata de representar de mane-
ra simple una sección transversal de un cerebro (veáse [Yu et al., 2012]). El maniquí
Shepp-Logan se ha mostrado anteriormente y en la figura 3.1 se muestra el For-
bild Head phantom junto con su gradiente para mostrar su escasez. Estos maniquís
50 100 150 200 250
50
100
150
200
250
(a)
50 100 150 200 250
50
100
150
200
250
(b)
Figura 3.1: (a) Maniquí Forbild Head . (b)Gradiende del maniquí Forbild.
se tomaron como imágenes de referencia para calcular el error cuadrático medio
introducido previamente en la ecuación 2.11:
MSE =
1
N2
‖~xref − ~x‖22.
3.2. Ruido en maniquís y perfiles de intensidad.
A los maniquís se les añadió ruido sintético blanco (es decir ruido con distribución
Gaussiana) de manera aditiva. Se tomó el valor máximo de pixel de las imágenes
(valor aproximado a 100 en la escala de grises mostrada en las figuras) y el 5% de
este valor se definió como la semianchura de la distribución del ruido (tal porcentaje
de ruido se eligió a modo de ensayo y error de tal forma que cualitativamente fuera
un nivel