Nuevo método de reconstrucción de mapa de dosis para verificación de tratamientos de radioterapia avanzados

FACULTAD DE F´ ISICA M´ ASTER INTERUNIVERSITARIO EN F´ ISICA NUCLEAR DEPARTAMENTO DE F´ ISICA AT´ OMICA, MOLECULAR Y NUCLEAR TRABAJO FIN DE M´ ASTER: NUEVO M´ ETODO DE RECONSTRUCCI´ ON DE MAPA DE DOSIS PARA VERIFICACI´ ON DE TRATAMIENTOS DE RADIOTERAPIA AVANZADOS Alumno: Antonio Dami´an Dom´ınguez Mu˜noz Tutora: D.aM. Isabel Gallardo Fuentes CURSO 2016-2017 RESUMEN Las t´ecnicas en radioterapia est´an en continua evoluci´on. Su objetivo es conseguir una distribuci´on de dosis tal que sea m´axima en la zona tumoral y m´ınima en el tejido sano para evitar efectos secundarios no deseados. El desarrollo de t´ecnicas avanzadas en radioterapia conlleva una mayor complejidad en los tratamientos, lo que hace necesario un desarrollo paralelo de t´ecnicas de verificaci´on dosim´etricas. En este trabajo se presenta un prototipo de sistema de detecci´on dise˜nado espec´ıficamente para la verificaci´on de tratamientos de radioterapia avanzados en el plano axial del paciente. El sistema est´a basado en dos detectores de tiras de silicio dispuestos en configuraci´on perpendicular y calibrados en dosis, que se insertan en el plano axial de un maniqu´ı cil´ındrico que puede girar sobre su eje. As´ı este prototipo mide la dosis media en agua depositada en cada tira en las distintas posiciones angulares. A partir de estos datos y mediante la descripci´on del mapa de dosis mediante p´ıxeles, puede desarrollarse un nuevo m´etodo de reconstrucci´on que presenta mayores ventajas que los m´etodos basados en la transformada de Radon que se ven´ıan usando hasta ahora. Teniendo en cuenta la naturaleza estad´ıstica del proceso de medida, se aplica el m´etodo de m´aximo a posteriori y se desarrolla el m´etodo de esperanza-maximizado (EM) para darle soluci´on. Los resultados obtenidos con este m´etodo y con una generalizaci´on del mismo (GEM) son satisfactorios, quedando futuras l´ıneas de investigaci´on pendientes para la obtenci´on de mapas de mejor calidad. ´ Indice general 1. Introducci´on 1 2. Material y M´etodo 3 2.1. Sistemadedetecci´on .............................. 3 2.2. Descripci´on matem´atica del problema . . . . . . . . . . . . . . . . . . . . . 5 2.2.1. Formalismo anal´ıtico . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2.2. Formalismo de discretizaci´on en p´ıxeles . . . . . . . . . . . . . . . . 7 2.3. Modelado geom´etrico del sistema de detecci´on . . . . . . . . . . . . . . . . 8 2.3.1. Descripci´on de los detectores . . . . . . . . . . . . . . . . . . . . . . 9 2.3.2. Descripci´on de los p´ıxeles . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.3. Estrategia a seguir . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4. M´etodos de reconstrucci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4.1. M´etodo iterativo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4.2. M´etodo estad´ıstico . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4.3. Reconstrucci´on de mapa de dosis usando el m´etodo estad´ıstico . . . 18 3. Resultados 29 3.1. Modeladodeldetector ............................. 29 3.2. Reconstrucci´on de mapa de dosis . . . . . . . . . . . . . . . . . . . . . . . 33 3.2.1. Reconstrucci´on de mapas artificiales sin informaci´on a priori usando elalgoritmoEM............................. 33 3.2.2. Reconstrucci´on de mapas artificiales con informaci´on a priori usandoelalgoritmoGEM.......................... 42 5 3.2.3. Reconstrucci´on del mapa de dosis de un tratamiento real sin informaci´on a priori usando el algoritmo EM . . . . . . . . . . . . . . . 42 4. Conclusiones 49 A. C´alculo de los coeficientes del modelado del detector 51 A.1. C´alculo del coeficiente Cijkαm para un pixel dado ij ............. 51 A.1.1. Condici´on de contribuci´on de los p´ıxeles . . . . . . . . . . . . . . . 51 A.1.2. C´alculo de los puntos de corte . . . . . . . . . . . . . . . . . . . . . 52 A.1.3. C´alculo de la contribuci´on del pixel por el que pasa la tira . . . . . 54 A.1.4. Apunte para los casos particulares . . . . . . . . . . . . . . . . . . . 55 A.2. C´alculo de los coeficientes Gijklm para los p´ıxeles en los que se produce el corte de dos l´ıneas del detector . . . . . . . . . . . . . . . . . . . . . . . . 56 B. Desigualdad de Jensen. Aplicaci´on a un caso caracter´ıstico. 59 C. Distribuci´on de probabilidad a priori de Gibbs 61 Cap´ıtulo 1 Introducci´on Los tratamientos de radioterapia con fotones est´an en continua evoluci´on y se han ido volviendo m´as complejos en los ´ultimos a˜nos. En la actualidad los aceleradores cl´ınicos son capaces de irradiar en condiciones que se pueden considerar no est´andares seg´un los patrones tradicionales. Se pueden conseguir campos de radiaci´on muy peque˜nos para los que puede que se hayan perdido las condiciones de equilibrio electr´onico lateral. Adem´as, en los nuevos tratamientos el n´umero de campos (aperturas e incidencias) usados puede ser muy elevado y la dosis depositada por algunos de ellos muy peque˜na. Entre estos nuevos tratamientos est´a la denominada radioterapia de intensidad modulada (IMRT, por sus siglas en ingl´es). La finalidad de estos tratamientos es ajustar al m´aximo la dosis en el tumor y proteger el tejido sano adyacente permitiendo en algunos casos escalar la dosis en el primero. Por esos motivos, el planificador, que es un programa inform´atico encargado de realizar el c´alculo de la dosis que se depositar´ıa en el paciente, no siempre est´a preparado para las condiciones extremas de estos tipos de tratamientos. Es preciso, entonces, verificar que los campos de tratamiento depositan la dosis que ha sido calculada por el planificador. La verificaci´on de un tratamiento se realiza habitualmente mediante medidas experimentales de la dosis, que se comparan con las calculadas por el planificador. Las medidas experimentales se obtienen lanzando el tratamiento real sobre un maniqu´ı que simula la forma humana, en el cual est´a insertado el sistema detector. Como la planificaci´on de tratamientos suele realizarse sobre los cortes axiales del paciente, obtenidos con alguna t´ecnica de imagen m´edica (TAC, RMI, TAC+PET, RMI+PET...), es interesante realizar la verificaci´on en estos planos. La calibraci´on de un detector que funcione como dos´ımetro se suele hacer a dosis en agua, esto es, la dosis que se depositar´ıa en el volumen activo del detector si este se sustituye por agua. La raz´on de utilizar el agua como patr´on se basa en su f´acil uso y su equivalencia a tejido. Como la dosis es una magnitud puntual, la lectura del dos´ımetro, que es la dosis promediada en su volumen activo, se le asigna a un punto efectivo de medida. 1 2 CAP´ ITULO 1. INTRODUCCI´ ON Figura 1.1: Plano axial del cuerpo humano. Como se ha comentado antes, es interesante medir la distribuci´on de dosis en los planos axiales. Sin embargo, resulta complicado encontrar un detector que pueda calibrarse a dosis en agua recibiendo la radiaci´on a lo largo de su plano activo, es decir, cuando la direcci´on del haz es paralelo al plano del detector. De este tipo, se puede hablar hasta ahora ´unicamente de la pel´ıcula radiocr´omica. El problema de ´esta es que su lectura no es en tiempo real, sino que se analiza posteriormente, lo cual ralentiza el proceso de verificaci´on y por tanto del posterior tratamiento de los pacientes. Otros sistemas de detecci´on que existen en el mercado solo realizan medidas puntuales del mapa de dosis en la superficie del maniqu´ı o en su interior. Con la informaci´on de la dosis en distintos puntos corrigen el mapa calculado por el planificador para inducir el mapa completo en los planos axiales. En este trabajo se presenta un nuevo m´etodo de obtenci´on de mapa de dosis en el plano axial, de forma casi inmediata, a partir del desarrollo de numerosas t´ecnicas de reconstrucci´on del mapa de dosis usando medidas experimentales realizadas en dicho plano por un prototipo de sistema detector dise˜nado espec´ıficamente para este fin [1]. En el segundo cap´ıtulo, primero se presentan las caracter´ısticas principales del prototipo de sistema de detecci´on. Esto permite describir la magnitud que es objeto de medida y la forma en la que ´esta se mide. Se procede entonces a realizar la descripci´on matem´atica del problema de reconstruir el mapa de dosis a partir de las medidas experimentales. Para obtener una descripci´on completa, resulta necesario modelar la respuesta del sistema detector ante un determinado mapa de dosis. Por ello se presenta tambi´en el algoritmo usado para calcular los coeficientes necesarios para traducir un mapa de dosis a la lectura del sistema detector. Por ´ultimo, se desarrolla un primer m´etodo de soluci´on a dicho problema usando t´ecnicas similares a las que se emplean hoy en d´ıa en reconstrucci´on de imagen m´edica (SPECT, PET...), pero que nunca antes se han aplicado a la reconstrucci´on de mapa de dosis para la verificaci´on de tratamientos. En el tercer cap´ıtulo se muestran los resultados obtenidos, tanto del modelado del detector como de la reconstrucci´on del mapa de dosis. Para ´este ´ultimo caso se expone tanto la reconstrucci´on de mapas artificiales como la realizada a partir de datos experimentales obtenidos con el prototipo. Con el cuarto cap´ıtulo concluye el trabajo. En ´el se recogen los principales avances conseguidos, lo que a´un queda por hacer y lo que se espera obtener en un futuro pr´oximo. Cap´ıtulo 2 Material y M´etodo 2.1. Sistema de detecci´on El sistema de detecci´on, como ya se ha mencionado anteriormente, ha sido dise˜nado, desarrollado y construido espec´ıficamente para obtener experimentalmente mapas de dosis 2D axiales de tratamientos radioterap´euticos para su verificaci´on. Est´a basado en un detector que se detallar´a a continuaci´on insertado un plano axial de un maniqu´ı cil´ındrico capaz de rotar alrededor de su eje. En la fecha en que se ha realizado este trabajo el prototipo actual es un modelo m´as avanzado [2] del que se us´o en un principio [3]. El prototipo est´a formado por dos detectores de tiras de silicio de una sola cara (en ingl´es single-sided silicon strip detector, SSSSD) de 64x64 mm2de ´area activa, 500µm de grosor y separados por 500µm de kapton del mismo tama˜no que las obleas de silicio. Cada una de las obleas de silicio est´a dividida en 32 tiras de 2 mm de grosor cada una. Con las tiras colocadas en configuraci´on perpendicular, el ensamblaje de estos dos detectores ha recibido el nombre de dual SSSSD (o DSSSSD). Figura 2.1: Versi´on final del DSSSSD para verificaci´on de radioterapia. En la foto se observa la parte inferior de la caja en la que se inserta el detector para su protecci´on. 3 10 CAP´ ITULO 2. MATERIAL Y M´ ETODO 2.3.2. Descripci´on de los p´ıxeles En el caso de los p´ıxeles, lo ´unico que hay que tener en cuenta en su representaci´on es la numeraci´on de la esquinas, la cual se ve en la figura 2.4. Figura 2.4: Descripci´on de un pixel. Los p´ıxeles no tienen por qu´e ser cuadrados, por lo que es ´util conocer el ´angulo de la diagonal del pixel con, por ejemplo, la vertical, θ∗. Este viene dado por θ∗= arctan w h.(2.9) 2.3.3. Estrategia a seguir El primer paso es descomponer el c´alculo de Gen dos partes Gijklm =Cancho ijklm ·Calto ijkm,(2.10) donde Cancho ijklm se corresponde con la fracci´on del pixel ij que quedar´ıa dentro del ancho de la tira klm suponiendo que su altura fuera infinita y Calto ijkm se corresponde con la fracci´on del pixel ij que quedar´ıa dentro del alto del detector men la posici´on angular ksin tener en cuenta la anchura finita. As´ı solo cuando ambos coeficientes sean distintos de cero, es decir, que est´en dentro de una tira dada de anchura finita y altura finita, el coeficiente buscado es distinto de cero. Cada uno de estos coeficientes puede calcularse estudiando las l´ıneas por separado de la manera que se describe a continuaci´on Cancho ijklm =Cijklm −Cijk(l+1)m,(2.11) Calto ijk,m=1 =Cijk,l=1,m=2 −Cijk,l=33,m=2, Calto ijk,m=2 =Cijk,l=1,m=1 −Cijk,l=33,m=2.(2.12) En (2.12) puede verse que, como ya se ha dicho, solo es necesario estudiar las l´ıneas de los bordes de las tiras ya que los bordes superior e inferior de un detector se corresponden, respectivamente, con los bordes 33 y 1 de las tiras del otro detector. 2.4. M´ ETODOS DE RECONSTRUCCI´ ON 11 Hay que especificar entonces la forma en que se calcula Cijkαm con α= 1,2, ...33, y esta es Cijkαm                                = 0,si los valores de las proyecciones de todos los puntos del pixel sobre el eje OX’ o OY’, seg´un sea m= 1 o m= 2, son menores que el valor de la posici´on de la l´ınea que define el borde (x0 αoy0 α) ∈(0,1),si la l´ınea que define el borde cruza a trav´es del pixel = 1,si los valores de las proyecciones de todos los puntos del pixel sobre el eje OX’ o OY’, seg´un sea m= 1 o m= 2, son mayores que el valor de la posici´on de la l´ınea que define el borde (x0 αoy0 α) . Para el caso de que la l´ınea cruce a trav´es del pixel hay que calcular la fracci´on de ´area del pixel formada por los puntos cuya proyecci´on en el sistema rotado sean m´as grandes que el valor de la posici´on de la l´ınea que define el borde (x0 αoy0 α, seg´un sea m= 1 o m= 2). Existen ciertos p´ıxeles en cada posici´on angular del sistema detector para los que los coeficientes Gno est´an bien determinados con la ecuaci´on (2.10). Estos p´ıxeles son aquellos en los que se produce el corte del borde de una tira con uno de los bordes superior o inferior del detector. Hay que localizar cuales son estos p´ıxeles y dise˜nar un algoritmo espec´ıfico para el c´alculo del coeficiente correcto. Por ´ultimo, el sistema detector en su giro barre un ´area del mapa de dosis. Hay que calcular los coeficientes de los p´ıxeles que est´an dentro de dicho ´area. El mapa completo tiene forma rectangular por lo que se calculan m´as coeficientes de los que interesa y se consume tiempo de c´alculo en coeficientes innecesarios. Esto es por lo que hay que adaptar el tama˜no del mapa rectangular al ´area que barre el detector. Lo m´as sencillo es calcular las posiciones de las esquinas del sistema detector en las distintas posiciones angulares que ser´an los conjuntos {Xk}y{Yk}, y buscar los valores m´aximo y m´ınimo de cada conjunto. 2.4. M´etodos de reconstrucci´on En esta secci´on se estudia la forma de obtener el mapa de dosis dado por el conjunto DP ={DPij}, a partir del conjunto de datos experimentales DT ={DTklm}, teniendo en cuenta la relaci´on (2.8). Para simplificar la notaci´on, de aqu´ı en adelante se agrupar´an los sub´ındices de la siguiente forma {i, j}→{u} {k, l, m}→{v}, quedando la relaci´on dosis promedio en el pixel ua dosis promedio en la tira vcomo sigue DTv=AP AT X u Guv ·DPu,∀v. (2.13) 12 CAP´ ITULO 2. MATERIAL Y M´ ETODO En los dos primeros apartados se muestran dos m´etodos distintos para resolver el problema, el iterativo y el estad´ıstico. El primero se describe en l´ıneas generales y se exponen algunos algoritmos para su consulta. Del segundo se explica su planteamiento y se desarrolla un tipo de algoritmo para su resoluci´on en un caso gen´erico. En el tercer apartado se aplica el m´etodo estad´ıstico al caso espec´ıfico que ocupa este trabajo. 2.4.1. M´etodo iterativo Existe una ecuaci´on (2.13) para cada v, por lo que se tiene un sistema de ecuaciones lineal. En la pr´actica, por lo general, la matriz de coeficientes de este sistema de ecuaciones no tiene inversa y, adem´as, el sistema suele ser compatible indeterminado, por lo que existen m´ultiples combinaciones {DPu}que reproduzcan los mismos valores {DTv}. As´ı resulta impracticable obtener una soluci´on directa. Los algoritmos iterativos obtienen una soluci´on al sistema lineal de ecuaciones mediante la repetici´on en bucle de un conjunto de pasos, partiendo de una estimaci´on inicial de los valores de los par´ametros desconocidos. Todo m´etodo iterativo contiene tres pasos principales. El primero es el paso de proyecci´on, que a partir de los valores de los par´ametros obtenidos en la iteraci´on anterior, obtiene artificialmente los datos experimentales correspondientes a dichos valores. En el caso estudiado ser´ıa aplicar (2.13). El segundo paso es el paso de comparaci´on, que consiste en comparar los datos experimentales creados artificialmente con los datos reales medidos, obteniendo as´ı un t´ermino de correcci´on. El ´ultimo es el paso de actualizaci´on, que consiste en obtener nuevos valores de los par´ametros aplicando la correcci´on obtenida en el paso anterior a los valores de partida en esa iteraci´on. Estos nuevos valores ser´an los valores de partida de la iteraci´on siguiente. Mientras que el paso de proyecci´on viene definido por el problema que es objeto de estudio, los otros dos pasos son de libre elecci´on, y, por lo general, suelen aparecer como uno solo (paso de comparaci´on + actualizaci´on). Adem´as, para que el proceso descrito acabe es necesario definir previamente una regla de detenci´on, que puede ser un n´umero m´aximo de iteraciones o que la variaci´on total de los valores de los par´ametros de una iteraci´on a otra queden por debajo de un umbral predefinido. En la figura 2.5 se esquematiza la estructura m´as simple de los algoritmos iterativos para su mejor entendimiento. Algunos de los algoritmos iterativos son ART (Algebraic Reconstruction Technique) [6], SIRT (Simultaneous Iterative Reconstruction Technique) [7] o MART (Multiplicative Algebraic Reconstruction Technique) [8], entre otros. En este trabajo no se aplicar´a este tipo de m´etodos, la raz´on se da en el siguiente apartado. 2.4.2. M´etodo estad´ıstico Este trabajo est´a enmarcado dentro de las aplicaciones m´edicas de la f´ısica nuclear. La base te´orica es la interacci´on radiaci´on materia, por lo que la naturaleza estad´ıstica de este fen´omeno debe aparecer en la formulaci´on del problema. Una soluci´on obtenida a partir de un sistema de ecuaciones lineales ser´ıa una soluci´on determinista que estar´ıa, entonces, omitiendo informaci´on inherente al problema. 2.4. M´ ETODOS DE RECONSTRUCCI´ ON 13 ENTRADA SALIDA BUCLE VALORES DE PARTIDA DATOS MEDIDOS PROYECCIÓN DATOS ARTIFICIALES COMPARACIÓN + ACTUALIZACIÓN NUEVOS VALORES VALORES INICIALES STOP? INICIO NO VALORES FINALES SI Figura 2.5: Esquema general de los algoritmos iterativos Los algoritmos estad´ısticos tratan las medidas experimentales como un conjunto de variables aleatorias, en las que las inc´ognitas son los par´ametros de la funci´on de distribuci´on de dicho conjunto. Los distintos m´etodos variar´an dependiendo del estad´ıstico escogido para estimar el valor de los par´ametros. El m´etodo que se emplear´a en esta secci´on es el de m´axima verosimilitud o maximum likelihood estimation (MLE) [9], en ingl´es. Es un procedimiento para estimar los valores de los par´ametros de una distribuci´on de probabilidad a partir de los valores que toman las variables aleatorias. Sea un conjunto de variables aleatorias y={yv}y un conjunto de par´ametros λ={λu}, si la distribuci´on de probabilidad de yviene dada por f(y|λ), la funci´on de verosimilitud, L, es la misma funci´on que da la distribuci´on de probabilidad, pero cambiados los roles de los par´ametros y de las variables aleatorias, es decir L(λ|y) = f(y|λ).(2.14) El estad´ıstico usado en la estimaci´on por m´axima verosimilitud, ˆ λ, es el valor del par´ametro que hace m´axima la funci´on de verosimilitud, esto es L(ˆ λ|y)>L(λ|y),∀λ6=ˆ λ. (2.15) En la mayor´ıa de los casos resulta conveniente trabajar con el logaritmo neperiano de la funci´on de verosimilitud `(λ|y) = ln L(λ|y).(2.16) Si se toma la primera derivada respecto a cualquier par´ametro ∂`(λ|y) ∂λu =1 L(λ|y) ∂L(λ|y) ∂λu ,(2.17) 14 CAP´ ITULO 2. MATERIAL Y M´ ETODO como L(λ|y)>0, entonces signo ∂`(λ|y) ∂λu!= signo ∂L(λ|y) ∂λu!, por lo que ambas funciones crecen y decrecen en los mismos intervalos, de manera que los m´aximos y m´ınimos de ambas funciones se producen para los mismos valores de λ. Puede as´ı trabajarse con el logaritmo de la funci´on verosimilitud, que resulta especialmente ´util cuando el conjunto yest´a constituido por variables aleatorias estad´ısticamente independientes debido a que `(λ|y) = ln "Y vL(λ|yv)#=X v ln L(λ|yv) = X v `(λ|yv).(2.18) Los estimadores de m´axima verosimilitud se obtienen a partir del sistema de ecuaciones que viene dado por ∂`(λ|y) ∂λuˆ λu = 0,∀u, (2.19) con la condici´on adicional det[H(`)(w)] <0, donde H(`)(w) es la matriz hessiana de `. Cuando se conoce informaci´on a priori de los par´ametros, f(λ), entonces la funci´on a posteriori, f(λ|y), est´a relacionada con la distribuci´on de probabilidad f(y|λ) mediante el Teorema de Bayes f(λ|y) = f(y|λ)·f(λ) f(y),(2.20) donde f(y) realiza el papel de una constante normalizadora, ya que los valores de yson conocidos. En este caso el m´etodo recibe el nombre de estimaci´on por m´aximo a posteriori (MAP). Al igual que antes puede definirse el logaritmo neperiano de la funci´on a posteri. Los estimadores de m´aximo a posteriori se obtiene a partir del sistema de ecuaciones ∂P(λ|y) ∂λuˆ λu = 0,∀u, (2.21) siendo P(λ|y) = ln f(y|λ) + ln f(λ).(2.22) Y la condici´on adicional det[H(P)] >0, donde H(P)es la matriz hessiana de P. Ambos m´etodos acaban planteando un sistema de ecuaciones ligadas para el conjunto λ, el cual puede ser no lineal. Esto complica su resoluci´on, por lo que hay que buscar alg´un m´etodo alternativo. Claramente, en el caso de ausencia de informaci´on a priori (f(λ)=cte) el m´etodo MAP se reduce al m´etodo MLE. Por ello, en lo que sigue, se abordar´a otra forma de darle soluci´on al m´etodo MAP que es el m´as general. 2.4. M´ ETODOS DE RECONSTRUCCI´ ON 15 Algoritmo de esperanza-maximizaci´on (EM) Este algoritmo, cuyo nombre en ingl´es es expectation-maximization algorithm, fue expuesto por primera vez en [10]. En este art´ıculo se generaliza el procedimiento a seguir aplicado al m´etodo ML. A continuaci´on se aplica el mismo desarrollo te´orico al caso de MAP. Para explicar este algoritmo hay que definir un entorno de trabajo basado en dos conjuntos de datos: el incompleto y el completo. El primero recibe su nombre debido a que o bien en el conjunto hay informaci´on perdida o bien la informaci´on que aporta no permite conocer ciertos aspectos del problema que se trata. Este conjunto, en la notaci´on usada, es el conjunto y, el cual puede ser medido. El segundo es un conjunto que, aunque no puede ser medido, puede definirse de manera que bien no tenga p´erdida de informaci´on o en el que el problema tenga una formulaci´on m´as natural y pueda obtenerse la informaci´on que con el incompleto no. Este nuevo conjunto se denota como x. Una vez definido el conjunto completo, se asume que existe la aplicaci´on x→y(x). As´ı, dado xse conoce y; pero dado ylo que se sabe es que hay un subconjunto xy⊂x con probabilidad distinta de cero de ocurrir. Dado que la distribuci´on de probabilidad de ydepende del conjunto de par´ametros λ, entonces la distribuci´on de probabilidad del conjunto x, de forma natural, depende tambi´en de los mismos. La estrategia es incluir en la ecuaci´on (2.22) la dependencia en x. Usando el teorema de Bayes, la distribuci´on de probabilidad f(y|λ) puede escribirse f(y|λ) = f(y|x, λ)·f(x|λ) f(x|y, λ).(2.23) Teniendo en cuenta que dado x, el conjunto yqueda bien definido f(y|x, λ) = f(y|x),(2.24) por lo que el logaritmo neperiano de (2.23) queda ln f(y|λ) = ln f(y|x) + ln f(x|λ)−ln f(x|y, λ).(2.25) Los distintos t´erminos dependen de x, pero la funci´on resultante tras realizar las operaciones no. Dado que el conjunto completo no es observable, no se pueden tener datos del mismo, por lo que incluirlo directamente har´ıa el problema intratable posteriormente. Lo que si puede conocerse es su esperanza matem´atica o la de cualquier funci´on del mismo dado el par de conjuntos y,˜ λ. Debido a que se busca estimar el valor de λ, los valores particulares del conjunto de par´ametros que se usar´a para determinar la esperanza se ha denotado como ˜ λ. Teniendo esto en cuenta y tomando la esperanza condicional respecto a y,˜ λde (2.25) como funci´on de x, se obtiene por un lado Ex[ln f(y|λ)|y, ˜ λ] = Ex[ln f(y|x)|y, ˜ λ] + Ex[ln f(x|λ)|y, ˜ λ]−Ex[ln f(x|y, λ)|y, ˜ λ],(2.26) por otro Ex[ln f(y|λ)|y, ˜ λ] = Zxy [ln f(y|λ)] ·f(x|y, ˜ λ)dx = ln f(y|λ),(2.27) 16 CAP´ ITULO 2. MATERIAL Y M´ ETODO y juntando ambas ln f(y|λ) = Ex[ln f(y|x)|y, ˜ λ] + Ex[ln f(x|λ)|y, ˜ λ]−Ex[ln f(x|y, λ)|y, ˜ λ].(2.28) La funci´on a maximizar en el m´etodo MAP (2.22) se puede expresar, incluyendo el conjunto completo xde forma tratable, como sigue P(λ|y) = Ex[ln f(y|x)|y, ˜ λ] + Ex[ln f(x|λ)|y, ˜ λ]−Ex[ln f(x|y, λ)|y, ˜ λ] + ln f(λ).(2.29) El t´ermino E[ln f(y|x)|y, ˜ λ] no depende de λ, por lo que la funci´on que hay que maximizar es ˆ P(λ|y)≡ Q(λ|˜ λ)−Ex[ln f(x|y, λ)|y, ˜ λ],(2.30) en la cual Q(λ|˜ λ) = Ex[ln f(x|λ)|y, ˜ λ] + ln f(λ).(2.31) Es importante resaltar que, mientras cada uno de los t´erminos a la derecha del igual de (2.30) dependen del conjunto ˜ λusado, su diferencia no. Se busca entonces el conjunto de valores de los par´ametros que maximice ˆ P, es decir el conjunto ˆ λque haga ˆ P(ˆ λ|y)>ˆ P(λ|y),∀λ6=ˆ λ, (2.32) lo que se traduce, teniendo en cuenta (2.30), en que el conjunto ˆ λes el que hace que Q(ˆ λ|˜ λ)−Ex[ln f(x|y, ˆ λ)|y, ˜ λ]>Q(λ|˜ λ)−Ex[ln f(x|y, λ)|y, ˜ λ],∀λ6=ˆ λ, ∀˜ λ. (2.33) o lo que es lo mismo Q(ˆ λ|˜ λ)−Q(λ|˜ λ)>Ex[ln f(x|y, ˆ λ)|y, ˜ λ]−Ex[ln f(x|y, λ)|y, ˜ λ],∀λ6=ˆ λ, ∀˜ λ. (2.34) Adem´as, si se toma ˜ λ6=ˆ λentonces puede escribirse sin p´erdida de generalidad Q(ˆ λ|˜ λ)−Q(˜ λ|˜ λ)>Ex[ln f(x|y, ˆ λ)|y, ˜ λ]−Ex[ln f(x|y, ˜ λ)|y, ˜ λ],∀˜ λ6=ˆ λ. (2.35) En el ap´endice B, identificando f(x) = f(x|y, ˜ λ) y g(x) = f(x|y, ˆ λ), est´a demostrado que Ex[ln f(x|y, ˆ λ)|y, ˜ λ]−Ex[ln f(x|y, ˜ λ)|y, ˜ λ]≥0,(2.36) para cualquier par de conjuntos ˜ λ, ˆ λ. La igualdad solo se alcanza cuando f(x|y, ˜ λ) = f(x|y, ˆ λ), por lo que en este caso nunca se dar´a ya que ˜ λ6=ˆ λ. As´ı para el conjunto ˆ λque maximiza ˆ P(λ|y) se cumple, combinando (2.35) y (2.36), que Q(ˆ λ|˜ λ)>Q(˜ λ|˜ λ),∀˜ λ6=ˆ λ. (2.37) Es decir, el conjunto ˆ λes un m´aximo global de Q(λ|˜ λ) para cualquier conjunto ˜ λ. Entonces, encontrar el m´aximo global de ˆ P(λ|y) es equivalente a encontrar el m´aximo global de Q(λ|˜ λ) para cualquier conjunto ˜ λ. La dificultad reside en la forma de encontrar dicho m´aximo, porque hay que tener en cuenta todos los posibles ˜ λ. 2.4. M´ ETODOS DE RECONSTRUCCI´ ON 17 El algoritmo EM propone una soluci´on iterativa a este problema. Dado el conjunto correspondiente a los valores de los par´ametros en la n-´esima iteraci´on λ(n), los valores de los par´ametros en la siguiente iteraci´on, λ(n+1), se obtienen en dos pasos Paso E: Obtener Ex[ln f(x|λ)|y, λ(n)], Paso M: Resolver m´ax λ(n+1) Q(λ|λ(n)) = Ex[ln f(x|λ)|y, λ(n)] + ln f(λ).(2.38) Como el m´aximo global no puede encontrarse de forma directa, en cada iteraci´on lo que se hace es obtener, a partir de unos valores dados de los par´ametros, λ(n), unos nuevos valores λ(n+1) que cumplen que maximizan Q(λ|λ(n)), por lo que cumplen tambi´en Q(λ(n+1)|λ(n))>Q(λ(n)|λ(n)).(2.39) Teniendo en cuenta la relaci´on (2.36) y la definici´on (2.30) en cada iteraci´on lo que se consigue entonces es que ˆ P(λ(n+1)|y)>ˆ P(λ(n)|y),(2.40) por lo que realizando el proceso repetidas veces al final se puede encontrar el conjunto de valores que maximice ˆ P(λ|y). En [11] puede encontrarse un estudio de dos aspectos esenciales de la convergencia del algoritmo EM. Este algoritmo es especialmente ´util cuando los par´ametros son independientes a priori, esto es f(λ) = Y u f(λu),(2.41) y el conjunto completo xse elige de forma que Q(λ|˜ λ) = X uQ(λu|˜ λ),(2.42) por lo que el paso M de (2.38) resulta en una ecuaci´on cerrada para cada λu. Algoritmo de esperanza-maximizaci´on generalizado Si bien no se cumple (2.42) , entonces el paso de maximizaci´on devuelve un sistema de ecuaciones para {λu}. Como se ha mostrado anteriormente, resolver un sistema de ecuaciones conlleva aplicar un algoritmo iterativo. En definitiva, el algoritmo EM ser´ıa un algoritmo iterativo de optimizaci´on, donde cada iteraci´on contiene otro algoritmo iterativo de optimizaci´on, lo cual es costoso computacionalmente. La clave del algoritmo EM es que, tras cada iteraci´on, se cumple (2.39), que es lo que lleva a (2.40). De esta manera es posible dise˜nar otros algoritmos que aseguren que tras cada iteraci´on se cumple (2.39) sin tener que maximizar Q(λ|˜ λ), evitando tener que resolver un sistema de ecuaciones. Este tipo de algoritmos se denotan, en ingl´es, como generalized expectation-maximization algorithms (GEM). 18 CAP´ ITULO 2. MATERIAL Y M´ ETODO 2.4.3. Reconstrucci´on de mapa de dosis usando el m´etodo estad´ıstico El conjunto de medidas experimentales DT ={DTv}posee una incertidumbre procedente del propio proceso de medida. La incertidumbre hace que dichas medidas deban modelarse como variables aleatorias con una cierta distribuci´on estad´ıstica, la cual depender´a de la naturaleza estad´ıstica de la fuente de incertidumbre. Las tiras de silicio del detector est´an calibradas a dosis en agua con factores de calibraci´on W={Wv}, esto es, experimentalmente lo que se mide es la carga colectada por cada tira en cada posici´on angular Q={Qv}de modo que DTv=Wv·Qv=q·Wv·Nv,∀v(2.43) donde qes la carga del electr´on y Nves el n´umero de pares electr´on-hueco creados en la tira v. Suponiendo que los factores de calibraci´on est´an perfectamente definidos, la fuente de incertidumbre en DT proviene de la incertidumbre en N={Nv}. La creaci´on de pares electr´on-hueco en un semiconductor es un proceso estad´ıstico. Cuando un electr´on llega al semiconductor interacciona tanto con los ´atomos para ionizarlos y producir pares electr´on-hueco, como con la propia red haci´endola vibrar y produciendo fonones. Por ello la estad´ıstica que sigue la producci´on de pares no es poissoniana, sino que hay que introducir el factor de Fano, F. Cuando un cierto electr´on µllega al detector y deposita una energ´ıa Eµ, si la energ´ıa necesaria para crear un par electr´on-hueco en el semiconductor es E0, el n´umero de pares que crea es nµ=Eµ E0 ,(2.44) con una varianza σ2(nµ) = F·nµ.(2.45) Si a la tira vllegan una serie de electrones con distintas energ´ıas, el n´umero total de pares que se crean en dicha tira es Nv=X µ nµ.(2.46) Como el n´umero de pares que crea cada electr´on es estad´ısticamente independiente del resto, la varianza del n´umero total de pares que se crean en la tira ves σ2(Nv) = σ2 X µ nµ!=X µ σ2(nµ) = X µ F·nµ=FX µ nµ=F·Nv.(2.47) En el silicio, F∼0.1, por lo que la varianza en el n´umero total de pares creados en una tira es menor que si la distribuci´on que siguiera fuese poissoniana, en la cual σ2(Nv) = Nv. Ahora bien si Nv>> F, la incertidumbre relativa a la medida (Nv) = σ(Nv) Nv =sF Nv << 1,(2.48) por lo que el factor de Fano puede omitirse y tratar Nvcomo una variable aleatoria con distribuci´on poissoniana de media Nv=1 q·Kv AP AT X u Guv ·DPu,(2.49) 2.4. M´ ETODOS DE RECONSTRUCCI´ ON 19 por lo que DTv=AP AT X v Guv ·DPu.(2.50) Ahora bien, el conjunto de factores de calibraci´on Wse obtiene tambi´en de forma experimental, por lo que est´a dotado de incertidumbre. As´ı Wes otra fuente de incertidumbre. La varianza de la dosis media en la tira vse calcula entonces teniendo en cuenta (2.43), obteniendo σ2(DTv) = DT2 v· σ2(Kv) K2 v +σ2(Nv) N2 v!.(2.51) Si la incertidumbre relativa de Kves mayor que la de Nvo ambas son parecidas, la dosis media en la tira vsigue la relaci´on DTv=AP AT X u Guv ·DPu+δv,(2.52) donde δves una variable aleatoria que sigue una distribuci´on normal de media nula y varianza dada por (2.51). As´ı DTves una variable aleatoria con distribuci´on normal de media dada por (2.50) y varianza dada por (2.51). Una vez descritas las distribuciones de probabilidad de las medidas experimentales, se aplica el m´etodo desarrollado en el apartado anterior para obtener el conjunto de dosis media en los distintos p´ıxeles DP ={DPu}a partir del conjunto de medidas experimentales DT. Para la reconstrucci´on de mapa de dosis se identifican respecto al caso gen´erico y→DT, λ→DP. (2.53) En primer lugar se aplicar´a el algoritmo EM a las dos distribuciones de probabilidad descritas sin incluir distribuci´on de probabilidad a priori. Hay que aplicar dicho algoritmo ya que en la expresi´on (2.50) aparecen los par´ametros en un sumatorio, por lo que el m´etodo de m´axima verosimilitud desembocar´a en un sistema de ecuaciones. Finalmente se incluye informaci´on a priori que reproduzca correlaciones locales entre pixeles y se desarrolla un algoritmo GEM para reconstruir el mapa de dosis sin necesidad de resolver ning´un sistema de ecuaciones. Aplicaci´on del algoritmo EM a la distribuci´on poissoniana El primer paso es definir el conjunto completo, denotado de forma gen´erica por x. En este caso el conjunto completo ser´a n={nuv}, donde nuv es, del n´umero total de pares creados en la tira v, el n´umero de pares electr´on-hueco que se crean en el pixel u. Esto es Nv=X u nuv.(2.54) Al igual que para Nv, se asume que nuv >> F, por lo que es una variable aleatoria con distribuci´on poissoniana de media ¯nuv =1 q·Wv·AP AT ·Guv ·DPu.(2.55) 26 CAP´ ITULO 2. MATERIAL Y M´ ETODO Los siguientes pasos son α2=1 α1        α2<0,si α1<0 α2≥1,si 0 < α1≤1 0< α2<1,si α1>1 (2.97) α3=α2−1       α3<0,si α1<0 α3≥0,si 0 < α1≤1 α3<0,si α1>1 (2.98) α4=α3+|α3| 2       α4= 0,si α1<0 α4=α3≥0,si 0 < α1≤1 α4= 0,si α1>1 (2.99) α5=α4+ 1        α5= 1,si α1<0 α5= 1/α1≥1,si 0 < α1≤1 α5= 1,si α1>1 (2.100) El m´aximo valor permitido se calcula entonces αm´ax =1 α5        αm´ax = 1,si α1<0 αm´ax =α1,si 0 < α1≤1 αm´ax = 1,si α1>1 ,(2.101) que es lo que se quer´ıa conseguir en (2.96). Una vez conocida la forma de calcular los distintos coeficientes, solo queda describir los pasos que forman este algoritmo GEM. En cada iteraci´on, a partir de DP(n)yDPEM(n+1), se visita pixel a pixel y se siguen los siguientes pasos PASO 1. Calcular los coeficientes Π1, Π2y Π3. PASO 2. Calcular αm´ax. PASO 3. Inicializar Π4= Π1−1 y α=αm´ax. PASO 4. Mientras que Π4<Π1 PASO 4.1. Calcular DP(n+1) u= (1 −α)DP(n) u+α·Π3. PASO 4.2. Calcular Π4con el valor DP(n+1) ucalculado en el paso anterior PASO 4.3. Hacer α=α/2. PASO 5. Hacer DPu=DP(n+1) u El paso 4 es el que asegura que se cumpla (2.76). Para ello, partiendo de α=αm´ax, actualiza el valor del pixel. Si este valor no lleva a que Π4>Π1, entonces toma un nuevo valor de αque es la mitad que el anterior y vuelve a actualizar el pixel para dicho valor. El proceso se repite hasta hallar el valor de αadecuado. 2.4. M´ ETODOS DE RECONSTRUCCI´ ON 27 Por ´ultimo, solo queda especificar los potenciales que se usar´an para realizar la reconstrucci´on en este trabajo. Solo tendremos en cuenta interacci´on de pares en un vecindario de orden dos. La funci´on energ´ıa ser´a de la forma U(DP) = X c∈C Vc(DP) = X uX u0>u u0∈Ωu κuu0·V(DPu, DPu0),(2.102) donde Ωues el vecindario del pixel uy el coeficiente κuu0toma los valores κuu0=   1,si u0est´a en la vertical o la horizontal de u 1/√2,si u0est´a en la diagonal de u.(2.103) En [13] se encuentra una breve discusi´on sobre las propiedades que deben de cumplir las funciones de potencial. En este trabajo se ver´a el efecto que producen cuatro potenciales distintos V1(DPu, DPu0) = (DPu−DPu0)2,(2.104) V2(DPu, DPu0) =    (DPu−DPu0)2/2δ, si |DPu−DPu0| ≤ δ |DPu−DPu0|−δ/2,otro caso ,(2.105) V3(DPu, DPu0) = ln "1 + (DPu−DPu0)2 δ2#,(2.106) V4(DPu, DPu0) = (DPu−DPu0)2 δ2+ (DPu−DPu0)2.(2.107) El primero de ellos, V1(·), penaliza la diferencia entre valores de los pixeles de manera que cuanto mayor es la diferencia, mayor es la penalizaci´on, por lo que favorecer´a que la imagen sea suave en todos sus puntos. La funci´on de Hubert, V2(·), se comporta como una funci´on cuadr´atica cuando la diferencia es peque˜na y de forma lineal cuando la diferencia sobrepasa el umbral determinado por δ. El comportamiento lineal hace que la tasa de penalizaci´on tras el umbral no cambie, por lo que no diferencia entre variaciones suaves y abruptas, no penalizando as´ı la existencia de bordes. En las dos ´ultimas, V3(·) y V4(·), la tasa de penalizaci´on disminuye una vez pasado el umbral dado por δ, favoreciendo as´ı la presencia de bordes. 28 CAP´ ITULO 2. MATERIAL Y M´ ETODO Cap´ıtulo 3 Resultados En este cap´ıtulo se muestran los resultados del modelado del detector as´ı como la reconstrucci´on de varios mapas de dosis usando los m´etodos EM y GEM. 3.1. Modelado del detector Se muestra un ejemplo del c´alculo del valor de los coeficientes Gijklm en la figura 3.1. Se comprueba que el algoritmo propuesto calcula dichos coeficientes correctamente. Los p´ıxeles que est´an contenidos completamente en la tira tienen Gijklm = 1, mientras que los que est´an al completo fuera tienen Gijklm = 0. Para los p´ıxeles que tienen parte de su ´area dentro y parte fuera de la tira, el valor del coeficiente se acerca m´as a la unidad cuanto mayor sea el ´area del mismo dentro de la tira. Para comprobar que la respuesta del detector modelada es correcta, se muestra la lectura del mismo para un mapa de dosis sencillo en distintas posiciones angulares en la figura 3.2. En el caso θ= 0◦la lectura de todas las tiras del detector con m= 1 es la misma, puesto que, como se ve en la figura 3.2(a), todas tienen la mitad del ´area la parte del mapa de dosis m´as baja y la otra mitad en la parte del mapa de dosis m´as alta, por lo que la dosis media en cada tira ser´a la mitad de esos dos valores. En el caso del detector con m= 2, ver figura 3.2(b), la primera mitad de las tiras est´an en la zona de dosis m´as baja y la segunda mitad en la zona de dosis m´as alta por lo que se tendr´an dos respuestas distintas. Como el mapa de dosis a lo largo de las tiras es homog´eneo, la dosis media en cada tira ser´a igual al valor de la dosis de la zona en la que est´e. Todo esto puede comprobarse en la figura 3.2(c). En el caso θ= 50◦las lecturas tienen un comportamiento diferente. Para el detector con m= 1, se observa en la figura 3.2(d), que tanto las tres primeras tiras como las tres ´ultimas tienen su ´area completa en una zona de dosis homog´enea, por lo que su medida ser´a igual al valor de la dosis de la zona en la que est´en. La cuarta tira tiene casi todo su ´area en la zona de dosis m´as alta, pero conforme se avanza en n´umero de tira, cada vez el ´area de la tira en la zona de dosis m´as alta se va reduciendo. Esto conlleva a que la dosis 29 30 CAP´ ITULO 3. RESULTADOS media, seg´un se avanza en el n´umero de tira, vaya decreciendo. Esto pasa tambi´en para todas las tiras en el detector con m= 2, que conforme se avanza en el n´umero de la tira hay m´as ´area de la tira en la zona de dosis m´as baja como se aprecia en la figura 3.2(e), por lo que la dosis media en la tira decrece con el n´umero de tira. Efectivamente esto es lo que se ve en la figura 3.2(f), donde adem´as se aprecia que la dosis en las tiras en el detector con m= 1 decrecen mucho m´as r´apido que en las tiras del detector con m= 2. (a) G(θ= 35◦, l = 8, m = 1) (b) G(θ= 35◦, l = 8, m = 1). Ampliado. -50 -25 0 25 50 X(mm) -50 -40 -30 -20 -10 0 10 20 30 40 50 Y(mm) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 G(θ=70º, l=18, m=2) (c) G(θ= 70◦, l = 18, m = 2) -15 -10 -5 0 5 X(mm) -10 -8 -6 -4 -2 0 2 4 6 8 10 Y(mm) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 G(θ=70º, l=18, m=2) (d) G(θ= 70◦, l = 18, m = 2). Ampliado. Figura 3.1: Valores de los coeficientes geom´etricos para una tira de cada uno de los dos detectores en una posici´on angular determinada. 3.1. MODELADO DEL DETECTOR 31 -50 0 50 X(mm) -50 -40 -30 -20 -10 0 10 20 30 40 50 Y(mm) 0 10 20 30 40 50 60 70 80 90 100 Dosis (a.u.) 132 (a) Detector con m= 1 en posici´on θ= 0◦ -50 0 50 X(mm) -50 -40 -30 -20 -10 0 10 20 30 40 50 Y(mm) 0 10 20 30 40 50 60 70 80 90 100 Dosis (a.u.) 1 32 (b) Detector con m= 2 en posici´on θ= 0◦ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 Nº Tira 50 60 70 80 90 100 Dosis(a.u.) Detector 1 Detector 2 (c) Respuesta de los detectores al mapa de dosis para θ= 0◦ 32 CAP´ ITULO 3. RESULTADOS (d) Detector con m= 1 en posici´on θ= 50◦ (e) Detector con m= 2 en posici´on θ= 50◦ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 Nº Tira 40 50 60 70 80 90 100 110 Dosis(a.u.) Detector 1 Detector 2 (f) Respuesta de los detectores al mapa de dosis para θ= 50◦ Figura 3.2: Lectura de los detectores para un caso sencillo. 3.2. RECONSTRUCCI´ ON DE MAPA DE DOSIS 33 3.2. Reconstrucci´on de mapa de dosis Con los algoritmos EM y GEM se pueden realizar una gran cantidad de pruebas distintas para conocer las ventajas y las limitaciones de dichos m´etodos de reconstrucci´on. En esta secci´on se recogen las reconstrucciones de distintos mapas para poner de manifiesto las caracter´ısticas de los mismos que se consiguen reproducir, as´ı como aquellas que necesitan de nuevos m´etodos en los que hay que investigar para obtener reconstrucciones de mayor calidad en un futuro. 3.2.1. Reconstrucci´on de mapas artificiales sin informaci´on a priori usando el algoritmo EM Para el caso de la reconstrucci´on con distribuci´on poissoniana, como experimentalmente los coeficientes Wuson pr´acticamente iguales, se tomar´a Wv=Wv0=cte ∀v, v0. Para el caso de la distribuci´on normal, la desviaci´on est´andar ser´a σv=DTv·v, donde ves el error relativo. Se tomar´an v=v0=cte ∀v, v0. Todos los mapas se reconstruyen partiendo de un mapa homog´eneo. Se reconstruye, en primer lugar, un mapa con las condiciones m´as ideales posibles. Estas son que sea nulo fuera del c´ırculo inscrito en el detector, que tenga una variaci´on suave y sin incertidumbre estad´ıstica. El mapa original se muestra en la figura 3.3. La reconstrucci´on se realiza con la misma cuadr´ıcula (tama˜no y posici´on de los p´ıxeles) que el mapa reconstruido. La reconstrucci´on se muestra en la figura 3.4 y las caracter´ısticas de la misma se muestran en el pie de figura. Se obtiene una reconstrucci´on pr´acticamente id´entica, con ligeras diferencias en los valores de los p´ıxeles. Se estudia que ocurre si se toma un tama˜no del pixel m´as peque˜no que el mapa original. El resultado se muestra en la figura 3.5. A menor tama˜no del pixel, mayor es la resoluci´on, pero tambi´en se hace mayor el ruido. En este caso, dividiendo cada pixel del mapa original en cuatro, el mapa reconstruido se difumina. Esto se debe a que pueden existir distintas combinaciones de valores medios en las cuatro ´areas de los nuevos p´ıxeles que mantengan el valor medio en el ´area total del pixel original. Tomando el valor medio de cuatro en cuatro p´ıxeles en la figura 3.5, se obtiene la figura 3.6, donde puede verse que se recupera la forma de la figura 3.4. Por tanto una de las ventajas de este m´etodo es que, con los mismos datos experimentales, pueden reconstruirse mapas con distinto tama˜no de pixel. Si en alguna zona del mapa la resoluci´on elegida es mayor que la original, esa zona aparecer´a difuminada en la reconstrucci´on, pero puede recuperarse su resoluci´on original promediando en los p´ıxeles de la zona del mapa difuminado. Se realiza ahora el estudio de la reconstrucci´on de un mapa tambi´en nulo en el exterior del c´ırculo circunscrito en el detector y sin incertidumbre estad´ıstica, pero en este caso est´a formado por una zona homog´enea que acaba en un gradiente alto de dosis, es decir, con bordes muy definidos. Este mapa se recoge en la figura 3.7 y la reconstrucci´on en la figura 3.8. Por un lado el algoritmo EM sin informaci´on a priori no reproduce bien las zonas homog´eneas, sino que aparecen numerosas irregularidades en ellas. Por otro el 34 CAP´ ITULO 3. RESULTADOS gradiente es mucho menos abrupto que el original. Entonces es necesario aplicar otros m´etodos para obtener mapas con zonas homog´eneas y bordes bien definidos. Una de las caracter´ısticas ya mencionadas que diferencia la reconstrucci´on de mapa de dosis de la reconstrucci´on de imagen m´edica y por la que no se puede aplicar los m´etodos basados en la transformada de Radon es que el mapa de dosis es no nulo en el exterior del c´ırculo inscrito al detector. El mapa de la figura 3.9 posee esta caracter´ıstica y en la figura 3.10 se presenta la reconstrucci´on para las dos distribuciones. Aunque el mapa sea no nulo en el exterior del c´ırculo inscrito en el detector es posible realizar la reconstrucci´on. Se distinguen dos zonas, la interior a la circunferencia inscrita en el detector y la zona que queda entre la circunferencia inscrita y la circunscrita. La reconstrucci´on en la primera zona es de mayor calidad que en la segunda. Esto ocurre debido a que el c´ırculo inscrito se mide por completo en todas las posiciones angulares, mientras que en una posici´on angular determinada solo se mide una parte de la zona que queda entre el c´ırculo inscrito y el circunscrito. Comparando las figuras 3.8 y 3.10(a), se aprecia que en ambas existen irregularidades parecidas en la zona homog´enea. En el segundo caso, se obtiene un valor de dosis en la zona homog´enea m´as peque˜no que en el primer caso y respecto al gradiente de la dosis, en el segundo caso es m´as suave que en el primero. Existen diferencias tambi´en en la reconstrucci´on usando el m´etodo EM para la distribuci´on de Poisson , figura 3.10(a), y para la distribuci´on normal, figura 3.10(b). En el caso de la distribuci´on normal, las irregularidades son m´as acusadas, pero el valor de la dosis y el gradiente se parecen m´as al mapa real que en el caso de la distribuci´on poissoniana. -50 -25 0 25 50 X (mm) -50 -40 -30 -20 -10 0 10 20 30 40 50 Y (mm) 0 10 20 30 40 50 60 70 80 90 100 Dosis (a.u.) Figura 3.3: Mapa nulo en el exterior del c´ırculo inscrito en el detector, con variaci´on suave y sin incertidumbre estad´ıstica. Tama˜no: 50x50 p´ıxeles. ´ Area del pixel: 2x2 mm2. 3.2. RECONSTRUCCI´ ON DE MAPA DE DOSIS 35 -50 -25 0 25 50 X (mm) -50 -40 -30 -20 -10 0 10 20 30 40 50 Y (mm) 0 10 20 30 40 50 60 70 80 90 100 Dosis (a.u) Figura 3.4: Reconstrucci´on del mapa de la figura 3.3 usando el M´etodo EM con distribuci´on poissoniana. Medidas experimentales tomadas entre 0◦y 85◦con paso de 5◦. Tama˜no: 50x50 p´ıxeles. ´ Area del pixel: 2x2mm2. N´umero de iteraciones: 15. Tiempo de reconstrucci´on: 3.29s. -50 -25 0 25 50 X (mm) -50 -40 -30 -20 -10 0 10 20 30 40 50 Y (mm) 0 10 20 30 40 50 60 70 80 90 100 Dosis (a.u) Figura 3.5: Reconstrucci´on del mapa de la figura 3.3 usando el M´etodo EM con distribuci´on poissoniana. Medidas experimentales tomadas entre 0◦y 85◦con paso de 5◦. Tama˜no: 100x100 p´ıxeles. ´ Area del pixel: 1x1mm2. N´umero de iteraciones: 15. Tiempo de reconstrucci´on: 12.81s. 42 CAP´ ITULO 3. RESULTADOS 3.2.2. Reconstrucci´on de mapas artificiales con informaci´on a priori usando el algoritmo GEM El algoritmo GEM desarrollado en este trabajo es un algoritmo primario para resolver problemas de m´aximo a posteriori con informaci´on a priori que modele las correlaciones locales entre p´ıxeles. Con los resultados de esta secci´on se pretende comprobar que la inclu- si´on de informaci´on a priori puede ayudar a solventar los problemas en la reconstrucci´on planteados en el apartado anterior. El primer paso es comprobar la dependencia de la reconstrucci´on con el par´ametro β. Para ello se usar´a el potencial V1(·), que favorece el suavizado, para el caso del mapa calculado por el planificador. Es de especial inter´es estudiar el caso β→ ∞, ya que este se corresponde con el caso sin informaci´on a priori, y as´ı se puede comprobar si el algoritmo GEM tiende a la misma soluci´on que el algoritmo EM. Este es el caso de la figura 3.14(a), donde se comprueba que se obtiene una soluci´on con las mismas caracter´ısticas de ruido que en la figura 3.13(c). En las figuras 3.14(b), (c) y (d) se observa el efecto de βen la reconstrucci´on. Cuanto m´as peque˜no se hace, mayor es la importancia de la informaci´on a priori en la reconstrucci´on y por eso en este caso m´as se suaviza la imagen. Es necesario entonces estimar el valor de βque mejor de adapte al problema, porque puede ocurrir que se tomen valores que acaben deformando el mapa como en 3.14(d). Por ´ultimo se aplica la informaci´on a priori con los potenciales V2(·), V3(·) y V4(·) al caso de zona homog´enea con bordes definidos. El resultado se presenta en la figura 3.15. En los tres casos se obtienen mapas con menos irregularidades en la zona homog´enea que en el caso EM, y respetando los bordes. Esto pone de manifiesto que usar informaci´on a priori en la reconstrucci´on conduce a mapas que se asimilan m´as al original. 3.2.3. Reconstrucci´on del mapa de dosis de un tratamiento real sin informaci´on a priori usando el algoritmo EM Por ´ultimo se muestra la reconstrucci´on de un mapa de dosis a partir de medidas experimentales realizadas en el Hospital Universitario Virgen Macarena. Las medidas se corresponden al mismo tratamiento calculado por el planificador 3.11, por lo que se espera que tenga una forma similar (no tienen por qu´e coincidir las dimensiones). El resultado se muestra en la figura 3.16. Adem´as se a˜nade la reconstrucci´on en dos casos adversos que se pueden dar en la pr´actica. Uno es que una o m´as tiras dejen de funcionar correctamente y haya que prescindir de ellas en la reconstrucci´on 3.17, de hecho las tiras 12, 13 y 14 se comportan de forma an´omala para este caso y son la que se eliminan. Y el otro es que se necesite una reconstrucci´on r´apida y se tomen solo la mitad de las medidas que en el caso anterior 3.18. Mientras que en la primera situaci´on lo que se tiene es una p´erdida de informaci´on en una zona localizada del mapa marcada con el rect´angulo fucsia, en la segunda lo que se obtiene es un mapa con menor resoluci´on. Pero en ambos es posible realizar la reconstrucci´on, lo cual verifica la flexibilidad del m´etodo de reconstrucci´on ante casos adversos. 3.2. RECONSTRUCCI´ ON DE MAPA DE DOSIS 43 -40 -20 0 20 40 X (mm) -40 -30 -20 -10 0 10 20 30 40 Y (mm) 0 10 20 30 40 50 60 70 80 90 100 Dosis (cGy) (a) β= 1020. Tiempo de reconstrucci´on: 144.08s. -40 -20 0 20 40 X (mm) -40 -30 -20 -10 0 10 20 30 40 Y (mm) 0 10 20 30 40 50 60 70 80 90 100 Dosis (cGy) (b) β= 5 ·105. Tiempo de reconstrucci´on: 337.80s. 44 CAP´ ITULO 3. RESULTADOS -40 -20 0 20 40 X (mm) -40 -30 -20 -10 0 10 20 30 40 Y (mm) 0 20 40 60 80 100 120 Dosis (cGy) (c) β= 105. Tiempo de reconstrucci´on: 476.23s. -40 -20 0 20 40 X (mm) -40 -30 -20 -10 0 10 20 30 40 Y (mm) 0 5 10 15 20 25 30 35 40 45 Dosis (cGy) (d) β= 500. Tiempo de reconstrucci´on: 653.47s. Figura 3.14: Reconstrucci´on del mapa de la figura 3.12 usando el M´etodo GEM con distribuci´on poissoniana y V1(·). Medidas experimentales tomadas entre 0◦y 85◦con paso de 5◦. Tama˜no: 95x95 p´ıxeles. ´ Area del pixel: 1x1mm2. N´umero de iteraciones: 30. 3.2. RECONSTRUCCI´ ON DE MAPA DE DOSIS 45 -50 -25 0 25 50 X (mm) -50 -40 -30 -20 -10 0 10 20 30 40 50 Y (mm) 0 10 20 30 40 50 60 70 80 90 100 Dosis (a.u.) (a) β= 100. V2(·). Tiempo de reconstrucci´on: 225.41s. -50 -25 0 25 50 X (mm) -50 -40 -30 -20 -10 0 10 20 30 40 50 Y (mm) 0 20 40 60 80 100 Dosis (a.u.) (b) β= 10. V3(·). Tiempo de reconstrucci´on: 145.46s. -50 -25 0 25 50 X (mm) -50 -40 -30 -20 -10 0 10 20 30 40 50 Y (mm) 0 20 40 60 80 100 Dosis (a.u.) (c) β= 10. V4(·). Tiempo de reconstrucci´on: 137.83s. Figura 3.15: Reconstrucci´on del mapa de la figura 3.7 usando el M´etodo GEM con distribuci´on poissoniana. Medidas experimentales tomadas entre 0◦y 85◦con paso de 5◦. Tama˜no: 100x100 p´ıxeles. ´ Area del pixel: 1x1mm2. N´umero de iteraciones: 15. 46 CAP´ ITULO 3. RESULTADOS -40 -20 0 20 40 X (mm) -40 -30 -20 -10 0 10 20 30 40 Y (mm) 0 10 20 30 40 50 60 70 80 90 Dosis (cGy) Figura 3.16: Reconstrucci´on un mapa de dosis de un tratamiento real usando el M´etodo EM con distribuci´on poissoniana. Medidas experimentales tomadas entre 0◦ y 85◦con paso de 5◦. Tama˜no: 184x184 p´ıxeles. ´ Area del pixel: 0.5x0.5mm2. N´umero de iteraciones: 15. Tiempo de reconstrucci´on: 60.20s. -40 -20 0 20 40 X (mm) -40 -30 -20 -10 0 10 20 30 40 Y (mm) 0 10 20 30 40 50 60 70 80 90 Dosis (cGy) Figura 3.17: Reconstrucci´on un mapa de dosis de un tratamiento real usando el M´etodo EM con distribuci´on poissoniana. Medidas experimentales tomadas entre 0◦ y 85◦con paso de 5◦. Tama˜no: 184x184 p´ıxeles. ´ Area del pixel: 0.5x0.5mm2. N´umero de iteraciones: 15. Tiempo de reconstrucci´on: 54.05s. La zona interior al recuadro fucsia es la que cambia al no tener en cuenta las tiras 12,13,14. 3.2. RECONSTRUCCI´ ON DE MAPA DE DOSIS 47 -40 -20 0 20 40 X (mm) -40 -30 -20 -10 0 10 20 30 40 Y (mm) 0 10 20 30 40 50 60 70 80 90 Dosis (cGy) Figura 3.18: Reconstrucci´on un mapa de dosis de un tratamiento real usando el M´etodo EM con distribuci´on poissoniana. N´umero de iteraciones: 15. Medidas experimentales tomadas entre 0◦y 80◦con paso de 10◦. Tama˜no: 184x184 p´ıxeles. ´ Area del pixel: 0.5x0.5mm2. Tiempo de reconstrucci´on: 36.89s 48 CAP´ ITULO 3. RESULTADOS Cap´ıtulo 4 Conclusiones Una vez explicada la motivaci´on del trabajo, desarrolladas las ideas principales y obtenidos los resultados, este cap´ıtulo cierra el trabajo recuperando los objetivos iniciales, viendo cuales se han conseguido e identificando cuales ser´an las futuras l´ıneas de investigaci´on en base a los resultados obtenidos. En los ´ultimos a˜nos los tratamientos de radioterapia se han vuelto m´as complejos. El planificador, que es el programa inform´atico encargado de calcular los mapas de dosis, no siempre est´a preparado para las condiciones extremas de los tratamientos m´as avanzados. Por ello es necesario la verificaci´on de los mismos, esto es, tomar medidas experimentales de la dosis y compararlas con las calculadas por el planificador. Como la planificaci´on de los tratamientos se realiza habitualmente en los planos axiales, es interesante que la verificaci´on se realice en estos planos. Resulta complicado encontrar en el mercado un detector calibrado a dosis que tome medidas cuando el haz es paralelo a su ´area activa. Solo puede hablarse de la pel´ıcula radiocr´omica, pero esta tiene una desventaja, que su lectura no es en tiempo real. Se ha dise˜nado y construido un prototipo de sistema de detecci´on espec´ıficamente para medir mapas 2D de dosis en el plano axial. El sistema est´a basado en dos SSSSD dispuestos en configuraci´on perpendicular y calibrados a dosis en agua, que se insertan el plano axial de un maniqu´ı cil´ındrico de dimensiones de una cabeza que es capaz de rotar alrededor de su eje. Lanzando el tratamiento que se quiere verificar sobre el detector en distintas posiciones angulares, se obtiene la dosis media en agua depositada en el volumen activo de cada tira en las distintas posiciones angulares. Debido a que el mapa de dosis no se anula dentro del ´area correspondiente al ´area activa del detector, los m´etodos basados en la transformada de Radon son dif´ıciles de aplicar en este caso. Para tratar los mapas de dosis se ha discretizado el mapa en p´ıxeles, de manera que se busca obtener es la dosis media en cada pixel. Existe una relaci´on que lleva de la dosis media en los p´ıxeles de un mapa a la dosis media en una tira colocada en una cierta posici´on angular. Esta relaci´on depende de un factor geom´etrico que es necesario calcular para modelar la respuesta del detector. 49 50 CAP´ ITULO 4. CONCLUSIONES Una vez calculado el factor geom´etrico hay que desarrollar un m´etodo que permita obtener la dosis media en cada pixel a partir de la dosis media en cada tira en las distintas posiciones angulares, todo ello en tiempo real. En el m´etodo hay que tener en cuenta la naturaleza estad´ıstica del proceso de medida, por lo que se modelan las medidas experimentales como variables aleatorias que siguen una determinada distribuci´on, donde los par´ametros de la misma son las dosis media en los p´ıxeles. Usando el m´etodo de m´aximo a posteriori puede estimarse los valores de las dosis media en los p´ıxeles como aquella que hace m´axima la funci´on a posteriori. En general este m´etodo conduce a un sistema de ecuaciones que no tiene por qu´e ser lineal. Por ello se desarrolla un algoritmo llamado de esperanza-maximizaci´on, que permite obtener una soluci´on al m´etodo MAP de forma iterativa. Se describe tambi´en la generalizaci´on del m´etodo EM, la cual permite estimar el valor de los par´ametros cuando la funci´on de distribuci´on a priori incluye correlaciones entre los valores de los p´ıxeles. Las reconstrucciones de mapas artificiales realizadas con el m´etodo EM sin tener en cuenta informaci´on a priori arroja muy buenos resultados para los casos en los que existe una variaci´on suave del mapa de dosis. A partir unas mismas medidas experimentales pueden reconstruirse mapas con distinto tama˜no de pixel, pudiendo siempre promediar sobre el mapa de pixeles m´as peque˜nos para obtener el mismo resultado que el mapa de p´ıxeles de mayor tama˜no. En el caso de reconstrucci´on de mapas con zonas homog´eneas y bordes muy definidos (alto gradiente), presenta irregularidades en la zona que deber´ıa ser homog´enea, y no reproduce con exactitud el gradiente. Es capaz de reconstruir im´agenes a las que se le ha aplicado un ruido estad´ıstico, pero es un algoritmo de convergencia lenta y la soluciones a las que converge son caracter´ısticas por la presencia de ruido. El m´etodo GEM incluye informaci´on a priori, a partir de la cual puede incluirse correlaciones locales entre los valores de los p´ıxeles, haciendo que la imagen que se obtenga sea m´as suave y permitiendo la existencia de bordes. Una de las tareas que queda por realizar es desarrollar un algoritmo de estimaci´on del par´ametro βque arroje el valor ´optimo para la reconstrucci´on, ya que una mala elecci´on del mismo puede da˜nar la reconstrucci´on. Cuando se reconstruye un mapa de dosis a partir de medidas experimentales reales usando el m´etodo EM sin informaci´on a priori se obtienen buenos resultados. Adem´as el m´etodo EM es flexible ante casos adversos, permitiendo obtener un mapa de dosis aunque alguna tira no funcione o no se disponga de tiempo para realizar todas las medidas experimentales deseadas. El siguiente paso a seguir ser´ıa incluir informaci´on a prori para este caso. El mapa obtenido por el planificador ofrece mucha informaci´on que puede usarse para desarrollar nuevas formas de los potenciales de la distribuci´on a priori para mejorar la reconstrucci´on. Lo ideal ser´ıa que se pudiese reconstruir los huecos producidos al fallar alguna tira, mejorar la calidad del mapa en el ´area existente entre el c´ırculo inscrito y el circunscrito y reconstruir el mapa m´as all´a del ´area barrida por el detector. Ap´endice A C´alculo de los coeficientes del modelado del detector A.1. C´alculo del coeficiente Cijkαm para un pixel dado ij Se realiza el c´alculo del coeficiente Cijklm para un pixel fijo ij. Para los diferentes c´alculos que se realizan en este apartado se usan las proyecciones de las esquinas del pixel que se denotan por x0 βpara el detector vertical y y0 βpara el horizontal, siendo β=1-4 seg´un la figura 2.4. Estos se calculan usando (2.1) a partir de las coordenadas de las esquinas en el sistema fijo, que se obtienen conociendo la posici´on del centro, la anchura y la altura del pixel. Adem´as, como se manejan l´ıneas verticales y horizontales, en cada una de las figuras ilustrativas que se muestran se especifica el intervalo angular correspondiente a esa figura usando θmcon m= 1 para el detector vertical (x0 α) y con m= 2 para el detector horizontal (y0 α). Lo ´ultimo a tener en cuenta es que se presenta la forma de calcular el coeficiente en el intervalo angular θ∈[0◦,90◦], siendo f´acilmente extensible al resto de cuadrantes. A.1.1. Condici´on de contribuci´on de los p´ıxeles El primer paso es saber si el pixel tienen contribuci´on nula C= 0, si contribuye en su totalidad C= 1 y si por el cruza el borde de una tira C∈(0,1). Para ello simplemente hay que tener en cuenta las proyecciones de las esquinas en el sistema rotado que depende de si la l´ınea es vertical u horizontal y del ´angulo en que est´e posicionado el sistema como se aprecia en la figura A.1. Matem´aticamente esto se puede escribir 51 58 Ap´endice B Desigualdad de Jensen. Aplicaci´on a un caso caracter´ıstico. Sea f(x) una funci´on real no negativa tal que Z∞ −∞ f(x)dx = 1,(B.1) h(x) una funci´on real cualquiera y φuna funci´on real c´oncava sobre el rango de h, la desigualdad de Jensen, v´ease [14], se formula φZ∞ −∞ h(x)f(x)dx≥Z∞ −∞ φ[h(x)]f(x)dx. (B.2) Si h(x) es una funci´on no negativa, entonces φ[h(x)] = ln[h(x)] es c´oncava en el rango de h. As´ı, desigualdad de Jensen establece que ln Z∞ −∞ h(x)f(x)dx≥Z∞ −∞ ln[h(x)]f(x)dx. (B.3) Ahora bien, si se toma h(x) = g(x)/f(x), donde g(x) es una funci´on real no negativa, entonces ln Z∞ −∞ g(x)dx≥Z∞ −∞ ln "g(x) f(x)#f(x)dx. (B.4) Si Z∞ −∞ g(x)dx = 1,(B.5) se sigue que 0≥Z∞ −∞ {ln[g(x)] −ln[f(x)]}f(x)dx, (B.6) o lo que es lo mismo Z∞ −∞ ln[f(x)]f(x)dx −Z∞ −∞ ln[g(x)]f(x)dx ≥0.(B.7) Como f(x) y g(x) son ambas funciones reales no negativas y normalizadas a la unidad, pueden ser identificadas como distribuciones de probabilidad. Usando notaci´on estad´ıstica puede escribirse entonces Ef(x){ln[f(x)]}−Ef(x){ln[g(x)]} ≥ 0,(B.8) donde la igualdad se alcanza cuando g(x) = f(x). 59 60 Ap´endice C Distribuci´on de probabilidad a priori de Gibbs La distribuci´on de probabilidad a priori de Gibbs permite incluir correlaciones locales entre pixeles en la reconstrucci´on. Para describir la forma de estas funciones hay que definir previamente algunos conceptos, todos recogidos en [15]. El primer concepto es el campo aleatorio de Markov (MRF, por sus siglas en ingl´es). Un MRF es un conjunto de variables aleatorias z={zw}asociadas a las posiciones de los elementos que forman un sistema, para las cuales la probabilidad de que se d´e un cierto valor de la variable del elemento u,zu, condicionada por los valores del resto de variables es igual a la probabilidad de se d´e dicho valor de la variable aleatoria del elemento u condicionada solo a los valores de las variables de un peque˜no subconjunto de elementos llamado vecindario. Matem´aticamente un MRF se define P(zu|zw;w6=u) = P(zu|zw;w∈Ωu),(C.1) donde Ωues el vecindario del elemento u. Todos los vecindarios cumplen que si wes vecino de u, entonces ues vecino de w, es decir u0∈Ωu↔u∈Ωu0.(C.2) La forma de los vecindarios depende del sistema tratado y generalmente se suelen clasificar por ´ordenes. En la figura C.1 se muestran distintos ´ordenes de vecindarios para una cuadr´ıcula 2D, que es el sistema estudiado en este trabajo. Una vez establecidos los vecindarios de los distintos elementos del sistema se define un nuevo tipo de subconjunto de elementos, el clique. Es un subconjunto formado por uno o m´as elementos de manera que cada una de ellos es vecino del resto de elementos del clique. Los distintos cliques que pueden darse dependen del sistema tratado y del orden de vecindario elegido. En la figura C.2 se muestran los distintos cliques posibles para vecindarios de distinto orden que se dan en una cuadr´ıcula 2D. 61 ORDEN 3 ORDEN 1ORDEN 0 ORDEN 2 Figura C.1: Vecindarios de distinto orden para una cuadr´ıcula 2D. Los p´ıxeles blancos son los vecinos del pixel negro. ORDEN 1ORDEN 0 ORDEN 2 Figura C.2: Cliques para vecindarios de distinto orden para una cuadr´ıcula 2D 62 El teorema de Hammersley–Clifford establece que un campo aleatorio definido en una red es un MRF si y solo si su distribuci´on de probabilidad se corresponde con una funci´on de Gibbs. Esta ´ultima tiene la forma f(z) = 1 Z(β)exp −H(z) β!,(C.3) donde βes una constante, H(z) se conoce como funci´on de energ´ıa y Z(β) es la funci´on de partici´on (constante de normalizaci´on) que viene dada por Z(β) = Z∞ −∞ exp −H(z) β!dz. (C.4) La funci´on energ´ıa tiene la forma H(z) = X c∈C Vc(z),(C.5) donde Ces el conjunto de los distintos cliques yVc(z) es la funci´on potencial de las variables aleatorias asociadas a los elementos del clique c. La funci´on de Gibbs es una clase de funciones muy interesante a la hora de usarlas como funciones de distribuci´on a priori para incluir correlaciones locales, debido a que la funci´on de partici´on permite elegir cualquier orden de vecindario y cualquier combinaci´on de funciones potenciales. La constante βcontrola la intensidad de la contribuci´on de la informaci´on a priori frente a la informaci´on que da la funci´on de distribuci´on de probabilidad del modelo experimental, por lo que es necesario dise˜nar alg´un algoritmo para su estimaci´on. 63 64 Bibliograf´ıa [1] M.I. Gallardo et al. Sistema y m´etodo de verificaci´on de tratamientos de radioterapia. Patente de invenci´on ES 2409760 B1 (2014). [2] M.A.G. Alvarez et al. Novel dual single sided silicon strip detector chip for radiotherapy verification. Proceedings of Science (X LASNPA) 049: 1-6 (2013). [3] A. Bocci et al. Silicon strip detector for a novel 2D dosimetric method for radiotherapy treatment verification. Nuclear Instruments and Methods in Physics Research A 673: 98–106 (2012). [4] C. Battaglia.Dosimetry studies for radiation therapy with photons and radiobiology using low-energy protons (Tesis doctoral). Universidad de Sevilla, Sevilla (2017). [5] Z. Abou-Haidar.A 2D Silicon Detector System for Complex Radiotherapy Treatment Verification (Tesis doctoral). Universidad de Sevilla, Sevilla (2012) [6] R. Gordon, R. Bender and G.T. Herman.Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography. J. Theor. Biol. 29: 471-481 (1970). [7] P. Gilbert.Iterative methods for the three-dimensional reconstruction of an object from projections. J. Theor. Biol. 36: 105-117 (1972). [8] C. Badea and R. Gordon.Experiments with the nonlinear and chaotic behaviour of the multiplicative algebraic reconstruction technique (MART) algorithm for computed tomography. Phys. Med. Biol. 49 : 1455-1474 (2004). [9] F.W. Scholz.Maximum Likelihood Estimation. Encyclopedia of Statistical Sciences. 7. (2006) [10] A. P. Dempster, N. M. Laird and D. B. Rubin.Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 39.1 : 1-38 (1977). [11] C.F. Jeff Wu.On the convergene properties of the EM algorithm. Ann. Statist. 11.1: 95-103 (1983). [12] T. Hebert and R. Leahy.A Generalized EM Algorithm for 3D Bayesian Reconstruction from Poisson Data Using Gibbs Priors. IEEE Transactions on Medical Imaging, 8.2: 194:202 (1989). 65 [13] K. Lange.Convergence of EM image reconstruction algorithms with Gibbs smoothing. IEEE Transactions on Medical Imaging, 9.4:439-446 (1990). [14] J.L.W.V. Jensen.Sur les fonctions convexes et les in´egalit´es entre les valeurs moyennes. Acta Math. 30: 175-193 (1906). [15] Markov Random Fields and Gibbs Distributions. Disponible online en: http://www.nlpr.ia.ac.cn/users/szli/mrfbook/Chapter1/node9.html 66